Domain decomposition for quasi-periodic scattering by layered media via robust boundary-integral equations at all frequencies
Carlos Pérez-Arancibia Stephen Shipman Catalin Turc Stephanos Venakides
math.NA physics.comp-ph 65N38, 35J05, 65T40, 65F08

Efficient Web Service Composition via Knapsack-Variant Algorithm
    Shiliang Fan^1
Yubin Yang^1,*
    
================================================================


  We develop a non-overlapping domain decomposition method (DDM) for scalar wave scattering by periodic layered media.  Our approach relies on robust boundary-integral equation formulations of Robin-to-Robin (RtR) maps throughout the frequency spectrum, including cutoff (or Wood) frequencies.   We overcome the obstacle of non-convergent quasi-periodic Green functions at  these frequencies by incorporating newly introduced shifted Green functions. Using the latter in the definition of quasi-periodic boundary-integral operators leads to rigorously stable computations of RtR operators. We develop Nyström discretizations of the RtR maps that rely on trigonometric interpolation, singularity resolution, and fast convergent windowed quasi-periodic Green functions. We solve the tridiagonal DDM system via recursive Schur complements and establish rigorously that this procedure is always completed successfully. We present a variety of numerical results concerning Wood frequencies in two and three dimensions as well as large numbers of layers.

 
  Keywords: Helmholtz transmission problem, domain decomposition, periodic layered media, lattice sum.

   
 AMS subject classifications: 
 65N38, 35J05, 65T40,65F08







§ INTRODUCTION




Simulation of electromagnetic wave propagation in periodic layered media has numerous applications in optics and photonics (photovoltaic devices, computation of plasmons, etc.). The use of periodic structures, such as diffraction gratings, which transmit and reflect waves along a discrete set of propagating directions,  opens up interesting possibilities to guide and direct waves in  unusual ways. Volumetric discretizations (finite-difference (FD) <cit.>, finite element (FE) <cit.>), that constitute the vast majority of numerical methods, require very large numbers of unknowns to suppress their inherent pollution effect, and thus produce very large linear systems requiring good preconditioners, which may not be readily available.  Furthermore, such methods must enforce radiation conditions in infinite domains by means of absorbing boundary conditions (ABC) or perfectly matched layers (PML) (see, for example, <cit.>), both of which meet difficulties in the treatment of surface waves and evanescent modes <cit.>. 


In the technologically relevant case of piecewise constant periodic layered media,  simulation methods based on boundary-integral equations (BIE) and  quasi-periodic Green functions are attractive candidates.  Radiation conditions are enforced automatically, and  discretizations of material interfaces are much smaller than volumetric discretizations and do not suffer from the pollution effect. 
Quasi-periodic Green functions are infinite sums of free-space Green functions with periodically distributed monopole singularities.  These double sums converge, although very slowly, for all but a discrete set of “cutoff" frequencies, for a given quasi-periodicity parameters (Bloch wavevector).  These are cutoff frequencies at which a Rayleigh diffraction mode transitions between propagating and evanescent and the number of propagating directions jumps.  Around these frequencies, the energy is rapidly redistributed along emerging new directions and is associated with anomalous scattering behavior.  These frequencies are often referred to as Wood frequencies (or Wood configurations of wavevector and frequency) because their problematic association in the literature to Wood's anomaly; see the works <cit.>,<cit.> and references therein for discussions on this phenomenon.
Popular methods for accelerating the slow convergence at non-Wood frequencies include Ewald summation <cit.> and lattice sums <cit.>.  At very high frequencies, asymptotic methods help to accelerate computation; see for example <cit.>.

While the underlying scattering problems are, with regard to the PDE, generically stable at Wood configurations of wavevector and frequency, the latter pose a challenge to BIE for quasi-periodic problems.  In three dimensions, they become increasingly close together at high frequency, and this puts the solution of quasi-periodic problems based on the quasi-periodic Green function out of reach.  For periodic layered media with large numbers of layers, such as thin films used in photovoltaic cells, the probability of encountering Wood frequencies is high.  Another difficulty is the need for an efficient algorithm for the evaluation of quasi-periodic Green functions and their integration into existing fast BIE solvers.  In the solution of the ensuing dense linear systems, the BIE formulations of periodic layered media give rise to tridiagonal solvers, whose structure can be exploited  to lead to efficient direct solvers <cit.>.


The above challenges faced by BIE-based quasi-periodic solvers were addressed
in two  recent computational methods. Alternative periodization schemes for boundary-integral formulations of quasi-periodic problems that do not rely on the classical quasi-periodic Green functions were proposed in <cit.>. These methods have ideas in common with the work presented in <cit.> as well as kernel-idependent FMM methods <cit.> and rely on representations of fields as sums of layer potentials and linear combinations of free-space fundamental solutions (radial basis functions) <cit.> whereby the quasi-periodicity and radiations conditions are enforced numerically and are not intrinsically satisfied.   This approach gives rise to efficient direct solvers for transmission problems in two-dimensional periodic layered media, and can yield results even at Wood frequencies <cit.>.  Its rigorous analysis appears to be absent in the literature, 
and we are not aware of evidence that these methods are capable of handling Wood frequencies in three dimensions.  

For the first time in this arena, a rigorous solution was provided to the problem of boundary-integral equation formulations of quasi-periodic problems at Wood frequencies in both two and three dimensions through a new method, in which the well-posedness of the formulation and the stability of the numerical scheme were proven, each in its own right <cit.>.   Smooth windowed truncations of the lattice sums for the Green functions were introduced and analyzed in <cit.>. It was first shown in the same reference <cit.> that the windowed Green functions (WGF) converge to their corresponding quasi-periodic Green functions superalgebraically away from Wood frequency/wavevector configurations as the radius of truncation increases.  The incorporation of WGF in existing fast boundary-integral solvers is relatively straightforward. Remarkably, the WGF method can be adapted to handle scattering problems in layered media whose infinite interfaces are no longer periodic <cit.>. Then, shifted Green functions that converge at and around Wood anomalies were used in a boundary-integral equation setting to provide accurate solutions of scattering problems for perfectly reflecting gratings throughout the frequency spectrum. The shifted Green functions converge algebraically fast at Wood frequencies, and the rate of convergence grows with the number of shifts. However, the shifted Green functions introduce new singularities (poles) in addition to those already present in the quasi-periodic Green functions. Remarkably, these additional singularities turn out to be benign in the case of perfectly reflecting periodic gratings as they can be arranged to be outside the computational domain if one uses indirect formulations <cit.>.


In this article, we extend the shifted Green function method to the case of scalar transmission problems in periodic layered media.  By using a domain decomposition method (DDM), we overcome the difficulty of poles of the shifted Green function inside the computational domain, and we establish the well-posedness of the ensuing system of boundary-integral equations.  There is a vast literature on DDM; the reader is referred to the seminal works of B. Després <cit.> and the expository books <cit.>.
DDM is well suited to the Helmholtz/Maxwell equations in periodic layered media because of the robustness of the Robin-to-Robin (RtR) operator for each layer <cit.>.  For a given periodic layer, it maps incoming (interior boundary) Robin data to outgoing (exterior boundary) Robin data on the interfaces that bound that layer.  In this way, Robin data are matched on each interface of material discontinuity.  This procedure produces a tridiagonal system whose unknowns are the Robin data on interfaces, and whose non-zero blocks consist of RtR operators.  If a particular layer has constant material properties, the RtR operators can be computed robustly in terms of boundary-integral operators that use  the ordinary quasi-periodic Green function for frequencies that are not Wood frequencies, and shifted Green functions for wavenumbers that are near or at Wood frequencies. Interestingly, the computations of RtR do not require use of hypersingular boundary-integral operators. We establish rigorously in this work two important facts.





  * The computations of RtR maps via boundary-integral operators are robust throughout the frequency spectrum if shifted Green functions are employed at Wood frequencies.
  
  * The DDM for solution of scalar transmission problems in periodic layered media with piecewise constant material properties presented in this paper is equivalent to the original PDE, assuming that the PDE problem is well-posed.
  

We develop a high-order discretization of the tridiagonal DDM system based  on Nyström discretizations of periodic boundary-integral operators. The latter, in turn, rely on trigonometric interpolation, logarithmic singularity extraction in two dimensions and analytic resolution of singularity in three dimensions, and the windowed Green function method <cit.>. We solve the DDM system using recursive Schur complements to eliminate sequentially the discretized Robin data corresponding to each layer in a top-down sweep, a procedure that leads to a computational cost that is linear in the number of layers. We also present theoretical arguments to explain why the Schur complement elimination procedure can be always completed successfully. The variety of two- and three-dimensional numerical results presented in this paper showcase the capability of our DDM solver to handle large numbers of layers, challenging Wood configurations, and inclusions in a periodic layered medium.  The DDM solvers presented in this paper, being built on quasi-periodic Green functions, must be re-assembled when the quasi-periodic parameter changes. Also, the computations of RtR operators require inversions of boundary-integral operators. In summary, the DDM solvers developed in this paper enjoy the following attractive features.


  * The computations of the RtR maps are stable across the frequency spectrum.

  * The DDM system can be solved via recursive Schur complements, leading to a computational cost and memory usage that are linear in the number of layers; and it can be shown rigorously that this procedure does not break down.

  * The DDM approach, being modular, allows for use of heterogeneous discretizations such as FE and BIE and use of non-conforming discretizations on interfaces pertaining to layers with different material properties.
  
  * DDM are easily parallelizable.
  
 The integration of WGF and shifted Green functions  in existing three-dimensional boundary-integral operator discretizations presented in this contribution is relatively seamless and results in a rigorous treatment of Wood configurations in three dimensions.
 
DDM approaches will be feasible for the solution of three-dimensional electromagnetic transmission problems in periodic layered media based  on quasi-optimal transmission conditions <cit.> that renders them amenable to Krylov subspace iterative solvers.  Quasi-optimal transmission conditions arise from a judicious choice of the complex wavenumber in the transmission operator that gives rise to a DDM whose rate of convergence is practically independent of frequency <cit.>.

The RtR DDM that employs a shifted Green-function scheme can handle interfaces between layers that are very general (including those that are not the graph of a function) and general frequencies.
There are of course situations in which other methods would be superior or should be used in combination with the RtR DDM.  In the case of small, smooth perturbations of flat interfaces, the method of variation of boundaries would provide increased acceleration <cit.>, even if the perturbations are not that small <cit.>.
And as noted above, at high frequencies, asymptotic methods should be used to accelerate the computation <cit.>.

The paper is organized as follows. In Section <ref> we present the scalar scattering problem in two-dimensional layered media and we review the main results about the well-posedness of these problems.  In Section <ref> we present a DDM formulation of the transmission problems that uses matching of classical Robin boundary conditions of the material interfaces, and we present computations of ensuing RtR maps that are shown to be stable throughout the frequency spectrum.  We continue in Section <ref> with a description of the Nyström discretization of the RtR maps and we provide and analyze a recursive Schur complement elimination algorithm for the direct solution of the discrete DDM system.  Finally, we present in Section <ref> a variety of numerical results of wave scattering at mostly Wood frequency configurations in periodic layered media. 



2pt plus2pt minus1pt



§ SCALAR TRANSMISSION PROBLEMS 

We consider the problem of quasi-periodic scattering by penetrable homogeneous periodic layers. For the sake of simpler notations, we present the two-dimensional case. We mention that all the derivations that we present are easily translatable to three-dimensional configurations. The periodicity of the layers is taken to be in the horizontal x_1 direction, that is the layers are given by Ω_j={(x_1,x_2)∈ℝ^2: F_j(x_1)≤ x_2≤ F_j-1(x_1)} for 0<j<N and Ω_0={(x_1,x_2)∈ℝ^2:F_0(x_1)≤ x_2} and Ω_N+1={(x_1,x_2)∈ℝ^2:x_2≤ F_N(x_1)}, and all the functions F_j are periodic with principal period d, that is F_j(x_1+d)=F_j(x_1) for all 0≤ j≤ N. We assume that the medium occupying the layer Ω_j is homogeneous and its permitivity is ε_j; the wavenumber k_j in the layer Ω_j is given by k_j=ω√(ε_j). A plane wave u^ inc(𝐱)=exp(i(α x_1+iβ x_2)) where α^2+β^2=k_0^2 impinges on the layered structure. We seek α-quasi-periodic fields u_j (i.e. u_j(x_1+d,x_2)=e^iα du(x_1,x_2) for all (x_1,x_2)∈ℝ^2) that satisfy the following system of equations:

    [                         Δ u_j +k_j^2 u_j                                        =                                        0                                       in     Ω_j^per:={(x_1,x_2)∈Ω_j: 0≤ x_1≤ d},;                            u_j+δ_0u^ inc                                        =                                    u_j+1                                       on Γ_j={(x_1,x_2):0≤ x_1≤ d, x_2=F_j(x_1)},;             γ_j(∂_n_ju_j+δ_0∂_n_ju^ inc)                                        =                       -γ_j+1∂_n_j+1u_j+1                                       on                                     Γ_j, ]

where δ_0 is the Dirac distribution supported on Γ_0 and n_j denote the unit normals to the boundary ∂Ω_j pointing to the exterior of the subdomain Ω_j. Note that we assigned to the partial derivatives on a given interface the index of the domain on whose side the partial derivative is taken; thus, on the interface Γ_j we have n_j=-n_j+1. We also assume that u_0 and u_N in equations (<ref>) are radiative in Ω_0 and Ω_N+1 respectively. The latter requirement amounts to expressing the solutions u_0 and u_N+1 in terms of Rayleigh series

    u_0(x_1,x_2)=∑_r∈ℤ C_r^+e^iα_r x_1+iβ_0,r x_2,   x_2>maxF_0

