Reconstruction of a Local Perturbation in Inhomogeneous Periodic Layers from Partial Near Field Measurements

TL;DR

This paper develops a method to reconstruct local perturbations in inhomogeneous periodic layers from near field measurements, overcoming challenges posed by the perturbation preventing reduction to a single periodic cell.

Contribution

It introduces a reformulation of the inverse scattering problem into quasi-periodic problems on bounded domains, enabling unique solution existence and numerical reconstruction of local perturbations.

Findings

Unique reconstruction of perturbations demonstrated in 2D and 3D.

Numerical algorithms show effective approximation of the perturbed medium.

Analysis of measurement operators highlights ill-posedness of the inverse problem.

Abstract

We consider the inverse scattering problem to reconstruct a local perturbation of a given inhomogeneous periodic layer in $\mathbb{R}^d$, $d=2,3$, using near field measurements of the scattered wave on an open set of the boundary above the medium, or, the measurements of the full wave in some area. The appearance of the perturbation prevents the reduction of the problem to one periodic cell, such that classical methods are not applicable and the problem becomes more challenging. We first show the equivalence of the direct scattering problem, modeled by the Helmholtz equation formulated on an unbounded domain, to a family of quasi-periodic problems on a bounded domain, for which we can apply some classical results to provide unique existence of the solution to the scattering problem. The reformulation of the problem is also the key idea for the numerical algorithm to approximate the solution, which we will describe in more detail. Moreover, we characterize the smoothness of the Bloch-Floquet transformed solution of the perturbed problem w.r.t. the quasi-periodicity to improve the convergence rate of the numerical approximation. Afterward, we define two measurement operators, which map the perturbation to some measurement data, and show uniqueness results for the inverse problems, and the ill-posedness of these. Finally, we provide numerical examples for the direct problem solver as well as examples of the reconstruction in 2D and 3D.

Figures

Tables

Table 1. Table 1: Relative L 2 ​ ( Ω 0 R ) superscript 𝐿 2 subscript superscript Ω 𝑅 0 L^{2}(\Omega^{R}_{0}) -error for Example 1.

		$N=\text{8}$	$N=\text{16}$	$N=\text{32}$	$N=\text{64}$	$N=\text{128}$	$N=\text{256}$
$2^{2\times M}=$	256	6.210e-02	1.741e-02	1.964e-02	1.825e-02	1.826e-02	1.826e-02
$2^{2\times M}=$	1 024	5.862e-02	8.930e-03	6.546e-03	4.473e-03	4.520e-03	4.519e-03
$2^{2\times M}=$	4 096	5.901e-02	9.473e-03	3.871e-03	1.088e-03	1.128e-03	1.127e-03
$2^{2\times M}=$	16 384	5.901e-02	9.843e-03	3.413e-03	2.788e-04	2.786e-04	2.776e-04
$2^{2\times M}=$	65 536	5.894e-02	9.910e-03	3.309e-03	1.662e-04	6.866e-05	6.847e-05

Table 2. Table 2: Time in seconds for Example 1.

		$N=\text{8}$	$N=\text{16}$	$N=\text{32}$	$N=\text{64}$	$N=\text{128}$	$N=\text{256}$
$2^{2\times M}=$	256	$1$	$2$	$3$	$6$	$12$	$26$
$2^{2\times M}=$	1 024	$1$	$2$	$4$	$7$	$14$	$28$
$2^{2\times M}=$	4 096	$4$	$7$	$13$	$25$	$47$	$93$
$2^{2\times M}=$	16 384	$36$	$70$	$139$	$265$	$539$	$1026$
$2^{2\times M}=$	65 536	$643$	$1786$	$3337$	$6737$	$\mathrm{13\hspace{0.17em}724}$	$\mathrm{27\hspace{0.17em}689}$

Table 3. Table 3: Relative L 2 ​ ( Ω 0 R ) superscript 𝐿 2 subscript superscript Ω 𝑅 0 L^{2}(\Omega^{R}_{0}) -error for Example 2.

		$N=\text{8}$	$N=\text{16}$	$N=\text{32}$	$N=\text{64}$	$N=\text{128}$	$N=\text{256}$
$2^{2\times M}=$	256	3.558e-01	5.420e-02	2.157e-02	8.410e-03	8.729e-03	8.722e-03
$2^{2\times M}=$	1 024	3.651e-01	5.247e-02	1.429e-02	1.900e-03	2.200e-03	2.195e-03
$2^{2\times M}=$	4 096	3.631e-01	5.221e-02	1.257e-02	4.486e-04	5.352e-04	5.305e-04
$2^{2\times M}=$	16 384	3.630e-01	5.213e-02	1.213e-02	4.461e-04	1.257e-04	1.236e-04
$2^{2\times M}=$	65 536	3.644e-01	5.195e-02	1.203e-02	5.003e-04	4.726e-05	5.279e-05

Table 4. Table 4: Relative L 2 ​ ( Ω 0 R ) superscript 𝐿 2 subscript superscript Ω 𝑅 0 L^{2}(\Omega^{R}_{0}) -error for Example 3.

		$N^2=\text{16}$	$N^2=\text{64}$	$N^2=\text{256}$
$4^{3\times M}=$	8	7.534e-01	7.591e-01	7.590e-01
$4^{3\times M}=$	64	4.807e-01	5.256e-01	5.260e-01
$4^{3\times M}=$	512	1.216e-01	8.489e-02	1.498e-01
$4^{3\times M}=$	4 096	2.005e-02	4.158e-02	2.779e-02
$4^{3\times M}=$	32 768	6.150e-03	9.793e-03	6.611e-03
$4^{3\times M}=$	262 144	1.577e-03	2.450e-03	-

Table 5. Table 5: Relative L 2 ​ ( Ω 0 R ) superscript 𝐿 2 subscript superscript Ω 𝑅 0 L^{2}(\Omega^{R}_{0}) -error for Example 4.

		$N^2=\text{16}$	$N^2=\text{64}$	$N^2=\text{256}$	$N^2=\text{1\hspace{0.17em}024}$	$N^2=\text{4\hspace{0.17em}096}$
$4^{3\times M}=$	8	3.656e-01	3.775e-01	3.775e-01	3.502e-01	3.490e-01
$4^{3\times M}=$	64	4.769e-01	1.017e-00	5.350e-01	4.858e-01	4.718e-01
$4^{3\times M}=$	512	6.182e-01	1.480e-01	4.499e-02	5.636e-02	4.387e-02
$4^{3\times M}=$	4 096	5.965e-01	1.018e-01	1.708e-02	1.926e-02	8.951e-03
$4^{3\times M}=$	32 768	5.959e-01	8.270e-02	1.716e-02	1.330e-02	-

Equations

$$\Delta u+k^2{n}^{2}u$$

$$u$$

$${\Omega }^R$$

$${\Gamma }_0^R$$

$$\Delta u+k^2{n}^{2}u=-f\text{ in }{\mathbb{R}}_{+}^d,u=0\text{ on }{\Gamma }^0.$$

$$u(x):=\frac{1}{(2\pi {)}^{\text{nicefrac}(d-1)2}}{\int }_{{\mathbb{R}}^{d-1}}{e}^{\mathrm{i}\underline{x}\cdot \xi +\mathrm{i}\sqrt{k^2-|\xi {|}^2}(x_d-R)}\hat{u}(\xi ,R)\mathrm{d}\xi \text{for }{x}_d>R,$$

$$\frac{\partial u}{\partial x_{d}}(\underline{x},R)=\frac{1}{(2\pi)^{\nicefrac{{(d-1)}}{{2}}}}\int_{\mathbb{R}^{d-1}}\mathrm{i}\sqrt{k^{2}-|\xi|^{2}}e^{\mathrm{i}\underline{x}\cdot\xi}\widehat{u}(\xi,R)\;\mathrm{d}\xi=:T(u\big{|}_{\Gamma^{R}})(\underline{x}),$$

$$a_{q}(u,v):=\int_{\Omega^{R}}\nabla u\cdot\nabla\overline{v}-k^{2}n^{2}u\overline{v}\;\mathrm{d}x-\int_{\Gamma^{R}}T(u\big{|}_{\Gamma^{R}})\overline{v}\;\mathrm{d}S=\int_{\Omega^{R}}f\overline{v}\;\mathrm{d}x$$

$$u(\underline{x}+2\pi j,x_d)=e^{-2\pi \mathrm{i}\alpha \cdot j}u(\underline{x},x_d)\text{for all }j\in {\mathbb{Z}}^{d-1}.$$

$${\mathcal{J}}_{{\mathbb{R}}^{d-1}}\varphi (\alpha ,\underline{x},x_d):=\sum _{j\in {\mathbb{Z}}^{d-1}}\varphi (\underline{x}+2\pi j,x_d)e^{2\pi \mathrm{i}\alpha \cdot j}.$$

$${\mathcal{J}}_{{\mathbb{R}}^{d-1}}^{-1}w(\underline{x}+2\pi j,x_d)={\int }_{I}w(\alpha ,\underline{x},x_d)e^{-2\pi \mathrm{i}\alpha \cdot j}\mathrm{d}\alpha ,x\in {\Omega }_0^R,j\in {\mathbb{Z}}^{d-1}.$$

$$u_{\alpha}(\underline{x},x_{d}):=\sum_{j\in\mathbb{Z}^{d-1}}\widehat{(u_{\alpha}\big{|}_{\Gamma_{0}^{R}})}_{j}e^{-\mathrm{i}\alpha_{j}\cdot\underline{x}+\mathrm{i}\beta_{j}(x_{d}-R)}\quad\text{for }x_{d}>R,$$

$${\hat{\varphi }}_j(\alpha ):={\int }_{{\Gamma }_0^R}\varphi (\underline{x})e^{\mathrm{i}\alpha \cdot \underline{x}}\overline{{\psi }^j(\underline{x})}\mathrm{d}\underline{x},\text{where}{\psi }^j(\underline{x}):=\frac{1}{(2\pi {)}^{\text{nicefrac}(d-1)2}}{e}^{-\mathrm{i}j\cdot \underline{x}}.$$

$$T_{\alpha }(\varphi )(\underline{x})=\frac{\mathrm{i}}{(2\pi {)}^{\text{nicefrac}(d-1)2}}\sum _{j\in {\mathbb{Z}}^{d-1}}\sqrt{k^2-|\alpha +j{|}^2}{\hat{\varphi }}_j(\alpha )e^{-\mathrm{i}(\alpha +j)\cdot \underline{x}}.$$

$$\displaystyle\int_{I}\left[\int_{\Omega_{0}^{R}}\nabla_{x}\mathcal{J}u\cdot\nabla_{x}\overline{v}-k^{2}n_{p}^{2}\mathcal{J}u\overline{v}\;\mathrm{d}x-\int_{\Gamma_{0}^{R}}T_{\alpha}(\mathcal{J}u\big{|}_{\Gamma_{0}^{R}})\overline{v}\big{|}_{\Gamma_{0}^{R}}\;\mathrm{d}S\right]\;\mathrm{d}\alpha$$

