Three-dimensional quasi-periodic shifted Green function throughout the spectrum--including Wood anomalies

TL;DR

This paper introduces a novel Green function method for efficiently solving wave scattering problems involving doubly periodic gratings at and near Wood anomaly frequencies, overcoming previous convergence issues.

Contribution

It presents a modified quasi-periodic Green function that ensures convergence and solvability at Wood anomalies, along with an efficient high-order numerical solver for these challenging problems.

Findings

Successfully solves scattering problems at Wood frequencies.

Achieves algebraic convergence with shifted Green functions.

Demonstrates applicability to acoustic scattering simulations.

Abstract

This work presents an efficient method for evaluation of wave scattering by doubly periodic diffraction gratings at or near "Wood anomaly frequencies". At these frequencies, one or more grazing Rayleigh waves exist, and the lattice sum for the quasi-periodic Green function ceases to exist. We present a modification of this sum by adding two types of terms to it. The first type adds linear combinations of "shifted" Green functions, ensuring that the spatial singularities introduced by these terms are located below the grating and therefore outside of the physical domain. With suitable coefficient choices these terms annihilate the growing contributions in the original lattice sum and yield algebraic convergence. Convergence of arbitrarily high order can be obtained by including sufficiently many shifts. The second type of added terms are quasi-periodic plane wave solutions of the Helmholtz equation which reinstate certain necessary grazing modes without leading to blow-up at Wood anomalies. Using the new quasi-periodic Green function, we establish, for the first time, that the Dirichlet problem of scattering by a smooth doubly periodic scattering surface at a Wood frequency is uniquely solvable. We also present an efficient high-order numerical method based on the this new Green function for the problem of scattering by doubly periodic three-dimensional surfaces at and around Wood frequencies. We believe this is the first solver in existence that is applicable to Wood-frequency doubly periodic scattering problems. We demonstrate the proposed approach by means of applications to problems of acoustic scattering by doubly periodic gratings at various frequencies, including frequencies away from, at, and near Wood anomalies.

Tables

Table 1. Table 1: Convergence of the p = 0 𝑝 0 p=0 Dirichlet solver (unshifted) as A 𝐴 A grows, for a configuration away from Wood anomalies. Normal incidence ( 𝜶 = 0 𝜶 0 \boldsymbol{\alpha}=0 ) was assumed for this example. The reference solution was produced using A = A r ​ e ​ f = 320 𝐴 subscript 𝐴 𝑟 𝑒 𝑓 320 A=A_{ref}=320 and 16 × 16 16 16 16\ \times\ 16 unknowns, for which ε = 1.0 × 10 − 6 𝜀 1.0 superscript 10 6 \varepsilon=1.0\times 10^{-6} .

$k$	Unknowns	$A$	$\max \left|G^A-G^{ref}\right|$	Iter	${\epsilon }_1$	$\epsilon$
$1$	$16\times 16$	30	4.4 $\times$ ${10}^{-2}$	13	1.0 $\times$ ${10}^{-1}$	1.8 $\times$ ${10}^{-1}$
$1$	$16\times 16$	60	3.1 $\times$ ${10}^{-3}$	13	6.3 $\times$ ${10}^{-3}$	6.6 $\times$ ${10}^{-3}$
$1$	$16\times 16$	120	4.3 $\times$ ${10}^{-4}$	13	4.8 $\times$ ${10}^{-4}$	4.3 $\times$ ${10}^{-4}$
$1$	$16\times 16$	160	3.7 $\times$ ${10}^{-5}$	13	1.3 $\times$ ${10}^{-4}$	2.3 $\times$ ${10}^{-4}$
$1$	$16\times 16$	240	2.4 $\times$ ${10}^{-6}$	13	4.9 $\times$ ${10}^{-6}$	6.7 $\times$ ${10}^{-6}$

Table 2. Table 2: Convergence of the p = 3 𝑝 3 p=3 Dirichlet solver with shift parameter d = 2.4 𝑑 2.4 d=2.4 , away from Wood frequencies, as A 𝐴 A grows. Normal incidence. The reference solution was produced using corresponds to A = A r ​ e ​ f = 120 𝐴 subscript 𝐴 𝑟 𝑒 𝑓 120 A=A_{ref}=120 , 32 × 32 32 32 32\ \times\ 32 unknowns, for which ε = 4.0 × 10 − 7 𝜀 4.0 superscript 10 7 \varepsilon=4.0\times 10^{-7} .

$k$	Unknowns	$A$	$G^A$			$G^{A,p},p=3$
			Iter	${\epsilon }_1$	$\epsilon$	Iter	${\epsilon }_1$	$\epsilon$
$6$	$16\times 16$	30	16	4.9 $\times$ ${10}^{-3}$	1.6 $\times$ ${10}^{-3}$	12	6.5 $\times$ ${10}^{-3}$	1.2 $\times$ ${10}^{-2}$
$6$	$16\times 16$	60	16	1.5 $\times$ ${10}^{-3}$	4.1 $\times$ ${10}^{-4}$	12	3.8 $\times$ ${10}^{-4}$	1.5 $\times$ ${10}^{-5}$
$6$	$16\times 16$	80	16	2.8 $\times$ ${10}^{-5}$	1.1 $\times$ ${10}^{-5}$	12	4.7 $\times$ ${10}^{-6}$	2.3 $\times$ ${10}^{-6}$

Table 3. Table 3: Convergence, as A 𝐴 A grows, of the p = 3 𝑝 3 p=3 Dirichlet solver at and around Wood frequencies. Shift parameter d = 1.4 𝑑 1.4 d=1.4 , GMRES residual tolerance equal to 10 − 6 superscript 10 6 10^{-6} . “Ref” refers to finely resolved solutions against which the error of the coarser solutions is evaluated.

$k$	Unknowns	$A$	Iter	${\epsilon }_1$	$\epsilon$
$2\pi$	$24\times 24$	20	19	3.2 $\times$ ${10}^{-2}$	1.7 $\times$ ${10}^{-2}$
$2\pi$	$24\times 24$	30	19	2.7 $\times$ ${10}^{-3}$	4.7 $\times$ ${10}^{-3}$
$2\pi$	$24\times 24$	40	19	6.9 $\times$ ${10}^{-4}$	4.0 $\times$ ${10}^{-4}$
$2\pi$	$24\times 24$	60	19	ref	2.4 $\times$ ${10}^{-6}$
$2\pi \pm {10}^{-6}$	$24\times 24$	40	19	6.8 $\times$ ${10}^{-4}$	4.3 $\times$ ${10}^{-4}$
$2\sqrt{2}\pi$	$24\times 24$	30	25	1.3 $\times$ ${10}^{-2}$	1.5 $\times$ ${10}^{-2}$
$2\sqrt{2}\pi$	$24\times 24$	40	25	6.8 $\times$ ${10}^{-3}$	5.5 $\times$ ${10}^{-3}$
$2\sqrt{2}\pi$	$24\times 24$	80	25	8.9 $\times$ ${10}^{-5}$	2.1 $\times$ ${10}^{-4}$
$2\sqrt{2}\pi$	$24\times 24$	120	25	ref	3.8 $\times$ ${10}^{-5}$
$2\sqrt{2}\pi \pm {10}^{-6}$	$24\times 24$	30	25	1.6 $\times$ ${10}^{-2}$	2.5 $\times$ ${10}^{-2}$
$2\sqrt{2}\pi \pm {10}^{-6}$	$24\times 24$	80	25	8.9 $\times$ ${10}^{-5}$	1.5 $\times$ ${10}^{-4}$
$4\pi$	$32\times 32$	30	28	4.6 $\times$ ${10}^{-2}$	4.9 $\times$ ${10}^{-2}$
$4\pi$	$32\times 32$	60	28	2.4 $\times$ ${10}^{-3}$	1.1 $\times$ ${10}^{-3}$
$4\pi$	$32\times 32$	180	28	ref	2.2 $\times$ ${10}^{-4}$