and

    u_N+1(x_1,x_2)=∑_r∈ℤ C_r^-e^iα_r x_1-iβ_N+1,r x_2,   x_2<minF_N

in which α_r=α+2π/dr and β_0,r=(k_0^2-α_r^2)^1/2 and β_N+1,r^2=(k_N+1^2-α_r^2)^1/2, where the square is root chosen such that √(1)=1 with branch cut along the negative imaginary axis. We assume that the wavenumbers k_j and the quantities γ_j in the subdomains Ω_j are positive real numbers.

Wood frequencies are those values of k for which there exist indices r_0 such that α_r_0^2=k^2. The well posedness of the equations (<ref>) was established in <cit.> in the case of two domains Ω_0 and Ω_1 separated by the periodic interface Γ_0 with γ_0=γ_1=1 and any real wavenumbers k_0 and k_1, including Wood frequencies. The techniques presented in <cit.> are easily applicable to periodic configurations with arbitrary number of layers.  To the best of our knowledge, any attempt at establishing uniqueness of solutions of equations (<ref>) was based on the aforementioned techniques. However, certain  requirements <cit.> must be imposed on the material parameters (k_j,γ_j), 0≤ j≤ N in order to establish rigorously the uniqueness of solutions of equations (<ref>) using those techniques. For the sake of completeness, we provide in Appendix <ref> a proof of uniqueness of solutions of equations (<ref>) under the assumption of monotonicity of the wavenumber k_j,0≤ j≤ N+1 and γ_j=1,0≤ j≤ N+1. In general, for a fixed periodic layered configuration with material properties ε_j, the transmission problem (<ref>) has a unique solution with the exception of a discrete set of frequencies ω whose only accumulation point is infinity <cit.>. The same comprehensive reference <cit.> contains a proof of existence of solutions for the transmission problem (<ref>) using both variational and boundary-integral equation arguments.





 

§ DOMAIN DECOMPOSITION APPROACH


 We present a domain decomposition method (DDM) based on boundary-integral equations (BIEs) for the numerical solution of transmission problems (<ref>). Just like BIE formulations, DDM formulations recast the original PDEs in terms of unknown quantities defined on the interfaces of material discontinuity. A non-overlapping domain decomposition approach for the solution of equations (<ref>) consists of solving Helmholtz subdomain problems in Ω_j,j=0,…,N+1 with matching Robin transmission boundary conditions on the common subdomain interfaces Γ_j for j=0,…,N. The main motivation for using DDM is the  seamless treatment of periodic configurations at Wood frequencies via BIE formulations, as well the ease with which it can handle inclusions in the periodic layers.  Specifically, DDM amount to computing α-quasi-periodic  subdomain solutions:

    Δ u_j +k_j^2 u_j    =   0     in  Ω_j^per,
    γ_0(∂_n_0u_0+∂_n_0u^ inc)-iη(u_0+u^ inc)   =   -γ_1∂_n_1u_1-iη u_1   on  Γ_0
    γ_1∂_n_1u_1-iη u_1   =   -γ_0(∂_n_0u_0+∂_n_0u^ inc)-iη(u_0+u^ inc)   on  Γ_0
    γ_j∂_n_ju_j-iη u_j   =   -γ_j+1∂_n_j+1u_j+1-iη u_j+1   on  Γ_j, 1≤ j≤ N
    γ_j+1∂_n_j+1u_j+1-iη u_j+1   =   -γ_j ∂_n_ju_j -iη u_j   on  Γ_j, 1≤ j≤ N.

In addition, we require that u_0 and u_N+1 be radiative and that η>0. The latter requirement ensures that the Robin problems in the semi-infinite domains Ω_0 and Ω_N+1 are well posed; see Theorem <ref>.

 The essence of the domain decomposition (<ref>)  is solving a Robin boundary-value problem in each layer subdomain and connecting the Robin boundary data across interfaces via the so-called Robin-to-Robin (RtR) maps <cit.>—also see below. For a given layer subdomain Ω_j with 1≤ j≤ N we seek w_j α-quasi-periodic solutions of the following Helmholtz boundary-value problem

    Δ w_j+k_j^2w_j   =   0   in Ω_j^per
    γ_j∂_n_jw_j-iη w_j   =   g_j-1,j   on Γ_j-1
    γ_j∂_n_jw_j-iη w_j   =   g_j,j   on Γ_j

where g_j-1,j and g_j,j are generic α-quasi-periodic functions defined on Γ_j-1 and Γ_j. The RtR map 𝒮^j is defined as
    𝒮^j[ g_j-1,j;   g_j,j ]=[ (γ_j∂_n_jw_j+iη w_j)|_Γ_j-1;   (γ_j∂_n_jw_j+iη w_j)|_Γ_j ].

The computation of the RtR maps 𝒮^j requires solving the Helmholtz boundary value problem (<ref>). Of the two indices of the boundary data g in equations (<ref>), the first index is the index of the interface and the second index is the index of the subdomain.  Thus, g_j-1,j refers to boundary data on the interface Γ_j-1 on the side of the subdomain Ω_j. The block structure of the RtR operators 𝒮^j defined in equation (<ref>) is 

    𝒮^j[ g_j-1,j;   g_j,j ]=[ 𝒮^j_j-1,j-1   𝒮^j_j-1,j;   𝒮^j_j,j-1     𝒮^j_j,j ][ g_j-1,j;   g_j,j ].

For the semi-infinite subdomain Ω_0 w_0 is the α-quasi-periodic outgoing solution of the Helmholtz boundary value problem

    Δ w_0+k_0^2w_0   =   0   in Ω_0^per
    γ_0∂_n_0w_0-iη w_0    =   g_0,0   on Γ_0 ,

in which g_0,0 is a α-quasi-periodic function defined on Γ_0, and we define the RtR map 𝒮^0 by

    𝒮^0g_0,0:=(γ_0∂_n_0w_0+iη w_0)|_Γ_0.

The RtR map 𝒮^N+1 corresponding to the semi-infinite subdomain Ω_N+1 is defined in a similar manner to 𝒮^0  but  for boundary data g_N-1,N defined on Γ_N. 

In DDM formulations (<ref>), the unknown Robin data associated with each interface Γ_j

    f_j=[   f_j,j; f_j,j+1 ]   :=   [         (γ_j∂_n_ju_j-iη u_j)|_Γ_j; (γ_j+1∂_n_j+1u_j+1-iη u_j+1)|_Γ_j ], 0≤ j≤ N

are matched via the subdomain RtR maps 𝒮^j,0≤ j≤ N+1 giving rise to a (2N+2)× (2N+2) operator linear system.  The unknown Robin data f=[f_0 f_1 … f_N]^⊤ are the solution of the following linear system

    𝒜f= ,

in which the DDM matrix 𝒜 is a tridiagonal block matrix whose first two rows, and the rows indexed by 2j+1 and 2j+2 (corresponding to the unknown Robin data f_j,j and f_j,j+1), and the last two rows, are given in explicit form

    𝒜=[           I     𝒮^1_0,0     𝒮^1_0,1           …           0           0           0           0           0           0           …           0           0;         𝒮^0           I           0           …           0           0           0           0           0           0           …           0           0;           …           …           …           …           …           …           …           …           …           …           …           …           …;           0           0           0           …           0           0           I   𝒮^j+1_j,j 𝒮^j+1_j,j+1           0           …           0           0;           0           0           0           …           0   𝒮^j_j,j-1     𝒮^j_j,j           I           0           0           …           0           0;           …           …           …           …           …           …           …           …           …           …           …           …           …;           …           …           …           …           …           …           …           …           …           …           …           I       𝒮^N+1;           …           …           …           …           …           …           …           …           …           …   𝒮^N_N,N-1     𝒮^N_N,N           I;             ]

 and in which the right-hand-side vector =[_0 _1 … _N]^⊤ has zero components _ℓ=[0 0]^⊤, 1≤ℓ≤ N, with the exception of the first component 
 
    _0=[ -(γ_0∂_n_0u^ inc-iη u^ inc)|_Γ_0; -(γ_0∂_n_0u^ inc+iη u^ inc)|_Γ_0 ].
 
 In what follows, we study spectral properties of the RtR operators 𝒮^j,0≤ j≤ N+1, properties that will shed light onto the solvability of the DDM system (<ref>).



 

 §.§ Spectral properties of the RtR operators 


 The first question that arises is whether the RtR operators are properly defined under the assumptions on wavenumbers k_j and coefficients γ_j>0, j=0,…,N+1. We establish the following result, whose proof is essentially a simple extension of arguments presented in <cit.>.
 
   Let w_0 be the α-quasi-periodic outgoing solution of the following Helmholtz equation 

    Δ w_0+k_0^2w_0   =   0   in Ω_0^per
    ∂_n_0w_0-iη γ_0^-1w_0    =   0   on Γ_0.

Then w_0 is identically zero in Ω_0^per.
 
  A similar uniqueness result holds for the homogeneous problem
 
    Δ w_N+1+k_N+1^2w_N+1   =   0   in Ω_N+1^per
    ∂_n_N+1w_N+1-iη γ_N+1^-1w_N+1   =   0   on Γ_N.

 
 
   Consider again h>maxF_0 and the domain Ω_0,h^per:={(x_1,x_2)∈Ω_0^per: F_0(x_1)≤ x_2≤ h}. A simple application of Green's identities leads to 

    ∫_Ω_0,h^per(|∇ w_0|^2-k_0^2|w_0|^2)dx   =   ∫_Γ_0∂_n_0w_0 w_0 ds + ∫_Γ_0,h∂_x_2w_0 w_0 dx_1
       =   iηγ_0^-1∫_Γ_0|w_0|^2 ds + ∫_Γ_0,h∂_x_2w_0 w_0 dx_1

where Γ_0,h:={(x_1,x_2): 0≤ x_1≤ d, x_2=h}. Taking into account the fact that w_0 is radiating, we can express w_0 on the line segment Γ_0,h in terms of the following Rayleigh series

    w_0(x_1,h)=∑_r∈ℤ C_r^+e^iα_r x_1+iβ_0,r h ,

from which it follows that

    ∫_Γ_0,h∂_x_2w_0 w_0 dx_1=id∑_r∈ℤ, β_0,r>0β_0,r|C_r^+|^2.

Consequently,

    ∫_Ω_0,h^per(|∇ w_0|^2-k_0^2|w_0|^2)dx=iηγ_0^-1∫_Γ_0|w_0|^2 ds+id∑_r∈ℤ, β_0,r>0β_0,r|C_r^+|^2.

The left-hand-side of this identity is real, whereas the right-hand side is a sum of non-negative imaginary terms, and thus each of these terms vanishes.  This implies that w_0=0 on Γ_0, and thus ∂_n_0w_0=0 on Γ_0 as well. The result now follows from Holmgren's uniqueness theorem <cit.>.
 

Consider now the following Helmholtz equation.   Let  w_0 be the α-quasi-periodic outgoing solution of

    Δ w_0+k_0^2w_0   =   0   in Ω_0^per
    ∂_n_0w_0-iη γ_0^-1w_0    =   g_0   on Γ_0

where g_0 is a α-quasi-periodic function defined on Γ_0. The matter of existence of such a solution will be settled in the next section through boundary-integral equation arguments. We are interested in estimating the norm of the RtR operator 𝒮^0 as a continuous operator from L^2_per(Γ_0) to itself.  We have

    g_0_2^2=∫_Γ_0(|∂_n_0w_0|^2+η^2γ_0^-2|w_0|^2)ds-2ηγ_0^-1∫_Γ_0∂_n_0w_0 w_0 ds

and

    𝒮^0g_0_2^2=∫_Γ_0(|∂_n_0w_0|^2+η^2γ_0^-2|w_0|^2)ds+2ηγ_0^-1∫_Γ_0∂_n_0w_0 w_0 ds.

Again, we have that

    ∫_Ω_0,h^per(|∇ w_0|^2-k_0^2|w_0|^2)dx=∫_Γ_0∂_n_0w_0 w_0 ds + ∫_Γ_0,h∂_x_2w_0 w_0 dx_1 .

Assuming the Rayleigh series expansion

    w_0(x_1,h)=∑_r∈ℤ C_r^+e^iα_r x_1+iβ_0,r h ,

we derive

    ∫_Γ_0∂_n_0w_0 w_0 ds=-d∑_r∈ℤ, β_0,r>0β_0,r|C_r^+|^2 ,

and hence

    𝒮^0g_0_2^2<g_0^2_2

for all g_0 is a α-quasi-periodic function defined on Γ_0.  It follows that 𝒮^0_L^2_per(Γ_0)→ L^2_per(Γ_0)≤ 1. Similar arguments lead to the estimate 𝒮^N+1_L^2_per(Γ_N+1)→ L^2_per(Γ_N+1)≤ 1. Green's identities establish the following theorem.


  The RtR operators 𝒮^j are unitary in the space L^2_per(Γ_j)× L^2_per(Γ_j+1) for all j : 1≤ j≤ N.


This unitarity can be used to establish the pointwise convergence of the Jacobi fixed-point iterations for the solution of the DDM formulation (<ref>) by a relatively straightforward adaptation of the arguments presented in <cit.> to the quasi-periodic setting.





 §.§ Calculations of RtR operators in terms of boundary-integral operators associated with quasi-periodic Green functions 