$$={\int }_I{\int }_{{\Omega }_0^R}(\mathcal{J}f)(\alpha ,\cdot )\overline{v}\mathrm{d}x\mathrm{d}\alpha$$

$${\int }_{{\Omega }^R}\nabla u\cdot \nabla \overline{v}-k^2{n}^{2}u\overline{v}\mathrm{d}x$$

$$={\int }_{{\Omega }^R}\nabla u\cdot \overline{{\mathcal{J}}^{\ast }({\nabla }_x\mathcal{J}v)}-k^2{n}^{2}u\overline{{\mathcal{J}}^{\ast }(\mathcal{J}v)}\mathrm{d}x$$

$$={\int }_I{\int }_{{\Omega }_0^R}{\nabla }_x(\mathcal{J}u)\cdot {\nabla }_x(\overline{\mathcal{J}v})-k^2{n}^2(\mathcal{J}u)(\overline{\mathcal{J}v})\mathrm{d}x\mathrm{d}\alpha$$

$$={\int }_I{\int }_{{\Omega }_0^R}\left({\nabla }_x{u}_{\alpha }\cdot {\nabla }_x{\overline{v}}_{\alpha }-k^2{n}^2{u}_{\alpha }{\overline{v}}_{\alpha }\right)\mathrm{d}x\mathrm{d}\alpha .$$

$${\int }_{{\Gamma }^R}T({\gamma }_{{\Gamma }^R}u){\gamma }_{{\Gamma }^R}\overline{v}\mathrm{d}S={\int }_I{\int }_{{\Gamma }_0^R}(\mathcal{J}T({\gamma }_{{\Gamma }^R}u))(\alpha ,\underline{x})(\overline{\mathcal{J}{\gamma }_{{\Gamma }^R}v})(\alpha ,\underline{x})\mathrm{d}S(\underline{x})\mathrm{d}\alpha .$$

$$T_{\alpha }(\mathcal{J}u)(\alpha ,\cdot )=(\mathcal{J}Tu)(\alpha ,\cdot )\text{for all }u\in H^{\text{nicefrac}12}({\mathbb{R}}^{d-1}).$$

$$(\tilde{J}\varphi )(\alpha ,\underline{x}):=\sum _{j\in {\mathbb{Z}}^{d-1}}\varphi (\alpha +j)e^{-\mathrm{i}\alpha \cdot \underline{x}}{\psi }^j(\underline{x}),(\alpha ,\underline{x},R)\in I\times {\Gamma }_0^R,$$

$$\mathcal{J}\circ Tu=\mathcal{J}\circ {\mathcal{F}}^{-1}\left(\xi \mapsto \mathrm{i}\sqrt{k^2-|\xi {|}^2}\mathcal{F}({\gamma }_{{\Gamma }^R}u)(\xi )\right)=\tilde{J}\left(\xi \mapsto \mathrm{i}\sqrt{k^2-|\xi {|}^2}\mathcal{F}({\gamma }_{{\Gamma }^R}u)(\xi )\right).$$

$$a_{\alpha }(w_{\alpha },v)$$

$$={\int }_{{\Omega }_0^R}({\mathcal{J}}_{{\mathbb{R}}^{d-1}}f)(\alpha ,\cdot )\overline{v}\mathrm{d}x$$

$$\displaystyle-\mathrm{Re}~{}\left(\int_{\Gamma_{0}^{R}}T_{\alpha}({w}_{\alpha}\big{|}_{\Gamma_{0}^{R}})\overline{{w}}_{\alpha}\big{|}_{\Gamma_{0}^{R}}\;\mathrm{d}S\right)=\sum_{\begin{subarray}{c}j\in\mathbb{Z}^{d-1}\\ |\alpha_{j}|>k\end{subarray}}|\beta_{j}|\ |\widehat{{w}}_{j}|^{2}\geq 0,$$

$$\mathrm{Re}{a}_{\alpha }(w_{\alpha },w_{\alpha })\ge ||w_{\alpha }|{|}_{H_{\alpha }^1({\Omega }_0^R)}^2-||1-k^2{n}^2|{|}_{L^2({\Omega }_0^{R_0})}||w_{\alpha }|{|}_{L^2({\Omega }_0^R)}^2.$$

$$\mathrm{Im}~{}\left(\int_{\Gamma_{0}^{R}}T_{\alpha}({w}_{\alpha}\big{|}_{\Gamma_{0}^{R}})\overline{w}_{\alpha}\big{|}_{\Gamma_{0}^{R}}\;\mathrm{d}S\right)=\sum_{\begin{subarray}{c}j\in\mathbb{Z}^{d-1}\\ |\alpha_{j}|<k\end{subarray}}|\beta_{j}|\ |\widehat{{w}}_{j}|^{2}\geq 0.$$

$$0=\mathrm{Im}\left({\int }_{{\Omega }_0^R}{k}^2{n}_p^2|w_{\alpha }{|}^2\mathrm{d}x+{\int }_{{\Gamma }_0^R}{T}_{\alpha }(w_{\alpha }){\overline{w}}_{\alpha }\mathrm{d}S\right)\ge {\int }_{{\Omega }_0^R}{k}^2\mathrm{Im}{n}_p^2|w_{\alpha }{|}^2\mathrm{d}x\ge 0.$$

$$0=\mathrm{Im}\left({\int }_I{\int }_{{\Omega }_0^R}{k}^2{n}_p^2|w_{\alpha }{|}^2\mathrm{d}x+{\int }_{{\Gamma }_0^R}{T}_{\alpha }(w_{\alpha }){\overline{w}}_{\alpha }\mathrm{d}S\mathrm{d}\alpha \right)\ge 0,$$

$$b_{\alpha }(u_{\alpha }^p,v_{\alpha }^p):=$$

$$\displaystyle\quad\quad\quad-\int_{\Gamma_{0}^{R}}S_{\alpha}({u}^{p}_{\alpha}\big{|}_{\Gamma_{0}^{R}})\overline{{v}}^{p}_{\alpha}\big{|}_{\Gamma_{0}^{R}}\;\mathrm{d}S,$$

$$\displaystyle\int_{\Gamma_{0}^{R}}S_{\alpha}({u}^{p}_{\alpha}\big{|}_{\Gamma_{0}^{R}})\overline{{v}}^{p}_{\alpha}\big{|}_{\Gamma_{0}^{R}}\;\mathrm{d}S:=\sum_{j\in\mathbb{Z}^{2}}\mathrm{i}\beta_{j}\widehat{({u}^{p}_{\alpha}\big{|}_{\Gamma_{0}^{R}})}_{j}(\alpha)\overline{\widehat{({v}^{p}_{\alpha}\big{|}_{\Gamma_{0}^{R}})}_{j}}(\alpha).$$

Full text

1 Introduction

The growing industrial interest for micro or nano-structured materials and the resulting challenge to construct an automated non-destructing testing method for the structures is one of the fundamental motivations to study perturbed periodic scattering problems. The direct and inverse scattering problems from unbounded periodic structures is a well-established topic in mathematics, especially if one considers quasi-periodic incident fields. This assumption allows to reduce the problem on the infinite periodic domain into one periodic cell, such that standard techniques for the existence theory and the standard numerical methods for bounded domains can be applied (see, e.g., [BS94], [DF92], [AN92], [BDC95], [Bao94], [Bao95], [Kir93], [Kir95]). If the periodicity is perturbed, or, one uses non-periodic incident fields, such as Gaussian beams, the reduction is typically impossible and one has to treat the problem as a scattering problem for an unbounded rough layer (see, e.g., [HL11], [Hu+15], [Mei+00]). The disadvantage is that for the existence theory one has to assume more regularity for the parameter, which we can avoid by considering the periodicity of the unperturbed parameter and applying the Bloch-Floquet transform to the variational problem to get an alternative problem. There are, however, some approaches for problems on locally perturbed periodic waveguides based on the Bloch-Floquet transform, see [JLF06], [FJ15], [ESZ09].

In this paper, we study the scattering problem formulated in the upper half space $\mathbb{R}_{+}^{d}:=\{x\in\mathbb{R}^{d}:x_{d}>0\}$, $d=2,3$,

[TABLE]

for a locally perturbed inhomogeneous layer, which is described by the refractive index $n^{2}\in L^{\infty}(\mathbb{R}_{+}^{d})$. Applying the Bloch-Floquet transform to decompose the (non-periodic) incident field into its quasi-periodic components, we can reformulate the scattering problem as a family of quasi-periodic scattering problems on a bounded domain. We show equivalence of the two problems and consider the latter to prove existence of the solution to the scattering problem by applying Fredholm theory for the reduced problem. Moreover, we stay in the framework of the equivalent formulation to introduce a numerical method to approximate the solution to the original problem, which is based on [LZ17] and [Zha18], where the algorithm for the sound-soft scattering layer is developed. Considering the regularity of the transformed solution w.r.t. the quasi-periodicity, we are able to improve the convergence rate of the inverse Bloch transform, approximated by the trapezoidal rule, and by choosing an adequate variable transform. For the implementation of the direct problem solver, we use the Finite-Element-Method library deal.II ([Arn+17]). The drawback of this method is that one needs to be able to compute analytically, or numerically, the transformed function of the incident wave. At least for incident point sources and Herglotz wave functions, which are models for Gaussian beams, some semi-analytic expressions are available in [LN15].

In the second part, we consider the inverse scattering problem to reconstruct the local perturbation by analyzing the measurement operator $\Lambda:\mathcal{D}(\Lambda)\to\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),L^{2}(\Omega^{R}_{0}))$, where $\Omega^{R}_{0}$ will be later defined as one periodic cell for $R>0$. The operator maps the perturbation $q$ to the solution operator dependent on $q$, which maps right hand sides in $L^{2}$ supported in $\Omega_{0}^{R_{0}}$ to scattered waves restricted to $\Omega^{R}_{0}$. Furthermore, we consider the second measurement operator $\mathcal{S}:\mathcal{D}(\mathcal{S})\to\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),L^{2}(\partial\Omega^{R}_{0}\cap\{x_{d}=R\}))$ mapping the perturbation to the operator, which maps $L^{2}(\Omega_{0}^{R_{0}})$ right hand sides to the upper trace of the scattered field, also restricted to one periodic cell. We show injectivity of $\mathcal{S}$ in the case that $d=3$, as long as the parameters are twice differentiable, and the whole trace on $\Gamma^{R}$ is given as data, considering the complex geometrical optics (see, e.g., [SU87]). In addition, we show the injectivity of $\Lambda$ (without these restrictions). Moreover, we compute the Fréchet derivative of these operators and show that the Fréchet derivative is a compact operator and the so-called tangential cone condition is satisfied by these operators, such that both inverse problems are locally ill-posed as well as the inverse problem for their linearizations. To show some numerical examples, we use the inexact Newton method CG-REGINN ([Rie05]), to reconstruct the perturbation from artificially generated noisy data.

