Fast Solver for Quasi-Periodic 2D-Helmholtz Scattering in Layered Media

We present a fast spectral Galerkin scheme for the discretization of boundary integral equations arising from two-dimensional Helmholtz transmission problems in multi-layered periodic structures or gratings. Employing suitably parametrized Fourier basis and excluding Rayleigh-Wood anomalies, we rigorously establish the well-posedness of both continuous and discrete problems, and prove super-algebraic error convergence rates for the proposed scheme. Through several numerical examples, we confirm our findings and show performances competitive to those attained via Nystr\"om methods.


Introduction
A vast number of scientific and engineering applications rely on harnessing acoustic and electromagnetic wave diffraction by periodic and/or multilayered domains. Current highly demanding operation conditions for such devices require solving thousands of specific settings for design optimization or the quantification of shape or parameter uncertainties in the relevant quantities of interest, challenging the scientific computing community to continuously develop ever more efficient, fast and robust solvers (cf. [8,18,31,40,39] and references therein). Assuming impinging time-harmonic plane waves, scattered and transmitted fields have been solved by a myriad of mathematical formulations and associated solution schemes. These range from volume variational formulations to various boundary integral representations and equations (cf. [2,4,9,11,21,33]), pure or coupled implementations of finite and boundary element methods (cf. [3,4,22,34,40] or [36,Chapter 5]) and Nyström methods [15,17,19,24,30].
In this work, we build upon our theoretical review given in [5] and present a spectral Galerkin method for solving second-kind direct boundary integral equations (BIEs) for the Helmholtz transmission problem for two-dimensional, periodic multi-layered gratings with smooth interfaces. Contrary to the low-order local basis functions used in the standard boundary element method, spectral bases are composed of high-order polynomials whose support lie on the whole scatterer boundary or on large portions of it. Successfully employed on twoand three-dimensional scattering problems [26,25,23], the main advantage of a spectral discretization is the ability to converge at a super-algebraic rate whenever solutions are smooth enough. Hence, our proposed method can in practice compete with Nyström methods while simultaneously inheriting all of the theoretical aspects of classical Galerkin methods.
In two dimensions, spectral methods are closely related to the theory of periodic pseudodifferential operators [37], since the discretization through spectral elements can be interpreted as a truncation of the associated Fourier series where the action of the operators is well understood. We show that wave scattering by periodic domains is closely connected to the bounded domain case, making it possible to reuse almost all the pseudo-differential operator theory for our analysis. Key to our analysis are the results in [34,41,22] regarding the unique solvability and eigenvalues of the associated volume problem. From here, we deduce that our BIE is uniquely solvable except at a countable set of wavenumbers composed of Rayleigh-Wood frequencies-wavenumbers for which the sum defining the quasi-periodic Green's function is not convergent-and of eigenvalues of the Helmholtz transmission problem. Mindless of the several remedies developed to tackle Rayleigh-Wood anomalies through BIEs [19,13,16,14], we choose to avoid them as they are not captured by our previous analysis in [5].
Our discretization method employs a quasi-periodic basis so that techniques forcing the quasi-periodicity of the discrete solutions are not necessary (cf. [24,46]). Instead, an accurate approximation of the quasi-periodic Green's function is required in order to extract its Fourier coefficients through the fast Fourier transform (FFT). Moreover, we prove that the chosen discretization basis enjoys a super-algebraic convergence rate on the degrees of freedom, which we then confirm through numerical experiments. In [35], a similar quasi-periodic exponential basis was employed to approximate solutions of a volume integral formulation.
The article is structured as follows. Section 2 presents the notation used throughout as well as the required quasi-periodic Sobolev spaces setting following [5]. In Section 3 we state the Helmholtz transmission problem for a multi-layered grating and study its solvability. Section 4 is concerned with the properties of quasi-periodic boundary integral operators (BIOs) along with an existence and uniqueness result for our BIEs. Section 5 provides rigorous error convergence rates of the spectral method and briefly describes the numerical algorithm used to compute the matrix entries associated with each integral operator. Numerical results are discussed in Section 6, followed by concluding remarks on Section 7.