Table 4. Table 4: Convergence, as A 𝐴 A grows, of the p = 3 𝑝 3 p=3 Neumann solver at and around Wood frequencies. Shift parameter d = 1.4 𝑑 1.4 d=1.4 , GMRES residual tolerance equal to 10 − 6 superscript 10 6 10^{-6} . “Ref” refers to finely resolved solutions against which the error of the coarser solutions is evaluated.

$k$	Unknowns	$A$	Iter	${\epsilon }_1$	$\epsilon$
$2\pi$	$24\times 24$	20	19	5.3 $\times$ ${10}^{-3}$	7.4 $\times$ ${10}^{-3}$
$2\pi$	$24\times 24$	30	19	5.7 $\times$ ${10}^{-4}$	1.7 $\times$ ${10}^{-3}$
$2\pi$	$24\times 24$	40	19	1.1 $\times$ ${10}^{-4}$	3.7 $\times$ ${10}^{-4}$
$2\pi$	$24\times 24$	80	19	ref	4.5 $\times$ ${10}^{-6}$
$2\pi \pm {10}^{-6}$	$24\times 24$	40	19	1.1 $\times$ ${10}^{-4}$	3.6 $\times$ ${10}^{-4}$
$2\sqrt{2}\pi$	$24\times 24$	30	25	7.3 $\times$ ${10}^{-2}$	5.1 $\times$ ${10}^{-2}$
$2\sqrt{2}\pi$	$24\times 24$	40	25	4.1 $\times$ ${10}^{-3}$	2.8 $\times$ ${10}^{-3}$
$2\sqrt{2}\pi$	$24\times 24$	80	25	2.6 $\times$ ${10}^{-4}$	3.4 $\times$ ${10}^{-4}$
$2\sqrt{2}\pi$	$24\times 24$	160	25	ref	4.2 $\times$ ${10}^{-5}$
$2\sqrt{2}\pi \pm {10}^{-6}$	$24\times 24$	80	25	2.5 $\times$ ${10}^{-4}$	3.5 $\times$ ${10}^{-4}$
$4\pi$	$32\times 32$	30	28	1.2 $\times$ ${10}^{-1}$	4.5 $\times$ ${10}^{-2}$
$4\pi$	$32\times 32$	60	28	2.7 $\times$ ${10}^{-3}$	1.6 $\times$ ${10}^{-3}$
$4\pi$	$32\times 32$	180	28	ref	1.1 $\times$ ${10}^{-4}$

Equations

$$G(\mathbf{x})=\frac{e^{ik|\mathbf{x}|}}{4\pi |\mathbf{x}|}$$

$$r_{mn}^2={\left|\mathbf{x}+m{\mathbf{v}}_1+n{\mathbf{v}}_2\right|}^2={\left|\tilde{\mathbf{x}}+m{\mathbf{v}}_1+n{\mathbf{v}}_2\right|}^2+z^2,$$

$$G^q(\mathbf{x})=\sum _{m,n\in \mathbb{Z}}G(\mathbf{x}+m{\mathbf{v}}_1+n{\mathbf{v}}_2)e^{-i\mathbf{\alpha }\cdot (m{\mathbf{v}}_1+n{\mathbf{v}}_2)}=\frac{1}{4\pi }\sum _{m,n\in \mathbb{Z}}\frac{e^{ik{r}_{mn}}}{r_{mn}}{e}^{-i\mathbf{\alpha }\cdot {\mathbf{v}}_{mn}},$$

$$G^q(\tilde{\mathbf{x}},z)=\frac{i}{2D}\sum _{j,\ell \in \mathbb{Z}}\frac{1}{{\gamma }_{j\ell }}{e}^{i{\mathbf{v}}_{j\ell }^{\ast }\cdot \tilde{\mathbf{x}}}{e}^{i{\gamma }_{j\ell }|z|},$$

$$\Lambda =\{{\mathbf{v}}_{mn}:m,n\in \mathbb{Z}\}\text{mbox}and{\Lambda }^{\ast }=\{{\mathbf{v}}_{j\ell }^{\ast }:j,\ell \in \mathbb{Z}\},$$

$$\Gamma =\{(\tilde{\mathbf{x}},z):\tilde{\mathbf{x}}\in {\mathbb{R}}^2\text{mbox}andz=f(\tilde{\mathbf{x}})\}$$

$$f(\tilde{\mathbf{x}}+m{\mathbf{v}}_1+n{\mathbf{v}}_2)=f(\tilde{\mathbf{x}})\text{for all}\tilde{\mathbf{x}}\in {\mathbb{R}}^2\text{and all}m,n\in \mathbb{Z}.$$

$${\Omega }^{+}=\{(\tilde{\mathbf{x}},z):\tilde{\mathbf{x}}\in {\mathbb{R}}^2\text{mbox}andz>f(\tilde{\mathbf{x}})\},$$

$$u^{inc}(\mathbf{x})=\exp [i(\mathbf{\alpha }\cdot \tilde{\mathbf{x}}-\gamma z)]$$

$$\{\begin{array}{l}\Delta u+k^{2}u=0\mathrm{in}{\Omega }^{+}\\ u=-u^{inc}\mathrm{on}\Gamma ,\end{array}$$

$$u(\tilde{\mathbf{x}},z)=\sum _{j,\ell \in \mathbb{Z}}{B}_{j\ell }^{+}\exp i[(2\pi j{\mathbf{v}}_1^{\ast }+2\pi \ell {\mathbf{v}}_2^{\ast })+\mathbf{\alpha }]\cdot \tilde{\mathbf{x}}\exp [i{\gamma }_{j\ell }z],z>z_{+}.$$

$${\Omega }^{-}=\{(\tilde{\mathbf{x}},z):\tilde{\mathbf{x}}\in {\mathbb{R}}^2\text{mbox}andz<f(\tilde{\mathbf{x}})\}$$

$$u(\tilde{\mathbf{x}},z)=\sum _{j,\ell \in \mathbb{Z}}{B}_{j\ell }^{-}\exp \{i[(2\pi j{\mathbf{v}}_1^{\ast }+2\pi \ell {\mathbf{v}}_2^{\ast })+\mathbf{\alpha }]\cdot \tilde{\mathbf{x}}\}\exp [-i{\gamma }_{j\ell }z],z<z_{-}.$$

$$\tilde{G}_{p}({\mathbf{x}})=\sum_{q=0}^{p}a_{pq}G\!\big{(}{\mathbf{x}}+(0,0,qd)\big{)}=\sum_{q=0}^{p}a_{pq}G\!\big{(}\tilde{\mathbf{x}}+(0,0,z+qd)\big{)},$$

$$a_{pq}=(-1{)}^q(\genfrac{}{}{0ex}{}{p}{q}),0\le q\le p.$$

$${\tilde{G}}_p^q(\mathbf{x})=\sum _{m,n\in \mathbb{Z}}{\tilde{G}}_p(\mathbf{x}+m{\mathbf{v}}_1+n{\mathbf{v}}_2)e^{-i\mathbf{\alpha }\cdot (m{\mathbf{v}}_1+n{\mathbf{v}}_2)}.$$

$${\mathbf{x}}_{mnq}^{\prime }=(m{\mathbf{v}}_1+n{\mathbf{v}}_2,-qd),m,n,q\in \mathbb{Z},0\le q\le p.$$

$$g(\rho ,\epsilon ):=\frac{e^{ik\rho \sqrt{1+{\epsilon }^2}}}{\rho \sqrt{1+{\epsilon }^2}},{\rho }_{mn}=|\tilde{\mathbf{x}}+m{\mathbf{v}}_1+n{\mathbf{v}}_2|,{\epsilon }_{mn}=z/{\rho }_{mn}\text{mbox}and{\hat{\epsilon }}_{mn}=d/{\rho }_{mn},$$

$${\tilde{G}}_p(\mathbf{x}+m{\mathbf{v}}_1+n{\mathbf{v}}_2)=\sum _{q=0}^p{a}_{pq}g({\rho }_{mn},{\epsilon }_{mn}+q{\hat{\epsilon }}_{mn}).$$

$$h(\rho ,\epsilon ,\hat{\epsilon }):=\sum _{q=0}^p{a}_{pq}g(\rho ,\epsilon +q\hat{\epsilon }),$$

$${\tilde{G}}_p(\mathbf{x}+m{\mathbf{v}}_1+n{\mathbf{v}}_2)=h({\rho }_{mn},{\epsilon }_{mn},{\hat{\epsilon }}_{mn}).$$

$$\sum _{q=0}^p{a}_{pq}f(\epsilon +q\hat{\epsilon }))=(-1{)}^p{\hat{\epsilon }}^p{f}^{(p)}(\epsilon +\xi )\xi \in [0,p\hat{\epsilon }],$$