The Bloch-Floquet transform is a well-known approach in electrical engineering, which is called the array scanning method, see, e.g., [MB79], [Val+08]. Nevertheless, the consideration of applying the transform to scattering problems was given just recently by constructing a numerical scheme and analyzing error bounds for the acoustic and electromagnetic scattering problem in the case of sound-soft boundary conditions (see [LZ17], [Zha18], [LZ17a]). Moreover, in [HN17] the acoustic scattering problem for an inhomogeneous layer was studied by applying the Bloch-Floquet transform and considering integral equations. The setting of the direct problem is close to the one in this paper, but with the somewhat easier assumption of a wave number with a positive imaginary part.

The remainder of this paper is structured as follows. In Section 2.1 we consider the direct problem, for which we present the setting of the scattering problem corresponding to the locally perturbed periodic layer. We use the Bloch-Floquet transform to show unique existence of the solution for the unperturbed case in Section 2.2 and consider the perturbed layer problem in Section 2.3. In Section 3 we analyze the inverse problem by defining a suitable parameter space, defining the parameter-to-state map, calculating the Fréchet derivative and show the ill-posedness as well as the uniqueness for the inverse problems. In the last two sections, we introduce the numerical method for the direct and inverse problem in Section 4 and show some numerical examples in Section 5.

2 Direct Scattering Problem

In this section we formulate the scattering problem for a perturbed periodic layer and prove unique existence of the scattered field. For that, we use the Bloch-Floquet transform to reduce the problem to a family of quasi-periodic problems on a bounded domain.

2.1 Formulation of the problem

Suppose $n_{p}^{2}\in L^{\infty}(\mathbb{R}_{+}^{d})$, $d=2,3$, is a $L$-periodic refractive index in $\underline{x}:=(x_{1},\ldots,x_{d-1})$, which satisfies $n_{p}^{2}=1$ for $x_{d}>R_{0}>0$ and characterizes the unperturbed scattering layer. To simplify the notation, we assume that $L$ equals to the scaled identity matrix $2\pi I_{d-1}\in\mathbb{R}^{(d-1)^{2}}$ and that the local perturbation $q\in L^{\infty}(\mathbb{R}_{+}^{d})$ has the support in $\Omega_{0}^{R}:=(-\pi,\pi)^{d-1}\times(0,R)$ for $R>R_{0}$, such that we consider the perturbed refractive index $n^{2}:=n^{2}_{p}+q$. Define for $R\geq 0$ the sets

[TABLE]

The scattering problem is to find the scattered field $u\in H_{0,\mathrm{loc}}^{1}(\mathbb{R}_{+}^{d})\cap H^{1}(\Omega^{R})$ for every $R>R_{0}$, such that

[TABLE]

Moreover, the scattering field is assumed to satisfy the so-called angular spectrum representation

[TABLE]

where $\widehat{u}$ is the Fourier transform of u\big{|}_{\Gamma^{R}} and the square root is extend by a branch cut at the negative imaginary axis. As a consequence, we can define the exterior Dirichlet-to-Neumann map $T$ as

[TABLE]

which is a bounded linear operator from $H^{\nicefrac{{1}}{{2}}}(\Gamma^{R})$ to $H^{\nicefrac{{-1}}{{2}}}(\Gamma^{R})$.

The analysis is easily extendable to the setting of free space scattering problem, assuming that the scattered field satisfies the angular spectrum representation in both directions. From now on, we call the space of $H^{1}(\Omega^{R})$-functions with vanishing trace on $\Gamma^{0}$ as $\widetilde{H}^{1}(\Omega^{R})$ and we consider an arbitrary function $f\in L^{2}(\Omega^{R})$, thus, the variational formulation is to

Problem 1**.**

Find a function $u\in\widetilde{H}^{1}(\Omega^{R})$, such that

[TABLE]

for all $v\in\widetilde{H}^{1}(\Omega^{R})$, where $n^{2}=n_{p}^{2}+q\in L^{\infty}(\Omega^{R})$.

Since for real wave numbers $k$ and for a real refractive index some surface waves can exist, we assume a small area of absorption.

Assumption 1**.**

The set $\{\mathrm{Im}~{}n_{p}^{2}>0\}$ is not empty and contains an open subset. Moreover, it holds $\mathrm{Im}~{}n_{p}^{2}\geq 0$ and $\mathrm{Im}~{}q\geq 0$.

The main result for this section is to prove unique existence of the scattered field.

Theorem 1**.**

If the 1 holds, the variational problem 1 has a unique solution.

To prove the theorem, we consider the quasi-periodic problem first.

2.2 Quasi-periodic inhomogeneous layer scattering

In this subsection we will be concerned with the quasi-periodic scattering problem and show the equivalence of the variational problem 1 to a family of quasi-periodic problems applying the Bloch-Floquet transform. For that, we treat the case that there is no perturbation at first, that means that $q=0$ and $n^{2}=n^{2}_{p}$. A function is called $\alpha$-quasi-periodic with $\alpha\in\mathbb{R}^{d-1}$ and period $2\pi$, if

[TABLE]

For smooth functions $\phi\in C^{\infty}_{0}(\overline{\Omega^{R}})$, the horizontal Bloch-Floquet transform $\mathcal{J}_{\mathbb{R}^{d-1}}$ is defined by

[TABLE]

Recall the spaces $H^{s}_{\alpha}(\Omega_{0}^{R})$ and $H^{s}_{\alpha}(\Gamma_{0}^{R})$ of $\alpha$-quasi-periodic Sobolev functions, and set $\widetilde{H}^{s}_{\alpha}(\Omega_{0}^{R})$ as the subspace of functions $u\in H^{s}_{\alpha}(\Omega_{0}^{R})$, such that u\big{|}_{\Gamma_{0}^{0}}=0. The Bloch-Floquet transform extends for $s\in\mathbb{R}$ to an isomorphism between $\widetilde{H}^{s}(\Omega^{R})$ and $L^{2}(I;\widetilde{H}^{s}_{\alpha}(\Omega_{0}^{R}))$ as well as between ${H}^{s}(\Gamma^{R})$ and $L^{2}(I;{H}^{s}_{\alpha}(\Gamma_{0}^{R}))$, where the index $\alpha$ indicates that the space depends on $\alpha\in I$ (see [Lec17]). The inverse of the transform is given by

[TABLE]

The scattered field $u_{\alpha}\in\widetilde{H}^{s}_{\alpha}(\Omega_{0}^{R})$ of the quasi-periodic scattering problem should satisfy the Rayleigh radiation condition

[TABLE]

where \widehat{(u_{\alpha}\big{|}_{\Gamma_{0}^{R}})}_{j} is the $j$-th Fourier coefficient of the trace. For $\phi\in H^{s}_{\alpha}(\Gamma_{0}^{R})$, $s\in\mathbb{R}$, $j\in\mathbb{Z}^{d-1}$ and $\alpha\in I$, the $j$-th Fourier coefficient $\widehat{\phi}_{j}(\alpha)$ of $\phi$ is defined by

[TABLE]

From the radiation condition, we derive the bounded quasi-periodic Dirichlet-to-Neumann operator $T_{\alpha}:H_{\alpha}^{\nicefrac{{1}}{{2}}}(\Gamma_{0}^{R})\to H_{\alpha}^{\nicefrac{{-1}}{{2}}}(\Gamma_{0}^{R})$ for $\phi=\sum_{j\in\mathbb{Z}^{d-1}}\widehat{\phi}_{j}(\alpha)e^{-\mathrm{i}\alpha\cdot\underline{x}}\psi^{j}(\underline{x})$ by

[TABLE]

Theorem 2**.**

Set $\mathcal{J}:=\mathcal{J}_{\mathbb{R}^{d-1}}$. A function $u\in\widetilde{H}^{1}(\Omega^{R})$ solves 1 for $q=0$ if and only if $\mathcal{J}u\in L^{2}(I;\widetilde{H}^{1}_{\alpha}(\Omega_{0}^{R}))$ solves the transformed variational problem

[TABLE]

for all $v\in L^{2}(I;\widetilde{H}^{1}_{\alpha}(\Omega_{0}^{R}))$. Furthermore, the radiation conditions (1) and (4) are equivalent for the corresponding problem.

Proof.

Set additionally $u_{\alpha}:=(\mathcal{J}_{\mathbb{R}^{d-1}}u)(\alpha,\cdot)\in\widetilde{H}^{1}_{\alpha}(\Omega_{0}^{R})$ for $u\in\widetilde{H}^{1}(\Omega^{R})$. From [Lec17] we know that the transform is an isomorphism between $\widetilde{H}^{s}(\Omega^{R})$ and $L^{2}(I;\widetilde{H}^{s}_{\alpha}(\Omega_{0}^{R}))$ for $s\in\mathbb{R}$, that the adjoint operator $\mathcal{J}^{*}$ can be identified with the inverse operator $\mathcal{J}^{-1}$, that one can interchange the transform with weak derivation and that the identity $\mathcal{J}(n^{2}w)=n^{2}\mathcal{J}w$ holds for every $w\in L^{2}(\Omega^{R})$. Applying these properties, we derive the equivalent sesquilinear form for the volume part as follows:

[TABLE]

The right hand side can be treated analogously. Now we have to show the equivalence on the boundary.

Calling $\gamma_{\Gamma^{R}}:\widetilde{H}^{1}(\Omega^{R})\to H^{\nicefrac{{1}}{{2}}}(\Gamma^{R})$ the trace operator on $\widetilde{H}^{1}(\Omega^{R})$ and $\gamma_{\Gamma_{0}^{R}}:\widetilde{H}_{\alpha}^{1}(\Omega_{0}^{R})\to{H}_{\alpha}^{\nicefrac{{1}}{{2}}}(\Gamma_{0}^{R})$ the trace operator on $\widetilde{H}^{1}_{\alpha}(\Omega_{0}^{R})$, we use the identification of the inverse Bloch-Floquet transform with its adjoint operator to get

[TABLE]

It holds the identity $\gamma_{\Gamma_{0}^{R}}\mathcal{J}u=\mathcal{J}\gamma_{\Gamma^{R}}u$, such that it remains to show that

[TABLE]

We define for every smooth function with compact support $\phi\in C_{0}^{\infty}(\mathbb{R}^{d-1})$ the operator

[TABLE]