General Notation
We denote the imaginary unit ı. Boldface symbols will denote vectorial quantities and will use greek and roman letters for data over boundaries and volume, respectively. Canonical vectors in R 2 are denoted e 1 , e 2 respectively. Also, we make use of the symbols À, Á andto avoid specifying constants irrelevant for the corresponding analysis.
Let H be a given Banach space. We shall denote its norm as ¨ H and its dual space by H 1 (set of antilinear functionals over H) with dual product denoted by x¨,¨y. If H is a Hilbert space, the inner product between two of its elements, x and y, is denoted as px, yq H . Moreover, if H is a Hilbert space over the complex field, the inner product will be understood in the anti-linear sense.
For an open domain ΩĂR 2 , its boundary shall be denoted as BΩ. Moreover, for any OĂR 2 such that Ω Ď O, we introduce the closure of Ω relative to O as Ω O :" ΩXO and the boundary of Ω relative to O as B O Ω :" Ω O zΩ. For n P N 0 :" N Y t0u, we denote by C n pΩq the set of scalar functions over Ω with complex values and continuous derivatives up to order n. C 8 pΩq refers to the space of functions with infinite continuous derivatives over Ω. We shall also make use of the following subset of C 8 pΩq: DpΩq :" tu P C 8 pΩq : supp u ĂĂ Ωu.
The space of p-integrable functions (for p ě 1) with complex values over Ω is denoted as L p pΩq.
We say that a one-dimensional Jordan curve Γ is of class C r,1 , for r P N 0 , if it may be parametrized by a function z : p0, 2πq Ñ Γ which has r Lipschitz-continuous derivatives and a non-vanishing tangential vector. The first derivative of the parametrization is denoted as 9 z. Moreover, we say Γ is of class C 8 if it is of class C r,1 for every r P N 0 (we will also use the notation C 8,1 to refer to the same class).
Throughout the following sections, we will consider periodic geometries along e 1 with a fixed period of 2π. Moreover, we say that a continuous function f is a θ-quasi-periodic function if, where the quasi-periodic shift θ is always assumed to be in r0, 1q. Finally, we define the canonic periodic cell on R 2 as G :" p0, 2πqˆR.

Quasi-periodic Sobolev Spaces
We denote by D θ pR 2 q the space of θ-quasi-periodic functions in C 8 pR 2 q that vanish for large |x 2 |, and denote by D 1 θ pR 2 q the space of θ-quasi-periodic distributions, which can be seen as the dual space of D θ pR 2 q (cf. [5,Proposition 2.4]). For G as before, we introduce D θ pGq the space of restrictions to G of elements in D θ pR 2 q. Moreover, for any open domain Ω Ă G we define D θ pΩq as the set of elements of D θ pGq with compact support on Ω and D 1 θ pΩq as the space of elements of D 1 θ pGq restricted to D θ pΩq. In what follows, for all j P Z we define j θ :" j`θ. Proposition 2.1 (Proposition 2.6 in [5]). Every u P D θ pR 2 q can be represented as a Fourier series, i.e.
upxq " ÿ jPZ u j px 2 qe ıj θ x 1 with u j px 2 q :" 1 2π so that u j P DpRq. On the other hand, every element F P D 1 θ pR 2 q can be identified with a formal Fourier series given by Let s P R. We define the θ-quasi-periodic Sobolev space of order s on G as follows, Figure 1: Example of a multi-layered grating. G is highlighted and the dotted lines represent its boundaries at 0 and 2π. Figure 1). Moreover, while all domains t r Ω i u M i"0 are unbounded along e 1 -due to their periodicity-only two of them, namely r Ω 0 and r Ω M , are unbounded in the second spatial dimension (along e 2 ). The restrictions of the aforementioned domains and surfaces to the periodic cell G are denoted by: Additionally, we fix H ą 0 so that holds. We will assume that the interfaces Γ i , i P t1, . . . , M u are all Jordan curves of class C 8 . Furthermore, for each i P t1, . . . , M u, the exterior and interior trace operators on Γ i are understood as and the normal vector on Γ i is chosen to point towards Ω i´1 .

Helmholtz transmission problem on periodic media
For a time-dependence e´ı ωt for some frequency ω ą 0, let the grating described in the previous subsection be illuminated by an incident plane wave, u (inc) pxq :" e ık 0¨x " e ıpk 0,1 x 1`k0,2 x 2 q , where k 0 " pk 0,1 , k 0,2 q P R 2 . Furthermore, we denote k 0 :" |k 0 |. For i " 0, . . . , M , the material filling each domain Ω i is assumed to be homogeneous and isotropic with refraction index η i -we assume η 0 " 1-and wavenumber k i :" where c i is the wave speed in Ω i . Throughout this section, we fix θ as the unique real in r0, 1q such that θ " k 0,1`n for some integer n and, for all j P Z, we define in Ω for i P t1, . . . , M u, " Equation (3b) represents the continuity of Dirichlet and Neumann traces across each interface. This condition can be generalized to include different transmission coefficients without much effort. The last two conditions, namely (3c) and (3d), correspond to radiation conditions for u 0 and u m , also known as the Rayleigh-Bloch expansions (cf. [34] for a detailed discussion), where tu p0q j u jPZ and tu pM q j u jPZ are the corresponding Rayleigh coefficients. Through an analogous analysis to that presented in [22,Section 3], one finds that-for a fixed choice of geometries tΓ i u M i"1 and refraction indices tη i u M i"1 -Problem 3.1 has a unique solution for all but a countable number of wavenumbers k 0 as all wavenumbers k i for i P t1, . . . , M u depend on k 0 . Assumption 3.2. The wavenumber k 0 is such that Problem 3.1 has a unique solution.
We shall make no further analysis of the volume problem as stated above, and limit ourselves to [7,10,27,34,41,45,22] and references therein for more detailed analyses of the radiation condition of similar problems.