Implementation of DDM requires computation of RtR maps. We present in this section explicit representations of RtR maps in terms of boundary-integral operators associated with quasi-periodic Green functions that will serve as the basis of the implementation of the DDM algorithm. 


  §.§.§ Quasi-periodic Green functions, layer potentials and integral operators

 For a given free-space wavenumber (normalized frequency) k, define the α-quasi-periodic Green function

    G^q_k(x_,x_2)=∑_n∈ℤ e^-iα ndG_k(x_1+nd,x_2)

where G_k(x_1,x_2)=i/4H_0^(1)(k|𝐱|), 𝐱=(x_1,x_2).   Define α_r:=α+2π/dr and β_r=β_r(k):=(k^2-α_r^2)^1/2, where the branch of the square roots in the definition of β_r is chosen in such a way that √(1)=1, and that the branch cut coincides with the negative imaginary axis. It can be shown that the series in the definition of the Green function G^q_k in equation (<ref>) converge for wavenumbers k for which none of the coefficients β_r is equal to zero. In such cases it can be shown that G^q_k can be expressed in the frequency domain in the form

    G^q_k(x_,x_2)=i/2d∑_r∈ℤe^iα_r x_1+iβ_r|x_2|/β_r.

  
  
When the wavenumber k is a Wood frequency, the set W=W(k):={r_0∈ℤ:β_r_0(k)=0} is nonempty.

For wavenumbers that are Wood frequencies, the series in the definition of the Green function G^q_k in equation (<ref>) does not converge. In the case when k is a Wood frequency, we introduce the following shifted Green functions <cit.>

    G^q,j_k,h(x_,x_2)=∑_n∈ℤ e^-iα nd∑_ℓ=0^j(-1)^ℓjℓG_k(x_1+nd,x_2+ℓ h)+∑_r∈ W c_r e^iα_r x_1+i(h)β_r x_2

for shifts h≠ 0, integers j>0, and non-zero coefficients c_r∈ℂ. The functions G^q,j_k,h are radiating α-quasi-periodic Green function in the halfplane x_2>0 for h>0 and respectively in the halfplane x_2<0 for h<0; these functions have poles at x_1=0 and x_2=-ℓ h, 0<ℓ≤ j. We note that the quantities G^q,j_k,h defined in equation (<ref>) still make sense when k is not a Wood frequency, in which case the set W can be defined as W:={r_0∈ℤ:|β_r_0|<ε}, where ε is chosen to be sufficiently small; obviously, the set W can be empty in some cases.

Assume now that the interface Γ^per is defined as Γ^per:={(x_1,F(x_1)):0≤ x_1≤ d} where F is a C^2 periodic function of principal period equal to d. Given a density φ defined on Γ^per (which can be extended by α-quasi-periodicity to arguments (x_1,F(x_1)), x_1∈ℝ) we define the single-layer potentials corresponding to a wavenumber k

    [SL_k^qφ](𝐱):=∫_Γ^perG^q_k(𝐱,𝐲)φ(𝐲)ds(𝐲)     [SL_k,h^q,jφ](𝐱):=∫_Γ^perG^q,j_k,h(𝐱,𝐲)φ(𝐲)ds(𝐲)

for 𝐱∉Γ^per and 𝐱=(x_1,x_2) such that 0≤ x_1≤ d.  The quantities SL_k^qφ can be extended by α-quasi-periodicity to define α-quasi-periodic outgoing solutions of the Helmholtz equation corresponding to wavenumber k in the domains {𝐱:x_2>F(x_1)} and {𝐱:x_2<F(x_1)}. Similarly, the quantities SL_k,h^q,jφ can be extended by α-quasi-periodicity to define α-quasi-periodic outgoing solutions of the Helmholtz equation corresponding to wavenumber k in the domains {𝐱:x_2>F(x_1)} for h>0 and respectively in the domain {𝐱:x_2<F(x_1)} for h<0. Assuming that 𝐧 is the unit normal to Γ^per pointing into the domain {𝐱:x_2>F(x_1)} one obtains the single layer potential on the interface Γ^per,

    [S_k^q(φ)](𝐱) := lim_ε→ 0 [SL_k^qφ](𝐱±ε𝐧(𝐱))=∫_Γ^perG_k^q(𝐱,𝐲)φ(𝐲)ds(𝐲),  𝐱∈Γ^per

and

    [S_k,h^q,j(φ)](𝐱):=lim_ε→ 0 [SL_k,h^q,jφ](𝐱+ ε𝐧(𝐱))=∫_Γ^perG_k,h^q,j(𝐱,𝐲)φ(𝐲)ds(𝐲),  𝐱∈Γ^per,     h>0

as well as

    [S_k,h^q,j(φ)](𝐱):=lim_ε→ 0 [SL_k,h^q,jφ](𝐱- ε𝐧(𝐱))=∫_Γ^perG_k,h^q,j(𝐱,𝐲)φ(𝐲)ds(𝐲),  𝐱∈Γ^per,     h<0.

Also, we have

    lim_ε→ 0∇[SL_k^qφ](𝐱±ε𝐧(𝐱))·𝐧(𝐱)=∓1/2φ(𝐱)+[(K_k^q)^⊤(φ)](𝐱),  𝐱∈Γ^per

and

    lim_ε→ 0∇[SL_k,h^q,jφ](𝐱+ ε𝐧(𝐱))·𝐧(𝐱)=-1/2φ(𝐱)+[(K_k,h^q,j)^⊤(φ)](𝐱),  𝐱∈Γ^per,     h>0

as well as

    lim_ε→ 0∇[SL_k,h^q,jφ](𝐱- ε𝐧(𝐱))·𝐧(𝐱)=1/2φ(𝐱)+[(K_k,h^q,j)^⊤(φ)](𝐱),  𝐱∈Γ^per,     h<0.

In equations (<ref>), the adjoint double-layer operators can be defined explicitly as

    [(K_k^q)^⊤(φ)](𝐱)=∫_Γ^per∂ G_k^q(𝐱,𝐲)/∂𝐧(𝐱)φ(𝐲)ds(y),  𝐱∈Γ^per

with similar definitions for the operators defined in equations (<ref>) and (<ref>) respectively.


  §.§.§ Boundary-integral  representation of RtR maps

Having defined the α-quasi-periodic boundary-integral operators above, we are now in a position to compute the various RtR operators 𝒮^j. We start with the RtR operator 𝒮^0 corresponding to problem (<ref>). We define Z_0=iηγ_0^-1 and seek w_0 in the form

    w_0:=SL^q_k_0φ_0,

in which the single-layer potential SL^q_k_0 is defined in (<ref>) integrating on the curve Γ_0.  From the equation


    𝒮^0 [(1/2) I + K - Z_0 S_Γ_0,k_0^q]
       = [(1/2) I + K - Z_0 S_Γ_0,k_0^q] ,


we obtain an explicit formula for the RtR operator 𝒮^0 defined in equation (<ref>),


    𝒮^0=I+2Z_0S_Γ_0,k_0^q(1/2I+(K_Γ_0,k_0^q)^⊤-Z_0S_Γ_0,k_0^q)^-1,

in which the operators (K_Γ_0,k_0^q)^⊤ are defined just as in equations (<ref>) but with unit normal n_0 pointing into Ω^-_0 (exterior of Ω_0). Here and in what follows we introduce an additional subscript to make explicit the curve that is the domain of integration of the boundary-integral operators. 

The invertibility of the operator featured in equations (<ref>) can be established in a straightforward manner.

  Under the assumptions that F_0 is C^2, and k_0 is not a Wood frequency, the operator
  
    𝒜_0:=1/2I+(K_Γ_0,k_0^q)^⊤-Z_0S_Γ_0,k_0^q,  𝒜_0:L^2_per(Γ_0)→ L^2_per(Γ_0)

  is invertible with continuous inverse. 


Because of the regularity of the boundary Γ_0, both operators (K_Γ_0,k_0^q)^⊤:L^2_per(Γ_0)→ L^2_per(Γ_0) and S_Γ_0,k_0^q:L^2_per(Γ_0)→ L^2_per(Γ_0) are compact. Thus, the conclusion of the Theorem follows once we establish the injectivity of the operator 𝒜_0. Let φ_0∈ Ker(𝒜_0) and define v_0:=SL_k_0^qφ_0 in ℝ^2∖{(x_1,F_0(x_1)), x_1∈ℝ)}.  The function v_0 is a radiating α-quasi-periodic solution of the Helmholtz equation in Ω_0 with impedance boundary conditions ∂_n_0v_0-Z_0v_0=0 on Γ_0, and thus, by Theorem <ref>, we have that v_0 is identically zero in Ω_0. In particular, it follows that v_0=0 on Γ_0. Thus, v_0 is a radiating, α-quasi-periodic solution of the Helmholtz equation with wavenumber k_0 in the domain Ω_0^-:={(x_1,x_2):x_2<F_0(x_1)} with zero Dirichlet boundary conditions on Γ_0. This implies that v_0 is identically zero in the domain Ω_0^- by uniqueness the Dirichlet problem. Finally, the jump conditions of the normal derivatives of single-layer potentials imply that φ_0=0 on Γ_0.



Alternatively, we can seek w_0 in the form

    w_0:=SL^q,j_k_0,hφ ,

from which we obtain a representation of the RtR operators 𝒮^0 in the form

    𝒮^0=I+2Z_0S_Γ_0,k_0,h^q,j(1/2I+(K_Γ_0,k_0,h^q,j)^⊤-Z_0S_Γ_0,k_0,h^q,j)^-1,   h>0.

The invertibility of the operators that feature in equation (<ref>) is much more subtle.  It can be established by modification of arguments presented in a recent paper of some of the authors <cit.>.  The proof is presented there for doubly periodic layered media in three dimensions.  The theorem below is also valid in three dimensions; its proof in the two-dimensional reduction is given in Appendix <ref>. 


  Under the assumption that F_0 is C^2 and that k_0 is a Wood frequency, the operator
  
    𝒜_0,h:=1/2I+(K_Γ_0,k_0,h^q,j)^⊤-Z_0S_Γ_0,k_0,h^q,j, j≥ 1, 𝒜_0,h:L^2_per(Γ_0)→ L^2_per(Γ_0)

  is invertible with continuous inverse for all but a discrete set of values of the shift h>0.


Note that the calculations of the RtR maps 𝒮^N+1 can be performed similarly to arrive at

    𝒮^N+1=I+2Z_N+1S_Γ_N,k_N+^q(1/2I+(K_Γ_N,k_N+1^q)^⊤-Z_N+1S_Γ_N,k_N+1^q)^-1,


    𝒮^N+1=I+2Z_N+1S_Γ_N,k_N+1,h^q,j(1/2I+(K_Γ_N,k_N+1,h^q,j)^⊤-Z_N+1S_Γ_N,k_N+1,h^q,j)^-1,   h<0.

where Z_N+1=iηγ_N+1^-1 and the adjoint double-layer operators are defined with respect to the unit normal n_N+1 pointing outside of the domain Ω_N+1. The invertibility of the operators in equations (<ref>) and (<ref>) can be established analogously to the results in Theorems <ref> and <ref>.

Finally, the RtR maps for the domains Ω_j,1≤ j<N can be expressed in closed form via boundary-integral equations. Indeed, consider a generic domain Ω^per:={(x_1,x_2):0≤ x_1≤ d, F_b(x_1)≤ x_2≤ F_t(x_1)} where F_t and F_b are d-periodic C^2 functions. Let us denote by Γ_t:={(x_1,x_2):0≤ x_1≤ d, x_2=F_t(x_1)}, Γ_b:={(x_1,x_2):0≤ x_1≤ d, x_2=F_b(x_1)}, and let n denote the unit normal to Γ_t∪Γ_b pointing outside of the domain Ω^per—see Figure <ref>. Then the Helmholtz problems (<ref>) can be all expressed in the generic form

    Δ w+k^2w   =   0   in Ω^per
    ∂_nw-Z w   =   g_t   on Γ_t
    ∂_nw-Z w   =   g_b   on Γ_b

where g_t and g_b are α-quasi-periodic functions and Z>0.  The RtR operators 𝒮^j,1≤ j<N are related to the following RtR operator associated with the Helmholtz problems (<ref>):

    𝒮[ g_t; g_b ]:=[ (∂_nw+Z w)|_Γ_t; (∂_nw+Z w)|_Γ_b ].

Seeking the solution w of equations (<ref>) in the form

    w = SL_k,t^qφ_t+SL_k,b^qφ_b ,

in which SL_k,t^q (SL_k,b^q) denotes the quasi-periodic single-layer potential whose domain of integration in Γ_t (Γ_b), we arrive at the following expression for the RtR operator 𝒮: 

    𝒮=[ I 0; 0 I ]+2Z[ S^q_k,t,t S^q_k,b,t; S^q_k,t,b S^q_k,b,b ][ (1/2)I+(K^q_k,t,t)^⊤ +Z S^q_k,t,t         (K^q_k,b,t)^⊤+Z S^q_k,b,t;         (K^q_k,t,b)^⊤+Z S^q_k,t,b (1/2)I+(K^q_k,b,b)^⊤ +Z S^q_k,b,b ]^-1.

We note that in equation (<ref>), the subscripts in the notation S^q_k,b,t signify that in equation (<ref>) the target point 𝐱∈Γ_t and the integration point 𝐲∈Γ_t, whereas in the notation S^q_k,b,t signify that in equation (<ref>) the target point 𝐱∈Γ_t and the integration point 𝐲∈Γ_b.  All the other additional subscripts in equation (<ref>) have similar meanings related to the locations of target and integration points for single and adjoint double-layer boundary-integral operators. The invertibility of the operators featuring in equation (<ref>) can be established using similar reasoning to that in the proof of Theorem <ref>.