which can be written as $\widetilde{J}=\mathcal{J}\circ\mathcal{F}^{-1}$, where $\mathcal{F}$ is the Fourier transform (see [Lec17]). This implies, in particular, that $\widetilde{J}$ is an isomorphism between the spaces $L^{2}_{s}(\Gamma^{R})$ and $L^{2}(I;H^{s}_{\alpha}(\Gamma_{0}^{R}))$ for $s\in\mathbb{R}$, where $L^{2}_{s}(\Gamma^{R})$ is the subspace of $L^{2}(\Gamma^{R})$ functions, for which the norm $||\xi\mapsto(1+|\xi|^{2})^{s/2}w(\xi)||_{L^{2}(\mathbb{R}^{d-1})}$ is finite. Putting the operator $\widetilde{J}$ into the definition of the Dirichlet-to-Neumann operator $T$, we conclude that

[TABLE]

Since it holds $(\xi\mapsto\mathrm{i}\sqrt{k^{2}-|\xi|^{2}}\mathcal{F}{(\gamma_{\Gamma^{R}}u)}(\xi))\in L^{2}_{s-1}(\mathbb{R}^{d-1})$ and $\mathcal{F}{(\gamma_{\Gamma^{R}}u)}(\alpha+j)=\widehat{(\gamma_{\Gamma_{0}^{R}}{u}_{\alpha})}_{j}(\alpha)$, we finally obtain the claimed identification.

For the radiation condition, one can use the same identity $\widehat{(\gamma_{\Gamma_{0}^{R}}{u}_{\alpha})}_{j}(\alpha)=\mathcal{F}(\gamma_{\Gamma^{R}}u)(\alpha+j)$ to directly calculate the equivalence of the radiation conditions. ∎

Theorem 3**.**

If the 1 holds, then the variational problem (6) is uniquely solvable.

We split the proof into three lemmas.

Lemma 4**.**

For all $\alpha\in\overline{I}$, there exists a unique solution ${w}_{\alpha}\in\widetilde{H}_{\alpha}^{1}(\Omega_{0}^{R})$ to the variational problem

[TABLE]

for every $\overline{v}\in\widetilde{H}_{\alpha}^{1}(\Omega_{0}^{R})$.

Proof.

Let $w_{\alpha}$ be in $\widetilde{H}_{\alpha}^{1}(\Omega_{0}^{R})$ for a fixed $\alpha\in\overline{I}$ and set \widehat{{w}}_{j}:=\widehat{({w}_{\alpha}\big{|}_{\Gamma_{0}^{R}})}_{j}(\alpha), then it holds for $\beta_{j}:=\sqrt{k^{2}-|\alpha+j|^{2}}$

[TABLE]

which implies

[TABLE]

Thus, the sesquilinear form fulfills the Gårding inequality. In the case of $||1-k^{2}n^{2}||_{L^{2}(\Omega_{0}^{R_{0}})}=0$, the problem is solvable by the theorem of Lax and Milgram. Because of the compact embedding of $\widetilde{H}_{\alpha}^{1}(\Omega_{0}^{R})$ into $L^{2}(\Omega_{0}^{R})$, the equation corresponds to a Fredholm operator of index zero. Consequently, by showing the injectivity, we obtain the unique existence of the solution.

For the boundary integral, it holds the inequality

[TABLE]

Since we assume $\mathrm{Im}~{}(n_{p}^{2})\geq 0$ in $\Omega_{0}^{R}$ and $\mathrm{Im}~{}(n_{p}^{2})>0$ on an open ball of $\Omega_{0}^{R}$, we derive for $(\mathcal{J}_{\mathbb{R}^{d-1}}f)(\alpha,\cdot)=0$

[TABLE]

We conclude that $w_{\alpha}$ vanishes on the open set, where $\mathrm{Im}~{}n_{p}^{2}>0$, and the theorem of unique continuation implies that $w_{\alpha}$ is equal to zero everywhere in $\Omega_{0}^{R}$. ∎

Using the same argumentation of the second part, we also get uniqueness for the integrated form (6).

Corollary 5**.**

Every solution to the variational problem (6) is unique.

Proof.

Because of the inequality (8), we have for $(\mathcal{J}_{\mathbb{R}^{d-1}}f)=0$ and the corresponding solution $w\in L^{2}(I;\widetilde{H}^{1}_{\alpha}(\Omega_{0}^{R}))$ the estimation

[TABLE]

where $w_{\alpha}:=w(\alpha,\cdot)$. This implies that $w_{\alpha}$ vanishes on an open ball for almost every $\alpha\in I$. Since $w_{\alpha}$ solves the Helmholtz equation almost everywhere in $I$, the unique continuation property implies that $w_{\alpha}$ vanishes everywhere w.r.t. to $x$ and almost everywhere in $I$. ∎

Now, we prove the connection between the pointwise variational problem and the integrated form.

Lemma 6**.**

The variational problem (6) is uniquely solvable.

Proof.

If we define the function $w$, $w(\alpha,\cdot):={u}_{\alpha}$, where $u_{\alpha}$ solve (4) for all $\alpha\in\overline{I}$, Lemma 4 implies that $w$ solves the problem (6). What still needs to be checked, is that $w$ lies in $L^{2}(I;\widetilde{H}^{1}_{\alpha}(\Omega_{0}^{R}))$. For that, we show that the solution operator $L_{\alpha}$ for the problem (4) is uniformly continuous in $\overline{I}$.

At first, we consider the continuity of the sesquilinear form (4). For every function ${v}_{\alpha}\in\widetilde{H}^{1}_{\alpha}(\Omega^{R})$, there exists a function ${v}^{p}_{\alpha}\in\widetilde{H}^{1}_{p}(\Omega_{0}^{R})$, where $\widetilde{H}^{1}_{p}(\Omega_{0}^{R})$ is the space with $\alpha=0$, such that ${v}_{\alpha}(\underline{x},x_{d})=e^{i\alpha\cdot\underline{x}}{v}^{p}_{\alpha}(\underline{x},x_{d})$. Moreover, the norms of the two functions are equal: $||{v}_{\alpha}||_{\widetilde{H}^{1}_{\alpha}(\Omega_{0}^{R})}=||{v}^{p}_{\alpha}||_{\widetilde{H}^{1}_{p}(\Omega_{0}^{R})}$. Now, we choose ${v}_{\alpha}$, ${u}_{\alpha}\in\widetilde{H}^{1}_{\alpha}(\Omega^{R})$ and ${v}^{p}_{\alpha}$, ${u}^{p}_{\alpha}\in\widetilde{H}^{1}_{p}(\Omega^{R})$ as described, and plugging them into the sesquilinear form (4) yields

[TABLE]

where

[TABLE]

In contrary to $T_{\alpha}$, the operator $S_{\alpha}$ only depends on $\alpha$ by the coefficients $\beta_{j}(\alpha):=\sqrt{k^{2}-|\alpha+j|^{2}}$.

Fix $\varepsilon>0$, $\alpha\in\mathbb{R}^{d-1}$ and $\alpha_{\varepsilon}\in\mathbb{R}^{d-1}$, where $|\alpha_{\varepsilon}-\alpha|<\varepsilon$, such that for every $u$, $v\in\widetilde{H}^{1}_{p}(\Omega_{0}^{R})$ it holds

[TABLE]

where

[TABLE]

For $j=0$ and for $j\in\mathbb{Z}^{d-1}$ with $k^{2}=|j+\alpha|^{2}$, the fraction $C(\alpha_{\varepsilon},\alpha)$ is continuous in $I$. For other $j\in\mathbb{Z}^{d-1}$, it holds

[TABLE]

For every $j\in\mathbb{Z}^{d-1}$ with $k^{2}\neq|j+\alpha|^{2}$, the value $\beta_{j}(\alpha)$ is contained either in $(-\infty,0)$, or $\mathrm{i}(0,\infty)$, and fulfills $|\beta_{j}(\alpha)|>\delta$ for a small constant $\delta>0$ independent of $j$. It follows

[TABLE]

for all $\alpha_{\varepsilon}\in\overline{I}$. Thus, it holds the estimation

[TABLE]

which implicates that the operator $\alpha\mapsto{b}_{\alpha}$ is continuous from $\overline{I}$ into $\mathcal{L}(\widetilde{H}^{1}_{\alpha}(\Omega_{0}^{R});\mathcal{L}(\widetilde{H}^{1}_{\alpha}(\Omega_{0}^{R});\mathbb{C}))$.

Since the sesquilinear form $a_{\alpha}$ is equivalent to $b_{\alpha}$, and since the norms of the spaces are equivalent, the sesquilinear form $a_{\alpha}$ is also continuous. Applying the Neumann series argument, we obtain that the solution operator $\alpha\mapsto L_{\alpha}$ is continuous on the compact set $\overline{I}$, and thus, bounded by a constant $C$ independent of $\alpha$. In particular, the function $w$ lies in ${L^{2}(I;\widetilde{H}^{1}_{\alpha}(\Omega_{0}^{R}))}$, since

[TABLE]

∎

2.3 Locally perturbed periodic inhomogeneous layer scattering

Combining Theorem 3 and Theorem 2, we obtain the unique existence of a solution for the unperturbed scattering problem. Now we consider the case that the perturbation $q\in L^{\infty}(\Omega_{0}^{R_{0}})$ is not vanishing. With the results from the subsection above, we are able to prove Theorem 1.

Proof of Theorem 1.

The sesquilinear form ${l}:\widetilde{H}^{1}(\Omega^{R})\times\widetilde{H}^{1}(\Omega^{R})\to\mathbb{C}$,

[TABLE]

is a compact perturbation, since $q$ vanishes outside of $\Omega_{0}^{R_{0}}$. As we showed earlier, the unperturbed problem is uniquely solvable, such that the variational formulation 1 corresponds to a Fredholm operator of index zero. Thus, we have to show uniqueness, which can be proven by using the same argumentation as in Lemma 4, if 1 holds, since for the solution $w=\mathcal{J}_{\mathbb{R}^{d-1}}u$ to $f=0$ it holds

[TABLE]

∎

As the last point, we show the regularity of the quasi-periodic solutions w.r.t. parameter $\alpha$. We will use this result for the implementation of the algorithm, since we can improve the convergence rate of the inverse transform with it. As the first step and defining the set

[TABLE]

one can show the regularity result for the unperturbed case by applying the Neumann series argument.

Theorem 7**.**

If the right hand side $\mathcal{J}_{\mathbb{R}^{d-1}}f$ is analytical in $\alpha\in\overline{I}$, then the map $\alpha\mapsto u_{\alpha}$, where $u_{\alpha}$ solves the quasiperiodic problem (4), is analytically in $\overline{I}\setminus\mathcal{A}$, and for any $\widehat{\alpha}\in\mathcal{A}$, there exists a $j_{0}\in\mathbb{Z}^{d-1}$ and a neighborhood $U(\widehat{\alpha})\subseteq\mathbb{R}^{d-1}$ of $\widehat{\alpha}$, such that the function can be decomposed into two analytical functions $u^{(1)}_{\alpha}$ and $u^{(2)}_{\alpha}$ in the form

[TABLE]

Proof.

