Oscar P Bruno

Bérangère Delourme

Rapidly convergent two-dimensional quasi-periodic Green function throughout the spectrum-including Wood anomalies

We introduce a new methodology, based on new quasi-periodic Green functions which converge rapidly even at and around Wood-anomaly configurations, for the numerical solution of problems of scattering by periodic rough surfaces in two-dimensional space. As is well known the classical quasi-periodic Green function ceases to exist at Wood anomalies. The approach introduced in this text produces fast Green function convergence throughout the spectrum on the basis of a certain "finite-differencing" approach and smooth windowing of the classical Green function lattice sum. The resulting Green-function convergence is super-algebraically fast away from Wood anomalies, and it reduces to an arbitrarily-high (user-prescribed) algebraic order of convergence at Wood anomalies.

Introduction

We consider the problem of evaluation of the fields scattered by a periodic perfectly conducting surface under planewave illumination. This problem has been extensively studied as it impacts upon a wide range of areas of science and engineering, including optics, photonics, communications, and stealth, and, through them, many fields of physics, astronomy, chemistry, biology and metallurgy [START_REF] Loewen | Diffraction gratings and applications[END_REF][START_REF] Valenzuela | Theories for the interaction of electromagnetic and oceanic waves-a review[END_REF]. A variety of numerical approaches have been used to tackle this important problem [START_REF] Bruno | Numerical solution of diffraction problems: a method of variation of boundaries III. doubly periodic gratings[END_REF][START_REF] Loewen | Diffraction gratings and applications[END_REF][START_REF] Nevière | Light propagation in periodic media: differential theory and design[END_REF][START_REF] Petit | Electromagnetic Theory Of Gratings[END_REF] including, notably, methods based on use of integral equations [START_REF] Maystre | Rigorous vector theories of diffraction gratings[END_REF]. Recent integral equation methods [START_REF] Bruno | Efficient high-order evaluation of scattering by periodic surfaces: deep gratings, high frequencies, and glancing incidences[END_REF][START_REF] Bruno | Efficient high-order evaluation of scattering by periodic surfaces: vectorparametric gratings and geometric singularities[END_REF], in particular, have made it possible to obtain, in reasonable computing times, highly accurate solutions for very challenging problems of scattering by periodic surfaces. The success of this methodology lies in part on its inherent dimensionality reduction (only the scattering surface needs to be discretized, not the surrounding volume) and associated automatic enforcement of radiation conditions; mathematical analyses of the integral equation method in various contexts, including scattering by periodic surfaces and bounded obstacles, can be found in [START_REF] Arens | On integral equation and least squares methods for scattering by diffraction gratings[END_REF][START_REF] Bruno | Solution of a boundary value problem for the helmholtz equation via variation of the boundary into the complex domain[END_REF][START_REF] Colton | Inverse Acoustic And Electromagnetic Scattering Theory[END_REF][START_REF] Desanto | Theoretical and computational aspects of scattering from rough surfaces: One-dimensional perfectly reflecting surfaces[END_REF][START_REF] Kress | Linear Integral Equations[END_REF][START_REF] Maystre | Rigorous vector theories of diffraction gratings[END_REF][START_REF] Petit | Electromagnetic Theory Of Gratings[END_REF] and references therein.

The properties of integral equations for periodic surfaces under plane-wave incidence are closely related to the character of the corresponding quasi-periodic Green functions used. As is well known, classical expressions for quasi-periodic Green functions converge extremely slowly, and a number of methods have therefore been introduced to produce rapidly convergent Green-function algorithms, including the well known Ewald summation method [START_REF] Arens | Analysing Ewald's method for the evaluation of Green's functions for periodic media[END_REF][START_REF] Capolino | Efficient computation of the 3d Green's function for the Helmholtz operator for a linear array of point sources using the Ewald method[END_REF][START_REF] Linton | The Green's function for the two-dimensional Helmholtz equation in periodic domains[END_REF] for two-and three-dimensional problems and, for the two-dimensional case, the highly efficient algorithm [START_REF] Veysoglu | Polarimetric passive remote sensing of periodic surfaces[END_REF]. Many other contributions have in fact been put forward over the years to facilitate evaluation of quasi-periodic Green functions; in addition to those mentioned above here we mention [START_REF] Chandler-Wilde | Efficient calculation of the Green function for acoustic propagation above a homogeneous impedance plane[END_REF][START_REF] Dienstfrey | Lattice sums and the two-dimensional, periodic Green's function for the Helmholtz equation[END_REF][START_REF] Kurkcu | Stable and efficient evaluation of periodized green's functions for the helmholtz equation at high frequencies[END_REF][START_REF] Mathis | A Comparison Of Acceleration Procedures For The Two-Dimensional Periodic Green's Function[END_REF][START_REF] Moroz | Exponentially Convergent Lattice Sums[END_REF][START_REF] Nicorovici | Lattice sums for off-axis electromagnetic scattering by gratings[END_REF][START_REF] Nicorovici | Efficient calculation of the Greens function for electromagnetic scattering by gratings[END_REF]; a recent survey can be found in [START_REF] Linton | Lattice sums for the Helmholtz equation[END_REF]. A combined approach which takes advantage of various methods, applying each algorithm for configurations for which it is most efficient (for the challenging three-dimensional Green function problem), was put forth in [START_REF] Guérin | Combined method for the computation of the doubly periodic Green's function[END_REF].

As is well known, none of these methods for evaluation of the quasi-periodic Green function can be applied to problems of scattering by periodic surfaces at Wood-anomaly configurations [START_REF] Rayleigh | Note on the remarkable case of diffraction spectra described by prof. wood[END_REF][START_REF] Wood | On a remarkable case of uneven distribution of light in a diffraction grating spectrum[END_REF] (at which one or more scattered waves propagate in a direction parallel to the scattering surface): for Wood-anomaly configurations the classical periodic Green function is not even defined (see Remark 2.2 for details on nomenclature concerning Wood anomalies). To address this difficulty we propose new quasi-periodic Green functions and associated series representation which converge rapidly even at and around Wood anomalies. More precisely, we present a set of rapidly convergent quasi-periodic Green functions G q j whose N -term truncated series converge, at Wood anomalies, at least as fast as (1/N ) (j-1)/2 for j even (resp. (1/N ) j/2 for j odd) as N → ∞; in view of the fact that this approach also incorporates the smooth windowing function methodology [START_REF] Monro | Windowing-Based Approach to the Evaluation of Scattering from Periodic Rough Surfaces[END_REF], the new Green functions also enjoy super-algebraically fast convergence (faster than any power of N ) away from Wood anomalies.

The approach introduced in this text produces fast Green-function convergence at and around Wood anomalies on the basis of a certain order-j "finite-differencing" method (for positive integer values of j). To our knowledge, this is the first approach ever presented that is applicable to problems of scattering by diffraction gratings at Wood anomalies on the basis of quasi-periodic Green functions.

It should be noted that a "method-of-images Green-function", which is related to our j = 1 Green function approach, was used in [START_REF] Chandler-Wilde | Scattering by rough surfaces: the Dirichlet problem for the Helmholtz equation in a non-locally perturbed half-plane[END_REF][START_REF] Zhang | Integral equation methods for scattering by infinite rough surfaces[END_REF] to treat problems of scattering by nonlocal, non-periodic perturbations of a line in two dimensions. The j = 1 method, which suffices to yield (slow) convergence in the two dimensional case, does not give rise to convergence in three dimensions: for three dimensional configurations convergence only results for j ≥ 2 [START_REF] Bruno | Efficient evaluation of doubly periodic Green functions in 3d scattering, including Wood anomaly frequencies[END_REF]. In this context we also mention the recent work [START_REF] Barnett | A new integral representation for quasi-periodic fields and its application to two-dimensional band structure calculations[END_REF][START_REF] Barnett | A new integral representation for quasi-periodic scattering problems in two dimensions[END_REF] which, for two-dimensional problems, introduces an alternative integral equation which does not utilize a quasi-periodic Green function, and which is also applicable at Wood anomalies: in that approach quasi-periodicity is enforced through use of auxiliary layer potentials on the periodic cell boundaries. The practical feasibility of an extension of this methodology to three dimensional problems has not as yet been established.

In order to demonstrate the character of the new approach we present efficient numerical methods, based on the new Green functions, for the solution of quasi-periodic scattering problems throughout the spectrum-even at and around Wood anomalies-; as shown in Section 4.4, certain slight modifications of the direct finite-differencing Green function expressions mentioned above need to be introduced to obtain uniquely solvable integral-equation problems. We further mention that, as indicated in [START_REF] Bruno | Efficient evaluation of doubly periodic Green functions in 3d scattering, including Wood anomaly frequencies[END_REF], additional acceleration can be induced in the Green-function convergence by means of the FFT-based equivalent-source methodology [START_REF] Bruno | A fast, high-order algorithm for the solution of surface scattering problems: basic implementation, tests, and applications[END_REF]; see Remark 6.2. Our numerical results demonstrate the capabilities of the new methodology: even in absence of the acceleration method [START_REF] Bruno | A fast, high-order algorithm for the solution of surface scattering problems: basic implementation, tests, and applications[END_REF][START_REF] Bruno | Efficient evaluation of doubly periodic Green functions in 3d scattering, including Wood anomaly frequencies[END_REF], the present approach can solve the complete scattering problem for rather challenging periodic surfaces (steep gratings) in the resonance regime, including Green function computations and matrix inversion by means of Gaussian elimination, in total computing times of a few tens to a few hundreds of milliseconds-depending on the complexity of the problem.

The remainder of this paper is organized as follows. Section 2 describes the scattering problem under consideration and it presents some background on quasi-periodic function and integral equations. Section 3 describes the smooth windowing method that gives rise to super-algebraically converging Green function away from Wood anomalies. Section 4 then presents the new rapidly convergent Green functions together with necessary theoretical discussions involving the Green function itself and associated integral equations. Section 5 describes the numerical implementation of the new Green functions and integral equations, and section 6, finally, presents a variety of numerical results demonstrating the properties of the overall proposed approach.

Preliminaries

Scattering problem

We consider the problem of scattering of a plane wave by a perfectly reflecting periodic surface

Γ = {(x, f (x)), x ∈ R} (1)
with f ∈ C r per (R), r ≥ 2, where, for any non-negative integer r, C r per (R) denotes the set of L-periodic r-times continuously differentiable functions defined in the real line. The propagation domain is thus the region Ω = {(x, y) ∈ R 2 , such that y > f (x)}.

(

) 2
Letting k be a positive wavenumber and, further, letting θ ∈ (-π/2, π/2), α = k sin(θ) and β = k cos(θ), we assume the periodic surface is illuminated by the incident plane wave u inc (x, y) = e i(α x-βy)

(3) so that θ is the angle of propagation of the incident field measured counterclockwise from the negative y-axis. The scattered field u s ∈ H 1 loc (Ω) is a quasi-periodic solution of the homogeneous Helmholtz equation that satisfies a radiation condition at infinity. More precisely, u s satisfies the Partial Differential Equation (PDE)

∆u s + k 2 u s = 0 in Ω (4)
as well as the quasi-periodicity condition

u s (x + L, y) = u s (x, y)e iαL (5)
together with either the Dirichlet boundary conditions

u s = -u inc on Γ (6)
or the Neumann boundary conditions

∂u s ∂ν = - ∂u inc ∂ν on Γ. ( 7
)
Here ν is the unit normal

ν(x) = (-f ′ (x), 1) 1 + f ′ (x) 2 . ( 8
)
The aforementioned condition of radiation at infinity results from consideration of the Rayleigh expansion [START_REF] Petit | Electromagnetic Theory Of Gratings[END_REF] u s (x, y) = n∈N a n e i(αnx+βny) + b n e i(αnx-βny) , y > H = max

x∈R f (x), (9)
for the solution u s . Here

α n = α + n 2π L , β n =    k 2 -α 2 n if α 2 n ≤ k 2 , i α 2 n -k 2 otherwise. (10)
We say that u s satisfies the condition of radiation at infinity (or that u s is outgoing) if b n = 0 for all n ∈ Z. In other words, u s is outgoing if and only if u s is given by a Rayleigh expansion of the form

u s (x, y) = n∈N a n e i(αnx+βny) , y > H, (11)
see e.g. [START_REF] Petit | Electromagnetic Theory Of Gratings[END_REF].