Boundary integral equations
Following our previous work [5], we introduce the quasi-periodic Green's function and recall some relevant properties. We then define the quasi-periodic single and double layer potentials and BIOs spanning from taking their respective traces on the periodic boundaries tΓ i u M i"1 . To conclude this section, we present an integral representation for the fields tu i u M i"0 and a proof of unisolvency for the corresponding BIE. As before, θ will denote the quasi-periodic shift, which is assumed to be in r0, 1q.

Quasi-Periodic Fundamental Solution
Consider a positive wavenumber k P R, we recall the definition of the Rayleigh-Wood frequencies.
Definition 4.1. We say k ą 0 is a Rayleigh-Wood frequency, if there is j P Z, such that where θ is the previously fixed quasi-periodic shift.
These frequencies correspond to values where the quasi-periodic Green's function can not be represented in a traditional manner. While a number of alternatives have been developed to circumvent this issue (e.g., [14,16,19]) their analysis is out of the scope of our current work.
Hence, in what follows, we will work under the following assumption over the wavenumber k.
Assumption 4.2. The wavenumber k ą 0 is not a Rayleigh-Wood frequency for the given θ P r0, 1q .
Under Assumption 4.2 we can define the θ-quasi-periodic Green's function as (cf. [34,29] and references therein) for all x, y in R 2 such that x´y ‰ 2πne 1 for all n P Z, wherein G k px, yq is the fundamental solution for the Helmholtz equation with wavenumber k, namely, where H p1q 0 p¨q denotes the zeroth-order first kind Hankel function. Moreover, the quasi-periodic Green's function is a fundamental solution of the Helmholtz equation in the following sense: for all x P R 2 and satisfies the radiation condition specified in the preceding section (cf. [34, Proposition 3.1]).

Layer Potentials and Boundary Integral Operators
On this section, we will assume a given boundary Γ satisfying the following assumption.
Assumption 4.4. Given r P r0, 8s, the interface Γ is a Jordan curve of class C r,1 .
Moreover, we denote by Ω the part of G below Γ (see Figure 1). For φ P D θ pΓq we define the single and double layer potentials as where γ n,y denotes the interior (with respect to Ω) Neumann trace operator acting on functions with argument y.
Lemma 4.5 (Theorems 4.7 and 4.10 in [5]). Let k and Γ be as in Assumptions 4.2 and 4.4 with r ě 0, respectively. Then, the single and double layer potentials can be extended as continuous operators acting on Sobolev spaces as follows, We then define BIOs by taking traces of the layer potentials as follows Moreover, due to the jump properties of the layer potentials [5,Lemma 4.11], the following relations hold: Remark 4.6. When considering interior and exterior traces acting on layer potentials, note that the normal vector on Γ is to be fixed so that the only difference between exterior and interior traces is the direction from which we approach Γ. Additionally, note that, having fixed the normal vector to Γ, the choice of trace taken in the definition of W k θ,Γ is arbitrary and makes no difference.