$$h(\rho ,\epsilon ,\hat{\epsilon })=(-1{)}^p{\hat{\epsilon }}^p\frac{{\partial }^{p}g}{\partial {\epsilon }^p}(\rho ,\epsilon +{\xi }_{\rho ,\hat{\epsilon }}),$$

$$h\left(\rho ,z/\rho ,d/\rho \right)=(-1{)}^p{\left(\frac{d}{\rho }\right)}^p\frac{{\partial }^{p}g}{\partial {\epsilon }^p}(\rho ,z/\rho +{\xi }_{\rho ,d/\rho }),{\xi }_{\rho ,d/\rho }\in [0,pd/\rho ].$$

$$\left|\frac{1}{{\rho }^p}\frac{{\partial }^{p}g}{\partial {\epsilon }^p}(\rho ,\epsilon )\right|\le \frac{C}{{\rho }^{\lceil \frac{p}{2}\rceil +1}}(\rho \ge 1,|\epsilon |<1),$$

$$h\big{(}\rho,\,z/\rho,\,d/\rho\big{)}\,=\,{\mathcal{O}}\left(\frac{1}{\rho^{\lceil\frac{p}{2}\rceil+1}}\right)\qquad(\rho\to\infty).$$

$$\frac{{\partial }^{p}g}{\partial {\epsilon }^p}(\rho ,\epsilon )=\sum _{(m,n,\ell )\in S_p}{C}_{m,n,\ell }^p\frac{{\epsilon }^m(ik\rho {)}^n}{{\sqrt{1+{\epsilon }^2}}^{\ell }}\frac{e^{ik\rho \sqrt{1+{\epsilon }^2}}}{\rho },$$

$$T_{m,n,\ell }(\rho ,\epsilon )=\frac{{\epsilon }^m(ik\rho {)}^n}{{\sqrt{1+{\epsilon }^2}}^{\ell }}\text{for }m,n,\ell \ge 0,$$

$$\frac{{\partial }^{p}g}{\partial {\epsilon }^p}(\rho ,\epsilon )=\sum _{(m,n,\ell )\in S_p}{C}_{m,n,\ell }^p{T}_{m,n,\ell }(\rho ,\epsilon )\frac{e^{ik\rho \sqrt{1+{\epsilon }^2}}}{\rho }.$$

$$g(\rho ,\epsilon )=T_{0,0,1}(\rho ,\epsilon )\frac{e^{ik\rho \sqrt{1+{\epsilon }^2}}}{\rho },$$

Full text

$\mathbf{Keywords}$: scattering, periodic Green function, lattice sum, smooth truncation, Wood frequency, Wood anomaly, boundary-integral equations, electromagnetic computation

1 Introduction

This work presents the second part of a two-part contribution. The first part [8], which will be referenced as Part I throughout this paper, introduced a “windowed Green function” method, which, utilizing a smooth cutoff, approximates the quasi-periodic Green function with super-algebraically small errors—that is, errors that admit upper bounds proportional to any negative power of the numbers of terms used—for configurations that are not close to a certain set of “Wood frequencies”. As discussed in Part I and references therein, at Wood frequencies the classical quasi-periodic Green function ceases to exist, and, therefore, integral equation methods based on such Green functions are inapplicable. Following upon the two-dimensional work [4], the present Part II introduces an additional element in the method: the “shifted” quasi-periodic Green function. By adding copies of the quasi-periodic Green function that are shifted perpendicular to the plane of periodicity, one obtains an algebraic convergence rate that increases with the number of shifts, even at and around Wood frequencies. Using the new quasi-periodic Green function ${G}^{q}_{p}({\mathbf{x}})$, we establish, for the first time, that the Dirichlet problem of scattering by a smooth doubly periodic surface, or diffraction grating, at a Wood frequency is uniquely solvable. We present an efficient high-order numerical method based on the Green function ${G}^{q}_{p}({\mathbf{x}})$ for the problem of scattering by doubly periodic three-dimensional gratings at and around Wood frequencies. We believe this is the first solver in existence that is applicable to doubly periodic scattering problems at Wood frequencies.

Wood frequencies depend on the wave vector parallel to the grating, and they occur when one of the Rayleigh waves, or diffraction orders, is at the transition from propagating to evanescent in the ambient medium and is thus exactly grazing the surface. Because of this, they are often called “cutoff frequencies”. Certain scattering anomalies at or near these frequencies, known as “Wood anomalies”, were first observed experimentally by Wood [25] and treated mathematically by Rayleigh [23] and Fano [14]. A brief discussion concerning historical aspects in these regards can be found in [4, Remark 2.2].

Computation of scattering by three-dimensional doubly periodic structures at or near Wood frequencies is particularly challenging. As mentioned above, in boundary integral methods the classical quasi-periodic Green function ceases to exist at Wood frequencies, even when the scattering problem admits a unique solution. A boundary-integral method that employs the free-space Green function and enforces periodicity through auxiliary layer potentials on the boundary of a period has been developed for two-dimensional problems [1, 15], but a three-dimensional version of this method does not as yet exist. Approaches based on finite-element methods [12] often rely on the classical quasi-periodic Green functions in order to enforce the radiation condition at infinity; such approaches must also necessarily fail at Wood anomalies. Finally, finite-element methods exist, such as that presented in the contribution [11], which enforce the radiation condition on the basis of sponge layers such as the perfectly-matched-layer technique. The results of that reference indicate that such approaches may become problematic for truly quasi-periodic problems and, we suggest, the difficulties are compounded at Wood frequencies, at which it would be necessary to damp waves which travel in directions parallel to the absorbing layer.

The classical quasi-periodic Green function $G^{q}({\mathbf{x}})$ (${\mathbf{x}}=(x,y,z)\in\mathbb{R}^{3}$) can be constructed as an infinite sum of translated copies of the free-space Helmholtz Green function

[TABLE]

with doubly periodically distributed monopole singularities. Indeed, let $\tilde{\mathbf{x}}=(x,y)$, $\boldsymbol{\alpha}:=(\alpha,\beta)$ (the Bloch wave-vector) and

[TABLE]

in which $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$ denote two independent vectors in $\mathbb{R}^{2}$ that characterize the periodicity. The quasi-periodic Green function can be expressed in the form

[TABLE]

in which $\mathbf{v}_{mn}=m\mathbf{v}_{1}+n\mathbf{v}_{2}$. As is well known, this expansion suffers from notoriously poor convergence properties. Various methods to accelerate its convergence, notably the Ewald method [13, 10, 22], have been proposed. Importantly, however, the sum does not converge at Wood frequencies. This is a difficulty that is not addressed by any of the aforementioned acceleration approaches.