Remark 2.1.
In what follows we denote by U the (finite) set of integers for which α 2 n < k 2 . For n ∈ U the function e i(αnx+βny) is a propagative (outgoing) plane wave. In the case α 2 n > k 2 , in contrast, the corresponding functions e i(αnx+βny) are evanescent: they decrease exponentially as y → ∞. In the limiting case α 2 n = k 2 we have a Wood anomaly frequency (cf. Remark 2.2 below). In this case the function e iαnx+iβny = e iαnx is a grazing plane wave, that is, it is a plane wave that propagates parallel to the grating. For given period L and angle θ we denote by K L,θ = {k : k = |α n |} the corresponding set of Wood frequencies. (We also write K instead of K L,θ when specification of the period and incidence angle is either not crucial or is otherwise clear from the context.) Remark 2.2. With reference to Remark 2.1, it is clear that continuous variations of one or more of the quantities k, L and θ can cause the n-th diffracted mode e iαnx+iβny to change from evanescent (β n imaginary) to propagative (β n real) or viceversa, passing by a wave with β n = 0 (k ∈ K L,θ ) at some "pass-off" intermediate configuration.

Rayleigh [START_REF] Rayleigh | Note on the remarkable case of diffraction spectra described by prof. wood[END_REF] connected such pass-off k ∈ K L,θ configurations, which are characterized by energy redistribution as cut-off of diffraction orders change, with previous experimental observations by Wood [47] of anomalous diffraction behavior. Eventually it was noticed that Rayleigh's pass-off configurations account for some but not all of Wood's observations of anomalous diffraction; see [31, p. 285] or, for a detailed discussion, [START_REF] Maystre | Theory of Wood's anomalies[END_REF]. It might therefore be appropriate to call pass-off configuration "Rayleigh-Wood configurations". For the sake of brevity, and in keeping with much of the existing literature, however, throughout this paper a configuration for which k ∈ K L,θ will be called a called a Wood configuration, and k will be called a Wood frequency, or a Wood anomaly frequency.

Remark 2.3. With reference to Remarks 2.1 and 2.2 we note that for a Wood frequency the Neumann problem for a flat interface (f = 0) is not uniquely solvable: for that case the Rayleigh mode e iαnx is a non-zero solution with zero Neumann boundary data.

For each n ∈ U we define the associated efficiency

e n = |a n | 2 β n β .
The (finite) set of all efficiencies satisfies the energy balance relation [START_REF] Petit | Electromagnetic Theory Of Gratings[END_REF] n∈U

e n = 1. ( 12
)
Remark 2.4. As is known [START_REF] Bonnet-Bendhia | Guided waves by electromagnetic gratings and nonuniqueness examples for the diffraction problem[END_REF][START_REF] Kirsch | Diffraction by periodic structures[END_REF][START_REF] Kirsch | Uniqueness theorems in inverse scattering theory for periodic structures[END_REF][START_REF] Nédélec | Integral equation methods in a quasi-periodic diffraction problem for the timeharmonic Maxwell's equations[END_REF][START_REF] Petit | Electromagnetic Theory Of Gratings[END_REF], the Dirichlet problem ( 4)-( 6) admits a unique outgoing solution u s ∈ H 1 loc (Ω) even at Wood anomalies (cf. Remark 2.3). A uniqueness result is not available under the Neumann boundary conditions [START_REF] Bruno | Efficient high-order evaluation of scattering by periodic surfaces: deep gratings, high frequencies, and glancing incidences[END_REF], although it is known that the Neumann problem ( 4), ( 5), ( 7), [START_REF] Bruno | Numerical solution of diffraction problems: a method of variation of boundaries III. doubly periodic gratings[END_REF] does admit a unique solution except possibly for a discrete set of frequencies k that can only accumulate at infinity. We note, however, that numerical evidence suggests that uniqueness (and thus existence) of solution for this problem does hold at least for all k ∈ K (cf. [23, p. 147]).

Classical quasi-periodic Green function

Let k / ∈ K (k is not a Wood anomaly) and let (X, Y ) ∈ R 2 . The quasi-Periodic Green function G q is defined by

G q (X, Y ) = n∈Z e -iαnL G(X + nL, Y ), (13)
where, denoting by H

0 the first Hankel function of order zero, G is the two-dimensional free space Green function

G(X, Y ) = i 4 H (1) 0 (k X 2 + Y 2 ). (14)
As is known, the series (13) converges for (X, Y ) = (rL, 0), r ∈ Z. Further, the truncated series

G q (X, Y ) = n∈Z,|n|≥2 e -iαnL G(X + nL, Y )
converges uniformly in any compact set not containing singularities of G q (see e.g. theorem 4.1 in [START_REF] Bruno | Solution of a boundary value problem for the helmholtz equation via variation of the boundary into the complex domain[END_REF]). For reference, finally, we note that the Green function also admits the Rayleigh expansion

G q (X, Y ) = n∈Z i 2L e iαnX+iβn|Y | β n (15)
(see e.g. theorem 4.4 in [START_REF] Bruno | Solution of a boundary value problem for the helmholtz equation via variation of the boundary into the complex domain[END_REF]).

Integral equation formulations

Here and in what follows we denote by Ω # and Γ # the intersection of Ω and Γ, respectively, with the set x ∈ (-L 2 , L 2 ) × R:

Ω # = (x, y) ∈ - L 2 , L 2 × R, such that y > f (x) , Γ # = (x, f (x)), x ∈ - L 2 , L 2 . ( 16
)
As is well known, the Dirichlet and Neumann scattering problems described in Section 2.1 can be reduced to second kind integral equations over the curve Γ # (cf. for instance [START_REF] Kirsch | Diffraction by periodic structures[END_REF][START_REF] Petit | Electromagnetic Theory Of Gratings[END_REF]). Indeed, a variety of such integral formulations exist. For definiteness, in what follows we present two types of integral equations, one for each one of the problems presented in Section 2.1, which we will consider throughout this text.

In the Dirichlet case the scattered field can be expressed in the form

u s D (x, y) = Γ # ∂ ν(x ′ ) G q x -x ′ , y -f (x ′ ) µ(x ′ ) ds(x ′ ),
where the density µ is the solution of the integral equation

-u inc |Γ # = Γ # ∂ ν(x ′ ) G q x -x ′ , f (x) -f (x ′ ) µ(x ′ ) ds(x ′ ) + 1 2 µ(x ′ ).
For the Neumann problem, in turn, the scattered field is given by the expression

u s N (x, y) = Γ # G q x -x ′ , y -f (x ′ ) η(x ′ ) ds(x ′ ).
where η is the solution of the integral equation

-∂ ν(x) u inc |Γ # = Γ # ∂ ν(x) G q x -x ′ , f (x) -f (x ′ ) η(x ′ ) ds(x ′ ) - 1 2 η(x ′ ).

Super-algebraically convergent representation away from Wood anomalies

In this section we present a novel methodology, put forth recently [START_REF] Monro | Windowing-Based Approach to the Evaluation of Scattering from Periodic Rough Surfaces[END_REF], for evaluation of the quasi-periodic Green function away from Wood anomalies. In fact, in view of its qualities at non-Wood anomaly frequencies and for ease of implementation of high-order integral Nyström solvers, the basic element in this method (smooth windowing) is incorporated as part of our all-frequency scheme-valid for both Wood and non-Wood frequencies-which we introduce in Sections 4 and 5. As established in [START_REF] Monro | Windowing-Based Approach to the Evaluation of Scattering from Periodic Rough Surfaces[END_REF] (Theorem (3.1) below) and as demonstrated by means of a simple example below in this section, integrals of products of smoothly windowed Green functions and quasiperiodic functions converge exponentially fast to the corresponding integrals involving the exact periodic Green function. Using an argument based on summation by parts (instead of the integration by parts calculation utilized below in this section) it can be shown that a corresponding convergence result holds for the windowed periodic Green function itself.

To proceed with our smooth-windowing Green-function algorithm we consider the smooth cut-off function

S(x, x 0 , x 1 ) =        1 if |x| ≤ x 0 , exp 2e -1/u u-1 if x 0 < |x| < x 1 , u = |x|-x 0 x 1 -x 0 , 0 if |x| ≥ x 1 , (17)
(see Figure 1), and we define the approximate periodic Green function and a corresponding approximation of its normal derivative

G q A by G q A (X, Y ) = i 4 n∈Z e -iαnL H (1) 0 (k (X + nL) 2 + (Y ) 2 ) S(X + nL, cA, A) (18)
H q A (x, x ′ ) = i 4 n∈Z e -iαnL ∂ ν(x ′ ) H (1) 0 (k (x -x ′ + nL) 2 + (y -y ′ ) 2 ) S(x -x ′ + nL, cA, A) , (19)
where A > L and 0 < c < 1. Then, as established in [START_REF] Monro | Windowing-Based Approach to the Evaluation of Scattering from Periodic Rough Surfaces[END_REF] the following result holds:

Theorem 3.1. Let k / ∈ K and let f be a smooth function. Then, for any α-quasi-periodic smooth function µ the integrals

Γ # H q A (x, x ′ ) µ(x ′ )ds(x ′ ) and Γ # G q A (x -x ′ , f (x) -f (x ′ )) µ(x ′ )ds(x ′ )
converge super-algebraically fast (that is, faster than any power to 1/A) to

Γ # ∂ ν(x ′ ) G q (x -x ′ , f (x) -f (x ′ )) µ(x ′ )ds(x ′ ) and Γ # G q (x -x ′ , f (x) -f (x ′ )) µ(x ′ )ds(x ′ ),
respectively, as A tends to infinity.

The main idea of the proof of this result can be conveyed by a simplified example that results as the Hankel function is substituted by its large-argument asymptotic expression (cf. reference [START_REF] Monro | Windowing-Based Approach to the Evaluation of Scattering from Periodic Rough Surfaces[END_REF] which also contains a complete proof of Theorem 3.1). For definiteness we consider the case of double-layer potential; the case of the single layer is even more direct. To construct our simplified example note that, since µ is α-quasiperiodic, the double layer potential can be expressed in the form

Γ # ∂ ν(x ′ ) G q x -x ′ , f (x) -f (x ′ ) µ(x ′ )ds(x ′ ) = ∞ -∞ ∂ ν(x ′ ) G x -x ′ , f (x) -f (x ′ ) µ(x ′ ) 1 + f ′ (x ′ ) 2 dx ′ . (20)
Taking into account the large z asymptotic expression H

∂ ν(x ′ ) G x -x ′ , f (x) -f (x ′ ) = C(x, x ′ ) e ikx ′ √ x ′ , ( 21
)
where C is a function which remains bounded from zero and ∞ as x-x ′ → ∞. On the other hand µ(x ′ ) 1 + f ′ (x ′ ) 2 is an α-quasiperiodic function, and it can therefore be expanded in a Rayleigh series of the form ∞ n=-∞ a n e iαnx ′ .

In our simplified example we thus replace the Green function by the asymptotic form [START_REF] Guérin | Combined method for the computation of the doubly periodic Green's function[END_REF] with C(x, x ′ ) = 1, so that the integrand equals an infinite sum of constant multiples of terms of the form e iknx √ x (where

k n = k -α - 2πn L
). And, in fact, for our example we consider just one such term, that is, we study the integration problem

I ex = +∞ 0 e iknx ′ √ x ′ dx ′ .
As stated above, throughout this section we assume k is not a Wood anomaly frequency (k / ∈ K), or, in other words, k n = 0 for all n. To illustrate Theorem 3.1 in the present context we investigate theoretically and numerically the convergence of the approximation

I S,A = +∞ 0 S(x ′ , cA, A) e iknx ′ √ x ′ dx ′
to the exact value I ex . For comparison purposes we also consider the classical approximation

I H,A = A 0 e iknx ′ √ x ′ dx ′ , (22)
which can be viewed as the result of substituting the smooth windowing function S by a suitable Heaviside function H, and which corresponds, in the context of this section, to the direct truncation of the series (13) keeping a number of the order of O(A/L) terms. The error I ex -I S,A in the windowed approximation is given by

I ex -I S,A = +∞ 0 (1 -S(x ′ , cA, A)) e iknx ′ √ x ′ dx ′ = +∞ cA P (x ′ , cA, A) e iknx ′ √ x ′ dx ′ ,
where P (x ′ , c, A) = 1 -S(x ′ , cA, A). Using the the change of variables x = x ′ ca we obtain

I ex -I S,A = √ cA +∞ 1 P (x, 1, 1 c ) √ x e ikncAx dx.
Given that P (1, 1, 1 c ) = 0, integration by parts yields

I ex -I S,A = - 1 ik n √ cA +∞ 1 ∂ ∂x P (x, 1, 1 c ) √ x e ikncAx dx,
while, more generally, p-times iterated integration by parts gives rise to the relation

I ex -I S,A = (-1) p (ik n ) p (cA) p-1/2 +∞ 1 ∂ p ∂x p P (x, 1, 1 c ) √ x e ikncAx dx.
Noting that the derivatives of

P (x,1, 1 c ) √ x
do not depend on A, and that they are bounded functions with support equal to the interval [1, 1/c], for all p ∈ N we obtain