Compacteness Properties
Until this point, we have established continuity properties of the four BIOs defined in (7). However, the BIEs we consider in the coming section require the subtraction of two instances of the same BIO with different wavenumbers. This will require a number of results from pseudodifferential operator theory [37] as well as a version of the Rellich theorem on quasi-periodic Sobolev spaces on boundaries. After our analysis, we will see that the difference between any two of the operators in (7)-with different wavenumbers-will result in a compact operator.
Proof. Follows directly from the definition of the quasi-periodic spaces and the result for standard Sobolev spaces (see [28,Theorem 8.3]).
Remark 4.9. No smoothness assumptions are needed for the proof of the previous theorem. Thus, it can be extended to Lipchitz boundaries for any pair of real numbers s 1 , s 2 ă 1, and potentially less regular cases if we restrict s 1 , s 2 to be non-negative.
Theorem 4.10 (Theorem 6.1.1 in [37]). Let a : RˆR Ñ C be a bi-periodic function of class C 8 and S be a 2π´periodic distribution in R. Consider the following formal operator acting on a periodic smooth function u P C 8 pRq: Aupsq " where integration is to be understood as a duality pairing. Furthermore, let us assume the Fourier coefficients of S to behave as for some p P R. Then, for any s P R, A in (9) may be continuously extended as an operator mapping from H s r0, 2πs to H s´p r0, 2πs, i.e., We also recall a classical result from Fourier analysis (c.f. [43]).
Lemma 4.11. Let m P N, f : R Ñ C be a periodic C m -class function such that its distributional derivative of order m`1 belongs to L 1 pp0, 2πqq. Then, its Fourier coefficients tf n u nPZ are such that |f n | À |n|´m´1.
In order to employ Theorem 4.10 we will need to express the quasi-periodic BIOs in a convenient way: with periodic functions as kernels. Let k and Γ be as in Assumptions 4.2 and 4.4, respectively. We begin by considering a periodic version of the fundamental solution in (5) and its derivatives on Γ as which may be expressed as with J k θ ps, tq :" e´ı θps´tq where J 0 p¨q is the zeroth-first kind Bessel function, P p0, 2πq and χ p¨q is a smooth function satisfying and R k θ ps, tq " p G k θ ps, tq´Spt´sqJ k θ ps, tq. Using known expansions of the Hankel functions (see [1, 9.1.12-9.1.13]) one can check that R k θ belongs to C 8 pRˆRq.
Before we proceed any further, it is necessary to introduce a second wavenumber. We will denote r k ą 0 a wavenumber (not necessarily different from k) that also satisfies Assumption 4.2.  (7) and where we have dropped the Γ subscript for brevity. Both operators may be considered as pseudo-differential operators of order´1, whence as a bounded linear operator for every s P R.
Proof. That V k θ (and V r k θ ) may be extended as claimed follows directly from Theorem 4.10, the kernel representation (11) and the decay of the Fourier coefficients of Sptq in (12) Employing Lemma 4.11, Theorem 4.10 and [1, Equation 9.1.13] we see that the second term of the right-hand side of (13) gives rise to a bounded operator from H s r0, 2πs to H s`p r0, 2πs for any p ą 0. On the other hand, the first term in the right-hand side of (13) may be decomposed as One can see (cf. [1,Equation 9.1.12]) that the term pJ k θ ps, tq´J r k θ ps, tqq |sinpt´sq|´2 belongs to C 8 pRˆRq, whereas the term |sinpt´sq| 2 Spt´sq give rise to an operator of order´3. In fact, its Fourier transform iś 1 2π where the last equality follows from [37, Example 5.6.1]. Finally, define We may now bound the last term in (14) by Theorem 4.10: The proof is completed by the density of D θ pΓq in the corresponding Sobolev space.
For the hyper-singular BIO, a similar result requires a technical lemma. To this end, let us define the tangential curl operator: for any ϕ P D θ pΓq and where z is a suitable (arbitrary) parametrization of Γ.
Lemma 4.13. Let k and Γ satisfy Assumptions 4.2 and 4.4 for r " 0, respectively, and let λ and ϕ belong to D θ pΓq. Then, where x¨,¨y Γ represents the duality product between H s θ pΓq and H´s θ pΓq for any s ą 0 and q V k θ is the extension by density of the operator given by Proof. Notice that for λ, ϕ in D θ pΓq, it holds that where the border terms cancel each other out due to the quasi-periodicity of λ and ϕ. Hence, the result for quasi-periodic functions follows verbatim from the standard case (see, for instance, [42,Theorem 6.15]). Proof. Let λ, ϕ in D θ pΓq. By Lemma 4.13, we have that Using Proposition 4.12, one obtainšˇˇA Where the inequality for the second term of the right-hand side is obtained using that both p q V k θ , q V r k θ q are operators of order´1 (this follow from Theorem 4.10 and [37, Example 5.6.1]. Then, since the curl Γ operator is a first-order differential operator, it holds that and the result follows by a duality argument and recalling the density of D θ pΓq in our quasiperiodic Sobolev spaces.
are compact for every ą 0.
Proof. The result is direct from the mapping properties shown and Theorem 4.8.
Lastly, we require the compactness of the operator resulting form taking traces of the single and double layer operators acting on densities lying on a boundary Γ 1 over another x 1 -periodic curve, say Γ 2 , that does not intersect with Γ 1 . Let us denote by γ 2 d , γ 2 n Dirichlet and Neumann traces over Γ 2 , respectively. Then, by an application of Lemma 4.11, Theorem 4.8 and Theorem 4.10, we obtain the following result.
Proposition 4.17. Let k satisfy Assumption 4.2. If Γ 1 and Γ 2 are x 1 -periodic C 8 -Jordan curves then the application of the following traces to the layer potentials: are compact operators for any choice of s 1 , s 2 P R. The result holds regardless of the direction from which the traces are taken.
Remark 4.18. For the main results in this section, we have assumed the interfaces to be of class C 8 . While this simplifies the analysis, we could obtain similar results with less stringent regularity requirements. Consider k and r k satisfying Assumption 4.2 and Γ as in Assumption 4.4 with r P r1, 8q, and the weakly-singular operator V k θ (where we have omitted the Γ subindex momentarily). The expression in (11) still holds for the kernel of V k θ , but R k θ and J k θ would be only of class C r,1 , instead of arbitrarily smooth.