The classical quasi-periodic Green function additionally admits a spectral representation that results from application of the Poisson Summation Formula to the series (3): letting $\mathbf{v}_{1}^{*}$ and $\mathbf{v}_{2}^{*}$ denote the dual vectors defined by $\mathbf{v}_{i}^{*}\cdot\mathbf{v}_{j}=\delta_{ij}$, and letting $D=\|\mathbf{v}_{1}\times\mathbf{v}_{2}\|$, $\mathbf{v}_{j\ell}^{*}=(2\pi j\ \mathbf{v}_{1}^{*}+2\pi\ell\ \mathbf{v}_{2}^{*})+\boldsymbol{\alpha}$ and $\gamma_{j\ell}=(k^{2}-\|\mathbf{v}_{j\ell}^{*}\|^{2})^{\frac{1}{2}}$, we have the alternative expansion

[TABLE]

which converges provided $\gamma_{j\ell}\neq 0$ for all $j,\ell$—that is, away from Wood configurations $(\alpha,\beta,k)$ where one of these exponents vanishes. The branch of the square root that defines $\gamma_{j\ell}$ is selected in such a way that $\sqrt{1}=1$, with a branch cut that coincides with the negative imaginary semiaxis.

The primary and dual periodicity lattices are denoted by

[TABLE]

respectively. The lattice sum (3) is only defined if $\gamma_{j\ell}\neq 0$ for all integer pairs $(j,\ell)$. A Wood frequency (for given $\boldsymbol{\alpha}$) is a value of $k$ for which at least one of the constants $\gamma_{j\ell}$ vanishes. In such cases both the spatial expansion (3) and the spectral representation (4) cease to exist.

This paper presents a Green function method that enables efficient and accurate evaluation of wave scattering by doubly periodic structures throughout the spectrum, including frequencies at and around Wood anomalies. The present contribution additionally incorporates the windowing approach introduced in Part I, which accelerates the Green-function convergence: a windowed version of the series (3) converges superalgebraically fast (i.e., faster than any power of the window size) for non-Wood configurations. The convergence of this windowed series deteriorates near Wood frequencies: the constants in the superalgebraic convergence estimates established in Part I grow without bound as Wood configurations are approached, and the windowed version of the lattice sum (3) once again fails to converge at Wood frequencies. When applied to the shifting method introduced in this work, however, the smooth windowing method increases the algebraic convergence rate at Wood anomalies by a factor equal to the truncation size raised to the power $-1/2$.

In order to re-establish convergence at Wood frequencies the proposed method replaces the free-space Green function term $\frac{e^{ikr_{mn}}}{r_{mn}}$ by a $p$-th order equispaced finite-difference for this function with respect to $z$, with $p\geq 3$ and with step (or “shift”) $d>0$. The combined effect of the windowing and shifting/finite-differencing procedure yields a Green function which converges rapidly at all frequencies. The approach is demonstrated in Section 5 via an application to the problems of sound-soft and sound-hard scattering by doubly periodic surfaces throughout the spectrum. Other scattering problems can be treated similarly—as demonstrated in the contributions [4] and [5] concerning two-dimensional diffraction gratings and periodic arrays of cylinders, respectively.

The remainder of this paper is organized as follows. After some preliminaries in Section 2, the shifted quasi-periodic Green function is introduced in Section 3. This section contains a proof of high-order algebraic convergence of the shifted Green function and an existence and uniqueness proof for configurations at and around Wood anomalies under Dirichlet boundary conditions. (Appendix A contains a lemma used in the aforementioned uniqueness proof, as well as a simplified existence and uniqueness proof which is only valid away from Wood anomalies.) Section 4 then outlines our numerical implementation. A variety of results presented in Section 5 demonstrate the character of the proposed solvers, for all configurations, far, near and at Wood anomalies.

2 Preliminaries

We consider the sound-soft problem of scattering by a doubly periodic scattering surface

[TABLE]

where $f$ is a smooth doubly periodic function of periodicity $\Lambda$:

[TABLE]

Letting

[TABLE]

and assuming an incident field

[TABLE]

impinges upon the surface from above (where $(\boldsymbol{\alpha},-\gamma)$ is the wavevector and $|\boldsymbol{\alpha}|^{2}+\gamma^{2}=k^{2}$), the scattered field $u$ under sound-soft conditions satisfies the equations

[TABLE]

together with the Sommerfeld radiation condition [24]: letting $z_{+}=\max{f}$ we have

[TABLE]

Although not physically relevant for the grating problems considered in this paper, the set

[TABLE]

below $\Gamma$ and the associated radiation condition

[TABLE]

($z_{-}=\min{f}$) will be used in the existence and uniqueness proofs in Section 3.4 and Appendix A.

3 Shifted quasi-periodic Green function

The finite-difference half-space shifted Green function we use can be viewed as a generalization of the Dirichlet half-space Green that results from the method of images. The method-of-images Green function decays more rapidly at infinity than the free-space Green function itself. In the three-dimensional case under consideration, such decay does not suffice to induce fast convergence, or even absolute convergence, in a corresponding series of the form (3). (In contrast, absolute convergence does result from use of the method-of-images Green function in the two-dimensional case [4, Lemma 4.3 and Th. 4.4]). But viewing the method-of-images Green function as a finite difference of the lowest order, a generalization of this idea emerges: as shown in what follows, a higher-order finite difference of the free-space Green function tends to zero more rapidly at infinity, and thus gives rise to an absolutely convergent quasi-periodic Green function series. In fact, the resulting quasi-periodic Green function can be made to converge with a prescribed order of accuracy provided finite-differences of sufficiently high order are used.

To pursue this idea, fix a “shift” value $d>0$ and define the shifted half-space Green function

[TABLE]

where $a_{pq}$ denote the finite-difference coefficients [4, 19]

[TABLE]

Clearly, the function $\tilde{G}_{p}({\mathbf{x}})$ has poles at the points $(0,0,-qd)$ for $0\leq q\leq p$. As shown in Lemma 3.1 below, the shifted half-space Green function $\tilde{G}_{p}({\mathbf{x}})$ tends to zero algebraically fast as $\tilde{\mathbf{x}}$ tends to infinity while $z$ remains bounded. The shifted quasi-periodic Green function is then defined by

[TABLE]

The function $\tilde{G}^{q}_{p}({\mathbf{x}})$ has poles at the points

[TABLE]

Theorem 3.2 below shows that the sum on the right-hand side of equation (14) converges algebraically fast: a truncation of this series to $|m\mathbf{v}_{1}+n\mathbf{v}_{2}|\leq A$ results in errors of order $\lceil p/2\rceil-1$ (where $\lceil r\rceil$ denotes the smallest integer greater than or equal to $r$) for all wavevectors $\mathbf{k}=(\boldsymbol{\alpha},-\gamma)$ and all periodicity lattices $\Lambda$, including Wood-anomaly configurations. That theorem also establishes that somewhat improved accuracies result when a smooth windowed truncation, as introduced in Part I, is additionally employed.

When evaluated at $(\tilde{\mathbf{x}}-\tilde{\mathbf{x}}^{\prime},z-z^{\prime})$ with $(\tilde{\mathbf{x}},z)\in\Gamma$ and $(\tilde{\mathbf{x}}^{\prime},z^{\prime})\in\Gamma$, the terms with $q=1,\dots,p$ in the series (13)-(14) are weighted copies of the free-space Green function with sources at the points $(\tilde{\mathbf{x}}^{\prime},z^{\prime}-qd)$ below the grating surface $\Gamma$, while the terms corresponding to $q=0$ produce the necessary sources on $\Gamma$. It follows that the layer potentials associated with the shifted Green function (cf. equation (56) below) are well defined for all points ${\mathbf{x}}=(\tilde{\mathbf{x}},z)$ on and above the grating surface, and satisfy the Helmholtz equation for ${\mathbf{x}}$ above $\Gamma$.

It is important to note that, as shown in Section 3.2, the addition of the shifted terms effectively suppresses all contributions in equation (4) that contain vanishing denominators at Wood frequencies, and thereby reinstates convergence of the quasi-periodic Green function even at such frequencies. But such modes are required in the complete solution of the grating scattering problem. Therefore corresponding quasi-periodic plane-wave terms need to be added to the Green function to incorporate (now with controlled coefficients) all the necessary grazing modes, as detailed in Section 3.2. The corresponding treatments for simply periodic scattering surfaces and arrays of cylinders in two dimensions is presented in references [4] and [5], respectively.

