Features in the diffraction of a scalar plane wave from doubly-periodic Dirichlet and Neumann surfaces

The diffraction of a scalar plane wave from a doubly-periodic surface on which either the Dirichlet or Neumann boundary condition is imposed is studied by means of a rigorous numerical solution of the Rayleigh equation for the amplitudes of the diffracted Bragg beams. From the results of these calculations the diffraction efficiencies of several of the lowest order diffracted beams are calculated as functions of the polar and azimuthal angles of incidence. The angular dependencies of the diffraction efficiencies display features that can be identified as Rayleigh anomalies for both types of surfaces. In the case of a Neumann surface additional features are present that can be attributed to the existence of surface waves on such surfaces. Some of the results obtained through the use of the Rayleigh equation are validated by comparing them with results of a rigorous Green's function numerical calculation.


I. INTRODUCTION
In several resent papers the present authors have studied theoretically and computationally the diffraction of ppolarized light from a perfectly conducting grating [1] and from a high-index dielectric grating [2], and the diffraction of a shear horizontally polarized acoustic wave from a gating ruled on the surface of an isotropic elastic medium [3]. Each of these media does not support a surface wave when the surface bounding it is planar, but supports one when it is periodically corrugated. The dependence of the diffraction efficiencies of some of the lowest-order Bragg beams on the angle of incidence was found to posses two types of anomalies. The first type of anomaly occurred at angles of incidence at which a diffractive order begins to propagate or ceases to propagate. They were first observed by Wood in 1902 [4] in the diffraction of light from a metallic grating. Their origin was explained by Lord Rayleigh [5], and they are now called Rayleigh anomalies. The second type of anomaly occurred at angles of incidence at which the incident field excites a surface wave supported by the grating, when the difference of the components of the wave vectors of the incident field parallel to the mean scattering plane and of the surface wave is made up by the addition of a grating vector. These anomalies were also first observed by Wood in 1902 [4], and were shown to be due to the grating-induced excitation of surface plasmon polaritons by the incident light by Fano [6]. We refer to them as Wood anomalies. Since surface waves don't exist on planar surfaces of the media studied in Refs. [1][2][3], our results obtained in these papers emphasized the necessity of surface waves for the existence of Wood anomalies, and that the surface waves need not be surface plasmon polaritons, but can be of a quite different nature.
In this paper we extend the work presented in Refs. [1][2][3] to the case of diffraction of a scalar plane wave from a doubly-periodic grating (often called a bigrating or a cross-grating), fabricated on a surface of the medium that when planar does not support a surface wave, but can support one when it is doubly periodic. Specifically, we consider the diffraction of a scalar plane wave from a doubly-periodic surface on which either the Dirichlet or the Neumann boundary condition is satisfied. For brevity, we will refer to these two types of surfaces as Dirichlet and Neumann surfaces, respectively. It is known that a doubly-periodic Neumann surface supports a surface wave while a doubly-periodic Dirichlet surface does not [7]. We calculate the dependence of several of the lowest-order diffraction efficiencies on the polar angle of incidence for a fixed azimuthal angle of incidence, and look for Rayleigh anomalies in these dependencies for both types of surfaces, and for Wood anomalies in the diffraction from a Neumann surface.
Calculations of the dependence of the efficiencies of diffraction of light from, or its transmission through, several types of doubly-periodic structures on the wavelength of the incident light, or on its polar angle of incidence for a fixed azimuthal angle of incidence, have been carried out by a variety of approaches. Some of them have been devoted to the wavelength dependence of the reflectivity and the dip it displays that arises from the excitation of a surface wave supported by the doubly-periodic structure [8][9][10][11], others to the phenomenon of total absorption [12], and still others to the wavelength dependence of the enhanced transmission of light through a doubly-periodic array of nanoscale holes piercing a thin metal film [13][14][15][16]. Similar calculations of higher-order diffraction efficiencies, such as the ones carried out in the present work, do not appear to have been carried out in these earlier studies.
The calculations in the body of this paper will be carried out on the basis of the Rayleigh hypothesis [17], perhaps the simplest approach to solving this scattering problem. The validity of this approach has been questioned on occasion [18][19][20], primarily on the grounds that if the indentations of the surface are sufficiently deep some of the scattered waves can be propagating downward toward the surface within them, before becoming upward outgoing propagating waves due to multiple scattering. Such waves are not taken into account by the Rayleigh hypothesis, which assumes only upward outgoing propagating waves everywhere above the surface. Nevertheless, in subsequent work limits of validity of this hypothesis have been determined for the scattering of a scalar plane wave from a singly-periodic surface [21][22][23][24][25][26] and from a doubly-periodic surface [27], defined by profile functions that are analytic functions of the coordinates in the mean scattering plane. It has recently been argued that the Rayleigh hypothesis is always valid [28].
With this background, in this paper we derive the Rayleigh equations for the diffraction of a scalar plane wave from doubly-periodic Dirichlet and Neumann surfaces, and solve them numerically. For greater generality, and for pedagogical reasons, we begin the derivation by first obtaining the Rayleigh equation for the scattering of a scalar plane wave from an arbitrary rough two-dimensional Dirichlet and Neumann surface, and then specialize this equation to the case of a doubly-periodic surface. From the solutions of these equations we will calculate the angular dependencies of several of the diffraction orders. The validity of the Rayleigh hypothesis in the context of the problem studied will be demonstrated by a comparison of some of the results obtained by its use with those obtained by a rigorous numerical method [29,30].