In the case when k is a Wood frequency, an equivalent representation of the RtR operator 𝒮 can be obtained if we replace the quasi-periodic boundary-integral operators in equation (<ref>) by shifted quasi-periodic boundary-integral operators, provided the shift h>0 is chosen larger than the width of the domain Ω^per; the latter requirement is needed to ensure that no poles of the shifted quasi-periodic functions are contained in the domain Ω^per. The invertibility of the ensuing matrix operator is the subject of the following theorem:

  Assume k is a Wood frequency. Then the operator
  
    𝒜_h:=[ (1/2)I+(K^q,j_k,h,t,t)^⊤ +Z S^q,j_k,h,t,t         (K^q,j_k,h,b,t)^⊤+Z S^q,j_k,h,b,t;         (K^q,j_k,h,t,b)^⊤+Z S^q,j_k,h,t,b (1/2)I+(K^q,j_k,h,b,b)^⊤ +Z S^q,j_k,h,b,b ]

  is invertible with continuous inverse in the space L^2_per(Γ_t)× L^2_per(Γ_b) for all but a discrete set of values of the shift h>0.

 Given that all the boundary-integral operators that enter the definition of the matrix operator 𝒜_h are compact in L^2_per(Γ_t)× L^2_per(Γ_b), the result follows once we establish the injectivity of the operator 𝒜_h. Let (φ_t,φ_b)∈ Ker(𝒜_h) and define
  
    w = SL_k,h,t^q,jφ_t+SL_k,h,b^q,jφ_b   in ℝ^2∖(Γ_t∪Γ_b).

  Clearly w is a α-quasi-periodic solution of equation (<ref>) with zero Robin boundary conditions on Γ_t and Γ_b, and as such w=0 in Ω^per. In particular, w vanishes on Γ_t. Also, given that the shift h is chosen so that the poles of G^q,j_k,h are in the domain Ω_b^-:={(x_1,x_2):0≤ x_1≤ d, x_2≤ F_b(x_1)}, w is a radiating α-quasi-periodic solution of the Helmholtz equation in the domain Ω_t^+={(x_1,x_2):0≤ x_1≤ d, F_t(x_1)≤ x_2}, which vanishes on Γ_t. This means that w=0 in Ω^+ <cit.>. Using the jump conditions of the normal derivatives of single-layer potentials across Γ_t, we get that φ_t=0 on Γ_t. Accordingly, we have that
  
    w = SL_k,h,b^q,jφ_b   in ℝ^2∖Γ_b

  vanishes in the domain Ω_b^+={(x_1,x_2):0≤ x_1≤ d, F_b(x_1)≤ x_2}. The arguments in the proof of Theorem <ref> can be repeated verbatim to conclude that φ_b=0 on Γ_b.
  
  



  The computation of the layer RtR maps described in equations (<ref>) can be extended in a straightforward manner to the case when impenetrable or penetrable inclusions are present in the domain Ω^per.  In this case, the matrix 𝒮 in equations (<ref>) needs be augmented by blocks that account for the interactions of the inclusions D with Γ_t and Γ_b, as well as its self-interactions that account for the boundary conditions to be imposed on ∂ D. 
  






§ DDM NYSTRÖM DISCRETIZATION


Our numerical solution of equations (<ref>)

relies on Nyström discretizations of the boundary-integral operators featured in the computation of the RtR operators given in Section <ref>. In order to speed up the notoriously slow convergence of the quasi-periodic Green function G^q_k defined in equation (<ref>) for frequencies that are away from Wood frequencies, we make use of the recently introduced windowed Green function Method <cit.>. Specifically, let χ(r) be a smooth cutoff function equal to 1 for r<r_1 and equal to 0 for r>r_2 (0<r_1<r_2) and define the windowed Green functions

    G^q,A_k(x_,x_2)=∑_n∈ℤ e^-iα ndG_k(x_1+nd,x_2)χ(r_n/A),   r_n=((x_1+nd)^2+x_2^2)^1/2

and

    G^q,j,A_k,h(x_,x_2)   =   ∑_n∈ℤ e^-iα nd∑_ℓ=0^j(-1)^ℓjℓG_k(x_1+nd,x_2+ℓ h)χ(r_n,ℓ/A)
       +   ∑_r∈ W c_r e^iα_r x_1+i(h)β_r x_2,    r_n,ℓ=((x_1+nd)^2+(x_2+ℓ h)^2)^1/2.

On account of the windowed function χ, the summations in equations (<ref>) and (<ref>) are over a finite range of indices n. The functions G^q,A_k were shown to converge superalgebraically fast to G^q_k as A→∞ when k is not a Wood frequency <cit.>, whereas the functions G^q,j,A_k,h were shown to converge algebraically fast to a α-quasi-periodic Green function as A→∞ (the rate increases as the number of shifts j grows) in the half-plane x_2>0 when h>0 and respectively x_2<0 when h<0 for all frequencies k, including at and around Wood frequencies <cit.>.

Our Nyström discretizations rely on trigonometric collocation in two dimensions. As such, we reformulate the DDM system in terms of periodic quantities by extracting the phase e^-iα x_1 from all Robin data, the right-hand side, as well as RtR maps. The calculation of the RtR maps is performed via boundary-integral operators acting on periodic densities φ defined as φ(x_1,x_2):=φ(x_1,x_2)e^-iα x_1 and periodic kernels e^iα(x_1-y_1)G_k^q(x_1,x_2;y_1,y_2). Furthermore, the discretization of the boundary-integral operators featured in Section <ref> is done by replacing the Green functions G^q_k and G^q,j_k in their definitions by the fast convergent windowed approximations G^q,A_k and G^q,j,A_k,h defined in equations (<ref>) and (<ref>) respectively.  Finally, our numerical scheme requires a simple modification of the Martensen-Kussmaul (MK) periodic logarithmic splitting Nyström approach <cit.> in order to enable high-order evaluations of boundary-integral operators whose kernels are windowed periodic Green functions—full details of this approach are given in <cit.>. In a nutshell, boundary-integral operators that feature the windowed Green functions G^q,A_k and G^q,j,A_k,h defined in equations (<ref>) and (<ref>) respectively are recast in a form that involves integration around target points 𝐱 but with domains of integration that span the whole real axis—the latter is achieved via the windowing functions χ and by periodic extensions of the densities φ. This setting allows for a direct extension of the periodic logarithmic splitting of the Green functions that is central to MK Nyström approach. In three dimensions, our Nyström discretizations also rely on global trigonometric interpolation, use of floating partitions of unity and analytic resolution of singularities, as well as the use of windowed Green functions <cit.>. Interestingly, using global trigonometric interpolation in conjunction with changes of variables to polar coordinates in the resolution of Green function singularities allows for straightforward constructions of Nyström collocation matrices of three-dimensional boundary-integral operators  via two-dimensional Discrete Fourier Transform matrices. The availability of such Nyström collocation matrices in three dimensions plays an important role in the efficient computations of RtR maps, as we explain next.

Following the prescriptions outlined above, a boundary-integral operator whose kernel is a windowed Green function (or its normal derivative) acting on a periodic density φ and whose domain of integration is a generic curve Γ^per per the definition given in Section <ref> (i.e. Γ^per:={(x_1,F(x_1)):0≤ x_1≤ d}, where F is a C^2 periodic function of principal period d) is Nyström discretized as a M× M matrix where the periodic density φ is trigonometrically collocated at the equi-spaced mesh {(t_ℓ,F(t_ℓ)): t_ℓ=ℓ d/M, 0≤ℓ<M=2m}. For a fixed M=2m, m>0, assuming that the Robin data f_j=[f_j,j f_j,j+1]^⊤ on each interface Γ_j,0≤ j≤ N is collocated at the mesh L_j:={(t_ℓ,F_j(t_ℓ)): t_ℓ=ℓ d/M, 0≤ℓ<M}, it follows that the RtR maps 𝒮^0 and 𝒮^N+1 are discretized as M× M Nyström matrices 𝒮^0_M and 𝒮^N+1_M via Nyström discretizations of the boundary-integral operators featured in equations (<ref>) and (<ref>) respectively in the case when neither k_0 nor k_N+1 are Wood frequencies or in equations (<ref>) and (<ref>) respectively in the case when k_0 and k_N+1 are Wood frequencies (same considerations apply in three dimensions). We note that according to equations (<ref>) and (<ref>) (and their analogues (<ref>) and (<ref>)), the calculation of the Nyström matrices 𝒮^0_M and 𝒮^N+1_M require inversions of M× M matrices, which is done using LU factorizations. Similarly, the RtR maps 𝒮^j,1≤ j≤ N are discretized as (2M)× (2M) Nyström matrices 𝒮^j_M  via Nyström discretizations of the boundary-integral operators featured in equations (<ref>), and their calculations require, in turn, inversions of (2M)×(2M) matrices; these inversions are also performed through LU factorizations. It is also possible to employ Schur complements to perform the inversion of the matrices needed in the calculations of the RtR maps 𝒮^j,1≤ j≤ N—see the proof of Theorem <ref>; in that case matrices of size M× M need be inverted. The Nyström discretization matrices 𝒮^j_M are further expressed in M× M block form

    𝒮^j_M=[ 𝒮^j_j-1,j-1,M   𝒮^j_j-1,j,M;   𝒮^j_j,j-1,M     𝒮^j_j,j,M ],

where each of the matrices above constitutes a Nyström discretization matrix of the operators on the right-hand-side of equation (<ref>).  

Based on these Nyström discretizations of RtR maps, the DDM algorithm computes matrix approximations of all the RtR maps needed. Clearly, the procedure outlined above allows us to assemble a block tridiagonal 2M(N+1)× 2M(N+1) Nyström discretization matrix A_M of the operator matrix A in equation (<ref>), where each operator block in equation (<ref>) is replaced by its corresponding Nyström discretization matrix. The computation of collocated Robin data f_j,M at the grids L_j for 0≤ j≤ N requires solution of a linear system featuring the matrix A_M. We present in what follows an efficient algorithm for the solution of this system that eliminates sequentially  the unknowns f_j,M, 0≤ j≤ N using N+1 stages of recursive Schur complements; storage of all of the non-zero blocks in the matrix A_M is not required by this algorithm. The key technical ingredient is that in the case when matrices

    𝒟:=[ I 𝒜; ℬ I ]

are invertible, then their inverses can be computed explicitly

    𝒟^-1=[ I+𝒜(I-ℬ𝒜)^-1ℬ   -𝒜(I-ℬ𝒜)^-1;   -(I-ℬ𝒜)^-1ℬ     (I-ℬ𝒜)^-1 ].


The Schur complement elimination algorithm begins with

Stage 1: elimination of the unknowns f_0,M. We express the discrete DDM system in the following block form that separates the contribution of the Robin data f_0,M from the rest of the Robin data. In detail,

    [ 𝒟_0,M 𝒜_0,M; ℬ_0,M 𝒞_0,M ][ f_0,M; f_0,M ]   =   [    _0,M; 0_2NM,1 ]
    𝒟_0,M   =   [       I_M 𝒮^1_0,0,M;     𝒮^0_M       I_M ]
    𝒜_0,M   =   [   𝒮^1_0,1,M         0_M 0_M,2(N-1)M;         0_M         0_M 0_M,2(N-1)M ]
    ℬ_0,M   =   [         0_M         0_M;         0_M   𝒮^1_1,0,M; 0_2(N-1)M,M 0_2(N-1)M,M ]
    𝒞_0,M   =   [       I_M 𝒮^2_1,1,M         ⋯; 𝒮^1_1,1,M       I_M         ⋯;         ⋯         ⋯         ⋯ ],

where f_0,M=[f_1,M f_2,M … f_N,M]^⊤. In the notations above and in what follows, we make explicit the matrix size of various zero matrices; for instance, the notation 0_p,q denotes a zero matrix with p rows and q columns, and 0_p denotes the zero p× p matrix.  We have

    f_0,M=-𝒟_0,M^-1𝒜_0,Mf_0,M+𝒟_0,M^-1_0,M,

and hence

    (𝒞_0,M-ℬ_0,M𝒟_0,M^-1𝒜_0,M)f_0,M=_1,M, _1,M:=-ℬ_0,M𝒟_0,M^-1_0,M,

which can be written in expanded form using formula (<ref>) to compute the inverse of 𝒟_0,M:

    [ 𝒟_1,M 𝒜_1,M; ℬ_1,M 𝒞_1,M ][ f_1,M; f_1,M ]   =   [        _1,M; 0_2(N-1)M,1 ]
    _1,M   =   -ℬ_0,M𝒟_0,M^-1_0,M
    𝒟_1,M   =   [       I_M 𝒮^2_1,1,M; 𝒮^top_1,M       I_M ]
    𝒮^top_1,M   =   𝒮^1_1,0,M(I_M-𝒮^0_M𝒮^1_0,0,M)^-1𝒮^0_M𝒮^1_0,1,M+𝒮^1_1,1,M
    𝒜_1,M   =   [   𝒮^2_1,2,M         0_M 0_M,2(N-2)M;         0_M         0_M 0_M,2(N-2)M ]
    ℬ_1,M   =   [         0_M         0_M;         0_M   𝒮^2_2,1,M; 0_2(N-2)M,M 0_2(N-2)M,M ]
    𝒞_1,M   =   [       I_M 𝒮^3_2,2,M         ⋯; 𝒮^2_2,2,M       I_M         ⋯;         ⋯         ⋯         ⋯ ],

where f_1,M=[f_2,M f_3,M … f_N,M]^⊤.