3.1 Convergence of the modified Green function at Wood frequencies

Lemma 3.1 below states that, as desired, the shifted half-space Green functions (13) enjoy enhanced degrees of decay as $\tilde{\mathbf{x}}\to\infty$. On the basis of this result, Theorem 3.2 then establishes the fast convergence of two different truncations (using discontinuous and and smooth window functions) for the spatial lattice sum (14). In preparation for the proofs, we examine the individual terms $\tilde{G}^{q}_{p}({\mathbf{x}}+m\mathbf{v}_{1}+n\mathbf{v}_{2})$. Considering the relations (2), (3) and (13) and using the notations

[TABLE]

the translated Green-function terms in the sum (14) can be expressed in the form

[TABLE]

Equivalently, letting

[TABLE]

we have

[TABLE]

In order to estimate the asymptotics of the function $h$ as $\rho\to\infty$ we use the finite-difference relation [17, p. 262, eq. 7]

[TABLE]

which is valid for every $p$-times continuously differentiable function $f$. Thus, for each pair of values of $\rho$ and $\hat{\varepsilon}$ there exists $\xi_{\rho,\hat{\varepsilon}}\in[0,p\hat{\varepsilon}]$ such that

[TABLE]

and, therefore

[TABLE]

In view of the asymptotic bounds on $\partial^{p}g/\partial\varepsilon^{p}$ in Lemma 3.1 below, one obtains h_{p}\big{(}\rho,\,z/\rho,\,d/\rho\big{)}\,=\,{\mathcal{O}}\left(1/\rho^{\lceil\frac{p}{2}\rceil+1}\right) as $(\rho\to\infty)$ where the constant in the $\cal{O}$-term, which depends on $z$, $d$ and $p$, can be taken to be fixed if $d$ and $p$ are given and $z$ is contained in a bounded subset of $\mathbb{R}$.

Lemma 3.1**.**

The $p$-th order derivative of the function $g$ with respect to $\varepsilon$ satisfies

[TABLE]

where the constant $C$ is independent of $\rho$ and $\varepsilon$ for all $\rho\geq 1$ and all $\varepsilon$ satisfying $|\varepsilon|<1$. Here, for real $x$, $\lceil x\rceil$ denotes the smallest integer larger than or equal to $x$. The estimate (22) is sharp: the left-hand side in that equation does not decay like $1/\rho^{t}$ for any $t>\lceil\frac{p}{2}\rceil+1$. Furthermore,

[TABLE]

Proof. It suffices to establish that (22) holds and is sharp; Equation (23) then follows directly from (21) and (22). In order to obtain the relation (22) we first note that the $p^{\text{th}}$ derivative $\frac{\partial^{p}g}{\partial\varepsilon^{p}}$ of $g$ with respect to $\varepsilon$ can be expressed as a finite linear combination of the form

[TABLE]

in which $S_{p}$ is a certain set of triples of non-negative integer indices and $C_{m,n,\ell}^{p}$ denote real valued coefficients. Defining

[TABLE]

we may thus write

[TABLE]

But, it is easy to check that

[TABLE]

and that, defining the unary operators

[TABLE]

we have

[TABLE]

We now may (and do) redefine the set $S_{p}$ to ensure that it contains exactly the triples $(m,n,\ell)$ corresponding to a sequence of $p$ applications of the unary operators (30). Thus calling $F=\{f^{-},f^{0},f^{+}\}$ we let

[TABLE]

Application of the operator $f^{+}$ results in a decrease by one in the power of $\varepsilon$ in the expression (25) for $T_{m,n,\ell}(\rho,\varepsilon)$, while application of either $f^{0}$ or $f^{-}$ results in an increase by one in that power. Susbstituting $\varepsilon=z/\rho$ in (25), on the other hand, we obtain

[TABLE]

and we thus define the $\rho$-order $Q$ of $T_{m,n,\ell}$ by

[TABLE]

Notice that an application of the operator $f^{+}$ (resp. $f^{0}$, $f^{-}$) to $T_{m,n,\ell}$ results in an increase (resp. no change, decrease) in the $\rho$-order $Q$.

Let $(m_{0},n_{0},\ell_{0})$ be such that $Q_{+}=Q(T_{m_{0},n_{0},\ell_{0}})$ is greater than or equal to $Q(T_{m,n,\ell})$ for all $(m,n,\ell)\in S_{p}$. We claim that $C_{m_{0},n_{0},\ell_{0}}^{p}\neq 0$. To establish this fact, first note that, in view of equations (24) and (31), the coefficient $C_{m,n,\ell}^{p}$ in the sum (26) equals a sum of several contributions, each one of which results from a sequence of $p$ operations from the set $\{f^{0},f^{+}\}$ applied to the root expression $T_{0,0,1}(\rho,\varepsilon)$. Indeed, while in principle all three elements in the set $\{f^{-},f^{0},f^{+}\}$ appear as a contribution to a coefficent $C_{m,n,\ell}^{p}$, it is easy to check that $f^{-}$ cannot appear as a contribution towards the maximum order coefficient $C_{m_{0},n_{0},\ell_{0}}$, since, as pointed out above, the $f^{-}$ operator decreases the $\rho$-order $Q$.

Since none of the $f^{0}$ and $f^{+}$ contributions are negative, it only remains to check that there is at least one positive contribution to the coefficient of $C(m_{0},n_{0},\ell_{0})$ of $T(m_{0},n_{0},\ell_{0})$ of $\rho$-order $Q_{+}$. But, in view of equation (30), a nonzero contribution of the form $f_{1}f_{2}\dots f_{p}T_{0,0,1}$ to the coefficient $C(m_{0},n_{0},\ell_{0})$ can only result provided no more than half of $p$ operators used equals $f^{+}$—which implies, in particular, that $Q_{+}\leq p/2$. In fact we have $Q_{+}=p/2$ (resp. $Q_{+}=(p-1)/2$) for $p$ even (resp. for $p$ odd), and a positive contribution to $C(m_{0},n_{0},\ell_{0})$ is provided by $T(m_{0},n_{0},\ell_{0})=(f^{+})^{\frac{p}{2}}(f^{0})^{\frac{p}{2}}T(0,0,1)$ for $p$ even and by $T(m_{0},n_{0},\ell_{0})=(f^{+})^{\frac{p-1}{2}}(f^{0})^{\frac{p+1}{2}}T(0,0,1)$ for $p$ odd. Taking into account the factor of $1/\rho$ in each one of the terms in equation (26), equation (22) follows and, thus, in view of (21), so does (23). The proof is now complete. $\blacksquare$

Now we are able to prove the algebraic convergence of the lattice sum for $\tilde{G}^{q}_{p}({\mathbf{x}})$.

Theorem 3.2** (Modified Green function for all frequencies; algebraic convergence).**

Let $\chi(r)$ be a smooth truncation function equal to $1$ for $r<r_{1}$ and equal to [math] for $r>r_{2}$ ($0<r_{1}<r_{2}$), and let $p$ denote an integer such that $p\geq 3$. Then, for all real triples $(k,\alpha,\beta)=(k,\boldsymbol{\alpha})$ ($k\not=0$) the sums

[TABLE]

where $(r_{mn}^{q})^{2}=|\tilde{\mathbf{x}}+m\mathbf{v}_{1}+n\mathbf{v}_{2}|^{2}+(z+qd)^{2}$  and  $\tilde{r}_{mn}=|\tilde{\mathbf{x}}+m\mathbf{v}_{1}+n\mathbf{v}_{2}|$ , converge to a radiating quasi-periodic modified Green function $\tilde{G}^{q}_{p}({\mathbf{x}})$ which satisfies the Partial Differential Equation

[TABLE]