II. SCATTERING THEORY
The system that we consider consists of a liquid in the region x 3 > ζ(x ) and an impenetrable medium in the region x 3 < ζ(x ) where x = (x 1 , x 2 , 0) is a position vector in the plane x 3 = 0. The surface profile function ζ(x ) is assumed to be a single valued function of x , and to be differentiable with respect to x 1 and x 2 .
In the scattering of a scalar plane wave from either type of surface, the incident field ψ(x|ω) inc can be written as where k = (k 1 , k 2 , 0) and α 0 (k , ω) = (ω 2 /c 2 ) − k 2 Similarly, the field scattered from either type of surface, ψ(x|ω) sc , can be written where α 0 (q , ω) = (ω 2 /c 2 ) − q 2 1/2 with Re α 0 (q , ω) > 0 and Im α 0 (q , ω) > 0. Of course the scattering amplitude R(q |k ) will be different for the scattering from a Dirichlet surface then it is for scattering from a Neumann surface due to the different boundary conditions satisfied on the two types of surfaces. We now substitute the sum of Eqs. (5) and (6) into the boundary conditions (2) and (3). This step represents our assumption of the validity of the Rayleigh hypothesis [17]. The resulting equations for the scattering amplitude can be written as We now introduce the function I(γ|Q ) by so that If we differentiate both sides of Eq. (8a) with respect to x α (α = 1, 2), we obtain the useful result When we substitute Eqs. (8) into Eqs. (7) and equate to zero the coefficient of exp[ip ·x ] in the resulting equations, the equations satisfied by the scattering amplitudes R(q |k ) become where M (p |q ) = 1, N (p |k ) = 1 (10) for a Dirichlet surface, and for a Neumann surface. Equations (9)-(11) constitute the Rayleigh equations for the scattering amplitude in the scattering of a scalar plane wave from a two-dimensional rough Dirichlet or Neumann surface.