This can be showed analogously to [Kir93a, Theorem a], which treats the case of the quasi-periodic scattering problem with sound-soft boundary conditions. Loosely speaking, one can split the differential operator $D_{\alpha}$ into $D^{(1)}_{\alpha}+\sqrt{k^{2}-|\alpha+j_{0}|^{2}}D^{(2)}_{\alpha}$, where both operators $D^{(1)}_{\alpha}$ and $D^{(2)}_{\alpha}$ are analytical in $\alpha$. Since $\sqrt{k^{2}-|\alpha+j_{0}|^{2}}\to 0$ for $|\alpha+j_{0}|^{2}\to k^{2}$, the Neumann series argument implies that the inverse of $D_{\alpha}$ can be decomposed in the same way. ∎

Since the compact perturbation of the sesquilinear form is independent of $\alpha$, one gets an analogous decomposition result to Theorem 7.

Theorem 8**.**

If the right hand side $\mathcal{J}_{\mathbb{R}^{d-1}}f$ is analytical in $\alpha\in\overline{I}$, the function $u_{\alpha}=\mathcal{J}_{\mathbb{R}^{d-1}}u(\alpha,\cdot)$, where $u\in\widetilde{H}^{1}(\Omega^{R})$ solves the (perturbed) variational problem 1, is analytically dependent on $\alpha\in\overline{I}\setminus\mathcal{A}$. For any $\widehat{\alpha}\in\mathcal{A}$, one can find a $j_{0}\in\mathbb{Z}^{d-1}$, a neighborhood $U(\widehat{\alpha})$ of $\widehat{\alpha}$, and two analytical functions $u^{(1)}_{\alpha}$ and $u^{(2)}_{\alpha}$, such that $u_{\alpha}$ can be written as

[TABLE]

Proof.

Let $K_{q}\in\mathcal{L}(L^{2}(I;\widetilde{H}^{1}_{\alpha}(\Omega_{0}^{R})))$ be the Riesz representation of $(w\mapsto k^{2}q\mathcal{J}_{\mathbb{R}^{d-1}}^{-1}w)\in\mathcal{L}(L^{2}(I;\widetilde{H}^{1}_{\alpha}(\Omega_{0}^{R})),L^{2}(I;\widetilde{H}^{1}_{\alpha}(\Omega_{0}^{R}))^{\prime})$,

[TABLE]

The operator $K_{q}$ maps functions from $L^{2}(I;\widetilde{H}^{1}_{\alpha}(\Omega_{0}^{R}))$ to functions, which are independent of $\alpha$, and thus, in particular, analytical in $\alpha$.

Let $w=\mathcal{J}_{\mathbb{R}^{d-1}}u\in L^{2}(I;\widetilde{H}^{1}_{\alpha}(\Omega_{0}^{R}))$ be the solution to the perturbed variational problem 1, and $A\in\mathcal{L}(L^{2}(I;\widetilde{H}^{1}_{\alpha}(\Omega_{0}^{R})))$ the Riesz representation of the unperturbed invertible differential operator for $q=0$. If we call the Riesz representation of the right hand side as $\widetilde{f}$, then it holds

[TABLE]

Since the right hand side $\widetilde{f}$ and the function ${K}_{q}w$ are analytical in $\alpha$, Theorem 7 implies that $w$ can be represented in the form of (9). ∎

In [Zha18] you can find comparable results for the sound-soft obstacle scattering problem and a detailed description, how to use the regularity to get a better convergence of the discretized inverse Bloch-Floquet transform.

Remark 9**.**

One can extend the regularity result in Theorem 8 easily for the case that the right hand site $f_{\alpha}$ can be decomposed in the same way as $f_{\alpha}=f_{\alpha}^{(1)}+\sqrt{k^{2}-|\alpha+j_{0}|^{2}}f_{\alpha}^{(2)}$, where $f_{\alpha}^{(1)}$ and $f_{\alpha}^{(2)}$ are analytical in $\alpha$.

3 The Inverse Problem

In this section we consider the inverse problem of reconstructing the perturbation. For that, we will consider the operator, which maps the perturbation to the solution operator for every right hand side $f\in L^{2}(\Omega^{R})$ with the support in one periodic cell $\Omega_{0}^{R_{0}}$.

At first, we will define the domain of definition for the measurement operators. For that, notice that in $d=2,3$ the space $\widetilde{H}^{1}(\Omega^{R}_{0})$ is continuously embedded in $L^{4}(\Omega_{0}^{R})$. Thus, for a $q_{0}\in L^{\infty}(\Omega^{R}_{0})$, $\mathrm{Im}~{}q_{0}\geq 0$, and $q\in L^{2}(\Omega_{0}^{R_{0}})$, such that $||q-q_{0}||_{L^{2}(\Omega_{0}^{R_{0}})}<\delta$ with $\delta>0$ small, for every $u$ and $v\in\widetilde{H}^{1}(\Omega^{R})$ it holds the estimation

[TABLE]

Consequently, for a small $\delta(q_{0})$, the sesquilinear form $a_{q}$ is a small perturbation of $a_{q_{0}}$, and the Neumann series argument guaranties the invertibility of the differential operator for perturbation $q$ of $n^{2}_{p}$. Since we need an open set as the domain of definition of the measurement operators, and the inversion methods for inverse problems depend on Hilbert spaces, we define the domain of definition $Q$ as

[TABLE]

where $B_{\delta(q_{0})}(q_{0})\subseteq L^{2}(\Omega_{0}^{R})$ is an open ball in $L^{2}(\Omega^{R}_{0})$ around a perturbation $q_{0}\in L^{\infty}(\Omega_{0}^{R})$ with the radius $\delta(q_{0})$ depending on $q_{0}$. Because of the Neumann series and the continuity of the sesquilinear form, the solution operator is well-defined for every $q\in Q$.

Definition 10**.**

Consider the linear and bounded operator $\Lambda_{q_{0}}:L^{2}(\Omega_{0}^{R_{0}})\to L^{2}(\Omega_{0}^{R})$, which maps a right hand side $f\in L^{2}(\Omega_{0}^{R_{0}})$ to the restriction of the solution $u_{q_{0}}\in\widetilde{H}^{1}(\Omega^{R})$ of 1 with $q=q_{0}$ to $\Omega_{0}^{R}$. We define the first measurement operator as

[TABLE]

mapping the perturbation $q\in Q$ to the operator $\Lambda_{q}$.

Definition 11**.**

Let $\widetilde{\Lambda}_{q_{0}}:L^{2}(\Omega^{R_{0}})\to\widetilde{H}^{1}(\Omega_{0}^{R})$ be the operator from above with codomain $\widetilde{H}^{1}(\Omega_{0}^{R})$ and let $\gamma_{\Gamma_{0}^{R}}:\widetilde{H}^{1}(\Omega^{R})\to H^{\nicefrac{{1}}{{2}}}(\Gamma_{0}^{R})$ be the trace operator, restricted to $\Gamma_{0}^{R}$. We define the second measurement operator

[TABLE]

which only measures the scattered field on one periodic cell of the upper boundary.

3.1 Uniqueness of the inverse problem

In this section, we will proof the injectivity of both operators $\Lambda$ and $\mathcal{S}$.

Theorem 12**.**

Consider two perturbations $q_{1}$ and $q_{2}\in Q$. Then it holds:

[TABLE]

Proof.

For a fixed right hand side $f\in L^{2}(\Omega_{0}^{R_{0}})$, we have two solutions $u_{1}$ for the variational problem 1 with $n^{2}=n^{2}_{p}+q_{1}$ and $u_{2}$ the solution to the 1 with $n^{2}=n^{2}_{p}+q_{2}$. Since $u_{1}$ equals $u_{2}$ on the set $\Omega_{0}^{R}$, the function $w:=u_{1}-u_{2}\in\widetilde{H}^{1}(\Omega^{R})$ solves the problem

[TABLE]

and vanishes, in especially, on $\Omega_{0}^{R}\setminus\Omega^{R_{0}}_{0}\neq\emptyset$. Applying the theorem of unique continuation, it follows that $u_{1}=u_{2}$ on $\Omega^{R}\setminus\Omega_{0}^{R}$, and consequently, the functions are identical on the whole domain $\Omega^{R}$. Thus, for every $v\in\widetilde{H}^{1}(\Omega^{R})$ it holds

[TABLE]

and the lemma of fundamental calculus implies $(q_{1}-q_{2})u_{2}=0$ almost everywhere. Since we can choose an arbitrary function $u_{2}\in C^{\infty}_{0}(\Omega_{0}^{R_{0}})$, we conclude the identity $q_{1}=q_{2}$ in $\mathbb{R}^{d}$. ∎

In the case of $d=3$, and additional regularity of the parameter $n_{p}^{2}$ and $q$, we can moreover prove injectivity of the operator $\mathcal{S}$, at least, if the whole trace on $\Gamma^{R}$ is given as data instead of data on one periodic cell. For that, we utilize the so called complex geometrical optics. The following proposition is adapted from [ILW16, Proposition 3.2] (see also [SU87]).

Proposition 13**.**

Let $D\subseteq\mathbb{R}^{3}$ be a bounded domain with Lipschitz boundary $\partial D$, $\xi\in\mathbb{C}^{3}$ satisfying $\xi\cdot\xi=0$ and $\rho\in H^{2}(D)$. Then there exist constants $C_{0}$ and $C_{1}$ depending on $D$, such that for $|\xi|>C_{0}||\rho||_{H^{2}(D)}$ there exists a solution $u$ of the form

[TABLE]

which solves the equation

[TABLE]

and satisfies

[TABLE]

Theorem 14**.**

Consider for $d=3$ two perturbations $q_{1}$ and $q_{2}\in Q\cap C^{2}(\mathbb{R}^{3}_{+})$ with compact support in $\mathbb{R}^{3}_{+}$, and assume $n_{p}^{2}\in C^{2}_{p}(\Omega_{0}^{R})$. If we call the solution operator $\widetilde{\Lambda}_{q}:L^{2}(\Omega_{0}^{R_{0}})\to\widetilde{H}^{1}(\Omega^{R})$ for $q\in Q$, and define $\widetilde{\mathcal{S}}:Q\to\mathcal{L}(L^{2}(\Omega^{R_{0}}),L^{2}(\Gamma^{R}))$, $q\mapsto\gamma_{\Gamma^{R}}\circ\widetilde{\Lambda}_{q}$, where $\gamma_{\Gamma^{R}}$ is the trace operator, then it holds:

[TABLE]

Proof.

For a fixed $f\in L^{2}(\Omega_{0}^{R_{0}})$, the operators $\widetilde{\mathcal{S}}(q_{1})$ and $\widetilde{\mathcal{S}}(q_{2})$ map the right hand side $f$ to the traces of the solutions $u_{1}$ and $u_{2}$ of 1 with $n^{2}=n^{2}_{p}+q_{1}$, or, $n^{2}=n^{2}_{p}+q_{2}$, respectively. The functions $u_{1}$ and $u_{2}$ coincide on $\Gamma^{R}$ and, of course, on $\Gamma^{0}$. The difference $w:=u_{1}-u_{2}$ satisfies the equation