as well as the quasi-periodicity condition  $\tilde{G}^{q}_{p}({\mathbf{x}}+(m\mathbf{v}_{2}+n\mathbf{v}_{2},0))=\tilde{G}^{q}_{p}({\mathbf{x}})\,e^{i\boldsymbol{\alpha}\cdot(m\mathbf{v}_{1}+n\mathbf{v}_{2})}$.  Further, there exists a constant $C_{p}=C_{p}(k,\alpha,\beta)$ for which

[TABLE]

for all sufficiently large values of $A$.

Proof. The sum (33) can be re-expressed in the form

[TABLE]

But, in view of (23) and letting $\nu=\lceil p/2\rceil+1$ we see that

[TABLE]

for some function $H$ which, for certain constants $C$ and $M$ satisfies $|H({\mathbf{x}};m,n)|<C$ as long as $|\mathbf{v}_{mn}|>M$. Thus, the sum in (38) converges to a limit $\tilde{G}^{q}_{p}({\mathbf{x}})$ as $A\to+\infty$, and for $A>M$ we have

[TABLE]

—a relation which establishes the desired result (36).

The proof of the bound (37) follows in part the proof of Theorem 2.1 in Part I [8]. For simplicity we assume $\tilde{\mathbf{x}}=(x,y)=0$ and $\boldsymbol{\alpha}=0$; the extension of the proof to nonzero values of these quantities is handled easily via consideration of standard properties of the Fourier transform, as described explicitly at the end of the proof in [8]. Let $U$ denote the finite set

[TABLE]

This set is nonempty exactly when $k$ is a Wood frequency. Since $k\not=0$, one has $(0,0)\not\in U$. The assumption $p\geq 3$ implies $\nu\geq 3$.

Let $z_{q}=z+dq$. We may re-express (34) in the form

[TABLE]

Using the Poisson Summation Formula, this sum is transformed into a lattice sum in the Fourier variable:

[TABLE]

in which

[TABLE]

The proof of in [8, Theorem 2.1] establishes that the sum $S_{1}$ ($U$ is empty in that work) converges superalgebraically to a limit $L$, that is, it is equal to $L+{\mathcal{O}}\left(A^{-n}\right)$ for each positive integer $n$. Thus, it only remains to show that the sum $S_{2}$ converges, with the difference from its limit being of order ${\mathcal{O}}\left(A^{-(\nu-3/2)}\right)$ (or, equivalently, ${\mathcal{O}}\left(A^{-(\lceil p/2\rceil-1/2)}\right)$) as $A\to+\infty$. Each fraction in $S_{2}$ can be expressed as an exponential in $r=|\boldsymbol{r}|$, multiplied by a Laurent expansion,

[TABLE]

The coefficients $a^{q}_{j}$ depend on $z$, and the expansion is convergent when $r>|z_{q}|$. In view of (18) and (23) we see that

[TABLE]

and clearly

[TABLE]

with uniform convergence over the set $r\geq r_{1}$.

In order to study the contribution by the sum $S_{2}$ we define the polar coordinates $\boldsymbol{r}=(r,\theta)$ and we note that, since $\boldsymbol{\xi}\in U$, we must necessarily have $\boldsymbol{\xi}=(k/2\pi,\gamma)$. Then, using the rescaling $\rho=r/A$ and the notation $\psi(\rho)=1-\chi(\rho)$, and in view of the fact that $\nu\geq 3$, we obtain

[TABLE]

where $J_{0}$ denotes the Bessel function of order [math]. Clearly, the number $\tilde{L}$ does not depend on $A$. In view of the well known Bessel function asymptotics $J_{0}(x)\sim(2/x\pi)^{1/2}\cos(x-\pi/4),\ x\to+\infty$, and since $k\not=0$, we obtain the relation

[TABLE]

It follows that the sum $S_{2}$ converges with an error of ${\mathcal{O}}\left(1/A^{\lceil p/2\rceil-1/2}\right)$ as $A\to\infty$ (since $\nu=\lceil p/2\rceil+1$). Together with the superalgebraic convergence of $S_{1}$ as $A\to\infty$, this fact establishes (37). The proof is now complete. $\blacksquare$

3.2 Complete Green function in Fourier space

In view of the Fourier expression (4) for the Green function away from Wood anomalies we obtain the corresponding expression

[TABLE]

for the shifted Green function. For $z>0$ this expression can be made to read

[TABLE]

In view of the limit $\lim_{\gamma_{j\ell}\to 0}\frac{(1-e^{i\gamma_{j\ell}d})^{p}}{\gamma_{j\ell}}=0$ for $p\geq 2$, we see that that (50) can be evaluated even at Wood anomalies: letting

[TABLE]

we may continuously extend the function $\tilde{G}^{q}_{p}(\tilde{\mathbf{x}},z)$ to all frequencies, including Wood configurations, by means of the expression

[TABLE]

Unfortunately, however, if $\gamma_{j\ell}=0$ for some $(j,\ell)$, the corresponding Fourier component is not present in the Green function (52) and, therefore, use of this Green function cannot give rise to a uniquely solvable system of integral equations. To tackle this difficulty we follow [4] and introduce a modified version $G^{q}_{p}$ of $\tilde{G}^{q}_{p}$, that is outgoing for $z\to\infty$ (but not for $z\to-\infty$) and which contains all necessary Fourier harmonics, even at Wood frequencies. The modified Green function is given by

[TABLE]

where $v(\tilde{\mathbf{x}},z)$ denotes a solution of the homogeneous Helmholtz equation of the form

[TABLE]

where $b_{j\ell}\neq 0$ are arbitrary non-zero complex constants. The function ${G}^{q}_{p}(\tilde{\mathbf{x}},z)$ is $\boldsymbol{\alpha}$-quasi-periodic in $\tilde{\mathbf{x}}$ with periods $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$, it satisfies (35) and, crucially, it contains all Fourier harmonics.

As shown in Section 3.4 and demonstrated numerically in Section 5, the complete Green function ${G}^{q}_{p}$ can be used to obtain uniquely-solvable integral-equation formulations around Wood frequencies. The analysis presented in Section 3.4 relies, in part, on use of yet another Green function, namely, a non-radiating Green function defined by a slowly-convergent series which, however, 1) is well defined at Wood frequencies; and 2) unlike the shifted Green function ${G}^{q}_{p}$, is a Helmholtz solution for $z<0$, and is therefore well suited for use as part of a proof that concerns the PDE domain and its complement at a Wood frequency. This rather peculiar Green function is introduced in the following section, and it is then used in the uniqueness proof presented in Section 3.4.

3.3 All space (non-radiating) Green function at and around Wood

frequencies

In view of the relation

[TABLE]

each term in the classical quasi-periodic Green function (4) equals the sum of two quantities, the first of which diverges and the second of which tends to $i|z|$ as $\gamma_{j\ell}\to 0$. In view of the relation (52), both terms can be made to vanish by using a $p$-th order finite-difference of shifted Green functions. As an alternative, a direct removal of the diverging term (which amounts to addition of a solution of Helmholtz’ equation) does produce a Helmholtz Green function for $(k,\alpha,\beta)$ in a neighborhood of a given Wood anomaly triple $(k_{0},\alpha_{0},\beta_{0})$ at which $\gamma_{j\ell}=0$. The resulting Green function at the Wood configuration $(k_{0},\alpha_{0},\beta_{0})$ is thus given by

[TABLE]

This is not an outgoing Green function, on account of the terms that contain $|z|$ as a factor. But as indicated in the previous section, the function $B^{q}(\tilde{\mathbf{x}},z)$ is a solution of the Helmholtz equation except at the periodically distributed singular points $(m\mathbf{v}_{1}+n\mathbf{v}_{2},0)$, and, unlike the shifted Green function ${G}^{q}_{p}$, it does not contain any additional singularities.

3.4 Uniquely solvable integral equations around Wood frequencies