|I ex -I S,A | ≤ C √ k n (k n c A) p-1/2 .
In other words, the error in the windowed approximation I S,A is super-algebraically small: it is a quantity of order A -p for all p ≥ 1.

The widowed integration method thus results in much closer approximations than the direct un-windowed approximation I H,A : a simple argument involving a single integration by parts shows that the error in the latter approximation only decays like A -1 2 . Table 3 illustrates the super-algebraically fast convergence of I S,A as well as the extremely slow convergence of I H,A .

4 Rapidly convergent quasi-periodic series at and around Wood anomalies

New quasi-periodic Green function: introduction

With reference to equation ( 13), we seek a quasi-periodic Green function Gq j that can be expressed by a rapidlyconvergent series

Gq j (X, Y ) = ∞ n=-∞ e -iαnL G j (X + nL, Y ), (23)
(j ∈ N), where G j (X, Y ) is a j-dependent half space Green function, defined in Section 4.2, which decays rapidly as X tends to infinity. (The super-index q on the left hand side of equation ( 23) denotes quasi-periodicity, and the tilde accent is used there to differentiate the Green function [START_REF] Kirsch | Uniqueness theorems in inverse scattering theory for periodic structures[END_REF] with a related but different version of it which is introduced in Section 4.4; see also Remark 4.1.) Because of this fast decay, a truncated sum of the form

GN,q j (X, Y ) = |n|≤N -1 e -iαnL G j (X + nL, Y ) (24)
converges more rapidly than does a corresponding truncation of the classical quasi-periodic Green function [START_REF] Capolino | Efficient computation of the 3d Green's function for the Helmholtz operator for a linear array of point sources using the Ewald method[END_REF]: as shown in Section 4.3 below, there exists a constant C M > 0 such that for all Y ∈ (-M, M ) and and for all X ∈ (-L, L) 2 we have

Gq j (X, Y ) -GN,q j (X, Y ) ≤        C M N j-1 2
for j even

C M N j 2
for j odd.

In particular, the series expansion [START_REF] Kirsch | Uniqueness theorems in inverse scattering theory for periodic structures[END_REF] for the new Green function converges for all frequencies k, even at and around Wood anomaly frequencies.

Remark 4.1. Our proposed quasi-periodic Green-function, which is denoted by G q j , results from a slight but important modification of Gq j that is necessary to induce uniqueness in associated integral equation formulations; see Section 4.4 for details.

Rapidly decaying (non-quasi-periodic) half-space Green function G j

Rapidly decaying half-space Green functions G j can be constructed on the basis of linear combinations of the regular free-space Green functions with arguments shifted by a number j of shifts. For example, letting j = 1 and

G 1 (X, Y ) = H (1) 0 k X 2 + Y 2 -H (1) 0 k X 2 + (Y + h) 2 (25)
in view of the mean value theorem we have

G 1 (X, Y ) = h k(Y + ξ) X 2 + (Y + ξ) 2 H (1) 1 k X 2 + (Y + ξ) 2 (26)
for some real value of ξ, Y < ξ < Y + h. But in view of the asymptotic formula

H (1) n (t) = 2 πt e i(t-n π 2 -π 4 ) 1 + O 1 t , t ∈ R, n ∈ N 0 , (27)
there exists a constant C > 0 such that, for large |X| we have

|G 1 (X, Y )| ≤ C |X| 3/2 . ( 28
)
In particular, the function G 1 , which is a Green function in the domain Y > -h, decays faster, as X → ∞, than the free space Green function [START_REF] Chandler-Wilde | Efficient calculation of the Green function for acoustic propagation above a homogeneous impedance plane[END_REF].

Green functions with arbitrarily fast algebraic decay can be obtained by a suitable generalization of these ideasusing, for a given positive integer j, a finite-difference operator that approximates a Y -derivative operator of order j. In this paper we use the finite-difference operator F j : C j+1 → C given by

F j (u 0 , . . . , u j ) = j ℓ=0 (-1) ℓ C j ℓ u ℓ , (29)
where

C j ℓ = j! ℓ!(j -ℓ)!
are the binomial coefficients. Denoting by I j : C j+1 → C j+1 , E j : C j+1 → C j+1 and P j : C j+1 → C the identity operator, the left-shift operator, and the first coordinate projection operator, respectively, I j (u 0 , u 1 , . . . , u j ) = (u 0 , u 1 , . . . , u j ), E j (u 0 , u 1 , . . . , u j ) = (u 1 , . . . , u j , 0),

P j (u 0 , u 1 , . . . , u j ) = u 0 , (30)
as a result of a binomial expansion we obtain

F j = P j • I j -E j j .
We thus see that F j results from j applications of the first-order difference operator I j -E j .

Remark 4.2. Since for each polynomial

P = P (Z) of degree k we have that Q(Z) = P (Z + ℓδ) -P (Z + (ℓ -1)δ) is a polynomial of degree (k -1)
, it follows that, as is known, the operator (29) produces exactly the j-th order derivatives of polynomials of degree strictly smaller than j: for all polynomials of the form

P (Z) = j-1 ℓ=0 a ℓ Z ℓ we have j ℓ=0 (-1) ℓ C j ℓ P (Z + δℓ) = 0, (31)
see e.g. [20, eqn. 5.42].

Using a Taylor expansion of order j -1 and including the j-th order remainder, it follows from Remark 4.2 that for a sufficiently smooth function v we have

j ℓ=0 (-1) ℓ C j ℓ v(Z + δℓ) = δ j j ℓ=1 (-1) ℓ C j ℓ ℓ j j! v (j) (Z + ξ ℓ ), 0 ≤ ξ ℓ ≤ δ ℓ. (32)
For given h > 0 and j ∈ N we thus define the j-th rapidly-decaying half-space Green function by

G j (X, Y ) = i 4 j ℓ=0 (-1) ℓ C j ℓ H (1) 0 k X 2 + (Y + ℓ h) 2 (33)
for (X, Y ) ∈ R 2 , (X, Y ) = (0, ℓh), and for all non-negative integers ℓ ≤ j. With this selection for the function G j , the corresponding quasi-periodic Green function ( 23) is given by

Gq j (X, Y ) = i 4 ∞ n=-∞ e -iαnL j ℓ=0 (-1) ℓ C j ℓ H (1) 0 k (X + nL) 2 + (Y + ℓ h) 2 . ( 34
)
In view of the Rayleigh expansion ( 15) together with equations ( 13) and ( 34), further, it is easy to check that for Y ≥ 0 Gq j can also be expressed in the form

Gq j (X, Y ) = n∈Z i 2Lβ n j ℓ=0 (-1) ℓ C j ℓ e iβnℓ h e iαnX e iβnY (35)
which plays a crucial role in the null-space study presented in Section 4.4.

As shown in Lemma 4.3 below, the Green function G j and its derivatives do indeed decay rapidly; the corresponding fast convergence of the left hand series expression in equation ( 23) is established in Section 4.3.

Lemma 4.3. Let j ∈ N, h > 0, ℓ = 0 or ℓ = 1, m = 1 or m = 2
and k > 0 be given. Then, for each M > 0 there exists a positive constant C M (that also depends on j, k and h) such that for all Y ∈ (-M, M ) and for all real numbers X with |X| > 1 we have

∂ ℓ m G j (X, Y ) ≤          C M |X| j+1 2 if j even, C M |X| j 2 +1 if j odd, ( 36
)
where ∂ ℓ m denotes differentiation of order ℓ in the m-th coordinate direction, that is, for ℓ = 0 ∂ ℓ m G j = G j , and for ℓ = 1 and m = 1 (resp. ℓ = 1 and m = 2), ∂ ℓ m G j denotes the first derivative of G j in the direction X (resp. in the direction Y ).

Proof. We consider the case ℓ = 0

(∂ ℓ m G j = G j ) first. Let f (X, Z) = i 4 H (1) 0 (k |X|u(Z)) , u(Z) = 1 + Z 2 . ( 37
)
Since, as it is easily seen, we have

G j (X, Y ) = j ℓ=0 C j ℓ (-1) ℓ f X, Y |X| + ℓh |X| ,
to study the decay rate of G j (X, Y ) as X → ∞ we may apply the relation [START_REF] Martensen | Über eine methode zum räumlichen neumannschen problem mit einer anwendung für torusartige berandungen[END_REF] to the function v(Z) = f (X, Z), for a fixed value of X, and with Z = Y X and δ = h |X| . We thus obtain

G j (X, Y ) = h |X| j j ℓ=1 (-1) ℓ C j ℓ ℓ j j! f [j] X, Y |X| + ξ ℓ , 0 ≤ ξ ℓ ≤ ℓh |X| , (38)
where f [j] denotes the Z-derivative of the composite function [START_REF] Moroz | Exponentially Convergent Lattice Sums[END_REF]:

f [j] (X, Z) = i 4 d j dZ j H (1) 0 k |X| u(Z) . (39)
The j-th derivative (39) can be obtained by means of an application of Faà di Bruno's formula [START_REF] Craik | Prehistory of Faà di Bruno's formula[END_REF] for differentiation of composite functions, which in the present context gives

f [j] (X, Z) = i 4 j! m 1 ! . . . m j ! (H (1) 0 ) (m) k|X| u(Z) j q=1 k |X| u (q) (Z) q! mq , (40)
where m = m 1 +• • •+m j , and where the sum is evaluated over all j-tuples of non-negative integers

(m 1 , m 2 , • • • , m j ) satisfying m 1 + 2 • m 2 • • • + j • m j = j.
In view of equations ( 38) and ( 40), to estimate the asymptotics of G j (X, Y ) as X → ∞ it suffices to estimate the corresponding asymptotics for the following two families of quantities:

(a) (H (1) 0 ) (m) k|X| u Y |X| + ξ ℓ , (1 ≤ ℓ, m ≤ j) , and (b) j q=1 k |X| u (q) Y |X| + ξ ℓ mq , (1 ≤ ℓ ≤ j) , j q=1 qm q = j.
It is easy to estimate the quantities (a). Indeed, in view of the well known expression [27, eqn. (5.

6.3)] 2d/dz H (1) n (z) = H (1) n-1 (z) -H (1) n+1 (z), we obtain (H (1) 0 ) (m) k |X| u(Z) = m q=0 c m q H (1) q k |X| u(Z)
and

|H (1) q (t)| ≤ L 1 t 1/2 , q ∈ Z , 0 ≤ q ≤ j , t ≥ 0 , for certain constants c m q and L 1 . Since Y |X| + ξ ℓ ≤ L 2
|X| for some constant L 2 , estimates on the quantities (a) follow: for some constant L 3 we have

(H (1) 0 ) (m) k|X| u Y |X| + ξ ℓ ≤ L 3 |X| 1/2 . (41)
Clearly, L 3 depends on k, L 2 , ℓ and M . In order to estimate the quantities (b), in turn, we first note that for small Z and for all q ≤ j we have

u (q) (Z) ≤ L 4 if q is even, L 4 Z if q is odd (42)
for some constant L 4 . It follows that

j q=1 k |X| u (q) Y |X| + ξ ℓ mq ≤    L 5 |X| j q=1 mq- j/2 q=1 m 2q-1 = L 5 |X| j/2
q=1 m 2q j is even,

L 5 |X| j q=1 mq- (j+1)/2 q=1 m 2q-1 = L 5 |X| (j-1)/2 q=1 m 2q j is odd. (43)
But, as is easily shown we have

j/2 q=1 m 2q ≤ j 2 for j even, (44)
and

(j-1)/2 q=1 m 2q ≤ j -1 2 for j odd. (45)
Thus

j q=1 k |X| u (q) Y |X| + ξ ℓ mq ≤    L 5 |X| j 2
for j even,

L 5 |X| j-1 2 for j odd. (46)
Combining ( 41) and ( 46) we obtain

f [j] X, Y |X| + ξ ℓ ≤    L 6 |X| j 2 -1 2 if j is even, L 6 |X| j 2 -1 if j is odd (47)
and, thus, in view of [START_REF] Nédélec | Integral equation methods in a quasi-periodic diffraction problem for the timeharmonic Maxwell's equations[END_REF], the estimates (36) for the case ℓ = 0 result. The estimates [START_REF] Monro | Windowing-Based Approach to the Evaluation of Scattering from Periodic Rough Surfaces[END_REF] for first derivatives (cases ℓ = 1, m = 1 and ℓ = 1, m = 2) can be obtained in an entirely analogous manner in view of the relation [START_REF] Lebedev | Special Functions And Their Applications[END_REF] d dz H

(1)

0 (z) = -H (1)
1 (z).

The proof is now complete.

Fast convergence of the quasi-periodic