Boundary Integral Formulation
We recall the notation and geometry configuration introduced in Section 3, that is: 1. u pincq denotes a plane incident wave with wavenumber k 0 , which is assumed to be quasiperiodic with shift θ P r0, 1q.
2. tΓ i u M i"1 denotes a set of M P N non-intersecting C r,1 -Jordan curves, with r P r1, 8s, ordered downwards.
4. tηu M i"1 denotes a parameter set such that the wavenumber in Ω i is given by k i " η i k 0 for i P t1, . . . , M u. Assumption 4.20. For the given shift, θ, the wavenumber k 0 and the parameters tη i u M i"1 are such that neither k 0 nor the wavenumbers k i " η i k 0 are Rayleigh-Wood frequencies.
Following the notation of Problem 3.1, the scattered field-defined as the total field u ptotq minus the incident field u pincq -is written as Under Assumption 4.20, we make the following representation Ansatz for the scattered field: in Ω i , for i P t1, . . . , M´1u where, for each i P t1, . . . , M u, the boundary data λ i and µ i are assumed to belong to H s θ pΓ i q for some possibly different values of s P R, i.e., s may be different for each boundary datum. SL k j θ,Γ i and DL k j θ,Γ i are, respectively, the single and double layer potentials of wavenumber k j on Γ i .
As shorthand, in what follows, we denote, for each i P t1, . . . , M u, where λ i and µ i are defined over Γ i . For s 1 , s 2 P R, we define the Cartesian product spaces: where all of these spaces are equipped with their natural graph inner products. For each i P t1, . . . , M u let us define the following operators: corresponding to self-interactions between the potentials defined over each Γ i with themselves. Analogously, for i, j P t1, . . . , M u, we define the following operators: corresponding to interactions between potentials defined over Γ i with those defined over Γ j .
are compact operators for any s 1 , s 2 P R with s 2 ă s 1 ă s 2`2 . Furthermore, the crossinteraction operators (17) are compact for any choice of s 1 , s 2 P R.
Proof. The first result is directly found using Proposition 4.16, whereas the second one follows from Proposition 4.17.
With the above definitions and using the jump properties of the BIOs, it holds that for each i P t1, . . . , M u, where I i corresponds to the identity map over V s 1 ,s 2 θ,Γ j , with s 1 , s 2 P R. We now introduce the following operator matrix over V s 1 ,s 2 θ , M :"¨A 1´I1´B1,2 Imposing the boundary conditions of Problem 3.1 to u (sc) leads to the following system of BIEs.

‹ ‹ ‚
where M corresponds to the operator matrix in (19) and γ e,1 corresponds to the exterior trace vector on Γ 1 .
In order to ensure the well-posedness of Problem 4.22, we introduce the following set of auxiliary problems.

Problem 4.23 (Auxiliary problems). We seek tv
for each i P t1, . . . , M u, where H ą 0 is as in Section 3.1, and tk i u M i"0 are the wavenumbers in each tΩ i u M i"0 , as introduced in Section 3. By the same analysis as that presented in [41,Section 3.4], each interface Γ i , i P t1, . . . , M u, potentially adds a countable set of wavenumbers, k 0 , such that Problem 4.23 is unsolvable. This justifies the following Assumption (recall k i " η i k 0 for all i P t1, . . . , M u).  Proof. Note that the operator matrix M may be written as Then, by the Fredhom alternative, we need only show uniqueness of Problem 4.22, as the above tridiagonal block is compact by Proposition 4.21. The proof is very similar to that for the classical scattering problem of a bounded object in free space (cf. [20,Theorem 3.41]). Let Λ P V s 1 ,s 2 θ , with s 1 " s`1 2 and s 2 " s´1 2 , be such that MΛ " 0. We define r u 0 pxq :"´L k 0 θ,Γ 1 pΛ 1 q¯pxq @ x P GzΓ 1 , and further define which is well defined in each Ω i , but could potentially have non-zero jumps across each interface Γ i . Moreover, r u solves the Helmholtz equation with wavenumber k i in each Ω i and satisfies the appropriate radiation conditions at infinity [5,Section 4]. Hence, r u solves Problem 3.1, and Assumption 3.2 implies r u " 0. We continue by defining the following auxiliary functions v i pxq :" It is clear from this definition that Furthermore, each v i satisfies the appropriate radiation conditions at infinity. Using the jump relationships of BIOs (see [5,Lemma 4.11]), we have that Since r u " 0, we have that