We seek a scattered field above $\Gamma$ in the form of a combined single- and double-layer potential

[TABLE]

in terms of a quasi-periodic density $\phi$ defined on $\Gamma^{per}$, where

[TABLE]

The domain $\Gamma^{per}$ is that part of $\Gamma$ that lies above the unit cell

[TABLE]

of the periodic lattice $\Lambda$.

The function $u$ defined in (56) is quasi-periodic and outgoing as $z\to\infty$, and it satisfies the Helmholtz equation in $\mathbb{R}^{3}\setminus\bigcup_{q=0}^{p}(\Gamma^{per}-(0,0,qd))$. This function is a solution of (9) if and only if $\phi$ solves the integral equation

[TABLE]

The following theorem establishes that this integral equation is uniquely solvable as long as the shift distance $d>0$ is selected in such a way that the restrictions

[TABLE]

are satisfied. It is important to note that, since $\gamma_{j\ell}$ is an imaginary quantity for $j,\ell$ large enough, equation (60) amounts to a finite number of constraints.

Theorem 3.3**.**

Let $p\geq 0$, let $\xi\neq 0$ and $\eta\neq 0$ denote real numbers satisfying $\eta/\xi<0$, and let $d>0$ be such that (60) holds. Then equation (59) admits a unique solution for all triples $(k,\alpha,\beta)$ of wavenumbers, including Wood anomalies.

Proof.

Given that the surface $\Gamma$ is smooth and that $G^{q}_{p}$ has the same singularity as $G$, it follows that the integral operators on the left-hand side of equation (59) are compact operators in the space $L^{2}(\Gamma^{per})$. Thus, by the Fredholm theory, the unique solvability of equation (59) is equivalent to the injectivity of the operator on the left-hand side of that equation. In order to establish injectivity, and thereby complete the proof of the theorem, in what follows we show that any solution $\phi\in L^{2}(\Gamma^{per})$ of the homogeneous equation

[TABLE]

must necessarily vanish.

Let $\phi$ satisfy (61) and let $u$ denote the corresponding Helmholtz solution (56). It follows that $u$ vanishes on $\Gamma^{per}$. Since $u$ satisfies the outgoing radiation condition (10) at $+\infty$ (because ${G}^{q}_{p}({\mathbf{x}})$ does), by uniqueness of solution of the Dirichlet problem (9), which holds even at Wood frequencies [24], it follows that $u({\mathbf{x}})=0$ for all ${\mathbf{x}}\in\Omega^{+}$. In particular, letting $u_{j\ell}(z)$ denote the Fourier coefficients of the doubly periodic function $u(\tilde{\mathbf{x}},z)e^{-i\boldsymbol{\alpha}\cdot\tilde{\mathbf{x}}}$ with respect to $\tilde{\mathbf{x}}$ for $z>z_{+}$, we have $u_{j\ell}(z)=0$ for all $j,\ell\in\mathbb{Z}$.

Given that for $z>0$

[TABLE]

(see Section 3.2) it follows that

[TABLE]

in which

[TABLE]

Under the assumptions of this theorem, $u$ may vanish only if $c_{j\ell}^{+}=0$ for all $j,\ell\in\mathbb{Z}$.

Using the non-radiating Green function introduced in Section 3.3 we then set

[TABLE]

In view of (55), for $z>z_{+}$ the function $v$ admits the Fourier expansion

[TABLE]

in which

[TABLE]

in which (recalling that $|z|=z-z^{\prime}$ when ${\mathbf{x}}$ is above $\Gamma$, or $z>f(\tilde{\mathbf{x}})$)

[TABLE]

Given that $c_{j\ell}^{+}=0$ for all $j,\ell\in\mathbb{Z}$, it follows that

[TABLE]

Therefore $v(\tilde{\mathbf{x}},z)\,=\,\sum_{(j,\ell)\in U}(-ic^{\prime}_{j\ell}+c^{\prime\prime}_{j\ell})\,e^{i\mathbf{v}_{j\ell}^{*}\cdot\tilde{\mathbf{x}}}$ for $z>z_{+}$. The field $v({\mathbf{x}})$ satisfies the Helmholtz equation for ${\mathbf{x}}\not\in\Gamma$, and $v({\mathbf{x}})$ is independent of $z$ for $z>z_{+}$. In view of the real-analyticity of $v({\mathbf{x}})$ for ${\mathbf{x}}\not\in\Gamma$ and the uniqueness of analytic continuation, it follows that $v({\mathbf{x}})$ is independent of $z$ everywhere above $\Gamma$ and we have

[TABLE]

Now we turn our attention to the Fourier expansion

[TABLE]

of the function $v$ in the region $z<z_{-}$. For $(j,\ell)\in U$, the Fourier coefficients in this expansion satisfy $v_{j\ell}^{-}(z)=-v_{j\ell}^{+}$, since the term $|z-z^{\prime}|$ in $B^{q}({\mathbf{x}}-{\mathbf{x}}^{\prime})$ equals $-(z-z^{\prime})$ for $z<z^{\prime}$. (There is no such relation for $(j,\ell)\not\in U$.) For $(j,\ell)\in U$ we thus have

[TABLE]

The function $v({\mathbf{x}})$ satisfies the radiation condition (12) for $z<z_{-}$ since, in spite of linear terms that are part of the Green function $B^{q}$, $v$ itself does not contain such linear growth for $z<z_{-}$.

The right-hand side expression in (69), considered as a function defined on all of $\mathbb{R}^{3}$, defines a Helmholtz field that satisfies the outgoing condition for $z\to\infty$ and $z\to-\infty$. Thus, by subtracting it from $v({\mathbf{x}})$, one obtains the field

[TABLE]

that satisfies the radiation conditions at both $+\infty$ and $-\infty$, and that vanishes for $z>f(\tilde{\mathbf{x}})$. The jump conditions of the single- and double-layer potentials that define $v({\mathbf{x}})$ in (64) imply that the limits of $\tilde{v}({\mathbf{x}})$ and its normal derivative from below $\Gamma$ are

[TABLE]

It follows that the function $\tilde{v}({\mathbf{x}})$ restricted to the domain $\{z\leq f(\tilde{\mathbf{x}})\}$ satisfies the homogeneous impedance boundary condition

[TABLE]

Since, per the discussion above, $\tilde{v}$ additionally satisfies the radiation condition (12), Lemma A.1 tells us that we must have $\tilde{v}({\mathbf{x}})=0$ for all ${\mathbf{x}}$ below $\Gamma$. In view of (70) and/or (71) and the assumption $\eta/\xi<0$ it follows that $\phi({\mathbf{x}})=0$ for all ${\mathbf{x}}\in\Gamma$. ∎

4 High-order numerical evaluation of the boundary-layer

potentials with quasi-periodic Green functions

For our numerical treatment we reformulate the quasi-periodic scattering integral equation (59) in terms of only periodic functions. We use the fact that the solution $\phi=\phi^{qper}$ of the integral equation (59) is $\boldsymbol{\alpha}$-quasi-periodic with respect to the lattice $\Lambda$, and that, therefore, the quantity

[TABLE]

is periodic with respect to $\Lambda$. This allows us to express the integral equation (59) in the form

[TABLE]

where the $\Lambda$-periodic Green function $G^{per}_{p}$ is defined by

[TABLE]

In practice, for any given $p=0,1,\dots$ the quantity $\hat{G}^{p,A}$ (equation (34)) is used as an approximation for $G^{q}_{p}$ which, when subsituted in (73), results in the needed numerical approximation of $G^{per}_{p}$. We solve equation (72) by means of the unaccelerated high-order Nyström procedure introduced in [6] and Part I. For all the numerical experiments presented in Section 5 a single patch was used to represent the biperiodic surfaces under consideration.

5 Numerical results