Green-function series [START_REF] Maystre | Rigorous vector theories of diffraction gratings[END_REF] In view of the relations [START_REF] Monro | Windowing-Based Approach to the Evaluation of Scattering from Periodic Rough Surfaces[END_REF] it is easy to estimate the rate of convergence of the series [START_REF] Maystre | Rigorous vector theories of diffraction gratings[END_REF] and its term-wise derivatives. Details in these regards are presented in the following theorem.

e -iαnL ∂ ℓ m G j (X + nL, Y ) ≤        D M N (j-1)/2 for j even, D M N j/2
for j odd.

(

) 48
It follows that for ℓ = 0 and ℓ = 1 with m = 1 or m = 2:

1. The truncated series N n=-N e -iαnL ∂ ℓ m G j (X + nL, Y ) (which, in the case ℓ = 1 result from truncation and term-wise differentiation of the Green-function series [START_REF] Maystre | Rigorous vector theories of diffraction gratings[END_REF]) converge as N → ∞ to the corresponding quantities ∂ ℓ m Gq j , and

The corresponding approximation errors

∂ ℓ m Gq j -N n=-N e -iαnL ∂ ℓ m G j (X + nL, Y
) decrease at least as fast as N -(j-1)/2 for j even and as fast as N -j/2 for j odd.

Proof. Follows directly from Lemma 4.3.

Remark 4.5. The convergence estimates (48) can be strengthened for configurations that do not correspond to Wood anomalies. Indeed, noting that the inequalities in Lemma (4.3) can be expressed as the sum of an asymptotic term (which, as is given by a product of 1/N s , for some s, and a smooth function, is a monotone function of 1/N for N large enough) plus a higher-order correction and mirroring the proof of [START_REF] Bruno | Solution of a boundary value problem for the helmholtz equation via variation of the boundary into the complex domain[END_REF]Thm. 4.1] to estimate a sum between N and P in the present context (for some P larger than our truncation integer N ), the error bound that results from use of [START_REF] Bruno | Solution of a boundary value problem for the helmholtz equation via variation of the boundary into the complex domain[END_REF]Lemma 4.2] as in [START_REF] Bruno | Solution of a boundary value problem for the helmholtz equation via variation of the boundary into the complex domain[END_REF]Thm. 4.1] can be used to obtain the improved error estimate

n∈Z,|n|≥N e -iαnL ∂ ℓ m G j (X + nL, Y ) ≤        D M N (j+1)/2 for j even, D M N j/2+1
for j odd.

(49)

The constant D M , which depends on the configuration under consideration, grows without bound when, as a result of variations in frequency, period or incidence angle, a Wood configuration is approached.

Gq j null-space and complete rapidly-convergent Green function G q j

Following the discussion in Section 2.3 and in view of the introduction of the rapidly convergent Green function Gq j in the previous sections it may be expected that the scattered fields u s D and u s N for the Dirichlet and Neumann problems could be expressed in terms of the rapidly converging Green functions-using, e.g., the representations

u s D (x, y) = Γ # ∂ Gq j ∂ν(x ′ ) x -x ′ , y -f (x ′ ) µ(x ′ ) ds(Γ # ), (50)
u s N (x, y) = Γ # Gq j x -x ′ , y -f (x ′ ) η(x ′ ) ds(Γ # ). (51)
Such representations give the desired solutions for the Dirichlet and Neumann problems presented in Section 2.1 as long as µ and η are solutions of the equations

Γ # ∂ ν(x ′ ) Gq j x -x ′ , f (x) -f (x ′ ) µ(x ′ ) ds(x ′ ) + 1 2 µ(x ′ ) = -u inc |Γ # , (52)
Γ # ∂ ν(x) Gq j x -x ′ , f (x) -f (x ′ ) η(x ′ ) ds(x ′ ) - 1 2 η(x ′ ) = -∂ ν(x) u inc |Γ # . ( 53
)
As it happens, however, in many cases the operators on the left hand sides of these equations are not invertible. More precisely, as shown in what follows, the operators on the left hand sides of Equations ( 52) and (53) are not invertible I. For certain "resonant" values of the shift parameter h (for which non-invertibility arises any frequency k, see Section 4.4.1); and II. At Wood anomalies (for which these operators are not invertible for any value of the shift parameter h, see Section 4.4.2).

Roughly speaking, the difficulties I and II are associated with absence of a (finite) number of Rayleigh modes in the spectral expansion of the quasi-periodic Green function Gq j . Fortunately, however, these difficulties are easily circumvented. We describe these issues in detail in Sections 4.4.1 and 4.4.2 and then, in Section 4.4.3, we put forth a modified version G q j of the Green function Gq j which does not suffer from either of the difficulties I or II.

Gq j null-space I: resonant shift parameter values h c

In view of equations ( 15), ( 23) and [START_REF] Mathis | A Comparison Of Acceleration Procedures For The Two-Dimensional Periodic Green's Function[END_REF], away from Wood anomalies we have

Gq j (X, Y ) = n∈Z i 2Lβ n j m=0 (-1) m C j m e iβnmh e iαnX e iβnY = n∈Z i 2Lβ n (1 -e iβnh ) j e iαnX e iβnY , (54)
so that, given an integer n 0 , for values of h = h c for which the resonance condition

(1 -e iβn 0 hc ) j = 0 (55) holds, the Green function Gq j (X, Y ) as a function of X (for any fixed value of Y ) does not contain the Rayleigh mode e iαn 0 X : the Rayleigh coefficient of order n 0 of the function Gq j (X, Y ) equals zero. In such cases the operators on the left hand sides of equations ( 52) and (53) are not invertible. Indeed, assuming equation (52) (resp. (53)) is invertible, let us consider the solution µ (resp. η) that results from a right-hand side given by -u inc (x, f (x)) (resp. given by -∂ ν(x) u inc (x, f (x))) with u inc (x, y) = e iαn 0 x+iβn 0 y (note that the exponent here differs from the one in equation ( 3)). By uniqueness of solution of the Helmholtz equation (see Remark 2.4) we would then have u s D = -e iαn 0 x+iβn 0 y (resp. u s N = -e iαn 0 x+iβn 0 y ) throughout Ω-which is clearly not possible since, as mentioned above, for any fixed Y the Rayleigh expansion of Gq j (X, Y ) does not contain the Rayleigh mode e iαnX for any value of the frequency k.

Gq

j null-space II: Wood anomalies Similar resonance conditions and associated non-uniqueness issues occur at Wood anomalies-as demonstrated in what follows. Throughout this section we use the notation α n = α n (k, θ) and β n = β n (k, θ) to display explicitly the (k, θ) dependence of the wave-numbers α n and β n , respectively. We assume k = k 0 and a given incidence angle θ = θ 0 gives rise to a Wood configuration for the n-th Rayleigh mode-that is, β n (k 0 , θ 0 ) = 0. For the present demonstration purposes, in this section we further assume that for the given configuration there is one and only one value of n ∈ Z, namely n = n 0 , for which β n (k 0 , θ 0 ) vanishes:

β n 0 (k 0 , θ 0 ) = 0 and β n (k 0 , θ 0 ) = 0 for n = n 0 .
Then, it is easy to check that, as (k, θ) approaches (k 0 , θ 0 ), the limit of the expansion (54) for G q j is given by

Gq j (X, Y ) = n∈Z,n =n 0 i 2Lβ n (k 0 , θ 0 ) j m=0 (-1) m C j m e imβn(k 0 ,θ 0 )h e iαn(k 0 ,θ 0 )X e iβn(k 0 ,θ 0 )Y + lim β→0 i 2Lβ j m=0 (-1) m C j m e imβh e iα n 0 (k 0 ,θ 0 )X e iβY .
But, using the Taylor expansion of e imβh around β = 0 we obtain

lim β→0 i 2Lβ j m=0 (-1) m C j m e imβh = lim β→0 i 2Lβ (1 -e iβh ) j = 0 for j ≥ 2, (56)
and, thus, at the Wood anomaly configuration, Gq j does not contain the Rayleigh mode e iα n 0 (k 0 ,θ 0 )X . Thus, using an argument similar to the one presented in Section 4.4.1 for resonant-shift cases, here we find that the operators on the right-hand sides of equations ( 52) and (53) are not invertible for any value of the shift parameter h. Remark 4.6. Figures 6 and7 demonstrate the lack of invertibility of integral equation formulations based on the quasi-periodic Green function Gq j at resonance values of the shift parameter and for Wood anomaly configurations. Note that equation (55) can only hold for propagative modes (n 0 ∈ U ) and, in particular, there can only be a finite number of values of n for which a resonance of type (55) may occur.

Remark 4.7. Not only do the two non-uniqueness issues discussed in Sections 4.4.1 and 4.4.2 translate into noninvertible operators at resonant-shift and Wood-anomaly cases, but they also give rise to ill-posed numerical solvers around both, Wood Anomalies and resonant-shift values. Note that the non-uniqueness problem at Wood anomalies (which, incidentally, occurs for j ≥ 2 but not for j = 1, cf. equation ( 56)) is in some sense more fundamental that the non-uniqueness issue arising from resonant-shift condition (55): the latter could be bypassed by adequate selection of the somewhat arbitrary shift parameter h; the former, however, presents an impediment for evaluation of scattering solutions at and around physically realizable Wood configurations.

Complete rapidly-convergent Green function G q j

The present section introduces slightly modified Green functions and corresponding integral equations that overcome the null-space difficulties described in the previous two sections. As detailed in what follows, the modified Green functions are obtained by simply adding to Gq j a multiple of each one of the Rayleigh modes e iαnX+iβnY that it does not possess. Clearly, such a procedure results in new outgoing Green functions, since the individual modes added are themselves outgoing solutions of the Helmholtz equation ( 4) in the domain Ω. In particular, our algorithm does not resort to selection of non-resonant shift parameters (cf. Remark 4.7).

To address the non-uniqueness issues arising from the resonance condition (55) we proceed as follows: letting η r > 0 and

U ηr r = n ∈ U such that (1 -e iβnh ) j < η r (57)
we define a smooth Kernel M r by

M r (X, Y ) = n∈U ηr r e iαnX+iβnY , X ∈ (-L, L), Y ∈ R, (58)
which will be added to Gq j to eliminate resonant-kernel non-uniqueness in equations ( 52) and (53). Similarly, to eliminate Wood anomaly non-uniqueness in the integral equations ( 52) and (53), for η w > 0 we define

U ηw w = n ∈ Z, such that |(1 -e iβnh ) j | ≤ η w ( 59
)
and for n ∈ U ηw w we let

c n = max x∈(-L 2 , L 2 ) |iα n f ′ (x) -iβ n |. (60)
Then, noting that c n is necessarily non-zero unless f ′ vanishes identically, for non-constant f we consider the smooth kernel M w ,

M w (X, Y ) = n∈U ηw w e iαnX+iβnY c n , X ∈ (-L, L), Y ∈ R. ( 61
)
Remark 4.8. For the (trivial) flat case (f ′ = 0, constant f ) at Wood anomalies we have c n = 0 for some n, and thus M w is not defined. The needed kernels ∂ ν(x ′ ) M w and ∂ ν(x) M w can nevertheless be extended to such cases using a limiting process, and it can thus be seen that integral equations and PDE solutions vary continuously as the flat-surface case is approached (except for the Neumann problem at horizontal incidence (|α| = k) which, for the flat surface, does not admit a unique PDE solution). Indeed for f ′ identically equal to zero we have c n = |β n |, and thus

lim βn→0 ∂ ν(x ′ ) e iαn(x-x ′ )+iβn(y-y ′ ) |β n | = -ie iαn(x-x ′ ) lim βn→0 β n |β n |
Noting that lim βn→0 βn |βn| equals either 1 or i (depending on the sign of (k 2 -α 2 n ) as the Wood anomaly is approached), for the flat case at Wood anomalies we replace ∂ ν(x ′ ) M w (x -x ′ , y -y ′ ) and ∂ ν(x) M w (x -x ′ , y -y ′ ) by one of the two possible associated limits, say,

∂ ν(x ′ ) M w (x -x ′ , y -0) → -i n∈U ηw w e iαn(x-x ′ ) and
(62)

∂ ν(x) M w (x -x ′ , y -0) → i n∈U ηw w e iαn(x-x ′ ) . ( 63
)
It is easy to check that these alternative definitions for the flat case at Wood anomalies (or, indeed, any definition as a sum of arbitrary non-zero multiples of all exponentials e iαn(x-x ′ ) for n ∈ U w ) all give rise to one and only one integral-equation solution for the flat-interface Wood-anomaly case-as long as the right-hand-side of the integral equation is orthogonal to e iαnx for all n ∈ U ηw w . Under the latter condition and using considerations concerning uniqueness properties (Remark 2.4), further, it can be checked [START_REF] Bruno | All-frequency Green's functions, Wood anomalies and uniqueness proofs in periodic rough-surface scattering[END_REF] that these definitions lead to solutions µ and η of the integral equations as well as scattered fields u s D and u s N that vary continuously as the flat-interface/Woodanomaly case is approached.