[TABLE]

with homogeneous Dirichlet boundary conditions. Since the upper trace determines the extension by the radiation condition, we can conclude that $w$ vanishes on an open set for some $R^{\prime}>R>R_{0}>0$. The unique continuation theorem implies that the function $w$ vanishes on the biggest connected subset $D^{c}\subseteq\{q_{1}=q_{2}\}\subseteq\mathbb{R}^{3}$, which includes the boundary $\Gamma^{R}\cup\Gamma^{0}$. Hence $u_{1}=u_{2}$ on $D^{c}$. Putting the functions into the sesquilinear form, we obtain

[TABLE]

Now, we choose two arbitrary right hand sides $f^{\prime}$ and $g^{\prime}\in L^{2}(\Omega_{0}^{R_{0}})$ with support in $D^{c}$, and define, for $\rho(t,x):=t(n_{p}^{2}(x)+q_{2}(x))+(1-t)(n_{p}^{2}(x)+q_{1}(x))$, $\rho\in C^{\infty}([0,1];L^{2}(\Omega_{0}^{R_{0}}))$, the two solutions to 1 with $n^{2}=\rho(t,\cdot)$ as $u^{f^{\prime}}_{\rho(t,\cdot)}$, or, $v^{g^{\prime}}_{\rho(t,\cdot)}$, respectively. Then it holds

[TABLE]

Since $u^{f^{\prime}}_{\rho(0,\cdot)}=u^{f^{\prime}}_{\rho(1,\cdot)}$ and $v^{g^{\prime}}_{\rho(0,\cdot)}=v^{g^{\prime}}_{\rho(1,\cdot)}$ on $D^{c}$, as we showed earlier, and both function $f^{\prime}$ and $g^{\prime}$ are chosen to be zero on $D$, the complement of $D^{c}$, we have

[TABLE]

Set $u:=u^{f^{\prime}}_{\rho(t,\cdot)}$, $v:=v^{g^{\prime}}_{\rho(t,\cdot)}$, $u^{\prime}:=\frac{\partial}{\partial t}u$ and $v^{\prime}:=\frac{\partial}{\partial t}v$ for now, then, the upper equation implies

[TABLE]

Differentiating the variational problem 1 for $n^{2}=\rho$ w.r.t. $t$, we conclude that $u^{\prime}$ and $v^{\prime}$ solve the problem

[TABLE]

Since the derivative of $\rho$ w.r.t. $t$ is given by $\frac{\partial}{\partial t}\rho=(q_{2}-q_{1})$, we obtain

[TABLE]

for every $f^{\prime}$ and $g^{\prime}\in L^{2}(\Omega_{0}^{R_{0}})$ with f^{\prime}\big{|}_{D}=g^{\prime}\big{|}_{D}=0.

The assumption that $d=3$, allows us to choose vectors $\xi^{i}\in\mathbb{C}^{3}$, $i=1,2$ such that the norms $|\xi^{i}|^{2}$ are large for $i=1,2$, and both can be decomposed into

[TABLE]

with pairwise orthogonal real vectors $l$, $m$ and $p$, such that $|l|^{2}=|m|^{2}+|p|^{2}$. Choosing some Lipschitz domain $\widetilde{D}\supseteq D$, Proposition 13 gives us two functions $u$ and $v\in H^{2}(\widetilde{D})$ of the form (10). Multiplying a cut-off function $\chi\in C^{\infty}(\mathbb{R}^{3})$ to the functions $u$ and $v$, which fulfills \chi\big{|}_{D}=1 and \chi\big{|}_{\mathbb{R}^{3}\setminus\widetilde{D}}=0, one can see that these functions are solutions to 1 with suitable right hand sides $f^{\prime}$ and $g^{\prime}$ supported in $\Omega_{0}^{R_{0}}\setminus D$. Inserting these two functions into (11), we obtain

[TABLE]

Letting $|p|$ go to infinity, we deduce that the Fourier transform of the function $(q_{1}-q_{2})$ equals to zero. Consequently, the identity $q_{1}=q_{2}$ holds everywhere in $\mathbb{R}^{3}$. ∎

3.2 Fréchet differentiability and ill-posedness of the inverse problem

In the following, we will apply an inexact Newton-method, called CG-REGINN ([Rie05]), to reconstruct the shape of the perturbation. For that, we prove differentiability of the measurement operators $\Lambda$ and $\mathcal{S}$ as well as the ill-posedness of the inverse problems.

Theorem 15**.**

Fix $q\in Q$ and let $u_{f}$ be the solution of 1 for the right hand side $f\in L^{2}(\Omega_{0}^{R_{0}})$. Furthermore, for $h\in L^{2}(\Omega_{0}^{R_{0}})$ let $W_{h}:L^{2}(\Omega_{0}^{R_{0}})\to L^{2}(\Omega_{0}^{R})$, f\mapsto w_{h,f}\big{|}_{\Omega_{0}^{R}}, be the operator mapping $f$ to the solution $w_{h,f}$ for the 1 with the right hand side $k^{2}hu_{f}$ replacing $f$, i.e., the function $w_{h,f}$ solves

[TABLE]

for all $v\in\widetilde{H}^{1}(\Omega^{R})$.

Then, the derivative of $\Lambda$ is given by

[TABLE]

Proof.

Applying the Riesz theorem, we can reformulate the variational problem 1 as

[TABLE]

where $B_{q}\in\mathcal{L}(\widetilde{H}^{1}(\Omega_{0}^{R}))$ is the Riesz representation of the differential operator of 1 and $g$ the Riesz representation of $f$. One can check easily that the sesquilinear form is Fréchet differentiable w.r.t.  the perturbation $q$. It follows that the operator $q\mapsto B_{q}$ has a Fréchet derivative $q\mapsto B_{q}^{\prime}\in\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),\mathcal{L}(\widetilde{H}^{1}(\Omega_{0}^{R})))$.

Since the operator $B_{q}$ is invertible for every $q\in Q$, a corollary of the Neumann series argument implies that the operator $q\mapsto\Lambda_{q}=B_{q}^{-1}$ is also Fréchet differentiable and the linearization can be written as $\Lambda^{\prime}(q)[h]=-\Lambda_{q}B^{\prime}_{q}h\Lambda_{q}$, which corresponds to the claiming representation. ∎

As a consequence, we obtain the differentiability of the measurement operator $\mathcal{S}$.

Corollary 16**.**

The forward operator $\mathcal{S}$ is Fréchet differentiable in $q\in Q$. The derivative is given by $\mathcal{S}^{\prime}q\in\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),L^{2}(\Gamma_{0}^{R})))$, which maps a function $h\in L^{2}(\Omega_{0}^{R_{0}})$ to the operator $\gamma_{\Gamma_{0}^{R}}\circ\widetilde{W}_{h}:L^{2}(\Omega^{R_{0}})\to L^{2}(\Gamma_{0}^{R})$, where $\widetilde{W}_{h}:L^{2}(\Omega_{0}^{R_{0}})\to\widetilde{H}^{1}(\Omega_{0}^{R})$ does the same as $W_{h}$, just mapping to $\widetilde{H}^{1}(\Omega_{0}^{R})$.

In the rest of the section, we show that the measurement operator $\Lambda$, and its derivative $\Lambda^{\prime}$, yields an locally ill-posed inverse problem by proving that the operator satisfies the tangential cone condition. The local ill-posedness of the inverse problem for $\mathcal{S}$ and its derivative $\mathcal{S}^{\prime}$ can be showed analogously, which we will sketch afterwards.

For a general (non-linear) operator $\Phi:\mathcal{D}(\Phi)\subseteq X\to Y$ between Banach spaces $X$ and $Y$, the operator $\Phi$ is called locally ill-posed in $x^{*}\in\mathcal{D}(\Phi)$, if for all $r>0$, there exists a sequence $\{x_{n}\}_{n\in\mathbb{N}}\subseteq B_{r}(x^{*})\cap\mathcal{D}(\Phi)$, such that $||\Phi(x_{n})-\Phi(x^{*})||_{Y}\to 0$ but $||x_{n}-x^{*}||_{X}\not\to 0$ for $n\to\infty$ ([Sch+12, Definition 3.15]).

To prove ill-posedness, we will show that $\Lambda:Q\to\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),L^{4}(\Omega_{0}^{R}))$ is locally ill-posed, which implicates that also $\Lambda:Q\to\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),L^{2}(\Omega_{0}^{R}))$ is locally ill-posed. For that, we first show ill-posedness of the linear operator $\Lambda^{\prime}q_{0}:L^{2}(\Omega_{0}^{R_{0}})\to\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),L^{4}(\Omega_{0}^{R}))$ for a $q_{0}\in Q$, and conclude afterwards that the inverse problem for $\Lambda$ is locally ill-posed by proving the tangential cone condition.

Lemma 17**.**

The operator $\Lambda^{\prime}(q):L^{2}(\Omega_{0}^{R_{0}})\to\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),L^{4}(\Omega_{0}^{R}))$ is a compact operator for all $q\in Q$. In particular, the linearized operator equation is locally ill-posed in $\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),L^{4}(\Omega_{0}^{R}))$.

Proof.

Let $\{h_{j}\}_{j\in\mathbb{N}}\subseteq L^{2}(\Omega_{0}^{R_{0}})$ be a weakly convergent sequence, which means that for every functional $\psi\in L^{2}(\Omega_{0}^{R_{0}})^{\prime}$, it holds $\psi(h_{j})\to 0$ for $j\to\infty$. Thus, the right hand side $k^{2}h_{j}\Lambda_{q}f$ converges weakly in $H^{-1}(\Omega^{R})$ to zero for every $f\in L^{2}(\Omega_{0}^{R_{0}})$. The Sobolev space $H_{0}^{1}(\Omega_{0}^{R})$ is compactly embedded in $L^{4}(\Omega_{0}^{R})$, such that for every $f\in L^{2}(\Omega_{0}^{R_{0}})$ the sequence of solutions $\{W_{h_{j}}f\}$ converges to zero in $L^{4}(\Omega_{0}^{R})$ for $j\to\infty$. Applying the theorem of Banach-Steinhaus, we conclude that the sequence of operators $\{W_{h_{j}}\}$ converges to zero. Thus, $\Lambda^{\prime}(q)$ is a compact operator. ∎

Theorem 18**.**

The inverse problem related to the operator $\Lambda:Q\to\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),L^{2}(\Omega_{0}^{R}))$ is locally ill-posed.

Proof.

We show that the operator $\Lambda$ satisfies the tangential cone condition, which means that for some $q_{0}\in Q$ there exist a constant $0\leq\omega<1$ and an $r>0$, such that

[TABLE]