Stage j: elimination of the unknowns f_j-1,M. Repeating the same steps outlined above we continue the elimination process until we arrive at the following linear system 

    [ 𝒟_j-1,M 𝒜_j-1,M; ℬ_j-1,M 𝒞_j-1,M ][ f_j-1,M; f_j-1,M ]   =   [        _j-1,M; 0_2(N-j+1)M,1 ]
    𝒟_j-1,M   =   [           I_M 𝒮^j_j-1,j-1,M;   𝒮^top_j-1,M           I_M ]
    𝒜_j-1,M   =   [ 𝒮^j_j-1,j,M         0_M 0_M,2(N-j)M;         0_M         0_M 0_M,2(N-j)M ]
    ℬ_j-1,M   =   [         0_M         0_M;         0_M 𝒮^j_j,j-1,M; 0_2(N-j)M,M 0_2(N-j)M,M ]
    𝒞_j-1,M   =   [         I_M 𝒮^j+1_j,j,M           ⋯;   𝒮^j_j,j,M         I_M           ⋯;           ⋯           ⋯           ⋯ ],

where f_j-1,M=[f_j,M f_j+1,M … f_N,M]^⊤. We proceed by eliminating the unknowns f_j-1,M from the linear system above. First, we have that

    f_j-1,M=-𝒟_j-1,M^-1A_j-1,Mf_j-1,M+𝒟_j-1,M^-1_j-1,M

and thus

    (𝒞_j-1,M-ℬ_j-1,M𝒟_j-1,M^-1A_j-1,M)f_j-1,M=_j,M, _j,M:=-ℬ_j-1,M𝒟_j-1,M^-1_j-1,M.

The last linear system can be in turn written in expanded form if we make use of formula (<ref>) to compute the inverse of 𝒟_j-1,M,

    [ 𝒟_j,M 𝒜_j,M; ℬ_j,M 𝒞_j,M ][ f_j,M; f_j,M ]   =   [        _j,M; 0_2(N-j)M,1 ]
    _j,M   =   -ℬ_j-1,M𝒟_j-1,M^-1_j-1,M
    𝒟_j,M   =   [         I_M 𝒮^j+1_j,j,M;   𝒮^top_j,M         I_M ]
    𝒮^top_j,M   =   𝒮^j_j,j-1,M(I_M-𝒮^top_j-1,M𝒮^j_j-1,j-1,M)^-1𝒮^top_j-1,M𝒮^j_j-1,j,M+𝒮^j_j,j,M
    𝒜_j,M   =   [ 𝒮^j+1_j,j+1,M           0_M 0_M,2(N-j-1)M;           0_M           0_M 0_M,2(N-j-1)M ]
    ℬ_j,M   =   [           0_M           0_M;           0_M 𝒮^j+1_j+1,j,M; 0_2(N-j-1)M,M 0_2(N-j-1)M,M ]
    𝒞_j,M   =   [             I_M 𝒮^j+2_j+1,j+1,M               ⋯; 𝒮^j+1_j+1,j+1,M             I_M               ⋯;               ⋯               ⋯               ⋯ ],

where f_j,M=[f_j+1,M f_j+2,M … f_N,M]^⊤.

The algorithm ends with a linear system involving only the Robin data f_N,M corresponding to the last interface Γ_N,

    [       I_M   𝒮^N+1_M; 𝒮^top_N,M       I_M ]f_N,M=_N,M,  _N,M:=-ℬ_N-1,M𝒟_N-1,M^-1_N-1,M,

which is solved using again formula (<ref>) and LU factorizations. Once the discretized Robin data f_N,M is computed, all the other discretized Robin data corresponding to the interfaces Γ_j,0≤ j<N are computed using backward substitution via the recursions

    f_j-1,M=-𝒟_j-1,M^-1A_j-1,Mf_j-1,M+𝒟_j-1,M^-1_j-1,M,1≤ j≤ N,

where we recall that f_j-1,M=[f_j,M f_j+1,M … f_N,M]^⊤. In order to streamline the recursions above, we store in the elimination process the quantities 𝒟_j-1,M^-1_j-1,M and the non-zero blocks of the matrices 𝒟_j-1,M^-1A_j-1,M (for each index j there are only two M× M such blocks). Thus, the storage require to perform the elimination algorithm followed by the backward substitution process is 𝒪(2NM^2).

We first note that the elimination process presented above consists of N+1 steps, each step requiring (a) the calculation of the Nyström matrices 𝒮^j_M whose cost is 𝒪(M^3), as well as (b) the inversion of a M×M matrix for the calculation of 𝒟^-1_j,M per formula (<ref>), which leads to a total computational cost proportional to (N+1)M^3. One drawback of the elimination algorithm presented above is that it is sequential: essentially the Robin data is peeled off layer by layer from the DDM linear system. We are exploring alternative strategies based on hierarchical RtR matrices mergings that could lead to efficient parallelization strategies <cit.>.  Of course, for a large numbers of layers, iterative solvers such as GMRES could also provide an alternative strategy for the solution of the DDM linear system. However, large numbers of layers/subdomains produce large numbers of DDM iterations <cit.> if classical Robin conditions are used on subdomain interfaces.



§ ANALYSIS OF THE SCHUR COMPLEMENT ELIMINATION ALGORITHM


A natural question is why the elimination procedure described in Section <ref> does not break down. This issue has been explored in a different context in <cit.>, where the elimination process was shown to be equivalent to merging of RtR maps. Indeed, the matrices 𝒮^top_j,M are themselves Nyström  discretization matrices of the RtR operators

    𝒮^top,j(ψ_j):=(γ_j∂_n_j u+ iη  u)|_Γ_j,

 where u is the solution of the following well posed problem:
 
    Δ u+k(x)^2u   =   0  in ∪_ℓ=0^jΩ_ℓ,   k(x)=k_ℓ, x∈Ω_ℓ, 0≤ l≤ j
    γ_j∂_n_j u+ iη  u    =   ψ_j  on Γ_j,

 with (i) u and γ_ℓ∂_n_ℓu continuous across Γ_ℓ for 0≤ l<j; and (ii) u radiating in Ω_0. Thus, the recurrence formula
 
    𝒮^top_j,M=𝒮^j_j,j-1,M(I_M-𝒮^top_j-1,M𝒮^j_j-1,j-1,M)^-1𝒮^top_j-1,M𝒮^j_j-1,j,M+𝒮^j_j,j,M,   for   1≤ j,

 where 𝒮^top_0,M:=𝒮^0_M, can be viewed as a means to compute Nyström discretization of the RtR map 𝒮^top,j defined in equation (<ref>) via recursive merging of
 the Nyström discretization matrices of the RtR maps 𝒮^ℓ. We will establish in this section the invertibility of the operators I-𝒮^top_j-1𝒮^j_j-1,j-1 in appropriate functional spaces. 
 
The matrices I_M-𝒮^top_j-1,M𝒮^j_j-1,j-1,M in (<ref>) are high-order Nyström discretizations of the operators I-𝒮^top_j-1𝒮^j_j-1,j-1. As such, the invertibility of the former matrices is a consequence of the invertibility of the latter operators owing to the fact that  𝒮^top_j-1,M and 𝒮^j_j-1,j-1,M can be shown to converge in strong operator norms to 𝒮^top_j-1 and 𝒮^j_j-1,j-1 as M→∞ <cit.>.

 We start by analyzing the invertibility of the operators I-𝒮^0𝒮^1_0,0, on which hinges Stage 1 of the elimination algorithm. We will make use of the quasi-periodic Sobolev spaces H^s_per(Γ) of α-quasi-periodic distributions defined on a generic interface Γ that is the graph of a periodic function of period d. These spaces can be defined in terms of Fourier series. We use the mapping properties of boundary-integral operators associated with quasi-periodic Green functions G^q_k and shifted quasi-periodic Green functions G_k,h^q,j. Since both of these functions have the same singularity as the free space Green function G_k corresponding to the same wavenumber k, the mapping properties of the boundary-integral operators associated with quasi-periodic Green functions can be easily derived by simply translating to the periodic setting the classical mapping properties of boundary-integral operators whose domain of integration is a closed curve in ℝ^2. The invertibility of the operator  I-𝒮^0𝒮^1_0,0 is established via Fredholm arguments, and it relies on the explicit representations of the operators 𝒮^0 and 𝒮^1 derived in equations (<ref>) and (<ref>) in the case when neither k_0 nor k_1 are Wood frequencies. In the cases when k_0 or k_1 are Wood frequencies, representations (<ref>) and that given in Theorem <ref> ought to be used.  Nevertheless, given that the shifted quasi-periodic Green functions have the same singularities as the quasi-periodic Green functions, in the computational domains under considerations, the arguments in the proof of the result below essentially do not change in the Wood-frequency case.  In what follows, we suppose that (1) k_j<k_j+1 for all j and (2) γ_j=1 for all j (we will show in Appendix <ref> that these two assumptions guarantee the well-posedness of the transmission problem under consideration). We establish
 
   The operator I-𝒮^0𝒮^1_0,0:H^-1/2_per(Γ_0)→ H^1/2_per(Γ_0) is Fredholm of index 0 under the assumption that Γ_0 is C^2. 
   
 
   Let us assume that k_0 is not a Wood frequency. We first establish that, given the representation
   
    𝒮^0=I+2Z_0S_Γ_0,k_0^q(1/2I+(K_Γ_0,k_0^q)^⊤-Z_0S_Γ_0,k_0^q)^-1,

 the operator 𝒮^0 can be expressed in the form
   
    𝒮^0=I+4Z_0S_Γ_0,k_0+iε_0^q+𝒯_0,  𝒯_0:H^-1/2_per(Γ_0)→ H_per^3/2(Γ_0),  ε_0>0,   Z_0=iη.

   The decomposition in equation (<ref>) can be achieved if we choose the operator 𝒯_0 in the following manner:
   
    𝒯_0=4Z_0(S_Γ_0,k_0^q-S_Γ_0,k_0+iε_0^q)+8Z_0S_Γ_0,k_0^q(1/2I+2(K_Γ_0,k_0^q)^⊤-2Z_0S_Γ_0,k_0^q)^-1((K_Γ_0,k_0^q)^⊤-Z_0S_Γ_0,k_0^q).

   Now, given that S_Γ_0,k_0^q-S_Γ_0,k_0+iε_0:H^-1/2_per(Γ_0)→ H_per^3/2(Γ_0), S_Γ_0,k_0^q:H^s_per(Γ_0)→ H_per^s+1(Γ_0) for -1/2≤ s≤ 1, and (K_Γ_0,k_0^q)^⊤:H^s_per(Γ_0)→ H_per^s+2(Γ_0) for -1/2≤ s≤ 0, it follows  that 𝒯_0:H^-1/2_per(Γ_0)→ H_per^3/2(Γ_0). The role of the complex wavenumber k_0+iε_0 in the definition of the quasi-periodic single-layer operators in equation (<ref>) will be made clear in what follows.

   We establish next a decomposition of the operator 𝒮^1_0,0 similar in spirit to that in equation (<ref>). To this end, we revisit representation (<ref>), valid for a general periodic layer and a wavenumber k that is not a Wood frequency. Let us define 
   
    𝒜:= [ (1/2)I+(K^q_k,t,t)^⊤ +Z S^q_k,t,t         (K^q_k,b,t)^⊤+Z S^q_k,b,t;         (K^q_k,t,b)^⊤+Z S^q_k,t,b (1/2)I+(K^q_k,b,b)^⊤ +Z S^q_k,b,b ]=[ 𝒜_t,t 𝒜_t,b; 𝒜_b,t 𝒜_b,b ].

  We have
  
    𝒜^-1=[ 𝒜_t,t^-1+𝒜_t,t^-1𝒜_t,b𝒟𝒜_b,t𝒜_t,t^-1                      -𝒜_t,t^-1𝒜_t_b𝒟;                      -𝒟𝒜_b,t𝒜_t,t^-1                                    𝒟 ]=[ 𝒜_t,t 𝒜_t,b; 𝒜_b,t 𝒜_b,b ]

  where
  
    𝒟=(𝒜_b,b-𝒜_b,t𝒜_t,t^-1𝒜_t,b)^-1.

  The invertibility of the operator 𝒜_t,t needed in formulas above can be established similarly to the result in Theorem <ref>.  The invertibility of the operator 𝒟, in turn, can be seen to be equivalent to the invertibility of the matrix operator 𝒜. Given that the kernels of the boundary-integral operators that enter in the definition of 𝒜_t,b and 𝒜_b,t are regular, we have 𝒜_t,b:H^s_per(Γ_t)→ H_per^s+2(Γ_b) and 𝒜_b,t:H^s_per(Γ_b)→ H_per^s+2(Γ_t).   Since
  
    𝒮=[ I 0; 0 I ]+2Z[ S^q_k,t,t S^q_k,b,t; S^q_k,t,b S^q_k,b,b ]𝒜^-1=[ 𝒮_t,t 𝒮_t,b; 𝒮_b,t 𝒮_b,b ],

  we have
  
    𝒮_t,t   =I+2ZS^q_k,t,t𝒜_t,t+2ZS^q_k,b,t𝒜_b,t=I+4ZS^q_k_0+iε_0,t,t+𝒮_t,t, 
    𝒮_t,t   :H^-1/2_per(Γ_t)→ H_per^3/2(Γ_t)