Spectral Galerkin Method
We now provide a numerical method to approximate solutions of Problem 4.22, along with its corresponding error estimates. We restrict ourselves to cases where the interfaces tΓ i u M i"1 are C 8 -Jordan curves. By Theorem 4.25, the solution is of arbitrary smoothness, and a spectral method should converge at a super-algebraic rate (cf. [37,Chapter 9] and [25,26]).

Discrete Spaces
Let us define a suitable family of finite dimensional subspaces of V s 1 ,s 2 θ . From the definition of quasi-periodic Sobolev spaces, it is natural to consider the following finite dimensional functional spaces over p0, 2πq p E N θ :" spantp e n θ ptq :" e ıpn`θqt : n P t´N, . . . , N uu.
We can see that r E N θ,Γ i is the space spanned by finite Fourier basis parametrized on Γ i and that E N θ,Γ i is constructed from the previous space by dividing the basis by the norm of the tangential vector of the corresponding interface. As before, it is clear that both Finally, we define the Cartesian product of discrete spaces whose infinite union on N forms a dense subspace of V s 1 ,s 2 θ,Γ i for any pair s 1 , s 2 P R.

Discrete Problem
We now consider the Galerkin discretization of Problem 4.22 on the finite dimensional product space Problem 5.1 (Discrete BIEs). Let the parameters k 0 and tη i u M i"1 satisfy Assumption 4.20 and let the interfaces tΓ i u M i"1 be of class where the duality product xΨ, Ξy Γ :" denotes the sum of two standard duality pairings in H 1 2 θ pΓ i q and H´1 2 θ pΓ i q, and accounts for the right-hand side of Problem 4.22.
Since this is a second-kind BIE, we can deduce a quasi-optimality approximation result for the Galerkin discretization (cf. [38,Theorem 4.2.9]), i.e. there exists N ‹ " tN ‹ i u M i"1 such that for all N " tN i u M i"1 such that N i ą N ‹ i for all i P t1, . . . , M u, it holds that From (25) we see that, in order to establish error convergence rates for the discrete solution, we need to bound those of the best approximation. From the definition of our discrete and continuous spaces, the problem of bounding the best approximation on V s 1 ,s 2 θ is equivalent to that of establishing bounds for the best approximation of an element of H s r0, 2πs when approximated by elements of p E N r θ with r θ " 0. This issue was already addressed, for example, in [37, Theorem 8.2.1]. Specifically, for any pair r 1 , r 2 P R with r 2 ą r 1 and f P H r 2 r0, 2πs, there holds Theorem 5.2. Let the parameters k 0 and tη i u M i"1 satisfy Assumption 4.20 and let the interfaces tΓ i u M i"1 be of class C 8 . Further assume Assumptions 3.2 and 4.24 to be satisfied and let s ě 0, s 1 " s`1 2 and s 2 " s´1 2 . Then, there exists Proof. For any Ξ N P E N θ , we denote Ξ N i i " pξ N i i , ζ N i i q t for all i P t1, . . . , M u, so that By definition of our continuous and discrete spaces together with (26), we see that for all i P t1, . . . , M u, one deduces where the unspecified constant depends only on Γ i . Hence, Since the problem is well posed, we obtain where the unspecified constant now also depends on the wavenumbers tk i u M i"0 .
Remark 5.3. Theorem 5.2 states that the proposed spectral Galerkin method has a similar performance to the Nyström method, since if interfaces belong to C 8 then one obtains superalgebraic convergence (commonly observed with the Nyström method [46]). The super-algebraic convergence rate of the Nyström method for the transmission problem on a bounded object in two dimension was rigorously proved in [12]. Similar convergence results for quasi-periodic problems using the Nyström scheme are, to the best of our knowledge, not available.
Remark 5.4. It follows from Remark 4.26 that we can obtain convergence of limited order if the interfaces are of class C r,1 with r P r1, 8q.