holds for all $q$ and $q^{*}\in B_{r}(q_{0})\cap Q$. Applying the triangle inequality, one deduces the relation

[TABLE]

for $q\neq q^{*}$. This, on the other hand, implies, together with [GL17, Theorem 4.5], that the local ill-posedness of the inverse problem $\Lambda(q)=\Lambda_{q}$ follows from the ill-posedness of the inverse problem for the Fréchet derivative $\Lambda^{\prime}(q)$, which we showed in Lemma 17.

Fix a right hand side $f\in L^{2}(\Omega_{0}^{R_{0}})$ and set $u_{q}\in\widetilde{H}^{1}(\Omega^{R})$ as well as $u_{q^{*}}\in\widetilde{H}^{1}(\Omega^{R})$ as the solutions to 1 for $n^{2}=n^{2}_{p}+q$, or, $n^{2}=n^{2}_{p}+q^{*}$, respectively. Moreover, let $w_{h}\in\widetilde{H}^{1}(\Omega^{R})$ be the solution to 1 for $n^{2}=n^{2}_{p}+q^{*}$ and right hand side $k^{2}hu_{q^{*}}$. If we define $w$ as $w:=u_{q}-u_{q^{*}}-w_{h}$, then it holds

[TABLE]

The function $w$ solves the variational problem

[TABLE]

for every $v\in\widetilde{H}^{1}(\Omega^{R})$. Consequently, it holds

[TABLE]

If the distance $||q-q^{*}||_{L^{2}(\Omega_{0}^{R_{0}})}$ is small enough, we can set $\omega:=Ck^{2}||q-q^{*}||_{{L^{2}(\Omega_{0}^{R_{0}})}}<1$, wherefrom the tangential cone condition follows, if we take the supremum on both sides:

[TABLE]

We conclude that the operator $\Lambda:Q\to\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),L^{4}(\Omega_{0}^{R}))$ is locally ill-posed, and thus, $\Lambda:Q\to\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),L^{2}(\Omega_{0}^{R}))$ is also locally ill-posed. ∎

The ill-posedness of the inverse problem related to the operator $\mathcal{S}$ can be shown analogously, which we will summarize in the next corollary.

Corollary 19**.**

The inverse problem related to the operator $\mathcal{S}:Q\to\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),L^{2}(\Gamma_{0}^{R}))$ is locally ill-posed.

Proof.

For $q\in Q$ and $h\in L^{2}(\Omega_{0}^{R_{0}})$ the definition space of the operator $\Lambda_{q}$ and its Fréchet derivative $W_{h}$ is $L^{2}(\Omega_{0}^{R_{0}})$, such that both operators map into $H^{2}(\Omega^{R}_{0})$. Thus, one can show analogously to Lemma 17 that the Fréchet derivative ${\Lambda}^{\prime}q$ is a compact operator mapping into $\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),H^{1}(\Omega_{0}^{R}))$, and further, one checks analogously to Theorem 18 that the tangential cone condition is satisfied for the image space $\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),H^{1}(\Omega_{0}^{R}))$. Consequently, the inverse problem related to the operator $\widetilde{\Lambda}:Q\to\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),H^{1}(\Omega_{0}^{R}))$ is ill-posed, which implicates that the inverse problem related to $\mathcal{S}$ is also ill-posed, since $\mathcal{S}q=\gamma_{\Gamma_{0}^{R}}\circ\widetilde{\Lambda}q$, where $\gamma_{\Gamma_{0}^{R}}$ is the trace operator. ∎

4 Numerical Solution Scheme and Reconstruction Method

In this section, we discuss the discretization of the unbounded locally perturbed variational problem 1, after applying the Bloch-Floquet transform to the variational formulation. To avoid having $\alpha$-dependent spaces $\widetilde{H}^{1}_{\alpha}(\Omega_{0}^{R})$, we will consider functions $w^{p}_{\alpha}\in\widetilde{H}^{1}_{p}(\Omega_{0}^{R})$, where $\widetilde{H}^{1}_{p}(\Omega_{0}^{R})$ is the space with $\alpha=0$, instead of $w_{\alpha}\in\widetilde{H}^{1}_{\alpha}(\Omega^{R})$, since they can be identified by $w_{\alpha}(\underline{x},x_{d})=e^{i\alpha\cdot\underline{x}}w^{p}_{\alpha}(\underline{x},x_{d})$. As the gradient $\nabla w_{\alpha}$ transforms to $(\nabla_{x}-\mathrm{i}\alpha)w^{p}_{\alpha}$, the $\alpha$-quasi-periodic variational problem (4) for $w_{\alpha}$ is equivalently reformulated for ${w}^{p}_{\alpha}$ as

[TABLE]

for every $v\in\widetilde{H}^{1}_{p}(\Omega_{0}^{R})$, where the Dirichlet-to-Neumann operator $T_{\alpha}$ is defined in the same way, since the Fourier coefficients do not change. We set

[TABLE]

for $w$, $v\in L^{2}(I;\widetilde{H}^{1}(\Omega_{0}^{R}))$, such that we can write the transformed (6) problem as

[TABLE]

Due to the perturbation, the sesquilinear form $b^{\prime}$ couples the $\alpha$-quasi-periodic components of the transformed solution.

4.1 Discretization of the scattering problem

In this section, we discretize the variational problem (13) as a family of problems, solved by finite elements method. Let $\mathcal{T}$ be the triangulation of $\overline{\Omega_{0}^{R}}=[-\pi,\pi]^{d-1}\times[0,R]$, consisting of $2^{d\times M}$ hypercubes that satisfy $\overline{\Omega_{0}^{R}}=\bigcup_{T\in\mathcal{T}}\overline{T}$, where $M\in\mathbb{N}$ stands for refinement cycles. Let ${\widetilde{M}}$ be the number of nodal points $\{x^{m}\}_{m=1}^{\widetilde{M}}\subset\overline{\Omega^{R}_{0}}$, which are equidistant in every direction, and $\{\phi_{\widetilde{M}}^{m}\}_{m=1}^{\widetilde{M}}$ the piecewise linear nodal functions, where $\phi_{\widetilde{M}}^{m}$ equals to one at the $m$-th nodal point $x^{m}$, and which equals to zero for other nodal points. Since the solution vanishes on the boundary $\Gamma_{0}^{0}$, we do not consider the nodal points there. Define the uniformly distributed grid points for $n=1,2,\ldots,N^{d-1}$ as $\alpha_{N}^{n}:=-\tfrac{1}{2}+\tfrac{1}{2N}+\tfrac{n}{N}$ in the case of $d=2$ and

[TABLE]

in the case of $d=3$ as well as the nodal basis of functions $\{\psi_{N}^{n}\}_{n=1}^{N^{d-1}}$, where $\psi_{N}^{n}$ equals to one on $I_{N}^{n}:=\alpha_{N}^{n}+[\nicefrac{{-1}}{{2N}},\nicefrac{{1}}{{2N}}]^{d-1}$ and zero, otherwise. The finite element space $V_{N,{\widetilde{M}}}$ is defined as

[TABLE]

and we seek for a finite element solution $\widetilde{w}\in V_{N,{\widetilde{M}}}$, which solves

[TABLE]

For a function $\widetilde{w}\in V_{N,{\widetilde{M}}}$, the inverse operator $\mathcal{J}_{\mathbb{R}^{d-1}}^{-1}$ of the Bloch-Floquet transform equals to the trapezoidal rule for integration, since

[TABLE]

For $n\in\{1,\ldots,N^{d-1}\}$, we approximate the value $\int_{I^{n}_{N}}a^{\prime}_{\alpha}(\phi_{\widetilde{M}}^{l},\phi_{\widetilde{M}}^{m})\;\mathrm{d}\alpha$ by $\tfrac{1}{N^{d-1}}a^{\prime}_{\alpha_{N}^{n}}(\phi_{\widetilde{M}}^{l},\phi_{\widetilde{M}}^{m})$ and $b^{\prime}(\widetilde{w},\psi_{N}^{n}\phi_{\widetilde{M}}^{m})$ by

[TABLE]

Thus, the discrete solution

[TABLE]

solves the linear system

[TABLE]

where for $n=1,\ldots,N^{d-1}$ and $m=1,\ldots,{\widetilde{M}}$ the discrete right hand side is defined by

[TABLE]

The matrix representation is given by

[TABLE]

where $F_{n}:=\left(F_{1,n},\ldots,F_{\widetilde{M},n}\right)$ and the matrices $A_{n}$, $B_{n}$ and $C_{n}$ are defined as

[TABLE]

The error analysis is out of scope of this paper, we refer the reader to [LZ17] and [LZ17a]. But we note that considering Theorem 8 we can actually improve the convergence rate of the discrete inverse of the Bloch-Floquet operator, if the right hand side is smooth enough. For that, one has to find a variable transform $g:I\to I$, such that the integrand of

[TABLE]

is a smooth and periodic function on $I$. It is well-known that the trapezoidal rule is converging very fast in the case of smooth periodic functions. For $d=2$, one can choose a function $g\in C^{\infty}(I;I)$, such that $g^{\prime}(\widehat{t})\geq 0$ and all of the derivatives $\frac{d}{dt^{m}}g(\widehat{t})$, $m\in\mathbb{N}$, vanish at $\widehat{t}$, where $g(\widehat{t})=\widehat{\alpha}$ holds. In this case, one gets convergence of order $\mathcal{O}(N^{-n}+2^{-2\times M})$ for some $n\in\mathbb{N}_{0}$ w.r.t. the $L^{2}(\Omega_{0}^{R})$-norm, if the right hand side is smooth enough (see [Zha18] for details).

The unperturbed system matrix $A$ is a block diagonal matrix consisting of the blocks $A_{i}$, $i=1,\ldots,N^{d-1}$. This emphasizes to invert the matrix block-wise using GMRES and the incomplete LU-decomposition for every block as the preconditioner. Furthermore, this allows the distribution of block-wise inversion tasks over a cluster of computers using Message Passing Interface (MPI). To utilize the special structure of A, we first solve the Schur complement for $U$

[TABLE]

and in a second step, we solve the equation $AW=F-BU$.

4.2 Regularization by inexact Newton method