taking into account the fact that S^q_k,b,t:H^s_per(Γ_b)→ H_per^s+2(Γ_t). In addition, we obtain
 
    𝒮_b,b   =I+2ZS^q_k,t,b𝒜_t,b+2ZS^q_k,b,b𝒜_b,b=I+4ZS^q_k_0+iε_0,b,b+𝒮_b,b, 
    𝒮_b,b   :H^-1/2_per(Γ_b)→ H_per^3/2(Γ_b)

 as well as the following smoothing  properties of the cross operators 𝒮_t,b and 𝒮_b,t,
 
    𝒮_t,b:H^-1/2_per(Γ_t)→ H_per^3/2(Γ_b),    𝒮_b,t:H^-1/2_per(Γ_b)→ H_per^3/2(Γ_t).

Combining the results in (<ref>) and (<ref>), we obtain
  
    I-𝒮^0𝒮^1_0,0=-8Z_0S_Γ_0,k_0+iε_0^q+𝒯_0,1,    𝒯_0,1:H^-1/2_per(Γ_0)→ H_per^3/2(Γ_0).

  Classical arguments <cit.> can be adapted to the periodic setting to establish
  
    (S_Γ_0,k_0+iε_0^qφ,φ)≥ c_0φ^2_H^-1/2_per(Γ_0), c_0>0,

  from which we obtain
  
    (-8Z_0S_Γ_0,k_0+iε_0^qφ,φ)≥  8η c_0φ^2_H^-1/2_per(Γ_0), c_0>0, η>0.

  Finally, given that 𝒯_0,1:H^-1/2_per(Γ_0)→ H_per^3/2(Γ_0), it follows that the operator 𝒯_0,1:H^-1/2_per(Γ_0)→ H_per^1/2(Γ_0) is compact, and thus the operator I-𝒮^0𝒮^1_0,0:H^-1/2_per(Γ_0)→ H_per^1/2(Γ_0) can be seen to satisfy a Gårding inequality. The result of the Theorem is thus established.  The case when k_0 and/or k_1 are Wood frequencies can be treated similarly. 
 

 We are now in the position to prove
 
   The operator I-𝒮^0𝒮^1_0,0:H^-1/2_per(Γ_0)→ H^1/2_per(Γ_0) is invertible with continuous inverse.
 
  Owing to the Fredholm alternative, the theorem follows once we establish the injectivity of the operator I-𝒮^0𝒮^1_0,0. Let φ∈ Ker(I-𝒮^0𝒮^1_0,0) and consider the following α-quasi-periodic Helmholtz equation
  
    Δ w_1+k_1^2w_1   =   0     in Ω_1^per
    ∂_n_1w_1-Z_0w_1   =   φ     on Γ_0,
    ∂_n_1w_1-Z_0w_1   =   0     on Γ_1,

  where Z_0=iη. Then
  
    𝒮^1_0,0φ=(∂_n_1w_1+ Z_0w_1)|_Γ_0.

  Consider also the α-quasi-periodic Helmholtz equation
  
    Δ w_0+k_0^2w_0   =   0     in Ω_0
    ∂_n_0w_0-Z_0w_0   =   𝒮^1_0,0φ     on Γ_0,

 with w_0 radiating.  Then using the fact that 𝒮^0𝒮^1_0,0φ=φ on Γ_0, we obtain
  
    𝒮^0𝒮^1_0,0φ=∂_n_0w_0+Z_0w_0=∂_n_1w_1-Z_0w_1   on Γ_0.

  Thus, we have derived the following system of equations on Γ_0
  
    ∂_n_0w_0-Z_0w_0   =   ∂_n_1w_1+Z_0w_1
    ∂_n_0w_0+Z_0w_0   =   ∂_n_1w_1-Z_0w_1.

  from which we obtain
  
    w_0|_Γ_0=-w_1|_Γ_0,    ∂_n_0w_0|_Γ_0=∂_n_1w_1|_Γ_0.

  Recall from the proof of Theorem <ref> the following identity for w_0:
 
    lim_h→∞∫_Ω_0,h^per(|∇ w_0|^2-k_0^2|w_0|^2)dx=∫_Γ_0∂_n_0w_0 w_0 ds+id∑_r∈ℤ, β_0,r>0β_0,r|C_r^+|^2 .

On the other hand,

    ∫_Ω_1^per(|∇ w_1|^2-k_1^2|w_1|^2)dx=∫_Γ_0∂_n_1w_1 w_1 ds+iη∫_Γ_1 |w_1|^2 ds.

Adding the last two identities and taking into account equation (<ref>), we derive

    lim_h→∞∫_Ω_0,h^per(|∇ w_0|^2-k_0^2|w_0|^2)dx+∫_Ω_1^per(|∇ w_1|^2-k_1^2|w_1|^2)dx=id∑_r∈ℤ, β_0,r>0β_0,r|C_r^+|^2+iη∫_Γ_1 |w_1|^2 ds.

This implies that w_1=0 on Γ_1, and hence w_1=0 in Ω_1 by Holmgren's theorem <cit.>. Using (<ref>) again we obtain that w_0=0 and ∂_n_0w_0=0 on Γ_0, which, in turn, implies that w_0=0 in Ω_0. From this we finally conclude that φ=0 on Γ_0.
 

 As a consequence of the results in Theorem <ref> and Theorem <ref> we obtain 

  The operator 𝒮^top_1 defined in equation (<ref>) can be expressed in the form
   
    𝒮^top_1=I+4Z_0S_Γ_1,k_1+iε_1^q+𝒯_1,  𝒯_1:H^-1/2_per(Γ_1)→ H_per^3/2(Γ_1), ε_1>0.

 
  Consider the following Helmholtz scattering problem with Robin boundary conditions on Γ_1:  Find α-quasi-periodic solutions u_0 and u_1 such that 
   
    Δ u_0+k_0^2u   =   0  in Ω_0
    Δ u_1+k_1^2u_1   =   0  in Ω_1
    
         u_0   =   u_1  on Γ_0
    ∂_n_0u_0   =   -∂_n_1u_1  on Γ_0
    γ_1∂_n_1 u+ iη  u_1    =   ψ_1  on Γ_1,

and u_0 is radiating in Ω_0. The arguments in the proof of Theorem <ref> can be applied to show the well-posedness of the this problem. Reformulating it in terms of matching Robin data on Γ_0 and applying the same arguments as in Section <ref> but at the operator (continuous) level, we obtain
   
    𝒮^top_1=𝒮^1_1,0(I-𝒮^0𝒮^1_0,0)^-1𝒮^0𝒮^1_0,1+𝒮^1_1,1.

Using the representation above together with the properties recounted in equations (<ref>) and (<ref>) with Γ_t=Γ_0, Γ_b=Γ_1, and Z=Z_0=iη, the result follows.
 

 The procedure presented above can be repeated inductively at the operator level, and can be viewed as means to recursively merge the RtR operators 𝒮^top_j-1 and 𝒮^j in order to obtain the RtR operator 𝒮^top_j according to the formula  
 
    𝒮^top_j=𝒮^j_j,j-1(I-𝒮^top_j-1𝒮^j_j-1,j-1)^-1𝒮^top_j-1𝒮^j_j-1,j+𝒮^j_j,j, 2≤ j.

 Equation (<ref>) constitutes the continuous analogue of equation (<ref>). The invertibility of the operators I-𝒮^top_j-1𝒮^j_j-1,j-1 for 2≤ j≤ N can be established similarly to the results in Theorem <ref> and Theorem <ref> using the link in Corollary <ref>. Indeed, it is straightforward to establish by induction that
 
    𝒮^top_j-1=I+4Z_0S_Γ_j-1,k_1+iε_j-1^q+𝒯_j-1,  𝒯_j-1:H^-1/2_per(Γ_j-1)→ H_per^3/2(Γ_j-1), ε_j-1>0, 3≤ j≤ N.

 We note that this latter representation suffices to establish the Fredholm property of the operator I-𝒮^top_j-1𝒮^j_j-1,j-1 for all j: 3≤ j≤ N; see Theorem <ref>. The arguments in the proof of Theorem <ref> also translate almost verbatim to obtain the invertibility of the operators I-𝒮^top_j-1𝒮^j_j-1,j-1 for all j:3≤ j≤ N. It is the very last step in the algorithm when we merge 𝒮^top_N and 𝒮^N+1 that is markedly different on account of the fact that the layer Ω_N+1 is semi-infinite and thus the arguments in the proof of Theorem <ref> have to be modified according to those in Theorem <ref> in Appendix <ref>.
 
   The operator I-𝒮^top_N𝒮^N+1:H^-1/2_per(Γ_N)→ H^1/2_per(Γ_N) is invertible with continuous inverse.
 
   Note that the Fredholm property of the operator I-𝒮^top_N𝒮^N+1 essentially follows via the same arguments as in Theorem <ref>. Owing to the Fredholm alternative, the result follows once we establish the injectivity of the operator I-𝒮^top_N𝒮^N+1. Let φ∈ Ker(I-𝒮^top_N𝒮^N+1) and consider the following α-quasi-periodic Helmholtz equation
  
    Δ w_N+1+k_N+1^2w_N+1   =   0     in Ω_N+1^per
    ∂_n_N+1w_N+1-Z_0w_N+1   =   φ     on Γ_N

  with w_N+1 radiating in Ω_N+1 and Z_0=iη. Then
  
    𝒮^N+1φ=(∂_n_N+1w_N+1+ Z_0w_N+1)|_Γ_N.

  Consider also the following α-quasi-periodic Helmholtz equation
  
    Δ w+k(x)^2w   =   0     in ∪_ℓ=0^NΩ_ℓ, 
    
        k(x)   =   k_ℓ     in Ω_ℓ
    ∂_n_Nw-Z_0w   =   𝒮^N+1φ     on Γ_N,

  where w and ∂_n_ℓw are continuous across the interfaces Γ_ℓ for 0≤ℓ≤ N-1 and w is radiating in the domain Ω_0. We have then
  
    𝒮^top_N𝒮^N+1φ=∂_n_Nw+Z_0w=∂_n_N+1w_N+1-Z_0w_N+1   on Γ_N

  using the fact that 𝒮^top_N𝒮^N+1φ=φ on Γ_N. Thus
  
    w|_Γ_N=-w_N+1|_Γ_N    ∂_n_Nw|_Γ_N=∂_n_N+1w_N+1|_Γ_N.

  Applying Green's identities in each domain Ω_j^per for 0≤ j≤ N+1 and taking into account the continuity conditions (<ref>) together with the continuity of w and its normal derivatives across interfaces Γ_j,0≤ j≤ N-1, we obtain that C^+_r=0 for all indices r such that β_0,r>0 and C^-_r=0 for all indices r such that β_N+1,r>0. The proof of the Theorem follows by applying analogous arguments as in the proof of Theorem <ref> to ∂_x_2w (which is continuous across the interfaces Γ_j for 0≤ j≤ N-1) and -∂_x_2w_N+1.
 

 We presented in this section an explanation of the fact that the Schur complement elimination process described in Section <ref> does not break down. Incidentally, we obtain as a byproduct of the results in this section a proof of the equivalence between the scattering PDE problem (<ref>) and its DDM formulation (<ref>) under the assumption that the former is well posed.  According to the Fredholm results established in Theorem <ref> and Theorem <ref>, the DDM formulation (<ref>) requires inversions of operators that are compact perturbations of single layer boundary integral operators. As such, the DDM formulation (<ref>) is not particularly suitable to Krylov subspace linear algebra solvers, especially for configurations that involve large numbers of layers. It is possible to derive DDM formulations that are more amenable to Krylov subspace linear algebra solvers if more general transmission operators are used instead of the multiplicative factors iη in the Robin conditions <cit.>. We are investigating such an approach in the context of periodic layered media. 



§ NUMERICAL RESULTS


We present a variety of numerical results regarding transmission scattering problems in periodic layered media. In all cases we assume that all the coefficients γ_j that appear in equations (<ref>) are equal to 1, and all the numerical results presented are at normal incidence. Qualitatively similar results are obtained in the case of general γ_j and oblique incidence. We present two error indicators, one concerning the energy balance and one concerning errors in the Rayleigh coefficient B_0^+ in expansion (<ref>). The energy conservation defect is defined as

    ε_en=|∑_r∈ U^+β_0,r/β_0|C_r^+|^2+∑_r∈ U^-β_N+1,r/β_0|C_r^-|^2-1|

where U^+:={r∈ℤ:β_0,r≥ 0} and U^-:={r∈ℤ:β_N+1,r≥ 0}. In all the numerical tests presented in this section, the energy conservation defect turned out to be indicative of the errors achieved in the Rayleigh coefficients of the scattered/transmitted fields. We also denote by ε_1 the relative error achieved in the Rayleigh coefficient C_0^+, measured against a reference solution that was produced with refined discretizations and large enough values of the parameter A in the windowed quasi-periodic functions (<ref>) and (<ref>) respectively. This parameter is chosen to be large enough so that very small energy conservation defects were achieved.  We do not know of a theoretical way to determine the parameter A in the various layers that optimizes the balance between accuracy and efficiency. In practice, this parameter is selected using information gathered from numerical experiments.  Particularly, we choose values of A large enough so that the RtR discretizations 𝒮^j_M have norms as close to 1 as possible.