Using these additional kernels we can now define a periodic Green function for all smooth periodic interfaces with f ′ = 0, all incidence angles and all frequencies, namely

G q j (X, Y ) = Gq j (X, Y ) + M r (X, Y ) + M w (X, Y ); (64)
with the special considerations specified in Remark 4.8, the methodology is extended naturally, with continuous solutions, to the trivial flat case f ′ = 0. The resulting second kind integral equations for the Dirichlet and Neumann are given by

1 2 µ(x) + Γ # ∂ ν(x ′ ) G q j (x -x ′ , f (x) -f (x ′ ))µ(x ′ ) ds(x ′ ) = -u inc |Γ # and (65)
- 1 2 η(x) + Γ # ∂ ν(x) G q j (x -x ′ , f (x) -f (x ′ ))η(x ′ ) ds(x ′ ) = -∂ ν(x) u inc |Γ # . (66)
The Dirichlet and Neumann scattering problems under consideration for non-flat rough surfaces thus reduce to finding density functions µ ∈ C α # (Γ # ) and η ∈ C α # (Γ # ) such that equations ( 65) and (66) are satisfied, where C α # (Γ # ) denotes the space of continuous α-quasi-periodic functions. A general discussion of existence and uniqueness for these problems that extends previous results (cf. Remark 2.4) to Wood anomaly cases will be presented elsewhere [START_REF] Bruno | All-frequency Green's functions, Wood anomalies and uniqueness proofs in periodic rough-surface scattering[END_REF]. Numerical algorithms based on these equations and corresponding numerical results are presented in Sections 5 and 6 below.

Numerical algorithm

Our numerical method for the solution of equations ( 65) and (66) incorporates the windowing and shifting methodologies introduced in Sections 3 and 4 for evaluation of integral operators at all frequencies together with a suitable modified version of a well known Nyström approach for high-order evaluation of logarithmic integral operators. This section presents the details of a numerical implementation of the combined windowing-shifting based Nyström approach.

Kernel decomposition in a periodic setting

The Martensen-Kussmaul (MK) Nyström approach [START_REF] Colton | Inverse Acoustic And Electromagnetic Scattering Theory[END_REF][START_REF] Kussmaul | Ein numerisches verfahren zur lösung des neumannschen außenraumproblems für die helmholtzsche schwingungsgleichung[END_REF][START_REF] Martensen | Über eine methode zum räumlichen neumannschen problem mit einer anwendung für torusartige berandungen[END_REF][START_REF] Nyström | Über die praktische auflösung von integralgleichungen mit anwendungen auf randwertaufgaben[END_REF] relies on decomposition of the integral kernel as a sum of a smooth kernel and a logarithmic term. We present our decomposition for the Dirichlet equation (65); the method for Neumann equation ( 66) is entirely analogous.

Like previous MK approaches, our MK Nyström method is based on use of the periodic unknowns: using the change of unknown μ(x) = µ(x)e -iαx we re-express the integral on the left-hand side of equation (65) in the form

L 2 -L 2 ∂ ν(x ′ ) G q j (x -x ′ , f (x) -f (x ′ )) φ(x, x ′ ) μ(x ′ ) 1 + f ′ (x ′ ) 2 dx ′ + 1 2 μ = -u inc |Γ # e -iαx , (67)
where μ is a periodic function of period L and where φ(x, x ′ ) = e -iα(x-x ′ ) . In view of (64) and expressing the integral involving the periodic Green function Gq j as an integral between -∞ and ∞ (cf. equation [START_REF] Graham | Concrete mathematics: A foundation for computer science[END_REF] where the corresponding identity for the kernel G q is displayed) the integral on the left-hand side of (67) can be expressed in the form

L 2 -L 2 ∂ ν(x ′ ) {M r + M w } (x -x ′ , f (x) -f (x ′ )) φ(x, x ′ ) μ(x ′ ) 1 + f ′ (x ′ ) 2 dx ′ + +∞ -∞ ∂ ν(x ′ ) G j (x -x ′ , f (x) -f (x ′ )) φ(x, x ′ ) μ(x ′ ) 1 + f ′ (x ′ ) 2 dx ′ . ( 68
)
The integrand in the first term of equation ( 68) is smooth and periodic, and, therefore, the corresponding integral operator can be evaluated with high-order accuracy by means of the trapezoidal rule [START_REF] Colton | Inverse Acoustic And Electromagnetic Scattering Theory[END_REF]. To compute the second term in (68), in turn, we first truncate the integral using the windowing technique introduced in Section 3. That is we approximate the second integral (68) by

+∞ -∞ S(x -x ′ , cA, A) ∂ ν(x ′ ) G j (x -x ′ , f (x) -f (x ′ )) φ(x, x ′ ) μ(x ′ ) 1 + f ′ (x ′ ) 2 dx ′ (69)
for adequately selected positive constants c and A, and we then evaluate the integral (69), with high-order accuracy, by accounting explicitly for the logarithmic singularity at (x -x ′ ) = 0. To do this we follow [START_REF] Colton | Inverse Acoustic And Electromagnetic Scattering Theory[END_REF] and note that, for |x -x ′ | ≤ 2π, the kernel ∂ ν(x ′ ) G j can be expressed in the form

∂ ν(x ′ ) G j (x -x ′ , f (x) -f (x ′ )) = K 1 (x, x ′ ) ln 4 sin 2 x -x ′ 2 + K 2 (x, x ′ ), (70)
where

K 1 (x, x ′ ) = k 4π f (x ′ )(x -x ′ ) -(f (x ′ ) -f (x)) 1 + f ′ (x ′ ) 2 J 1 k (x -x ′ ) 2 + (f (x) -f (x ′ )) 2 (x -x ′ ) 2 + (f (x) -f (x ′ )) 2 ,
and

K 2 (x, x ′ ) = ∂ ν(x ′ ) G j (x -x ′ , f (x) -f (x ′ )) -K 1 (x, x ′ ) ln 4 sin 2 x -x ′ 2 ,
are smooth functions.

In the context [START_REF] Colton | Inverse Acoustic And Electromagnetic Scattering Theory[END_REF] the condition |x -x ′ | ≤ 2π is always satisfied; unfortunately, however, this is not the case for the integration problem (69): under the present conditions, use of the decomposition (70) for -∞ < x -x ′ < ∞, would give rise to singularities whenever x-x ′ is an integer multiple of 2π. To overcome this difficulty we utilize an additional cut-off function that vanishes outside a sufficiently small local neighborhood of the singular point x ′ = x: taking 0 < A ℓ < L/2, 0 < d ℓ < 1 and recalling equation ( 17) we define the "local" windowing function

P ℓ (x -x ′ , d ℓ , A ℓ ) = S(x -x ′ , d ℓ A ℓ , A ℓ ), d ℓ , A ℓ ∈ R >0 , (71)
and we decompose the integrand in (69) in the form

S(x -x ′ , cA, A)∂ ν(x ′ ) G q j (x -x ′ , f (x) -f (x ′ )) φ(x, x ′ ) 1 + f ′ (x ′ ) 2 = K ℓ (x, x ′ ) ln 4 sin 2 π L (x -x ′ ) + K r (x, x ′ ) (72)
where

K ℓ (x, x ′ ) = k 4π P ℓ (x -x ′ , d ℓ , A ℓ ) f (x ′ )(x -x ′ ) -(f (x ′ ) -f (x)) (x -x ′ ) 2 + (f (x) -f (x ′ )) 2 J 1 k (x -x ′ ) 2 + (f (x) -f (x ′ )) 2 , ( 73
)
and where

K r (x, x ′ ) = S(x -x ′ , cA, A)∂ ν(x ′ ) G q j (x -x ′ , f (x) -f (x ′ )) φ(x, x ′ ) 1 + f ′ (x ′ ) 2 -K ℓ (x, x ′ ) ln 4 sin 2 π L (x -x ′ ) . (74)
It is easy to check that, given the assumptions 0 < A ℓ < L/2 and 0 < d ℓ < 1, K ℓ and K r are smooth functions of x and x ′ , and, thus, the MK methodology can be used for evaluation of the integral (68)-as detailed in the following section.

MK quadrature rules

To derive MK quadrature rules in the present context we consider the decomposition (72) and we re-express the integral (68) in the form

L 2 -L 2 K per r (x, x ′ )μ(x ′ ) dx ′ + +∞ -∞ K r (x, x ′ )μ(x ′ ) dx ′ + +∞ -∞ K ℓ (x, x ′ ) ln 4 sin 2 π L (x -x ′ ) μ(x ′ ) dx ′ , ( 75
)
where

K per r (x, x ′ ) = ∂ ν(x ′ ) {M r + M w } (x -x ′ , f (x) -f (x ′ )) φ(x, x ′ ) 1 + f ′ (x ′
) 2 is a smooth and L-periodic kernel; note that all integrands in the (-∞, ∞) integrals in (75) are compactly supported functions with supports contained in the interval (x -A, x + A). We additionally use an equispaced discretization mesh containing an even number n I of points per period of the scattering surface: denoting by N = ⌈ A L ⌉ the number of periodic intervals contained in the integration domain to the right of the point x = L 2 (where ⌈x⌉ denotes the smallest integer larger than or equal to x), the quadrature rule uses a total of (2N + 1)n I quadrature points x j given by

x j = (- L 2 -N L) + (j -1) L n I , j ∈ N, j ≤ (2N + 1)n I ,
over an integration domain consisting of 2N + 1 periodic intervals.

Our MK quadrature rule for |x| < L/2 then results as follows: The first integral in equation ( 75) is evaluated with high-order accuracy by means of the trapezoidal rule formula

L 2 -L 2 K per r (x, x ′ ) μ(x ′ ) dx ′ ≈ L n I -L 2 ≤|x-x j |≤ L 2 K add (x, x j ) μ(x j ). (76)
Similarly, the second term in (75) is integrated with high-order accuracy by means of the trapezoidal rule expression

+∞ -∞ K r (x ′ , x)μ(x ′ )dx ′ ≈ L n I |x-x j |≤A K r (x, x j ) μ(x j ). (77)
Finally, the third term in (75) (the singular integral) is approximated by means of the MK quadrature rule

+∞ -∞ K ℓ (x, x ′ ) ln 4 sin 2 π L (x -x ′ ) μ(x ′ )dx ′ = x+ L 2 x-L 2 K ℓ (x ′ , x) ln 4 sin 2 π L (x -x ′ ) μ(x ′ )dx ′ ≈ L 2π |x-x j |≤A ℓ R j,n I 2πx L K ℓ (x, x j ) μ(x j ), (78)
where

R j,n I (t) = 1 n I n I 2 q=- n I 2 +1 2π 0 e iq( 2πx j L +τ -t) ln 4 sin 2 τ 2 dτ. (79)
Formula ( 78) is obtained by introducing the change of variables t = 2π L x, τ = 2π L x ′ and expressing K ℓ (x, x ′ ) μ(x ′ ) as a Fourier series. Using the relation [24, ch. 12]

1 2π 2π 0 ln 4 sin 2 τ 2 e imτ dτ = 0 if m = 0 - 1 m m = 1, 2, • • • we obtain R j,n I 2π x i L = - 4π n I   n I 2 q=1 1 |q| cos q 2π L (x i -x j )   - 4π n 2 I cos n I 2 2π L (x i -x j ) (80)
and our basic MK quadrature rule for the integral (75) is complete. The next two sections describe certain accuracy and efficiency difficulties that arise from direct application of the MK algorithm described in this section, and they put forth modified procedures that eliminate such difficulties. Thus, Section 5.3 introduces an inexpensive local quadrature method for nearly-singular integrals that arise for small values of the shift parameter h (see equation ( 33)), and Section 5.4 presents an asymptotic method that can be used to accelerate the computation of the large numbers of large-argument Hankel functions that are needed in our Greenfunction procedure. As demonstrated in Section 6, once these elements are introduced a highly accurate and efficient all-frequency algorithm results.

Accurate integration around "nearly singular points" in shifted Green functions