We present numerical computations of scattering by the doubly periodic scattering surface $f(x,y)=\frac{1}{2}\cos(2\pi x)\cos(2\pi y)$ with periodicity lattice vectors $\mathbf{v}_{1}=(1,0,0)$ and $\mathbf{v}_{2}=(0,1,0)$ and under Dirichlet and Neumann boundary conditions. For non-Wood configurations we utilize the $p=0$ (unshifted) version of the algorithm described in Section 4. At and around Wood configurations, on the other hand, we use the $p=3$ version of that algorithm. A fully three-dimensional single-core Matlab implementation of these methods was used, which was neither accelerated nor optimized; accordingly, our numerical error studies at Wood anomalies do not go beyond relative errors of the order $10^{-4}$. Clear high-order convergence is observed in all cases.

We report the quality of the solutions on the basis of two error indicators. The first of these indicators is the energy-conservation defect

[TABLE]

which we have verified (by means of numerical resolution studies) to be an excellent error predictor for these solvers. An additional error estimator we present, $\varepsilon_{1}$, on the other hand, equals the absolute error in the Rayleigh coefficient $B^{+}_{0,0}$ (as estimated by comparison with a reference solution obtained by means of a highly-refined discretization, a large value of the window parameter $A$, and a sufficiently small GMRES tolerance). The word “ref” on a table entry indicates that the parameter values on that row were used to produce the reference solution necessary for evaluation of the errors $\varepsilon_{1}$ for the corresponding frequency $k$ on the that table. The numbers of iterations required by the GMRES solvers to reach specified tolerances are provided in each case.

Table 1 demonstrates the high-order character, as the window-size parameter $A$ grows, for the proposed $p=0$ (unshifted) Dirichlet solvers at frequencies $k$ away from Wood anomalies. Table 2 demonstrates the high-order character of the $p=3$ (shifted) solver, as $A$ grows, also under Dirichlet boundary conditions and for values of $k$ away from Wood anomalies. Tables 3 and 4, in turn, concern configurations at and near Wood anomalies; Dirichlet (resp. Neumann) boundary conditions are considered in the first (resp. second) of these tables. In the normal-incidence case considered in those tables, the first three Wood anomalies occur at $k=2\pi$, $k=2\sqrt{2}\pi$, and $k=4\pi$. Once again, fast convergence is observed as $A$ grows, even at and around Wood anomalies. We note that the number of iterations required by the GMRES solvers based on Combined Field Integral Equations remains small even for Wood and near-Wood parameters for both Dirichlet and Neumann problems.

Appendix A Appendix: Integral equations away from Wood frequencies

In order to establish the unique solvability of the integral equations (59), the proof of Theorem 3.3 relies on the following classical result on solutions to the homogeneous surface-impedance problem.

Lemma A.1**.**

Let $v({\mathbf{x}})$ is a quasi-periodic field that satisfies the Helmholtz equation $\Delta v+k^{2}v=0$ for $z<f(\tilde{\mathbf{x}})$ (resp. $z>f(\tilde{\mathbf{x}})$), the outgoing condition (12) (resp. (10)) and the impedance condition

[TABLE]

with $\zeta>0$ (resp. $\zeta<0$). Then $v({\mathbf{x}})=0$ for $z<f({\mathbf{x}})$ (resp. $z>f({\mathbf{x}})$).

Proof.

We establish the result in the case $z<f({\mathbf{x}})$ and $\zeta>0$; the complementary case is handled analogously. Consider the truncated period

[TABLE]

with lower boundary $S=\{(\tilde{\mathbf{x}},z):\tilde{\mathbf{x}}\in Q,\,z=z_{-}\}$ oriented downward. Using integration by parts we obtain

[TABLE]

(The integrals over the lateral sides of $\Omega$ add up to zero as a result of the assumed quasi-periodicity of $v$.) In view of the impedance condition on $\Gamma$ and the outgoing condition (12) we obtain

[TABLE]

The imaginary part of this quantity equals the imaginary part of (75) and it must therefore vanish. It follows that $v=0$ on $\Gamma$. The impedance condition then yields $\partial v/\partial{n}=0$ on $\Gamma$ as well. By Green’s identity it follows that $v=0$ in $\Omega$ and therefore for all ${\mathbf{x}}$ with $z<f(\tilde{\mathbf{x}})$ ($\tilde{\mathbf{x}}\in\mathbb{R}^{2}$). ∎

For completeness we now present a simpler alternative proof of Theorem 3.3, which, however, is restricted to the case $p=0$ (unshifted Green function) and to configurations away from Wood anomalies.

Theorem A.2**.**

Let $\frac{\eta}{\xi}<0$, let $p=0$, and let us assume that $k$ is a wavenumber for which the quasi-periodic Green function $G^{q}_{0}$ exists, that is, $\gamma_{j\ell}\not=0$ for all pairs $(j,\ell)\in\mathbb{Z}$. Then the integral equation (59) is uniquely solvable in $L^{2}(\Gamma^{per})$.

Proof.

Given that the surface $\Gamma$ is smooth and that $G^{q}_{0}$ has the same singularity as $G$, it follows that the integral operators on the left-hand side of equation (59) are compact operators in the space $L^{2}(\Gamma^{per})$. Thus, by Fredholm theory, the unique solvability of equation (59) is equivalent to the injectivity of the operator on the left-hand side of that equation. In order to establish injectivity, let $\phi_{0}\in L^{2}(\Gamma^{per})$ denote a solution of equation (59) with zero right hand-side, and let $u^{\pm}$ denote the restrictions to $\Omega^{\pm}$ of the potentials $u$ defined by equation (56) above and below $\Gamma$ with density $\phi=\phi_{0}$. It follows that $u^{+}$ is a radiating solution of the Helmholtz equation in $\Omega^{+}$ with zero Dirichlet boundary conditions, and hence $u^{+}=0$ in $\Omega^{+}$ [24]. Using the jump relations satisfied by the layer potentials in equation (56) we obtain $u^{-}|_{\Gamma}=-\xi\phi_{0}$ and $\left(\frac{\partial u^{-}}{\partial{n}}\right)|_{\Gamma}=i\eta\phi_{0}$. Thus $u^{-}$ is a quasi-periodic radiating solution of the Helmholtz equation in $\Omega^{-}$ with zero impedance boundary conditions $\frac{\partial v}{\partial{n}}({\mathbf{x}})-i\zeta\,v({\mathbf{x}})=0,\ \zeta=-\frac{\eta}{\xi}>0$ on $\Gamma$. By Lemma A.1 it follows that $v^{-}=0$ in $\Omega^{-}$, and thus $\phi_{0}=0$ in $L^{2}(\Gamma^{per})$, as desired. The proof of the theorem is now complete. ∎

No uniqueness results exist for the Helmholtz scattering problem (9) under Neumann boundary-value conditions and the radiation condition (10), even away from Wood anomalies—although it has been repeatedly conjectured (cf. [18, p. 147]), [4]) that such a uniquess result does hold. Assuming that the wavenumber $k$ is such that this scattering problem does admit a unique solution, however, we may seek the scattered field in the form

[TABLE]

in terms of the unknown surface density $\psi$. Using the jump condition for the normal derivatives of single-layer potentials and the sound-hard (Neumann) boundary condition, the unknown density $\psi$ is seen to be a solution of the integral equation

[TABLE]

Given that, addtionally, the scattering problems from doubly periodic surfaces and Dirichlet boundary conditions admit unique solutions for all wavenumbers, it follows that the integral equations (77) have themselves unique solutions.

Data accessibility. All data applicable to this paper is included in the article.

Competing interests. There are no competing interests relevant to this article.

Authors contributions. All authors are equally considered co-contributors in this article.

Funding statement. This effort was supported by AFOSR, NSF and a NSSEFF Vannevar Bush Fellowship under contracts FA9550-15-1-0043, DMS-1411876 and N00014-16-1-2808 (OB); NSF DMS-0807325 (SPS); NSF DMS-1312169 and DMS-1614270 (CT); and NSF DMS-1211638 (SV).

Ethics statement. No ethics statement applies to this article.