At the heart of our DDM algorithm are computations of RtR maps, which rely on evaluations of boundary-integral operators involving quasi-periodic Green functions. The quasi-periodic functions are approximated via windowing functions cf. (<ref>) and (<ref>). We discuss the selection of various parameters that enter the definition of the windowed Green function defined in equation (<ref>) and the shifted windowed Green function defined in equation (<ref>). In cases when the wavenumbers k_ℓ are not Wood frequencies we had to choose only the parameter A in the definition of the windowed Green function G_k_ℓ^q,A defined in equation (<ref>); the windowed Green function G_k_ℓ^q,A converge superalgebraically to the quasi-periodic Green function G_k_ℓ^q as A→∞ <cit.>. On the other hand, in cases when k_ℓ is a Wood frequency, we had to select  two additional parameters in the definition of G_k_ℓ,h_ℓ^q,j,A (<ref>): the number j of shifts and the value of the shift h_ℓ. The rate of convergence of the Green functions G_k_ℓ,h_ℓ^q,j,A is algebraic in ℓ as A→∞ <cit.>. Naively, the cost of evaluating G_k_ℓ,h_ℓ^q,j,A is j+1 times more expensive than that of evaluating G_k_ℓ^q,A for a fixed value of A.  However, the quantities G_k_ℓ,h_ℓ^q,j,A can be evaluated at considerably reduced costs via accurate asymptotic expansions <cit.>. The values of the shifts h_ℓ should be such that the poles of the Green function G_k_ℓ,h_ℓ^q,j,A are outside the computational domain Ω_ℓ; this requirement entails that h_0>0 and h_N+1<0, whereas the shifts h_ℓ,1≤ℓ≤ N should be positive and larger than the width of the corresponding layer domain Ω_ℓ measured in the x_2 direction. In addition to these requirements on the shifts h_ℓ, there are discrete sets of values of the shifts for which the shifted Green functions G_k_ℓ,h_ℓ^q,j,A do not converge, and those sets can be explicitly computed <cit.>. Indeed, as explained in the proof of Theorem <ref>, in the domain Ω_ℓ the forbidden set of shifts consists of values h_ℓ such that e^iβ_r h_ℓ= 1 for an index r∈ U_ℓ, where U_ℓ is the set of propagating modes corresponding to the wavenumber h_ℓ. In practice it is straightforward to choose the shifts so that all the requirements specified above are met.
Of course, one can use a fixed number of shifts in the shifted Green function for all layers at all frequencies, and this is guaranteed to work as long as the discrete set of shifts h is avoided.  This set is easy to compute for a given structure.   If one wishes to optimize the computational performance, the Wood frequencies for each layer can be computed a priori to decide whether the shifted Green function or simply the smoothly windowed one should be used.


After selection of the parameters that enter the various windowed Green functions, the DDM algorithm is implemented according to its description in Section <ref>. We detail in all the numerical experiments the size M of the discretization points used to approximate each Robin data f_j. The DDM linear system consists of 2(N+1)M unknowns. We mention that it is also possible to use non-conforming discretizations of Robin data, that is to use different M_j for each layer Ω_j; the values of M_j are chosen to resolve the wavenumbers k_j as well as the profiles Γ_j <cit.>–see Table <ref>.  The Schur complement elimination algorithm described in Section <ref> allows for solution of large DDM linear systems using only limited memory storage. For instance, the numerical experiments in Table <ref> and Table <ref> involving DDM linear systems with 81920 and respectively 163840 unknowns were run on a MacBookPro machine with 8Gb of memory. An important drawback of the elimination algorithm described in Section <ref> is its sequential nature.


We organize the presentation of the numerical experiments into four categories.  First we treat the case of one interface Γ_0 which is important in its own right for several applications; then we present results for large numbers of layers; we continue with results involving periodic layers that contain periodic inclusions, as such configurations are relevant to photonics applications; and we conclude with three-dimensional results. We emphasize Wood frequencies to demonstrate the versatility of the shifted Green function method, as such cases are computationally more challenging. We indicate in the headings of each Table whether the windowed Green function G_k^q,A or the shifted Green function G_k,h^q,j,A were used in the numerical experiments. 

Results for periodic transmission problems with one grating interface. We start with an illustration in Table <ref> of the high-order accuracy that can be achieved by the DDM solver in the case of one periodic interface/grating Γ_0 given by the graph of the 2π periodic function
x_2=H/2cosx_1 for two values of the height H=0.6 and H=2; in this example we took k_0=4.1 and k_1=16.1 and thus the period of the interface Γ_0 is respectively 4.1 and 16.1 wavelength across.




We continue in Table <ref> and Table <ref> with numerical results concerning the convergence of the DDM algorithm in the case of one grating profile Γ_0 given by the graph of the 2π periodic function x_2=H/2cosx_1, H=0.6 for two values of the wavenumbers k_0 and k_1 that are simultaneously Wood frequencies.  In Table <ref> we consider medium frequencies—the set U^+:={r∈ℤ:β_0,r≥ 0} consists of 17 propagating modes and the set U^-:={r∈ℤ:β_N+1,r≥ 0} consists of 65 propagating modes; in Table <ref> higher frequencies—the set U^+:={r∈ℤ:β_0,r≥ 0} consists of 31 propagating modes and the set U^-:={r∈ℤ:β_N+1,r≥ 0} consists of 121 propagating modes. In the configuration in Table <ref> the period of the interface Γ_0 is 8 and 32 wavelengths across on each side, whereas in the configuration in Table <ref> the period of the interface Γ_0 is 15 and 60 wavelengths across on each side.






In Table <ref> we present results concerning a deep grating profile Γ_0 given by the graph of the 2π periodic function x_2=H/2cosx_1, H=2 for two values of the wavenumbers k_0 and k_1 that are simultaneously Wood frequencies, that is k_0=4 and k_1=16. We note that high-accuracy results can be achieved in this case by increasing the size of discretization. 


We conclude the numerical results in this part with a case in Table <ref> with a transmission experiment involving one interface of material discontinuity Γ_0 and two wavenumbers such that one of them is not a Wood frequency while the other is a Wood frequency. 


Multiple layers. The next set of results concern transmission experiments involving multiple periodic layers. In the example that follows we consider the first profile Γ_0 described by (a) x_2=F_0(x_1), F_0(x_1):=H/2cosx_1 and (b) x_2=F_0(x_1), F_0(x_1):=π H(0.4cos(x_1)-0.2cos(2x_1)+0.4cos(3x_1)) and the subsequent profiles Γ_ℓ being simple down shifted versions of the first profile, that is the grating Γ_ℓ is given by x_2=F_ℓ(x_1), F_ℓ(x_1):=-ℓ L + F_0(x_1), 0≤ℓ≤ N. The first set of results in Table <ref> concerns the convergence of the DDM transmission solver in layered configurations consisting of 4 layers (that is N=2) separated by interfaces Γ_ℓ,0≤ℓ≤ 2, when each wavenumber k_ℓ, 0≤ℓ≤ 3 is a Wood frequency. 



In the next set of results in Table <ref> we present numerical experiments concerning periodic configurations that involve large numbers of layers (i.e 10, 20, and 40 layers) and associated wavenumbers that are all Wood frequencies. We used shifted Green functions with a number j=3 of shifts, as this choice leads to small energy conservation defects. In the examples presented in Table <ref>, the interface Γ_j is j+1 wavelengths across, leading thus to problems that overall are 55, 210, and respectively 820 wavelengths in size.




In applications that involve high-contrast layer media, using non-conforming DDM discretizations leads to more efficient solvers. We present experiments in Table <ref> concerning configurations consisting of layers with alternating high-contrast material properties. In such settings it is natural to use coarser discretizations to compute the RtR maps corresponding to layers with smaller wavenumbers as well as restriction/interpolation Fourier matrices to match non-conforming interface Robin data. 



In the last set of results in this part we present in Table <ref> numerical experiments concerning very large numbers of layers and associated wavenumbers that are not Wood frequencies. In such cases, the Schur complement elimination algorithm for the solution of the DDM algorithm reduces the memory requirements via the forward/backward domain sweep. In the examples presented in Table <ref>, the interface Γ_j is approximately j+1 wavelengths across, leading thus to problems that overall are about 820 and respectively 3240 wavelengths in size. We mention that all the matrices 𝒟_j,M (see Section <ref>) that need be inverted in the Schur complement solution of the problems presented in Table <ref> are well conditioned, with condition numbers in the interval [19,400]; also, the condition numbers of the matrices 𝒟_j,M grow with the layer index j. However, the condition numbers of the matrices 𝒟_j,M grow with the size M of the discretization cf. Theorem <ref> and Theorem <ref>, yet not drastically. On the other hand, if the solution of  the DDM corresponding to transmission problems with N=9,19,29 layers and wavenumbers k_ℓ=ℓ+1.2, 0≤ℓ≤ N+1 were attempted via iterative solvers, the numbers of GMRES iterations required to reach a relative residual of 10^-4 are 270,718,1292 for the grating interfaces F_ℓ(x_1)=-ℓ L+H/2cosx_1, H=2, L=0.3, 0≤ℓ≤ N, and respectively 318,814,1498 for the grating interfaces F_ℓ(x_1)=-ℓ L+π H(0.4cos(x_1)-0.2cos(2x_1)+0.4cos(3x_1)), H=1, L=0.3, 0≤ℓ≤ N when the discretized RtR matrices are of size M^2, M=192 and η=1. Despite the large numbers of GMRES iterations required for DDM convergence in the experiments above, the condition numbers of the DDM matrices (which can be built in these cases) are reasonable: 99.1, 317.6, 671.4 in the first case, and respectively 119.2, 427.1, 947.6 in the second case.






Inclusions in periodic layers. In the last part of the numerical results section we present numerical experiments concerning periodic layers with embedded perfectly reflecting inclusions as presented in Figure <ref>. Specifically, we consider perfectly reflecting inclusions D whose boundary ∂ D is a smooth closed curve given in parametric form ∂ D:={(x_1(t),x_2(t)): x_1(t)=3.3 + r(t)cost, x_2(t)=-1+r(t)sint, r(t)=0.8 +0.4cos3t, 0≤ t≤ 2π}. The inclusions D are embedded periodically in a layered structure whose top boundary is explicitly given by either Γ_t:={(x_1,F_0(x_1)): F_0(x_1)=0.6+H/2cosx_1, H=0.6} or the flat interface Γ_t:={(x_1,0.6)}, and the bottom interface is given by shifting the top interface 3 units down the x_2 axis. The first set of numerical results presented in Table <ref> concerns such configurations in the case when all three wavenumbers k_0,k_1 and k_2 are Wood frequencies.





Finally, we conclude with a numerical experiment in Table <ref> concerning scattering by a periodic arrays of perfectly reflecting cylinders D described above at Wood frequencies. We treat this case via fictitious periodic layers bounded by flat interfaces as in the right panel of Figure <ref>, for which we apply the DDM transmission algorithm with k_0=k_1=k_2.  The DDM approach for the solution of scattering by array of cylinders at Wood frequencies requires discretization of fictitious boundaries and as such is not as computationally efficient as alternative approaches that rely on Sherman-Morrison formula <cit.>. Nevertheless, we believe that the DDM approach is more straightforward and more modular in the sense that it consists of several black-box solvers that can be easily assembled to treat complex periodic-layered cases.   



3D results. We start our presentation of three-dimensional results with the case of scalar Helmholtz transmission problems featuring one interface of material discontinuity given by the doubly periodic grating surface x_3=f(x_1,x_2)=1/2cos(2 π x_1)cos(2 π x_2).  For all the numerical experiments presented here, a single patch was used to represent the doubly periodic surfaces. We illustrate in Table <ref> the high-order convergence achieved by our DDM solvers in the case of transmission problems involving Wood frequencies in both semi-infinite domains. For the grating considered, under normal incidence, the first three Wood frequencies occur at 2π, 2√(2)π, and 4π respectively.



We continue in Table <ref> with an illustration of the accuracy achieved by our DDM solvers in the case of three layers separated by two doubly periodic gratings and wavenumber configurations that involve Wood frequencies. Finally, we conclude with an illustration in Table <ref> of the ability of the Schur complement DDM solvers to handle very large numbers of layers in three dimensions that require large discretizations; for instance, the largest problem considered in Table <ref> involves a periodic layered configuration consisting of 80 doubly periodic interfaces of material discontinuity, spanning about 160 wavelengths, whose DDM discretization required 163840 unknowns.  We mention that the condition numbers of the matrices 𝒟_j,M (see Section <ref>) that need be inverted in the Schur complement solution of the problems presented in Table <ref> belong to the interval [10^2,3.3× 10^3].






While the grating profiles considered in this paper are relatively simple, qualitatively similar results can be obtained for more geometrically complex profiles. We mention that extensions to grating profiles that involve corners is straightforward in two dimensions; graded meshes and weighted versions of RtR maps are needed to treat those cases <cit.>. Extensions to three-dimensional gratings with edges and corners can be done by applying existing technology <cit.>.  The results presented in this section were produced on a MacBookPro with a 2.7 GHz Intel processor and 8Gb of RAM based on a MATLAB implementation of the DDM algorithm. We did not strive to optimize the code in order to harness the best sequential computational performance; this is certainly possible, and it has been done in <cit.> via fast methods based on equivalent sources. Instead, we wanted to illustrate that, within a DDM approach, the use of windowed Green function method combined with shifted Green functions leads to a computational method for scattering by periodic layered media that is accurate and robust at all frequencies, including the challenging Wood frequencies.

  


§ CONCLUSIONS


We presented analysis and numerical experiments concerning  boundary-integral operators-based DDM for two and three-dimensional periodic layered media scalar scattering problems. We have shown that the RtR maps that are needed by DDM can be computed in a robust manner at all frequencies, including Wood frequencies. The Wood frequencies configurations were treated via boundary-integral operators that incorporate shifted quasi-periodic Green functions that converge at Wood frequencies. The tridiagonal DDM linear system associated with transmission problems in periodic layered media was solved via recursive Schur complements resulting in a computational cost that is linear in the number of layers. Extensions to full three-dimensional electromagnetic configurations are straightforward. We are currently investigating the design of DDM with quasi-optimal transmission conditions for the solution of transmission problems in periodic layered media.