It is useful to note that, as far as the decay of the Green function [START_REF] Mathis | A Comparison Of Acceleration Procedures For The Two-Dimensional Periodic Green's Function[END_REF] and associated rapid convergence of the quasiperiodic Green function [START_REF] Kirsch | Uniqueness theorems in inverse scattering theory for periodic structures[END_REF] are concerned, it is advantageous to use small values of the shift parameter h: small shift values give rise to correspondingly smaller values of the Green-function expression [START_REF] Mathis | A Comparison Of Acceleration Procedures For The Two-Dimensional Periodic Green's Function[END_REF] for a given value of X (see e.g. equation [START_REF] Nédélec | Integral equation methods in a quasi-periodic diffraction problem for the timeharmonic Maxwell's equations[END_REF]) and, thus, to accurate approximation of infinite integrals via integration over reduced numbers of periods. Unfortunately, however, although the kernel K r defined in (74) is smooth and periodic for arbitrarily small values of h, a direct application of the trapezoidal quadrature rule for the second term in (75) with a mesh-size that is no finer than is needed for resolution of the wavelength may yield inaccurate results when the surface Γ and its downward h-shifted counterpart Γ -h (see equation [START_REF] Mathis | A Comparison Of Acceleration Procedures For The Two-Dimensional Periodic Green's Function[END_REF]) are close to each other. This accuracy loss is easily understood. For example, for small values of h and for a given x ∈ (0; L), the kernel in the shifted integral

x i +A x i -A ∂ ν(x ′ ) H (1) 0 k (x i -x ′ ) 2 + (f (x i ) -f (x ′ ) + h) 2 µ(x ′ ) dx ′ , A > 0, (81)
which is part of the integrand in (75), is nearly singular (with a logarithmic near-singularity) for values of x ′ in a neighborhood of x ′ = x. Clearly, further, the distance between Γ and Γ -h is smaller for steep surfaces Γ than it is for gently-sloped surfaces; the situation is illustrated on the left portion of Figure 5.3: the small distance between Γ and Γ -h at regions where Γ is steep can be easily appreciated.

To reestablish accuracy we smoothly partition the integration domain and we use a finer mesh in a region around the observation point x ′ = x. Thus, using a smooth cut-off function χ of the general form [START_REF] Craik | Prehistory of Faà di Bruno's formula[END_REF], for the test case Note that the two surfaces are very close to each other at x = π: the minimum distance between the two surfaces is 0.24. Right. Grid refinement: coarse grid x i and fine grid y j .

considered above in this section, for example, we write

x+A x-A ∂ ν(x ′ ) H (1) 0 k (x -x ′ ) 2 + (f (x) -f (x ′ ) + h) 2 µ(x ′ ) dx ′ = x+A x-A ∂ ν(x ′ ) H (1) 0 k (x -x ′ ) 2 + (f (x) -f (x ′ ) + h) 2 (1 -χ(x -x ′ )) µ(x ′ ) dx ′ + x+A x-A ∂ ν(x ′ ) H (1) 0 k (x -x ′ ) 2 + (f (x) -f (x ′ ) + h) 2 χ(x -x ′ ) µ(x ′ ) dx ′ . ( 82
)
The first integral on the right hand side of this equation, which is not nearly singular provided the function χ on the right-hand side of this equation is adequately selected, is accurately approximated by means of the classical trapezoidal rule on the coarse mesh x i . For the second right-hand integral, in turn, we use a fine equispaced quadrature mesh {y ℓ }. The necessary coarse-grid to fine-grid interpolations are produced, in our algorithm, by evaluation of the discrete Fourier transform for the discrete version of the periodic quantity µ(x ′ )e -iαx ′ on the coarse mesh x i and subsequent evaluation of the resulting finite Fourier series on the fine mesh y ℓ , as detailed below in this section.

A few comments are in order with regards to the selection of the windowing function χ, the fine-integration region and the associated fine mesh. We first note that while the (nearly singular) integral containing the factor χ needs to be approximated by means of a fine integration mesh, the integral containing the factor (1 -χ) can be produced accurately by means of the coarser background integration mesh-provided the windowing function χ has a sufficiently broad support and a sufficiently mild slope. It is therefore not advisable to use a very narrow fine-integration region around x ′ = x: doing so would lead to a steep function χ, and, thus, to low accuracies arising from the coarse mesh integration of the integral containing the factor (1 -χ). We have therefore found it convenient to utilize a windowing function χ whose support spans several periods, and which, in particular, can be selected independently of x. We thus take our fine integration mesh to be given by

y ℓ = -QL + (ℓ -1)∆ y , ∆ y = L n f , 1 ≤ ℓ ≤ 2 Q n f
(right portion of Figure 5.3), where n f and Q are positive integer numbers. (Clearly n f the number of points per period in the fine integration mesh and 2Q is the number of periods in the fine-integration region.) Of course, n f must be taken to be larger than n I in order to make the near-singular integration method accurate for a given profile Γ and a given shift parameter h. Additionally, Q can and should be taken to be much smaller than the number N in the overall integration region, to avoid excessive fine-integration costs. In practice we have found the value Q = 4 is small enough not to give rise to significant additional computing cost as a result of the fine-integration procedure, and sufficiently large to enable use of an adequately gently sloped windowing function χ.

The second integral on the right hand side of equation ( 82) is thus approximated by trapezoidal-rule integration based on the fine grid {y ℓ }:

x+A x-A ∂ ν(x ′ ) H (1) 0 k (x -x ′ ) 2 + (f (x) -f (x ′ ) + h) 2 χ(x -x ′ ) µ(x ′ ) dx ′ ≈ ∆y 2 Q n f ℓ=1 ∂ ν (x ′ )H (1) 0 k (x -y ℓ ) 2 + (f (x) -f (y ℓ ) + h) 2 χ(x -y ℓ )µ(y ℓ ).
Since µ is a smooth and periodic function, the needed fine-mesh values µ(y ℓ ) can be obtained via fine-mesh evaluation of a truncated Fourier series of µ:

µ(y ℓ ) ≈ n I /2 p=-n I /2+1
µ p e 2ipπy ℓ /L with µ p ≈ ∆x

n I /2 q=-n I /2+1 µ(x q ) e -2iqπxq/L . ( 83
)
Remark 5.1. In practice we have found that, even though the procedure presented in this section gives rise to accurate integration for arbitrarily small values of the shift parameter h, use of values of h that are smaller than certain thresholds make the discrete versions of the integral equations ( 65) and (66) ill-conditioned and thus result in accuracy limitations. The parameter h should therefore be selected by means of test runs designed to optimize efficiency for a prescribed error tolerance; see also Remark 5.2.

Remark 5.2. In view of the limited expense associated with the singularity treatment algorithm introduced in this section, a simple prescription can be given for selection of the fine mesh-size value ∆y introduced in this section, namely, to use ∆ y = ∆x R (or, equivalently, n f = Rn I ), for some moderate value of R, in the prescribed number Q of nearby periods. These parameters can be chosen by means of integration convergence analyses involving explicit densities µ (e.g. µ equal to the values on Γ of the incident field). In practice we have found that for h large enough no such refinement is necessary. In the examples considered in this paper we have used 0.8 < h < 3; for lower values the ill conditioning difficulty mentioned in Remark 5.1 becomes noticeable, specially for deep gratings, and makes the accuracy levels demonstrated in Section 6 difficult to obtain. We have found that, for h = 0.8, for example, the values R = 8 and Q = 4 adequately resolve the near singularity, and that higher values of R do not give rise to improved accuracies in the numerical integration and solution processes.

Efficient evaluation of G q j via accurate asymptotic expansions

In view of our reliance on the half-space Green function G j for evaluation of the complete quasi-periodic Green function G q j (equations ( 64), ( 23) and ( 33)), the computation of G q j and its derivatives in principle requires evaluation of Hankel functions for a number of arguments-which could require relatively long computing times. Additionally, since the rapid decay of the half-space Green function G j is derived from cancellation of Hankel-function values that decay less rapidly than G j itself, a direct evaluation of the G j gives rise to cancellation errors which can be noticeable in some cases (see Figures 4 and5). Both of these difficulties can be addressed effectively via consideration of asymptotic expressions for the Green function G j , as indicated in what follows.

We consider the case j = 1 first. Letting

z n = (X + nL) 2 + Y 2 we obtain G 1 (X + nL, Y ) = H (1) 0 k (X + nL) 2 + Y 2 -H (1) 0 k (X + nL) 2 + (Y + h) 2 . ( 84
)
Using the large z asymptotic expression [27, eqn. 5.11.15]

H (1) 0 (z) = e iz √ z h 0 + h 1 z + h 2 z 2 + • • • + h m z m + O(z -(m+1) ) (85)
(where h r are constant complex coefficients given by simple expressions) together with the relation

(X + nL) 2 + (Y + h) 2 = z n 1 + g z 2 n 1/2 , (g = 2Y h + h 2 ),
and expanding (84) in powers of 1/z n , we find that the term involving the constant h 0 cancels exactly (thus eliminating numerical cancellation errors) and we obtain

G 1 (X + nL, Y ) = e izn √ z n d 1 z n + d 2 z 2 n + d 3 z 3 n + ... + d m z m n + O(z m+1 n ) (86)
for certain explicitly computable coefficients d r which depend on g and the coefficients h r .

A similar procedure can be applied for the evaluation of the shifted Green function G j of any order j; for general j we obtain an expansion of the form

G j (X + nL, Y ) = e izn √ z n d j 1 z J n + d j 2 z J+1 n + d j 3 z J+2 n + ... + d j m z J+m n + O(z J+m+1 n ) (87)
where J = j/2 for even j and J = (j+1)/2 for odd j (cf. Lemma 4.3). As in the case j = 1, the coefficients d j r can be evaluated explicitly in terms of the parameter g. The algebra becomes more and more substantial as j grows, but the calculations can be greatly facilitated by application of a computer algebra system such as Maple or Mathematica: what is needed is, simply, evaluation of finite order Taylor expansions in powers of the small parameter γ = 1/z n for the quantity

√ z n e izn G asym j (X + nL, Y ),
where G asym j is the asymptotic expansion of G j that results from use of the expansion (85) in the expression [START_REF] Mathis | A Comparison Of Acceleration Procedures For The Two-Dimensional Periodic Green's Function[END_REF].

Our numerical examples utilize values d j r for r between r = 1 and r = 4. The evaluation of the periodic Green function (64) can thus be achieved accurately and at a substantially reduced cost by substituting G j (X + nL, Y ) in ( 23) by a sum of the first m terms in an expansion of the form (87) for sufficiently large values of n. It is interesting to note that the expansion (86) can be used for all n large enough, and, thus, the cost of evaluation of the expansion coefficients d j r is readily amortized and, indeed, it becomes insignificant when compared to the cost of other portions of the Green-function evaluation.

Remark 5.3 (Algorithm for selection of parameters in the expansion (87)). The expansion order m used in equation (87) for a given range of values of n should be selected to meet a prescribed error tolerance ε in the final solution while minimizing computational cost. The parameter m can thus be approximately selected, for each n, as the smallest integer for which the integral of the error term in (87) is of the order of the error tolerance,

+∞ 0 O(z J+m+1 n ) < ε,
for all relevant values of Y . Estimating this remainder by the term of order J + m + 1 in expansion (87) and letting A j m+1 = sup (x,x ′ )∈(0,L) 2 |d j m+1 |, the "negligible-remainder" condition may be expressed in the form

+∞ nL A j m+1 (x -x ′ ) 2 + (f (x) -f (x ′ )) 2 J+m+3/2 dx ′ < ε.
Approximating, further, (x

-x ′ ) 2 + (f (x) -f (x ′ )) 2 ≈ |(x ′ -L)|, the negligible-remainder condition becomes A j m+1 ((n -1)L) J+m+1/2 (J + m + 1/2) < ε or, for large n, n > 1 L A j m+1 ε(J + m + 1/2) 1 J +m+1/2 . (88)
In practice we have used m ≤ 3 in order to avoid complex closed-form expressions that arise for the coefficients d j r for r large. (Since the value d j m+1 is needed to evaluate n in the negligible-remainder condition (88), our numerical examples have thus utilized values of d j r with 1 ≤ r ≤ 4.) Under this restriction our parameter-selection algorithm uses three cutoff values n 3 , n 2 and n 1 , each one of which is defined as the maximum between a certain "asymptotic cutoff" integer n min and one plus the integer part of the value formula on the right hand side of (88) with m = 3, 2 and 1, respectively. (Noting that our asymptotic analysis is valid as long as the regime for which the asymptotic expressions used are accurate, the integer n min is used to insure that the asymptotic regime has indeed been reached.) In terms of these cutoff integers, the needed values of the function G j (X + nL, Y ) are obtained as follows:

-For n < n 3 , G j (X + nL, Y ) is evaluated by means of the exact formula [START_REF] Kirsch | Uniqueness theorems in inverse scattering theory for periodic structures[END_REF].

-For n 3 ≤ n < n 2 , G j (X + nL, Y ) is evaluated by means of the approximation (87) with m = 3.

-For n 2 ≤ n < n 1 , G j (X + nL, Y ) is evaluated by means of the approximation (87) with m = 2.