Implementation
We continue with an overview of the procedure employed to compute the approximation Λ N . For a given N P N and l, m P Z such that´N ď |l| , |m| ď N , integrals where f and F are smooth periodic and bi-periodic functions, respectively, can be computed to exponential accuracy through the FFT to construct trigonometric interpolations of the corresponding functions (cf. [37,Theorem 8.4.1]). Since the associated kernels correspond to smooth bi-periodic functions, the computation of the block matrices B i,j on (19) is performed in this way. In terms of computational cost, the set of integrals tI 1 l u N l"´N involves 2N`1 evaluations of the function f and one FFT aplication to a vector of length 2N`1, whence the total computational cost is O pp2N`1q logp2N`1qq 1 arithmetic operations -plus 2N`1 function evaluations-to compute the 2N`1 integrals. For the set of integrals tI 2 l,m u N l,m"´N , we require p2N`1q 2 evaluations of the function F , and 2p2N`1q FFTs for vectors of length 2N`1, yielding a cost of O p2p2N`1q logp2N`1qq arithmetic operations (plus p2N`1q 2 function evaluations).
On the other hand, the block matrices A i in (19) consist of differences of the self-interaction operators on Γ i for the four BIOs. While the difference of two operators is compact-the resulting kernel is smoother than that associated to a single evaluation of the same operatorthe kernel is not arbitrarily smooth, even if the geometry is. Consequently, a deeper analysis is required before applying classical algorithms for the computation of Fourier transforms.
Let us consider, as an illustrative example, the weakly singular operator. We are required to compute integrals such as where p G k θ is as in (10). Decomposing p G k θ as shown in (11), we obtain two integrals, Since J k θ ps, tq is smooth and periodic, each of the integrals of the right-hand side is easy to compute. Moreover, the terms in the series in (28) decay exponentially and the series may be truncated at the cost of a small approximation error. Furthermore, the sum in (28) may be understood as a discrete convolution, allowing it to be computed by multiplying the corresponding Fourier transforms (see [25] for details).
The computational cost of computing tI S l,m u N l,m"´N and tI R l,m u N l,m"´N is dominated by the latter set of integrals, since it involves 2pN`1q 2 evaluations of the quasi-periodic Green's function, which is done following [13]. The evaluation cost of the quasi-periodic Green's function corresponds to p2N`1q 2 p2N 1`1 q evaluations of the Hankel function, with N 1 ą N is a truncation parameter for the series in (5) 2 . Meanwhile, the total cost for I S l,m is proportional to p2N`1q logp2N`1q.
For the operators K k θ and K 1k θ , a similar technique can be applied using (15). The integrals corresponding to the hyper-singular BIO are approximated by first using the integration-byparts formula in Lemma 4.13, reducing it to two different integrals which are then approximated as those corresponding to the weakly-singular BIO.
Considering M ą 1 interfaces, 2N`1 degrees of freedom on each interface and N 1 proportional to N the total cost of the matrix assembly process can be estimated as OpN 3 M q Hankel function evaluations and OpM N 2 log N q arithmetic operations. We point out that the cost could be reduced drastically by constructing an accurate algorithm to approximate the Hankel functions by pre-computing some values.
Remark 5.5. We have restricted ourselves to the analysis of the semi-discrete case, that is, we do not take into consideration the error coming from the approximation of the integrals for the error bound in Theorem 5.2. However, it is not difficult to incorporate it. Assuming that the parametrizations tz i u M i"1 correspond to Jordan curves of class C 8 and using the aliasing proprieties of the Fourier basis [44,Chapter 4] and Lemma 4.11 we have that the approximation error for the computation of the required integrals is OpN´ ´1 q, with 2N`1 being the number of degrees of freedom per interface and an arbitrarily large integer. Then, the fully discrete error can be obtained by an application of Strang's lemma [38,Section 4.2.4], from where it follows that the behaviour of the fully discrete error, with respect to N , is the same as in Theorem 5.2.  norm with respect to the analytic solution. We have included results for different numbers ghost layers (1,2 and 3, respectively), i.e., the first experiment considers only the first 3 layers (counting downwards), the second one considers the first 4 layers and the third considers all 5 layers.