III. THE DIFFERENTIAL REFLECTION COEFFICIENT
The scattering amplitude R(q |k ) is of great importance in calculations of scattering from rough surfaces because an experimentally accessible quantity, the differential reflection coefficient, is expressed in terms of it. The differential reflection coefficient ∂R/∂Ω s is defined such that (∂R/∂Ω s )dΩ s is the fraction of the total time averaged incident flux that is scattered into an element of solid angle dΩ s around the direction defined by the polar and azimuthal angles of scattering (θ s , φ s ).
The magnitude of the total time-averaged flux incident on the surface is given by where L 1 and L 2 are the lengths of the scattering surface along the x 1 and x 2 axes, while A is a coefficient that drops out of the expression for the differential reflection coefficient. (For the scattering of a particle of mass m, A = /m where denotes Planck's constant.). The minus sign that appears in on the right-hand side of Eq. (12) compensates for the fact that the incident flux is negative. For the form of the incident field biven by Eq. (5) we find easily that where we have used the fact that α 0 (k , ω) is real. Similarly, the magnitude of the total time-averaged scattered flux is given by With the use of the expression for ψ(x|ω) sc given by Eq. (6), this expression becomes In obtaining this result we have used the fact that α 0 (q , ω) is real for 0 < q < ω/c, while it is imaginary for q > ω/c, to obtain the domain of integration indicated. We now introduce the polar and azimuthal angles of incidence (θ 0 , φ 0 ) and of scattering (θ s , φ s ), respectively, through the relations It follows that α 0 (k , ω) = (ω/c) cos θ 0 , α 0 (q , ω) = (ω/c) cos θ s , and d 2 q = (ω/c) 2 cos θ s dΩ s where dΩ s , the element of solid angle, is dΩ s = sin θ s dθ s dφ s . With the use of these results the incident flux can be written as while the scattered flux becomes where By definition the differential reflection coefficient is where k and q are defined by Eqs. (16a) and (16b), respectively.

IV. DOUBLY-PERIODIC SURFACE
The results obtained in the preceding sections of this paper apply to an arbitrary two-dimensional rough surface defined by the single-valued surface profile function ζ(x ) that is differentiable with respect to x 1 and x 2 . In this section we specialize these results to the case where the surface profile function is a doubly-periodic function of x , namely a bigrating.
Thus the surface profile function ζ(x ) is assumed to possess the property ζ x + x ( ) = ζ(x ), where the vector x ( ) is a translation vector of a two-dimensional Bravais lattice [31]. It is defined by where a 1 and a 2 are the two noncolinear primitive translation vectors of the Bravais lattice, while 1 and 2 are any positive or negative integers or zero, which we denote collectively by . The area of a primitive unit cell of this lattice is a c = |a 1 × a 2 |. We also introduce the lattice reciprocal to the one defined by Eq. (21). Its lattice sites are defined by the vectors where the primitive translation vectors b 1 and b 2 are related to the primitive translation vectors of the direct lattice, a 1 and a 2 , by a i ·b j = 2πδ ij , while h 1 and h 2 are any positive or negative integers or zero, which we denote collectively by h.
We now proceed to transform the Rayleigh equation (9) for the scattering amplitude R(q |k ) into the Rayleigh equation for the corresponding amplitude that arises in the diffraction of a scalar plane wave from an impenetrable bigrating.
Due to the periodicity of the surface profile function ζ(x ), the field in the region x 3 > ζ(x ) must satisfy the Floquet-Bloch condition [32,33] This condition is satisfied if we rewrite the scattering amplitude R(q |k ) in the form In writing this equation we have replaced summation over h by summation over G . A second consequence of the periodicity of the surface profile function ζ(x ) is that the function I(γ|Q ) defined by Eq. (8b) can now be written where a c ( ) is the area of the unit cell containing the translation vector x ( ). The change of variable x = x ( ) + u , and the relation ζ(x + x ( )) = ζ(x ), yield the result The use of the relation [34] exp in Eq. (26) yields the result where When the results given by Eqs. (24) and (28) are substituted into Eq. (9), and the integration over q is carried out, we obtain the equation In writing this equation, to simplify the notation we have defined the two wave vectors and have replaced summation over G and G by summation over K and K , respectively. On equating the coefficients of δ(p − K ) on both sides of Eq. (30) we obtain the Rayleigh equation satisfied by r(K |k ) Equation (32) holds for all possible values of K (or G ), and hence it represents a linear system of equations of infinite dimension. To be able to solve the system numerically, we need a system of finite dimension. This can be achieved by restricting the vectors G (h) and G (h) to a domain for which their lengths are no more than several times ω/c, but at the same time no shorter than ω/c. In this way a finite dimensional linear system is obtained that can be solved for r(k + G (h)|k ) by standard methods.

V. DIFFRACTION EFFICIENCIES
The total time-averaged flux scattered from our doubly-periodic surface is obtained by substituting Eq. (24) into Eq. (15): The only nonzero terms on the right-hand side of this equation are those for which G = G . Then, with the result that in two-dimensions Eq. (33) becomes The prime on the sum indicates that it extends over only those values of G for which |k + G | < ω/c. Equation (35) demonstrates that each diffracted beam contributes independently to the scattered flux. When the scattered flux is normalized by the total time-averaged flux of the incident field, Eq. (13), the result can be written where The quantity e k + G |k , called the diffraction efficiency, is the fraction of the total time-averaged incident flux that is diffracted into a Bragg beam defined by the wave vector k + G (when the incident beam is defined by k ). It has a physical meaning for only those values of G for which α 0 (|k + G |, ω) is real. The propagating diffracted beams defined by this condition are called the open channels.
Since there is no absorption in the scattering from an impenetrable surface, all the power incident on it must be scattered back into the medium of incidence. Hence, the conservation of energy in the scattering process requires that G e k + G |k = 1.
The closeness to unity of the sum on the left-hand side of Eq. (38) is a good test of the quality of the numerical simulation calculations of the diffraction efficiencies.
The reflectivity of the bigrating is obtained from the diffraction efficiency for the beam defined by G = 0:

VI. RESULTS
We will illustrate the proceeding results by presenting simulation results for the dependence of the reflectivity and several other diffraction efficiencies on the polar and azimuthal angles of incidence θ 0 and φ 0 , respectively, when the bigrating defined by the surface profile function is illuminated by a scalar plane wave of frequency ω. The primitive translation vectors of the square Bravais lattice underlying this surface profile function are Those of the corresponding reciprocal lattice are An attractive feature of the form of the surface profile function in Eq. (40) is that the I-integral defined in Eq. (29) can be obtained analytically, and takes the form where J n (·) represents the Bessel function of the first kind and order n.
In the first set of calculations of the dependence of the reflectivity of the bigrating defined by Eq. (40) on the polar angle of incidence θ 0 for a value of the azimuthal angle of incidence φ 0 = 0 • [k = (1, 0, 0)], we assumed that the lattice constant a had the value a = 3.5λ, where λ is the wavelength of the incident wave, while the amplitude ζ 0 took several values. The calculated reflectivities for Neumann surfaces characterized by the amplitudes ζ 0 = 0.3λ, 0.5λ, 0.7λ are presented in Fig. 1. These results show a complex dependence of the reflectivity on the polar angle of incidence in the form of the presence of many sharp peaks and dips. These features are Rayleigh anomalies, which occur at values of θ 0 for a given value of φ 0 at which diffractive orders start or cease to propagate.
To determine the angles of incidence at which the Rayleigh anomalies occur, we note that the lateral wave vector of the diffractive beam characterized by the index pair h 1 , h 2 is given by This diffractive beam goes from being a propagating one to an evanescent one when |q (h 1 , h 2 )| = ω/c, which is the condition for a potential Rayleigh anomaly to be associated with this wave. On squaring both sides of Eq. (44) and using Eqs. (16a) and (22), we obtain a quadratic equation for sin θ 0 withk = (cos φ 0 , sin φ 0 , 0).

Equation (45) determines for a general grating where Rayleigh anomalies can exist. From its solutions
under the condition | sin θ 0 | ≤ 1, as h 1 and h 2 each run over all positive and negative integers and zero, the polar angles of incidence θ 0 at which Rayleigh anomalies can exist for a specified azimuthal angle of incidence φ 0 are obtained. The values of θ 0 obtained in this way are indicated by gray vertical dashed lines in Fig. 1. From the results of this figure we see that the majority of the peaks and dips present in the refelctivity are Rayleigh anomalies. It should be noted that even if a Rayleigh anomaly is predicted to exist at a particular polar angle of incidence, it may not be observed in the reflectivity, because it is too weak to be seen. We see from Fig. 1 that as the amplitude of the surface profile function ζ 0 is increased, the polar angles of incidence at which the Rayleigh anomalies occur do not change, as must be the case, but the forms of the anomalies can change. Peaks and dips can change their magnitudes, and broaden, and dips can change into peaks, and peaks can change into dips.
The numerical calculations that produced the results presented in Fig. 1 were performed under the assumption that G (h) ≤ 4ω/c, and the linear system of equations satisfied by r(K |k ), Eq. (32), was solved by the routine la_gesv from LAPACK95 [35]. For this value of max G (h) the simulation time required per angle of incidence to produce the results in Fig. 1 was 1.5 s, or less on average, when the simulations were performed on a machine equipped with an Intel i7-5930K CPU running at 3.50 GHz. The energy conservation condition (38) was found to be satisfied with an error no greater than 10 −10 for all the values of θ 0 and ζ 0 that we considered.
In Fig. 2 we present the dependence of the reflectivity on the polar angle of incidence when the azimuthal angle of incidence is φ 0 = 45 • . In this case the unit vectork becomesk = (1/ √ 2, 1/ √ 2, 0), and the values of θ 0 at which Rayleigh anomalies are predicted to occur are different from the values at which they occur in Fig. 1. With an increase of the amplitude of the surface profile function, these anomalies undergo the same kinds of changes in their forms as they do in the case where φ 0 = 0 • .
We now turn to the diffraction of a scalar plane wave from a doubly-periodic Dirichlet surface defined by Eq. (40). In Fig. 3 we present the reflectivity as a function of θ 0 for the case where φ 0 = 0 • . The parameters defining the surface profile function are a = 3.5λ and ζ 0 = 0.3λ, 0.5λ, 0.7λ, namely the values assumed in obtaining the results presented in Fig. 1 and 2. The values of θ 0 at which Rayleigh anomalies are predicted to exist are the same as those at which they are predicted to exist in Fig. 1. However, these anomalies are significantly weaker than those occurring at the same values of θ 0 in Fig. 1. This difference demonstrates the important role played by the boundary condition on the surface of the bigrating satisfied by the field in the region x 3 > ζ(x ) in forming these anomalies.
A further comparison of Figs. 1 and 3 prompts the following observation. In Fig. 1(a) we see a dip in the reflectivity at a value of θ 0 ≈ 88.0 • , an angle at which no Rayleigh anomaly is predicted to exist. In Fig. 1(b), for a larger amplitude of the bigrating profile function, this dip has broadened and shifted to a smaller value of θ 0 , namely θ 0 ≈ 84.5 • . Again, no Rayleigh anomaly is predicted to occur at this angle. With a further increase of the amplitude of the bigrating profile function, we see in Fig. 1(c) a break in the slope of the reflectivity curve at a value of θ 0 ≈ 80 • . Such a dip is more clearly visible at these three values of θ 0 in Figs. 2(a)-(c). No such feature is present at these (or other) angles in Figs. 3(a)-(c). It is known [7] that a doubly-periodic Neumann surface supports a surface wave, while a doubly-periodic Dirichlet surface does not. It is also known that changing the amplitude of the surface profile function shifts the nonradiative and radiative branches of the dispersion relation (in the reduced zone scheme) of the surface wave on a Neumann bigrating [7]. Since a Wood anomaly arises due to the excitation of a surface wave on a periodically modulated surface by the incident field [1,6] the angles of incidence at which these anomalies occur will shift with changes in the surface profile function. These properties of the large angle dip suggest that it represents a Wood anomaly. However, confirmation of this conjecture has to await the determination of the branches of the dispersion curve of the surface wave supported by the doubly-periodic Neumann surface, in the radiative region of the (k , ω) plane as well as in the nonradiative region.  Figs. 1(b), 2(b), and 3(b). Here also the scan over the polar angle of incidence was done in steps of ∆θ0 = 0.025 • , and G (h) ≤ 4ω/c was assumed in performing the numerical calculations.
Features similar to those observed in Figs. 1-3 are present in the dependence of other diffraction efficiencies on θ 0 for a given value of φ 0 . In Fig. 4 we plot this dependence for the efficiencies of the {1, 0}, {−1, 0}, {0, ±1}, {1, ±1}, and {−1, ±1} beams diffracted from the Neumann surface defined by Eq. (40) with ζ 0 = 0.5λ and a = 3.5λ. The azimuthal angle of incidence is φ 0 = 0 • . The notation {h 1 , ±h 2 } indicates that the {h 1 , h 2 } and {h 1 , −h 2 } efficiencies are identical. This identity is a consequence of the symmetry of the scattering system under reflection in the x 1 axis when φ 0 = 0 • . The predicted angular positions of the Rayleigh anomalies are indicated by the gray vertical dashed lines. It is seen that all of the peaks and dips in these dependencies occur at these angles, but not every one of these angles has an anomaly associated with it.
It is apparent from the results presented in Fig. 1, for instance, that the reflectivity of the doubly-periodic cosine Neumann surface depends strongly on its amplitude ζ 0 . To further investigate this dependence, we present in Fig. 5 as a solid line the reflectivity of such a surface of period a = 3.5λ as a function of the amplitude ζ 0 for polar and azimuthal angles of incidence θ 0 = 0 • and φ 0 = 0 • , respectively. These results were obtained on the basis of a numerical solution of the Rayleigh equation (32) for the same values of the numerical parameters assumed in obtaining the results in Fig. 1. Figure 5 shows that the reflectivity of the doubly-periodic Neumann surface decreases monotonically from unity to approximately 3 × 10 −5 when its amplitude increases from zero (planar surface) to ζ 0 = 0.371λ. Increasing the amplitude beyond this value causes the reflectivity of the surface to increase monotonically, and it reaches the value 0.1357 when ζ 0 = 0.7λ. What happens to the reflectivity when ζ 0 /λ > 0.7, we have not investigated here.
To validate our use of the Rayleigh equation in obtaining the results presented in this work, we performed additional calculations obtained on the basis of a rigorous Green's function-based numerical approach [29,30]. To this end, the latter approach was used to calculate the reflectivity for normal incidence as a function of the corrugation strength ζ 0 . The results of such calculations are presented as open symbols in Fig. 5, and they show satisfactory agreement with the corresponding results obtained on the basis of the Rayleigh equation approach. In particular, the five orders of magnitude variation of the reflectivity is consistently predicted by both approaches. It is only for ζ 0 /λ > 0.5 that some minor discrepancy starts to develop. As we will comment below, it it not entirely clear if this should be interpreted an indication that the Rayleigh equation approach starts to become less accurate.
We now briefly detail how the rigorous Green's function-based numerical calculations were performed; reference 30 gives additional details. Since the Green's function-based approach as formulated in Ref. 30 does not explicitly use the fact that the surface is periodic, the first step of the calculation is to restrict the doubly-periodic surface (40) to a square region of the x 1 x 2 plane of edges L. Next, this surface profile, as well as its derivatives up to order two, are discretized on a square lattice of points of lattice parameters ∆x. To avoid diffraction artifacts from the edges of the surface, the incident beam is assumed to be a Gaussian beam of 1/e half-width W . In the numerical calculations using the Green's function approach that we report in Fig. 5 the values of the numerical parameters were L = 38λ,  (32) while the open symbols were obtained by a rigorous Green's function-based numerical approach [29]. In performing the latter calculations it was assumed the edge of the square region of the x1x2 plane covered by the doubly-periodic surface was L = 38λ, the width of the Gaussian incident beam was W = 15λ, and ∆x = 0.15λ was the discretization interval used. W = 15λ, and ∆x = 0.15λ. The reflectivity of the surface was calculated by integrating the differential reflection coefficient ∂R/∂Ω s over an angular region around the specular direction in such a way that only the contribution from the fundamental diffractive order was included. Since the period of the surface that we consider is sufficiently long [a = 3.5λ], the diffractive orders were well separated. For the calculation of the reflectivity we used a region defined by |q − k | < 0.1 ω/c. We also checked and found that minor adjustments of the size of this angular region did not affect the reflectivity values obtained in any significant way. It should be remarked that for the largest values of ζ 0 that we considered we found a weak but detectable dependence of the reflectivity on the parameters L and W when it was calculated by the Green's function approach as described above. Hence the discrepancy seen in Fig. 5 in the reflectivity obtained by the two approaches for the largest values of the corrugation strength is not necessarily due to the Rayleigh approach becoming inaccurate. The calculations based on the Green's function approach, whose results are reported in Fig. 5, took about 34 min (or about 2000 s) of cpu-time to complete per value of ζ 0 when the calculation was performed on an Intel i7 960 processor running at 3.20 GHz, and the memory footprint of the calculation was almost 19 Gb. A similar calculation using the Rayleigh equation approach took about 1 s of cpu-time when performed on the same computer for the numerical parameters that we assumed and it required only a faction of the computer memory needed by the Green's function calculation.
Finally we should remark that the rigorous Green's function approach [30] described and used above is not ideal for a doubly-periodic system. An approach of this kind that is adapted to doubly-periodic systems uses periodic Green's functions [36]. However, the usual expressions for doubly-periodic Green's functions contain slowly convergent series [36], and have to be subjected to accelerated transformations [36][37][38] to make them useful in calculations. We have therefore decided not to pursue it in this work.

VII. CONCLUSIONS
We have derived the Rayleigh equation for the amplitude of the scattered field when a scalar plane wave is incident on a two-dimensional rough surface on which the Neumann or Dirichlet boundary condition is imposed. From this equation we have obtained the equation for the amplitudes of the diffracted Bragg beams, when the rough surface is a doubly-periodic one. This equation has been solved by a rigorous numerical approach, and from the solution the dependence of the diffraction efficiencies of several of the of the lowest-order diffracted beams on the polar and azimuthal angles of incidence has been determined. These dependencies display a rich structure of peaks and dips as functions of the polar angle of incidence for a specified azimuthal angle of incidence. These features occur at the angles at which a diffracted beam starts or ceases to propagate. Hence they are the analogues for a doubly-periodic grating of the anomalies that were first observed by Wood [4] in the diffraction of light by a classical metal grating, and were subsequently explained by Lord Rayleigh [5]. They are now called Rayleigh anomalies. These anomalies are observed in the diffraction of a scalar wave from both a Neumann and a Dirichlet surface. In the case of diffraction from a Neumann surface an additional anomaly, a dip, is observed in the reflectivity at angles of incidence for which no Rayleigh anomaly is predicted to occur. No such anomaly is presented in diffraction from doubly-periodic Dirichlet surfaces. A doubly-periodic Neumann surface supports a surface wave, while a doubly-periodic Dirichlet surface does not. From this and responses of the dip to changes of the surface profile function of the Neumann surface, it is conjectured that it is a Wood anomaly that was first reported by Wood in Ref. 4, and was subsequently explained by Fano [6] as due to the excitation of the surface electromagnetic wave supported by the grating by the incident light through the periodic modulation of the surface. This conjecture can only be verified when the branches of the dispersion curve of the surface wave on the Neumann bigrating in the radiative region of the (k , ω) plane have been determined. That will be the subject of a separate work. It should be noted that neither the Neumann nor the Dirichlet surface supports a surface wave when it is planar. Finally, by comparing results obtained from solutions of the Rayleigh equation with results obtained by a rigorous Green's function-based numerical approach [29], we have validated the use of the Rayleigh equation in the calculations reported here.