-For n ≥ n 1 , G j (X + nL, Y ) is evaluated by means of the approximation (87) with m = 1.

For all of our experiments we used the asymptotic cutoff value n min = 10.

Pseudo-code

The following pseudo-code summarizes the algorithm described in sections 5.1 through 5.4.

Initialization of physical and algorithmic parameters.

1.1 Input of physical parameters: wave number k, incidence angle θ, grating profile function f and polarization (Dirichlet or Neumann boundary conditions).

1.2 Input of basic algorithmic parameters: j, h, n I , m and N : Green-function order, shift, number of discretization points per period, number of terms used in the asymptotic expansion (87), and number of periods used in the periodic Green-function summation, respectively. Parameters should be selected so as to ensure the overall solution meets a prescribed error tolerance ε; guidance in these regards is provided in Remarks 5.1, 5.2 and 5.3.

1.3 Input of other necessary algorithmic parameters mentioned in point 2. below.

2. Construction of the system matrix for the Dirichlet equation (65) or the Neumann equation (66), as required per Point 1.1 above.

2.1 Detection of near-Wood configurations and near-resonant shift parameter values per equations ( 57) and (59) for prescribed tolerance parameters η w and η r (for all of the numerical experiments considered in this paper the values η r = η w = 0.2 were used). Inclusion of corresponding Green-function correction terms M w and M r if near-Wood or near-resonant conditions occur.

2.2 Detection of near singularities in accordance with Remark 5.2 and, if existent, construction of the matrix associated with the nearly singular kernel, as described in Section 5.3.

2.3

Construction of the matrix associated with the singular kernel (73) using the MK quadrature rule (78).

2.4 Construction of the matrix associated with the regular kernel (74):

-Evaluation of the coefficients d j r of the asymptotic expansion of the Green function (87) for r ≤ m; in all the examples in Section 6 we have used m = 4.

-Evaluation of the needed values of the regular Kernel K r (defined in equation ( 74)) via the asymptotic expansion (87) or the exact expression (74), following the prescriptions in Remark 5.3.

3. Evaluation of discrete right-hand-side values by means of a standard discretization of source terms in (65)-(66).

4. Solution of the linear system resulting from points 2. and 3. above by means of Gaussian elimination (LU decomposition and back-substitution).

5. Post-processing. Evaluation of scattering efficiencies and defect in energy balance (by substitution of a spectral expansion of the Green function G q j arising from ( 54) and (64) into the integral expression of the scattered field) and, if required, the scattered field at prescribed points in space.

Numerical results

This section presents numerical tests and examples that demonstrate the properties of the rapidly convergent Green functions and scattering solvers introduced in previous sections-including at and near Wood-anomaly configurations. For definiteness the examples presented in this section concern problems of scattering by two families of

for various values of L and H; see Figure 3. The profile y = f 1 (x) is the prototypical example of a rough surface; with the present conventions these are sinusoids of period L and height-to-period ratio H with peak-to-trough height equal to HL. The multi-frequency profile y = f 2 (x), in turn, is used in this section to demonstrate the new solvers for a somewhat more generic rough surface. The examples in this section include solutions for a range of frequencies and incidence angles. All solutions presented in this section were produced by means of a Fortran implementation of the methods and algorithms described in Sections 4 and 5 together with the LAPACK Gaussian elimination routines for solution of matrix equations, on a single core of a Xeon X 5670 processor. Detailed numerical results presented in the following sections concern Green-function convergence (Section 6.2), effect of Green-function completion terms (Section 6.3), convergence of forward-maps, solutions and efficiencies (Section 6.4), computing costs (Section 6.5), and frequency-efficiency curves under Littrow mounts (Section 6.5). As a preface to these numerical tests, Section 6.1 summarizes theoretical results on convergence rates for Green functions, forward maps, solutions and efficiencies that were presented in (or follow from) the theoretical results presented in Section 4.

Convergence-rate summary: Green functions, forward maps, solutions and efficiencies

The following list summarizes the convergence properties of the various quantities of interest at and away from Wood-anomaly configurations.

1. Decay order exhibited by the half-plane Green function G j (X, Y ) and its first derivatives as X → ∞ (Lemma 4.3): j+1 2 for j even and j 2 + 1 for j odd.

2. Convergence order of the -N ≤ n ≤ N truncated series for the rapidly convergent quasi-periodic Green function Gq j (X, Y ) and its derivatives-2.1 -at Wood-anomaly configurations (Theorem 4.4):

j-1 2 for j even and j 2 for j odd;

2.2 -away from Wood-anomaly configurations (Remark 4.5): j+1 2 for j even and j 2 + 1 for j odd.

3. Convergence order as N → ∞ of the (continuous) forward map-3.1 -at Wood-anomaly configurations: the same as the convergence orders shown in item 2.1 above.

3.2 -away from Wood-anomaly configurations: super-algebraically fast, in spite of the result in item 2.2 above, and in view of the super-algebraic convergence induced in the integration process by the smooth windowing methodology described in Section 3.

4. Convergence order as N → ∞ of individual efficiencies-4.1 -at Wood-anomaly configurations: as shown in items 2.1 and 3.1 for j odd, and one order higher than the ones shown items 2.1 and 3.1 for j even (Remark 6.1), that is to say j+1 2 for j even and j 2 for j odd. Remark 6.1. The convergence rate for the j odd case in point 4.1 above is in accordance with the corresponding Green function convergence rate (point 2.1). In contrast, a theoretical basis has not been found as yet for the higher efficiency convergence rates which, as mentioned in point 4.1, are observed in practice for even values of j.

Green-function convergence rates for various differencing orders

Gq j (X, Y ) - N n=-N e -iαnL G j (X + nL, Y ) (90)
in the rapidly convergent Green functions of various orders j for two frequencies, namely, non-Wood-anomaly frequency k = 1.5 and the Wood-anomaly frequency k = 1. For these experiments the values L = 2π, θ = 0 and h = 2.3 were used for the period, incidence angle and shift parameter, respectively. The maximum in equation ( 90) was approximated as the corresponding maximum over all points (X, Y ) = (0, 0) in a 201 × 201 Cartesian grid in the square (-π, π) × (-π, π).

The left portions of Figures 4 and5 present results that arise as the direct differencing formula [START_REF] Mathis | A Comparison Of Acceleration Procedures For The Two-Dimensional Periodic Green's Function[END_REF] is used for all values of N ; a convergence stall is clearly noticeable at large values of N , which, as mentioned in Section 5.4, results from cancellation errors in the differencing process. The corresponding right-hand figures, in turn, display results produced by the asymptotic-approximation algorithm described in Section 5.4. The asymptotic approach clearly eliminates the large-N cancellation errors. The slopes of the curves displayed in Figures 4 and5 are in agreement with the convergence rates cited in point 2. of Section 6.1. : Absolute approximation error in the rapidly convergent Green function G q j as a function of the truncation parameter N for various differencing orders j at the frequency k = 1.5-which is away from the set of {1, 2, 3, . . . } of Wood anomaly frequencies for this setup. Left: Direct evaluation of the expression [START_REF] Mathis | A Comparison Of Acceleration Procedures For The Two-Dimensional Periodic Green's Function[END_REF] for all values of n. Right: Asymptotic method described in Section 5.4. Note that the curves for j = 1 and j = 2, both of which have a slope of 3/2, are virtually indistinguishable from each other.

Integral equations: null-spaces and effect of Green-function completion terms

To demonstrate the ability of the completion terms introduced in equation (64) to eliminate the null-space problems described in Section 4.4, in Sections 6.3.1 and 6.3.2 below we present solutions of problems of scattering produced by means of the basic and complete rapidly-convergent Green functions Gq j and G q j (equations ( 23) and (64), respectively). In those experiments scattering problems are considered for the surface f (x) = π 10 cos(x) (L = 2π, H = 1/10) under an incident angle θ = π/9.

Effect of completion terms at and around resonant shift parameter values h = h c

Using the wavenumber k = 2 we consider solvers that from use of the basic and complete rapidly-convergent Green functions Gq j and G q j with j = 5 and with values of h at and around h = h c , where shift value (see equation ( 55)). As shown in the left portion of Figure 6, the discrete problem associated with the integral equation ( 52), which is based on use of the basic (incomplete) Green function Gq j , is ill-conditioned for shift values h close to h c . We have checked that similar ill-conditioning occurs for the integral equation (53) as well as for both of these equations at other values of h in the resonance set H = h > 0 : h = 2πq βn , q ∈ N, n ∈ U . The right portion of Figure 6 demonstrates that, as expected, large errors (as evidenced by corresponding departures from zero of the defect in the energy-balance condition ( 12)) also occur in such cases. These graphs also display the corresponding condition numbers and energy-balance errors resulting from use of the complete Green function G q j ; clearly good conditioning and high-accuracies are restored by the Green-function completion process.

h c = 2π β 0 ≈ 3.3432 is a resonant ! "# $# %# &# "'# ($# &## "&## "# "# "# ! "# # ) *++,+ - - ---.-/-" ---.-/-$ ---.-/-( ---.-/-% ---.-/-! ---.-/-' ---.-/-0 ---.-/-& ! "# $# %# &# "'# ($# &## "&## "# "# "# ! "# # ) *++,+ - - ---.-/-" ---.-/-$ ---.-/-( ---.-/-% ---.-/-! ---.-/-' ---.-/-0 ---.-/-&

Effect of completion terms at and around Wood anomaly configurations

A situation similar to the one described in the previous section takes place at Wood anomalies: ill conditioning and large errors arise at and around Wood anomaly configurations if the basic Green function Gq j is used; as indicated in Remark 4.7, however, at Wood anomalies such difficulties only arise for j ≥ 2. To demonstrate the restorative character of the Green-function completion terms at Wood anomalies we consider once again the scattering surface and incidence angle mentioned at the beginning of Section 6.3. The set of Wood anomaly frequencies for this problem is given by W = {k ∈ R 2 , k 2 = (k sin(θ) + n) 2 , n ∈ Z}; for the present example we consider wavenumbers in a neighborhood of the n = -1 Wood anomaly value k = k -1 = 1/(1 + sin(θ)) ≈ 0.74. To eliminate resonance effects associated with the shift parameter we select the non-resonant shift value h = 2 (see equation ( 55)) and we consider basic and complete Green functions of order j = 7.

Figure 7 displays the k-dependence of the condition number (left) and energy balance defect (right) which result from discretization of the incomplete and complete integral equations ( 52) and (65). The condition number resulting from use of the incomplete equation becomes very large as k tends to the Wood anomaly; for the complete equation, in contrast, the condition number remains small and extremely small errors result. A small error increase around the Wood anomaly is attributed to the algebraic (finite) Green-function convergence rate at Wood anomalies (Theorem 4.4) as opposed to the super-algebraic convergence (Theorem 3.1) afforded by the smooth-windowing methodology away from Wood anomalies. frequency k -1 ≈ 0.74. Solid (resp. dashed) curves: quantities resulting from use of the incomplete (resp. complete) Green function Gq j (resp. G q j ). A small but noticeable error increase around k = k -1 is attributed to the algebraic (finite) Greenfunction convergence rate at Wood anomalies as opposed to the super-algebraic convergence afforded by the smooth-windowing methodology away from Wood anomalies.

Convergence of forward-maps, solutions and efficiencies

Focusing on integral equations associated with the Dirichlet problem (4)-( 6), this section demonstrates the convergence of the numerical approximation of the integral operator (68) for a given density function μ (the discrete "forward-map") as well as the convergence of the solution itself, as the truncation parameter A in equation ( 69) is increased, or, equivalently, as the truncation parameter N in equation ( 24) grows. The corresponding convergence properties for the Neumann problem are entirely analogous. For these examples the profile Γ with f = f 1 (equation (89)) and with H = 1/2 and L = 2π was used, together with incident waves with incidence angle θ = π/18. In all of the examples considered in this section a sufficiently large number n I of integration points per period was used so that the approximation errors are dominated by the truncation error that arises from use of finite values of the parameter N .

Figure 8 displays errors in the forward map, in the energy balance and in the efficiencies e -1 and e 0 for the n = 1 Wood anomaly frequency k = k 1 = (1 -sin θ) -1 . As expected, the semi-log error curves are straight lines of various slopes, indicating convergence of various algebraic orders; the convergence orders implicit in these curves are in agreement with those laid down in Section 6.1. Figure 9, in turn, demonstrates super-algebraic convergence for the non-Wood-anomaly wavenumber k = 1.4-in agreement, once again, with the theoretical predictions summarized in Section 6.1.

Computing cost