In this section, we summarize the regularization scheme for the problem. Since the image spaces of the operators $\Lambda$ and $\mathcal{S}$ are not Hilbert spaces, we adjust these operators first. For that, we discretize $L^{2}(\Omega_{0}^{R_{0}})$ by the linear span of $N_{f}$ nodal functions $\{f_{R,m},f_{I,m}\}_{m=1}^{N_{f}}$, which are locally constant with the value of either zero or $1$ in the case of $f_{R,m}$, or, $\mathrm{i}$ in the case of $f_{I,m}$. All of the functions for the real part $\{f_{R,m}\}_{m=1}^{N_{f}}$ and all of the functions imaginary part $\{f_{I,m}\}_{m=1}^{N_{f}}$ are chosen to have disjoint support, such that it holds $\sum_{m=1}^{N_{f}}f_{R,m}=1$ and $\sum_{m=1}^{N_{f}}f_{I,m}=\mathrm{i}$. We define a modified operator $\widetilde{\Lambda}:Q\to L^{2}(\Omega^{R}_{0})^{2\times N_{f}}$ of $\Lambda:Q\to\mathcal{L}(L^{2}(\Omega_{0}^{R_{0}}),L^{2}(\Omega^{R}_{0}))$, where a perturbation $q\in Q$ is mapped to the $2\times N_{f}$ solutions of the variational problem 1 for the corresponding right hand sides. Analogously, we define $\widetilde{\mathcal{S}}:Q\to L^{2}(\Gamma^{R}_{0})^{2\times N_{f}}$, which maps the perturbation to the traces of these solutions. Since both operators $\widetilde{\Lambda}$ and $\widetilde{\mathcal{S}}$ map between Hilbert spaces, we can apply the regularization method CG-REGINN.

Assume now that

[TABLE]

We briefly summarize the regularization scheme CG-REGINN (“REGularization based on INexact Newton iteration”) stated and analyzed by Rieder in [Rie05], which we propose for the inversion. We will only consider the first inverse problem in (15) for the summary.

We have given the noisy version $U^{\varepsilon}$ of the exact measurement $U^{+}$ with the relative noise level $\varepsilon\in(0,1)$, i.e., $||U^{\varepsilon}-U^{+}||\leq\varepsilon||U^{+}||\approx\varepsilon||U^{\varepsilon}||$, which we assume to know a-priori. The algorithm generates a sequence $\{q_{m}\}_{m\in\mathbb{N}_{0}}\subseteq Q$ of approximations of $q^{+}$, starting with the initial guess $q_{0}\in Q$. If we write $q^{+}=q_{m}+s_{m}^{+}$ for each $m\in\mathbb{N}_{0}$, the best update $s_{m}^{+}$ solves the linearized problem

[TABLE]

where $E(q^{+},q_{m})$ is the linearization error. Since we do not know the linearization error, we only know the perturbed right hand side $b_{m}^{\varepsilon}:=U^{\varepsilon}-\widetilde{\Lambda}(q_{m})$ with the upper bound $||b_{m}^{\varepsilon}-b_{m}^{+}||\leq\varepsilon||U^{\varepsilon}||+\mathcal{O}(||q^{+}-q_{m}||)$ for the noise level.

CG-REGINN applies the regularization method of conjugate gradients (CG) for the linearized problem and stops, when the relative linear residuum is smaller than a tolerance times the non-linear residuum. CG creates an inner iteration that computes a sequence of approximations $\{s_{m,i}\}_{i\in\mathbb{N}_{0}}$ of $s^{+}_{m}$. The inner loop is terminated, when $||\widetilde{\Lambda}^{\prime}(q_{m})[s_{m,i}]-b_{m}^{\varepsilon}||<\mu_{m}||b_{m}^{\varepsilon}||$ for a tolerance $\mu_{m}\in(0,1)$ is satisfied for the first time, which index we call $i_{m}$. Then, we use backtracking, to get $\widetilde{s}_{m,i_{m}}:=\beta\widetilde{s}_{m,i_{m}}+(1-\beta)\widetilde{s}_{m,i_{m}-1}$, where $\beta\in[0,1]$ is chosen, such that $||\widetilde{S}^{\prime}(q_{m})[\widetilde{s}_{m,i}]-b_{m}^{\varepsilon}||=\mu_{m}||b_{m}^{\varepsilon}||$. We define the update as $s_{m}:=\widetilde{s}_{m,i_{m}}$ and set $q_{m+1}=q_{m}+s_{m}$, until the discrepancy principle with $\tau>0$ is satisfied for the outer loop.

Considering the suggestion in [Rie05], we chose $\mu_{1}=\mu_{2}=\mu_{start}=0.55$, $\gamma=0.9$, $\mu_{max}=0.99$, and

[TABLE]

where

[TABLE]

Taking the results in [EH18] into account, we use the adjoint matrix of the discretized problem for the inner loop of the numerical reconstruction, instead of the discretization of the theoretical adjoint of the Fréchet derivative, to have a more stable inversion.

5 Numerical Examples

In this section, we present some numerical results for the Bloch transform based method and the two inverse problems. We note at this point that the deal.II library does not support complex numbers, such that the values of the functions are considered as elements of $\mathbb{R}^{2}$. In this case, we get a system of two partial differential equations with some couplings, which double the number of degrees of freedom.

5.1 Example for the Bloch transform based method

For the first and second example, we choose $d=2$, $R=5$, $k=\sqrt{0.4}$, the cut-off of the Fourier expansion of the boundary for $|j|\leq 300$, and

[TABLE]

as well as

[TABLE]

as the reference solutions. The Bloch-Floquet transformed function of the second solution $u_{2}$ can be approximated by

[TABLE]

where $\operatorname{sinc}(t):=\nicefrac{{\operatorname{sin}(t)}}{{t}}$ is a smooth function. The transformed function of $u_{1}$, we simply approximate by

[TABLE]

since this function is decaying fast. Because of the extra factor $\nicefrac{{x_{2}}}{{5}}$, both do not satisfy the Neumann boundary condition, such that we add some correction factors $r_{i}$, $i=1,2$,

[TABLE]

where $r_{2}$ can be simplified to

[TABLE]

For the unperturbed refractive index, we take the function

[TABLE]

and the perturbation $k^{2}q$ is given by

[TABLE]

In Figure 2 both parameter are visualized. We set $2^{2\times M}$ as the number of cubic cells the domain is discretized in and $N$ the number of points for the discretization of the interval $I$. Note that $2^{2\times M}=$65,536\text{,}$$ corresponds to $132\,098\text{\,}$ degrees of freedom, since we have a system of two partial differential equations. The relative tolerance for GMRES is chosen to be $10^{-10}$. In Table 1 and Table 3 one can see the relative $L^{2}(\Omega^{R}_{0})$-errors for the two examples and in Table 2 the computation time for Example 1 using three computers (Intel i7-4790, $8\times 3.6$GHz cores, $32$GB memory) in parallel. In both cases we use the variable transformation $g:I\to I$, which is defined as

[TABLE]

where $\widehat{k}:=|k-\lfloor k+0.5\rfloor|$,

[TABLE]

If we use the identity as variable transform in the case of Example 1, then the error would decrease faster for smaller $N$, since the function $\mathcal{J}_{\mathbb{R}^{d-1}}u_{1}(\cdot,x)$ is $C_{p}^{\infty}(I)$ for every $x\in\Omega^{R}_{0}$. For $N=8$ we would already see near as good error values as for $N=256$ in Table 1. But in the case of the second example, the variable transform lets the error decrease much faster w.r.t. $\alpha$, since the second part of the decomposition of $\mathcal{J}_{\mathbb{R}^{d-1}}u_{2}$ shown in Theorem 8 does not vanish and the function has the square-root-line behavior.

To show some three dimensional examples, we choose $d=3$, $R=5$, $k=\sqrt{0.4}$ and the Fourier expansion cut-off $|j|\leq 10$. For the reference solutions, we choose

[TABLE]

and

[TABLE]

We approximate the Bloch-Floquet transformed functions $\mathcal{J}_{\mathbb{R}^{2}}u_{3}$ and $\mathcal{J}_{\mathbb{R}^{2}}u_{4}$ by

[TABLE]

or,

[TABLE]

respectively, and we add a correction term for the Neumann boundary condition. For the unperturbed refractive index, we choose the function

[TABLE]

and the perturbation $k^{2}q$ is given by

[TABLE]

Both parameter are visualized in Figure 4. The relative tolerance of GMRES is still chosen as $10^{-10}$, and we took the identity for the variable transform $g$ in both cases.

5.2 Examples for the reconstruction of the perturbation

In this subsection, we give the results, if we reconstruct the perturbation of the periodic refractive index, both shown in Figure 2. To generate the data, we use the algorithm from above, and refine some of the parameter, such that, we have $16\,384\text{\,}$ cells for $\Omega^{R}_{0}$, $128$ subintervals of $I$ and with the cut-off of the Fourier expansion of the boundary is $|j|\leq 600$. After that, we interpolate the solution down to $4096\text{\,}$ cells for $\Omega^{R}_{0}$, put some unified distributed noise of $5\%$ on it, and use this as the given data. For the reconstruction, we use $4096\text{\,}$ cells for $\Omega^{R}_{0}$, $N=64$, and a cut-off of $300$. For the right hand sides, we choose $R_{0}=4.5$, and split the domain $\Omega_{0}^{R_{0}}$ into $16$ equal parts. We approximate the $L^{2}(\Omega_{0}^{R_{0}})$ space with $32$ local constant functions $\{f_{R,m},f_{I,m}\}_{m=1}^{16}$, which are locally constant on the every part with the value of either zero or $1$ in the case of $f_{R,m}$, or, $\mathrm{i}$ in the other case, and such that it holds $\sum_{m=1}^{16}f_{R,m}=1$, $\sum_{m=1}^{16}f_{I,m}=\mathrm{i}$.

We approximate the perturbation $k^{2}q$ as a function in the finite element space, which is spanned by the finite elements $\{\phi_{\widetilde{M}}^{m}\}_{m=1}^{\widetilde{M}}$ in (14), and stop the outer iteration of REGINN by the discrepancy principle, when the relative discrepancy is smaller than $1.2\times 0.05$. In Figure 3 one can see the result of the reconstruction, where relative $L^{2}(\Omega_{0}^{{}^{R}})$ reconstruction error is about $38.0087\%$ in the case of $\widetilde{\Lambda}$, and a reconstruction error of $57.1368\%$ in the case of $\widetilde{\mathcal{S}}$. The results for the inversion of $\widetilde{\Lambda}$ are much better, since it has more data given to work with. Furthermore, the quality of the reconstruction depends highly on the size and the value of the absorption area. The bigger the set $\{\mathrm{Im}~{}n^{2}_{p}>0\}$ and the value inside is, the more the error of the reconstruction decreases.

For the three dimensional example, we use $4096$ cells for $\Omega^{R}_{0}$, $256$ cells of $I$ and a cut-off for the Rayleigh boundary condition of $|j|\leq 30$, $j\in\mathbb{Z}^{2}$. For the right hand side, we split the domain $\Omega_{0}^{R_{0}}$ into $8$ cubes, and for the data, we added $5\%$ of unified distributed noise. The relative reconstruction error in the case of $\widetilde{\Lambda}$ is about $51.3783\%$ (compare Figure 4).

Acknowledgement

The first author is very grateful for the devoted and generous support of Armin Lechleiter during his master’s and PhD program, who, although no longer with us, continues to inspire by his example and dedication to mathematics and teaching.

This project was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Projektnummer 281474342/GRK2224/1.