Numerical Examples
We now showcase computational experiments verifying the convergence estimates found in Theorem 5.2. The implementation of the aforementioned algorithms was carried through a C++ cpu-only library. All the experiments ran on a Intel I7-4770@3.4GHZ processor with 8 threads. The code was compiled with gcc 4.9.4, openmp and O2 flags on.

Code Validation
We begin by considering the simple case of a grating with two media separated by a single horizontal line segment acting as its layer. Hence, using the following expansion of the Green's function [5,Proposition 4.2]: it is possible to assemble the matrix analytically. The matrix M is then composed of only block diagonal terms. Since the right-hand side only has two non-null components 3 , only the corresponding components for the solution are non-zero, yielding a closed form for the solution.
In order to test the implementation, we consider an artificial (harder) problem by including ghost domains, i.e., we add extra smooth (ghost) layers that separate domains with the same refraction index. Hence, the solution is the same as if these additional domains did not exist and has a closed form, as before. The results for different ghost layers are reported in Figure  2.  Table 1: Value of the refraction indices tη p1q i u 12 i"1 and tη p2q i u 12 i"1 (corresponding to the two considered cases) for the grating in Figure 5 (counting downwards).
We also display the convergence behaviour of the method for interfaces with limited regularity by repeating the previous experiment (same incident field) with one ghost domain and an interface given by where a, b are real numbers that scale the interface, and p is an odd integer that determines the smoothness degree of the interface. In particular, z 3 is in C p´2,1 or, more precisely, C p´1 with an integrable p-th derivative. Results are reported in Figure 3.
For all experiments in this section, the frequency is chosen as k 0 " 1 and the incidence angle is 0.47 radians.

Convergence results
We now consider a smooth geometry composed of the 12 layers and varying refraction indices. Two different scenarios for the choice of indices are employed, reported in Table 1 (η p1q i and η p2q i for the first and second cases, respectively). We also consider three different wavenumbers for the incident wave, k 0 " 2.8, 14 and 28. Convergence results in the energy norm for the solution of Problem 5.1 for the different cases of parameters and wavenumbers are reported in Figure 4, where exponential convergence is observed for all considered scenarios, as expected. All errors were computed with respect to an overkill solution, with approximately 50 more bases per interface than the last plotted point for each curve. The incidence angle is, again, 0.47 radians.
Finally, in Figure 5 we present an illustration of the total field corresponding to the refraction indices given in Table 1.
Remark 6.1. Though establishing the relation between the parameters-tη i u M i"0 and k 0 -and the number of basis elements required to attain a certain desired accuracy is not straightforward, our experiments suggest that N should be chosen proportional to the maximum wavenumber k max :" max iPt0,...M u k i .

Conclusions
We have proposed a fast spectral method for the efficient representation, through surface potentials based on the quasi-periodic Green's function, for the solution of the Helmholtz equation with transmission boundary conditions on a periodic domain. Theorem 5.2, we obtained convergence estimates for the discrete approximation of the corresponding boundary data, and found that discrete solution converge at a super-algebraic rate to continuous solutions of the considered boundary integral equation. Though, we focused on the Helmholtz transmission problem, our approximation results and convergence estimates can be easily extended to other boundary integral equations on quasi-periodic Sobolev spaces whenever the formulation is well posed. We avoided Rayleigh-Wood anomalies from our analysis since the series in (5) fails to converge for said frequencies and, for the same reason, our previous results from [5] exclude them as well.
Though similar numerical results are known for the Nyström Method, theoretical results confirming the observed convergence rates are scarce [12], which is an advantage of Galerkin discretizations such as that presented in this article. Moreover, the convergence rate for the proposed discretization is equal to that expected of Nyström methods, so that it is numerically competitive with them while inheriting the theoretical benefits of a Galerkin discretization.
Future work considers: (i) including Rayleigh-Wood anomalies to our analysis, (ii) extending our results to three dimensional Helmholtz equations and Maxwell's equations on periodic domains and (iii) applications in uncertainty quantification [40] and shape optimization [6]. norm with respect to the analytic solution. The legend indicates an estimate of the slope of the error convergence curves for different values of p (degrees of smoothness). Classically, error convergence estimates for spectral methods indicate the slope to be at least equal to p. We also consider the case p " 2, where the extra layer is C 8 and the super-algebraic convergence rate is observed.   Table 1. Notice that the curves in red-corresponding to parameters η p2q i in Table 1-display a longer preasymptotic regime before convergence is observed for all considered values of k 0 , seemingly due the presence of layers with higher wavenumbers (see Remark 6.1).