The unaccelerated Green-function representations and solvers presented in this paper give rise to fast numerical algorithms. As demonstrated in this section, the required computing times to produce full solutions of challenging rough-surface scattering problems with very high accuracies can be measured in terms of tens to hundreds of milliseconds, both away, at and around Wood anomalies-highly competitive, in fact, with the performance of other available solvers even for configurations away from Wood anomalies (cf. [START_REF] Barnett | A new integral representation for quasi-periodic scattering problems in two dimensions[END_REF][START_REF] Bruno | Efficient high-order evaluation of scattering by periodic surfaces: deep gratings, high frequencies, and glancing incidences[END_REF][START_REF] Linton | One-and two-dimensional lattice sums for the three-dimensional helmholtz equation[END_REF] and references therein).

In Sections 6.5.1, 6.5.2 and 6.5.3 we present a detailed account of the computing times required by the present solver for a case study concerning the profile Γ with f = f 1 (equation ( 89)), L = 2π and H varying form 0.25 to 1 for configurations far from Wood anomalies, at Wood anomalies and near Wood anomalies. For the sake of diversity, in this section incidence angles are taken in Littrow mount of order ℓ = -1 (for which the scattered wave of order -1 propagates in the backscattering direction, and under which k = 1.5 is Wood-anomaly wavenumber but k = 1 is not; see Section 6.6 and Remark 6.3). Throughout these sections the computing times cited sufficed to produce full scattering solutions with an energy balance error of the order of 10 -8 . Remark 6.2. Significant additional computational acceleration can be induced in the integral equation considered in this section by means of an FFT-based equivalent-source methodology [START_REF] Bruno | A fast, high-order algorithm for the solution of surface scattering problems: basic implementation, tests, and applications[END_REF]-a three-dimensional version of which was outlined and demonstrated in the preliminary text [START_REF] Bruno | Efficient evaluation of doubly periodic Green functions in 3d scattering, including Wood anomaly frequencies[END_REF]. The use of acceleration methods lies beyond the scope of this paper, however; accelerated implementations for two-and three-dimensional configurations will be presented elsewhere.

Computing cost I: wave-numbers away from Wood anomalies

Tables 2, 3 and 4 display the computational times required by our solver to produce a solution with an energy balance error of the order 10 -8 for the wavenumber k = 1.0 (away from the set of Wood anomaly frequencies). We see that in this case the low order solvers are most efficient. This is easily explained: away from Wood anomalies the superalgebraic convergence of the windowed Green-function methodology we use is highly effective. The incremental gains in convergence speed resulting from use higher-order methods does give rise to reductions in the numbers N of periodic intervals necessary to achieve convergence within a given tolerance, but these gains are out-weighted by the additional computing cost required for evaluation of the relatively large numbers of shifted free-space Greenfunction values required in the higher-order Green function.

Clearly, short computing times suffice for evaluation of solution of rather challenging problems. According to Table 4, for example, an 80-millisecond single-core run (j = 1) suffices to produce a complete scattering solution with a 10 -8 energy-balance error for a problem of scattering by a sinusoidal diffraction grating with height-to-period ratio equal to one (cf. [START_REF] Bruno | Efficient high-order evaluation of scattering by periodic surfaces: deep gratings, high frequencies, and glancing incidences[END_REF] and references therein). 2 ) and corresponding single-core computing times (in seconds) necessary to produce the complete scattering solutions with an energy-balance error of the order of 10 -8 on the basis of the complete j-th order Green function G q j for the profile described at the beginning of Section 6.5 with H = 0.25 and for the wavenumber k = 1.0 (far from the set of Wood anomaly frequencies). 2 but for H = 1 and for the frequency k = 1.49 (close to the Wood anomaly frequency k = 1.5).

The ℓ-th order Littrow configuration takes place for sin(θ ℓ ) = -sin(θ), or, in view of (91), for sin(θ) = -πℓ/(kL). Remark 6.3. Under ℓ-th order Littrow mount (k sin(θ) = -πℓ/L), the Wood anomaly condition is given by k 2 = ((2n -ℓ)π/L) 2 . In particular, for the ℓ = -1 Littrow configurations with periods L = 2π considered in this section, Wood anomalies occur for k = ±(2n + 1)/2, where n is an arbitrary integer number.

Figure 10 displays the efficiency e -1 as a function of the frequency k for various gratings of the form y = f 1 (x) (with L = 2π) and y = f 2 (x), where f 1 and f 2 are given by equation (89) for H ∈ {0.1, 0.2, 0.5, 1}. Wood anomaly frequencies are displayed as vertical dotted lines on the efficiency graphs in Figure 10 (cf. Remark 6.3). Corresponding results for the Neumann problem are displayed in Figure 12. Clearly increasingly more oscillatory efficiency curves result as the grating depth is increased. Figure 11 displays corresponding errors: energy-balance errors and "maximum" errors (that is, the maximum departure, over all efficiencies, of each efficiency from its true value, as evaluated by means of a convergence analysis) for the experiment considered in Figure 10(a). Corresponding results for the gratings y = f 1 (x) and y = f 2 (x) under Neumann boundary conditions are presented in Figure 12; errors similar to those reported in Figure 11 arise for all of the experiments considered in this section.

Conclusions

This paper introduces a new methodology for solution of problems of scattering by periodic surfaces. To our knowledge, these are the first Green-function methods ever presented for the solution of problems of scattering by periodic diffraction gratings that are applicable even at and around Wood anomaly configurations. The computing costs required by the proposed methods are very small, and, as pointed out in Remark 6.2, significant additional acceleration can be induced in the proposed methodology by means of an FFT-based equivalent-source acceleration strategy.

Figure 1 :

 1 Figure 1: Windowing function S = S(x, x 0 , x 1 ) for x 0 = 1 and x 1 = 2.

  ) ∼ C e iz √ z , for large positive values of x ′ we obtain

Theorem 4 . 4 .

 44 Let j ∈ N, h > 0, ℓ = 0 or ℓ = 1, m = 1 or m = 2 and k > 0 be given. Then, for each M > 0 there exists a constant D M > 0 (that also depends on j, k and h) such that, for all X, Y satisfying -L ≤ X ≤ L and -M < Y < M , ad for all integers N > 1, we have n∈Z,|n|>N

Figure 2 :

 2 Figure 2: Left. Sample scattering surface y = f (x) = π cos(x) (solid line) and shifted surface y = f (x) -h with h = 0.8.

Figure 3 :

 3 Figure 3: Grating profiles y = f 1 (x) = HL 2 cos 2πx L for L = 2π and H = 1 (left) and y = f 2 (x) = π H (0.4 cos x -0.2 cos(2x) + 0.4 cos(3x)) with H = 1 (right).

1 (and f 2

 12 profiles Γ (equation (1)), namely, the sinusoidal profile family f = f 1 and the multi-frequency profile family f = f 2 where f (x) = π H (0.4 cos x -0.2 cos(2x) + 0.4 cos(3x))

4. 2 - 2 . 5 .

 225 away from Wood-anomaly configurations: super-algebraically fast, per the rationale in point 3.Convergence to zero of the energy balance error-5.1 -at Wood-anomaly configurations: The same as in point 4.1 above. 5.2 -away from Wood-anomaly configurations: super-algebraically fast, per the rationale in point 3.2.

Figures 4 and 5

 5 Figures 4 and 5 display the maximum absolute error

Figure 4

 4 Figure 4: Absolute approximation error in the rapidly convergent Green function G q

Figure 5 :

 5 Figure 5: Same as Figure 4 but for the Wood anomaly frequency k = 1. Once again the j = 1 and j = 2 curves, which in this case have a slope of 1/2, are virtually indistinguishable from each other.

Figure 6 :

 6 Figure 6: Condition number (left) and energy balance error (right) as functions of the shift parameter h around the resonant value h c ≈ 3.3432. Solid (resp. dashed) curves: quantities resulting from use of the incomplete (resp. complete) Green function Gq j (resp. G q j ).

Figure 7 :

 7 Figure 7: Condition number (left) and energy balance error (right) as functions of the frequency k around the Wood anomaly

Table 2 :

 2 Numbers n I (number of discretization points per period) and N (number of periodic intervals contained in the integration domain to the right of the point x = L

  j

Figure 10 :

 10 Figure 10: Efficiency e -1 as a function of the spatial frequency k in the Littrow configuration of order -1 under Dirichlet boundary conditions, for the profiles y = f 1 (x) with L = 2π and y = f 2 (x), for various values of H.

Figure 11 :

 11 Figure 11: Energy-balance error and maximum efficiency error for the Littrow configuration of order -1 considered in Figure 10(a).

Figure 12 :

 12 Figure 12: Same as Figure 10 but with configurations under Neumann boundary conditions.

Table 1 :

 1 Approximation errors for various values of A (k n = 2πc = 0.1)

	A 10 20 25 50 75 100	|I ex -I H,A | 5.0 × 10 -2 3.6 × 10 -2 3.2 × 10 -2 2.3 × 10 -2 1.8 × 10 -2 1.6 × 10 -2	|I ex -I S,A | 8.5 × 10 -5 9.7 × 10 -7 1.9 × 10 -7 4.9 × 10 -10 4.7 × 10 -11 7.7 × 10 -14

Table 6 :

 6 Time (s) 1.34 0.55 0.35 0.22 0.26 0.21 Same as Table2but for H = 0.5 and for the Wood anomaly frequency k = 1.5.

		3	4	5	6	7	8
	n I	34	36	38	38	38	38
	N	2600 600 300 200 140 100
	j	4	5	6	7	8
	n I	48	48	50	50	54
	N	1200 300 260 140 100
	Time (s) 1.27 0.57 0.55 0.48 0.63

Table 7 :

 7 Same as Table2but for H = 1 and for the Wood anomaly frequency k = 1.5. Time (s) 2.58 0.29 0.16 0.11 0.06 0.05 0.06

	j	0	1	2	3	4	5	6
	n I	26	26	26	26	26	26	26
	N	3000 1000 600 350 250 120 100

Table 8 :

 8 Same as Table 2 (H = 0.25) but for the frequency k = 1.49 (close to the Wood anomaly frequency k = 1.5). Time (s) 3.35 0.57 0.34 0.26 0.21 0.15 0.24

	j	0	1	2	3	4	5	6
	n I	34	34	34	34	34	34	38
	N	2500 1200 800 450 300 220 140

Table 9 :

 9 Same as Table2but for H = 0.5 and for the frequency k = 1.49 (close to the Wood anomaly frequency k = 1.5). Time (s) 6.61 1.89 1.30 1.92 1.37 1.22 0.90

	j	0	1	2	3	4	5	6
	n I	48	48	50	50	54	50	50
	N	2500 1800 1400 1200 800 250 250

Table 10 :

 10 Same as Table

Acknowledgments

The authors thank Dr. Santiago Fortes for his assistance during the preliminary phases of this project. Support from NSF and AFOSR is gratefully acknowledged.

Computing cost II: a Wood anomaly frequency

Tables 5, 6 and 7 present statistics for the diffraction gratings considered in the previous section, but for the Woodanomaly wavenumber k = 1.5. As shown in Table 7, in this case a computational time of less than five-hundred milliseconds suffices to solve the Wood-anomaly scattering problem for the deep grating H = 1 (height-to-period ratio equal to one) with an energy-balance error of the order of 10 2 (H = 0.25), but for the Wood anomaly frequency k = 1.5.

Computing cost III: frequencies near Wood anomalies

Tables 8, 9 and 10, finally, present statistics for the diffraction gratings considered in the previous sections, but for the wavenumber k = 1.49, which is near but different from the Wood anomaly value k = 1.5. Interestingly, the computing times required for this near-anomalous case are actually somewhat larger (albeit by a relatively small factor) than those presented in the previous section for the actual Wood anomaly. In the present case a computing time of nine-hundred milliseconds is necessary to solve the near Wood-anomaly scattering problem for the deep grating H = 1 (height-to-period ratio equal to one) with an energy-balance error of the order of 10 -8 , as opposed to the five-hundred milliseconds required for the Wood anomaly case k = 1.5.

First order Littrow configuration

To conclude our numerical results section we present an application of our solvers to problems involving diffraction gratings under the classical Littrow mount [START_REF] Petit | Electromagnetic Theory Of Gratings[END_REF]. Under Littrow mount of order ℓ the ℓ-th diffracted wave, whose efficiency is e ℓ , propagates precisely in the backscattering direction. Calling θ ℓ the angle of propagation of the ℓ-th wave, the ℓ-th order Littrow configuration occurs when sin(θ ℓ ) = -sin(θ). Since the propagation angle of the ℓ-th order scattered mode is given by α ℓ = α + 2πℓ L , in terms of the scattering angle θ ℓ we have sin(θ ℓ ) = sin(θ) + 2πℓ kL , ℓ ∈ Z .

(91)