§ ACKNOWLEDGMENTS

Stephen Shipman acknowledges support from NSF through contract DMS-0807325. Catalin Turc acknowledges support from NSF through contract DMS-1614270. Stephanos Venakides acknowledges support from NSF through contract DMS-1211638.



§ APPENDIX



  Under the assumptions that (1) the wavenumbers k_j are such that 0≤ k_j≤ k_j+1 for all 0≤ j≤ N and (2) the coefficients γ_j=1 for all 0≤ j≤ N+1, the system of Helmholtz equations (<ref>) has a unique solution when the functions F_j are C^2.


  Clearly, the uniqueness of solutions amounts to showing that when the incident field is zero, the only solution of the transmission equations (<ref>) is the trivial solution. The key ingredient in the proof is the application of Green's identities. Let us choose h>maxF_0 and define the domain Ω_0,h^per:={(x_1,x_2)∈Ω_0^per: F_0(x_1)≤ x_2≤ h}. A simple application of Green's identities leads to 

    ∫_Ω_0,h^per(|∇ u_0|^2-k_0^2|u_0|^2)dx=∫_Γ_0∂_n_0u_0 u_0 ds + ∫_Γ_0,h∂_x_2u_0 u_0 dx_1

where Γ_0,h:={(x_1,x_2): 0≤ x_1≤ d, x_2=h}. Note that the integrals over the vertical lines vanish due to the quasi-periodicity of the field. Taking into account the fact that u_0 is radiating, we can express u_0 on the line segment Γ_0,h in terms of the following Rayleigh series

    u_0(x_1,h)=∑_r∈ℤ C_r^+e^iα_r x_1+iβ_0,r h

 from which it follows immediately that

    lim_h→∞∫_Γ_0,h∂_x_2u_0 u_0 dx_1=id∑_r∈ℤ, β_0,r>0β_0,r|C_r^+|^2.

Hence, we get

    ∫_Ω_0^per(|∇ u_0|^2-k_0^2|u_0|^2)dx=∫_Γ_0∂_n_0u_0 u_0 ds + id∑_r∈ℤ, β_0,r>0β_0,r|C_r^+|^2.

Taking the imaginary part of the equation above we arrive at

    ∫_Γ_0(∂_n_0u_0) u_0 ds =-d∑_r∈ℤ, β_0,r>0β_0,r|C_r^+|^2.

On the other hand, application of the Green identities in the layers Ω_j^per,1≤ j≤ N that have a finite width in the x_2 leads to

    ∫_Ω_j^per(|∇ u_j|^2-k_j^2|u_j|^2)dx=∫_Γ_j-1∂_n_ju_j u_j ds + ∫_Γ_j∂_n_ju_j u_j ds.

Taking the imaginary part in the equation above we obtain

    ∫_Γ_j-1(∂_n_ju_j) u_j ds = -∫_Γ_j(∂_n_ju_j) u_j ds.

Applying the same arguments that led to the derivation of equation (<ref>) in the case of the semi-infinite layer Ω_N+1 we obtain

    ∫_Γ_N(∂_n_N+1u_N+1) u_N+1 ds =-d∑_r∈ℤ, β_N+1,r>0β_N+1,r|C_r^-|^2.

Adding the left-hand sides as well as the right  hand sides of equations (<ref>), (<ref>), and (<ref>) and taking into account the continuity conditions in the transmission system (<ref>) we obtain

    ∑_r∈ℤ, β_0,r>0β_0,r|C_r^+|^2+ ∑_r∈ℤ, β_N+1,r>0β_N+1,r|C_r^-|^2=0.

The last equation implies that the Rayleigh coefficients of the propagating modes corresponding to u_0 and u_N+1 are all equal to zero, that is  C_r^+=0 for all r such that β_0,r>0 as well as C_r^-=0 for all r such that β_N+1,r>0. We have then

    u_0(x_1,h)=∑_r∈ℤ,β_0,r=0 C_r^+e^iα_r x_1+iβ_0,r h + ∑_r∈ℤ,β_0,r>0 C_r^+e^iα_r x_1+iβ_0,r h.

We define v_0:=∂_x_2u_0 in the domain Ω_0 and we apply Green's third identity to the functions u_0 and v_0 in the domain Ω_0,h^per and take h→∞ to obtain


    ∫_Γ_0∂_n_0u_0 v_0 ds=∫_Γ_0u_0 ∂_n_0v_0 ds.

We similarly define v_j:=∂_x_2u_j in the domains Ω_j for 1≤ j≤ N and we obtain in a similar manner

    ∫_Γ_j-1∂_n_ju_j v_j ds+∫_Γ_j∂_n_ju_j v_j ds=∫_Γ_j-1 u_j ∂_n_jv_j ds+∫_Γ_ju_j ∂_n_jv_j ds, 1≤ j≤ N

as well as

    ∫_Γ_N+1∂_n_N+1u_N+1 v_N+1 ds=∫_Γ_N+1 u_N+1 ∂_n_N+1v_N+1 ds.

Now, using the continuity of the normal derivatives and the tangential derivatives of the fields u_j across interfaces Γ_j, it follows immediately that the quantities v_j are continuous across the interfaces Γ_j. Also, as shown in <cit.>, we have that

    ∂_n_jv_j+∂_n_j+1v_j+1=(k_j+1^2-k_j^2)n_j,x_2u_j   on Γ_j, 1≤ j≤ N,

where n_j,x_2 denotes the component of the normal n_j along the x_2 axis. Taking these last two facts into account, we add the left-hand sides and right hand sides of equations (<ref>), (<ref>), and (<ref>) and we get

    ∑_j=0^N∫_Γ_j(∂_n_ju_j+∂_n_j+1u_j+1)v_j ds=∑_j=0^N∫_Γ_ju_j(∂_n_jv_j + ∂_n_j+1v_j+1) ds

Now, given that ∂_n_ju_j+∂_n_j+1u_j+1=0 on Γ_j for 0≤ j≤ N, and taking into account the continuity condition in equation (<ref>), we obtain

    ∑_j=0^N(k_j+1^2-k_j^2)∫_Γ_jn_j,x_2|u_j|^2ds=0.

Now, given the assumption (1) and the fact the normals n_j are chosen to point to the exterior of the domains Ω_j and hence n_j,x_2<0 for all 0≤ j≤ N, it follows in particular that u_0=0 on Γ_0. In the light of this fact, we revisit formula (<ref>) and we get that

    ∫_Γ_0∂_n_0u_0 v_0 ds =0.

Given that u_0=0 on Γ_0, it follows that its tangential derivative is also equal to zero on Γ_0. Denoting by w_0:=∂_x_1u_0 in Ω_0, the latter fact translates into n_0,x_2w_0=n_0,x_1v_0 on Γ_0. Since ∂_n_0u_0 v_0=n_0,x_1w_0 v_0+n_0,x_2|v_0|^2=n_0,x_2|w_0|^2+n_0,x_2|v_0|^2, we get 

    ∫_Γ_0n_0,x_2(|w_0|^2+|v_0|^2)ds=0.

Hence, u_0=0, ∂_n_0u_0=0 on Γ_0, which implies that u_0=0 in Ω_0 by Holmgren's theorem. Now the use of continuity conditions across interfaces Γ_j for 1≤ j≤ N and Holmgren's theorem leads to the conclusion of the theorem.




§ APPENDIX

We devote this Appendix to proving the following result.

Theorem 3.4 
  Under the assumption that F_0 is C^2 and that k_0 is a Wood frequency, the operator
  
    𝒜_0,h:=1/2I+(K_Γ_0,k_0,h^q,j)^⊤-Z_0S_Γ_0,k_0,h^q,j, j≥ 1, 𝒜_0,h:L^2_per(Γ_0)→ l^2_per(Γ_0)

  is invertible with continuous inverse for all but a discrete set of values of the shift h>0.



  Since for a given 𝐲∈Γ_0, the kernels G_k_0,h^q,j(·-𝐲) have the same singularity on Γ_0 (that is at 𝐱=𝐲) as the kernels G_k_0^q(·-𝐲), then both operators (K_Γ_0,k_0,h^q,j)^⊤:L^2_per(Γ_0)→ L^2_per(Γ_0) and S_Γ_0,k_0,h^q,j:L^2_per(Γ_0)→ L^2_per(Γ_0) are compact. Thus, the conclusion of the Theorem follows once we establish the injectivity of the operator 𝒜_0,h. Let φ∈ Ker(𝒜_0,h) and define
  
    w_0:=SL^q,j_k_0,hφ   in ℝ^2∖Γ_0.

The function w_0 is a radiating α-quasi-periodic solution of the Helmholtz equation in the domain Ω_0 with zero Robin boundary values on Γ_0, that is

    ∂_n_0w_0-Z_0w_0=0   on  Γ_0.

It follows from Theorem <ref> that w_0=0 in Ω_0. It is straightforward to see that the shifted function G_k_0,h^q,j has the following frequency domain representation <cit.>:

    G_k_0,h^q,j(x_1,x_2;y_1,y_2)=i/2d∑_r∉ U e^iα_r(x_1-y_1)(1-e^iβ_r h)^j/β_re^iβ_r|x_2-y_2|+∑_r∈ U c_r e^iα_r(x_1-y_1),

with (x_1,x_2)∈Ω_0, (y_1,y_2)∈Ω_0 where U:={r∈ℤ: β_r=0}; clearly U is not empty since k_0 is a Wood frequency. Hence, the solution w_0 admits the following representation

    w_0(x_1,x_2)=∑_rw_0,r^+(x_2)e^iα_r x_1, x_2>maxF_0

where

    w_0,r^+(x_2)=e^iβ_r x_2(1-e^iβ_r h)^j/β_rb_r^+    r∉ U
    
    c_rb_r^+    r∈ U

with

    b_r^+=i/2d∫_Γ_0e^-iα_ry_1e^-iβ_ry_2φ(y_1,y_2)ds(y_1,y_2).

Assuming that the shift h is chosen such that 1-e^iβ_r h≠ 0 for all r∈ U, it follows immediately that w_0=0 in Ω_0 implies that all the coefficients b_r^+=0 for all r.

The key insight in the proof <cit.> is to introduce the following α-quasi-periodic Green function that is defined even at Wood frequencies

    B^q(x_1,x_2)=i/2d∑_r∉ Ue^iα_r x_1e^iβ_r|x_2|/β_r+i/2d∑_r∈ Ue^iα_r x_1i|x_2|.

The function B^q is not an outgoing Green function on account of the linear term |x_2|. Now define

    v(𝐱):=∫_Γ_0B^q(𝐱-𝐲)φ(𝐲)ds(𝐲), 𝐱∉Γ_0.

Then the function v admits the representation

    v(x_1,x_2)=∑_rv_r^+(x_2)e^iα_r x_1, x_2>maxF_0

where

    v_r^+(x_2)=e^iβ_r x_2/β_rb_r^+    r∉ U
    
    ix_2b_r^+-ib_r'    r∈ U

with

    b_r'=i/2d∫_Γ_0e^-iα_ry_1y_2φ(y_1,y_2)ds(y_1,y_2), r∈ U.

Given that we established the fact that b_r^+=0 for all r, we obtain v(x_1,x_2)=∑_r∈ U(-ib_r')e^-iα_rx_1 for x_2>maxF_0. Certainly, v(𝐱) is a solution of the Helmholtz equation for 𝐱∉Γ_0, and v(𝐱) is independent of x_2 for x_2>maxF_0. But v(𝐱) is real analytic for 𝐱∉Γ_0, so it follows from analytic continuation that v(𝐱) is actually independent of x_2 everywhere in Ω_0. Thus, we have

    v(x_1,x_2)=∑_r∈ U(-ib_r')e^-iα_rx_1   for x_2≥ F_0(x_1).

We study next the behavior of the function v(x_1,x_2) for x_2<minF_0. We get immediately that v also admits the representation

    v(x_1,x_2)=∑_rv_r^-(x_2)e^iα_r x_1, x_2<minF_0

where

    v_r^-(x_2)=e^-iβ_r x_2/β_rb_r^-    r∉ U
    
    -ix_2b_r^-+ib_r'    r∈ U

with

    b_r^-=i/2d∫_Γ_0e^-iα_ry_1e^iβ_ry_2φ(y_1,y_2)ds(y_1,y_2), r∈ℤ.

Comparing the definitions of the coefficients b_r^+ and b_r^- we see that they differ in general for r∉ U. However, and most importantly, we have that b_r^-=b_r^+ for r∈ U since β_r=0 for r∈ U. In conclusion, we have

    v_r^-(x_2)=e^-iβ_r x_2/β_rb_r^-    r∉ U
    
    ib_r'    r∈ U.

The very last fact we established implies that v is actually a radiating α-quasi-periodic solution of the Helmholtz equation in the domain Ω_0^-:{(x_1,x_2):x_2≤ F_0(x_1)}. In the last step of the proof we define

    ṽ(x_1,x_2):=v(x_1,x_2)-∑_r∈ U(-ib_r')e^-iα_rx_1,   for (x_1,x_2)∈ℝ^2

which is a α-quasi-periodic solution of the Helmholtz equation that satisfies the radiation condition as x_2→∞ as well as x_2→-∞, while vanishing in Ω_0. Clearly, ṽ is a radiating α-quasi-periodic solution of the Helmholtz equation in the domain Ω_0^- that vanishes on Γ_0. Consequently, ṽ=0 in Ω_0 <cit.>. Using the jump conditions of the normal derivative of the normal derivative of the single-layer potentials, we get that φ=0 on Γ_0 which concludes the proof of the theorem. 
  
  














































































































































  





































