Numerical solutions of the Rayleigh equations for the scattering of light from a two-dimensional randomly rough perfectly conducting surface

We present rigorous, nonperturbative, purely numerical solutions of the Rayleigh equations for the scattering of p - and s -polarized light from a two-dimensional randomly rough perfectly conducting surface. The solutions are used to calculate the reflectivity of the surface, the mean differential reflection coefficients, and the full angular distribution of the intensity of the scattered field. These results are compared with previously published rigorous numerical solutions of the Stratton – Chu equations, and very good agreement is found. © 2014 Optical Society of America OCIS codes: (240.0240) Optics at surfaces; (240.5770) Roughness; (240.6680) Surface plasmons; (290.5880) Scattering, rough surfaces.


INTRODUCTION
A randomly rough perfectly conducting surface can serve as a simpler testing ground for new computational approaches to rough surface scattering calculations than a randomly rough dielectric or metallic surface. At the same time it can serve as a simple model for a randomly rough metallic surface in the infrared and longer wavelength regions of the electromagnetic spectrum.
Although there have been several numerical calculations of the scattering of light from a two-dimensional randomly rough perfectly conducting surface by one approximate approach or another [1][2][3][4][5][6], there have been few exact solutions of the integral equations by numerical methods [7][8][9][10][11]. This is due largely to the fact that such calculations are still computationally intensive.
In the past, nonperturbative purely numerical solutions of the one-dimensional reduced Rayleigh equations (RREs) for various one-dimensional, randomly rough surface geometries [12,13] have proven rather useful for the study and prediction of several rough surface scattering phenomena [14][15][16][17].
In this paper we present rigorous, nonperturbative, purely numerical solutions of the Rayleigh equations for the scattering of p-and s-polarized light from a two-dimensional randomly rough perfectly conducting surface. The solutions are used to calculate the reflectivity of the surface for incident light of both polarizations, the mean differential reflection coefficients, and the full angular distribution of the intensity of the scattered field.
One motivation for our undertaking of these calculations is that the integral equations that have to be solved are simpler than those that are obtained from the Stratton-Chu formula [18] for the magnetic field in the vacuum region, which underlie the calculations carried out in [7][8][9][10]. The rigorous numerical solution of RREs in recent work by the present authors [19][20][21][22] has been shown to be an effective approach to the calculation of various properties of light scattered from two-dimensional randomly rough penetrable surfaces. We anticipated that this would also be the case in the scattering of light from a perfectly conducting surface. A second motivation for this work is that it provides an opportunity to assess the degree of roughness of a two-dimensional randomly rough perfectly conducting surface up to which numerical solutions of the Rayleigh equations can produce reliable results. This is done by comparing the results of Rayleigh-based calculations to the results of rigorous calculations such as those of [10]. Additionally, the degree to which the results conserve energy is investigated.

RAYLEIGH EQUATIONS
The physical system we consider in this work consists of vacuum in the region x 3 > ζx ∥ , where x ∥ x 1 ; x 2 ; 0, and a perfect conductor in the region x 3 < ζx ∥ (Fig. 1). The surface profile function ζx ∥ is assumed to be a single-valued function of x ∥ that is at least once differentiable with respect to x 1 and x 2 and constitutes a stationary, zero-mean, isotropic, Gaussian random process defined by Here, δ hζ 2 x ∥ i 1∕2 is the root-mean-square (rms) height of the surface, and Wjx ∥ − x 0 ∥ j is the normalized surface height autocorrelation function of the surface profile. The correlation functions to be considered in this paper will be of Gaussian form, given by where a is the correlation length. The angle brackets, here and in all that follows, denote an average over an ensemble of realizations of the surface profile function. Each realization of the surface profile function with these properties is generated numerically by a two-dimensional version of the Fourier filtering method [23,24].
The surface x 3 ζx ∥ is illuminated from the vacuum by an electromagnetic field in the form of a plane wave of angular frequency ω. The total electric field in the vacuum region is the sum of the incident field and the scattered field, where In these equations, k ∥ k 1 ; k 2 ; 0, c is the speed of light in vacuum, α 0 k ∥ ω∕c 2 − k 2 ∥ 1∕2 , with Re α 0 k ∥ > 0 and Im α 0 k ∥ > 0, and the amplitudes E 0 k ∥ and Aq ∥ are given by and The subscripts p and s denote the p-polarized and s-polarized components of the incident and scattered fields with respect to the planes of incidence and scattering, respectively. A caret over a vector denotes that it is a unit vector. The boundary conditions at the surface x 3 ζx ∥ require the vanishing of the tangential component of the electric field, where with ζ μ x ∥ ∂ζx ∥ ∕∂x μ (μ 1, 2), is a vector normal to the surface x 3 ζx ∥ at each point of it, directed into the vacuum. Equation (7) is a vector equation, which consists of the three equations Because n · fn × Exjω inc Exjω sc g 0, only two of the three Eqs. (9) are independent. We choose them to be Eqs. (9a) and (9b). Written out explicitly these two equations become It is convenient at this point to introduce the Fourier representation so that (μ 1, 2) For calculating IγjQ ∥ we need the inverse to Eq. (11) In terms of this function, Eqs. (10) take the forms Upon substituting these relations into Eqs. (15) we find that the matrix of scattering coefficients Rq ∥ jk ∥ satisfies the matrix integral equation where and

REFLECTIVITY AND THE MEAN DIFFERENTIAL REFLECTION COEFFICIENT
The scattering amplitudes fR αβ q ∥ jk ∥ g play a central role in the present theory because the mean differential reflection coefficient can be expressed in terms of these amplitudes. The differential reflection coefficient (∂R∕∂Ω s ) is defined such that (∂R∕∂Ω s ) dΩ s is the fraction of the total time-averaged flux incident on the surface that is scattered into the element of solid angle dΩ s about the scattering direction θ s ; ϕ s . As we are studying the scattering of light from a randomly rough surface, it is the average of this function over the ensemble of realizations of the surface profile function that we need to calculate. The contribution to the mean differential reflection coefficient when incident light of polarization β and wave vector k ∥ is scattered into light of polarization α and wave vector q ∥ , denoted by h∂R αβ ∕∂Ω s i, is Here S is the area of the plane x 3 0 covered by the rough surface. The parallel components of the wave vectors k ∥ k 1 ; k 2 ; 0 and q ∥ q 1 ; q 2 ; 0 are defined in terms of the polar and azimuthal angles of incidence θ 0 ; ϕ 0 and scattering θ s ; ϕ s , respectively, by k ∥ ω∕c sin θ 0 cos ϕ 0 ; sin ϕ 0 ; 0 and q ∥ ω∕c sin θ s cos ϕ s ; sin ϕ s ; 0.
If we write the scattering amplitude R αβ q ∥ jk ∥ as the sum of its mean value and its fluctuation about this mean value, each of these two terms makes its own contribution to the mean differential reflection coefficient (20), which we write now in the form where the first term on the right-hand side is the contribution from the light scattered coherently (specularly), while the second is the contribution from the light scattered incoherently (diffusely). These contributions are given by The mean value of the scattering amplitude R αβ q ∥ jk ∥ can be written in the form The factor 2π 2 δq ∥ − k ∥ in this expression is due to the assumed stationarity of the surface profile function, the Kronecker symbol δ αβ is due to the conservation of angular momentum in the scattering process, and the result that the factor R α k ∥ depends on the wave vector k ∥ only through its magnitude is due to the restoration of isotropy by the averaging process. When Eq. (25) is substituted into Eq. (23), the latter becomes In obtaining this result we have used the relations and k ∥ ω∕c sin θ 0 ; q ∥ ω∕c sin θ s : Upon integrating Eq. (26) over all solid angles in the upper hemisphere (dΩ s sin θ s dθ s dϕ s ) we obtain the reflectivity R α θ 0 of the surface for incident light of α polarization, From Eqs. (25) and (27) we find that Equations (29)-(31) are the basis for the calculation of the reflectivity of the randomly rough surface. Equation (24) is the basis for the calculation of the contribution to the mean differential reflection coefficient from the light scattered incoherently.

ENERGY CONSERVATION AND ITS CONSEQUENCES
In scattering from a perfectly conducting surface the total time-averaged scattered flux must equal the total timeaveraged incident flux. We recall the definition of the mean differential reflection coefficient h∂R αβ ∕∂Ω s i, namely that h∂R αβ ∕∂Ω s idΩ s is the fraction of the incident flux of β polarization that is scattered into light of α polarization in the element of solid angle dΩ s about the scattering direction θ s ; ϕ s . Consequently, the energy conservation condition, or unitarity [Uθ 0 ; ϕ 0 ], can be expressed in terms of the mean differential reflection coefficient as (β p; s) where the integral is to be taken over the upper hemisphere. Equation (32) is useful in assessing the quality of the simulation results, in particular in making certain that the discretization is sufficient, and that the integration interval in q ∥ in Eq. (17) is large enough. However, we emphasize that Eq. (32) is only a necessary condition, and its satisfaction does not guarantee that the simulation results are correct. In the following section, we will also make use of the contribution to U β θ 0 ; ϕ 0 that comes from light that has been scattered incoherently by the surface. This contribution is defined as (β p; s) and it signifies the fraction of the incident power flux that was scattered incoherently by the surface into any scattering direction.

SOLUTION OF THE RAYLEIGH EQUATIONS
Equations (17)- (19) were solved by the method described in detail in [25]. A realization of the surface profile function was generated on a grid of N x × N x points within a square region of the x 1 x 2 plane of edge L. In evaluating the q ∥ integral in Eqs. (17), the infinite limits of integration were replaced by finite limits jq ∥ j < Q∕2, and the integration was carried out by a two-dimensional version of the extended midpoint rule [26, p. 161] applied to a grid in the q 1 q 2 plane that is determined by the Nyquist sampling theorem [26, p. 605] and the properties of the discrete Fourier transform [25]. The function IγjQ ∥ was evaluated by expanding the integrand in Eq. (13) in powers of ζx ∥ and calculating the Fourier transform of ζ n x ∥ by the fast Fourier transform [25]. The resulting matrix equations were solved by LU factorization. These calculations were carried out for a large number N p of realizations of the surface profile function ζx ∥ for an incident field of p or s polarization. For each realization the scattering amplitude R αβ q ∥ jk ∥ and its squared modulus jR αβ q ∥ jk ∥ j 2 were obtained. An arithmetic average of the N p results for these quantities yielded the mean values hR αβ q ∥ jk ∥ i and hjR αβ q ∥ jk ∥ j 2 i entering Eqs. (31) and (24) for the reflectivity and mean differential reflection coefficient, respectively.

RESULTS
The purpose of this section is to present some examples of numerical results that can be obtained by the formalism just introduced and to discuss their quality and how they compare to results obtained by other methods.
Calculations were carried out for two-dimensional randomly rough perfectly conducting surfaces defined by an isotropic Gaussian height distribution of rms height δ λ∕10 and an isotropic Gaussian correlation function of lateral correlation length a 3λ∕5. The incident light consisted of a p-or s-polarized plane wave of wavelength λ (in vacuum) and well-defined angles of incidence θ 0 ; ϕ 0 . In these calculations, it was assumed that the surfaces covered a region of the x 1 x 2 plane of area L 2 15λ × 15λ. Realizations of the surface profile function were generated [23,24] in this region on a grid of N x × N x 321 × 321 points. With these spatial parameters, the corresponding momentum space parameters used in the simulations were Δq 2π∕L for the discretization intervals in momentum space, and the largest momentum value that can be resolved is given by Q Δq⌊N x ∕2⌋, where ⌊ · ⌋ denotes the floor function [26]. In the simulations, we used Q 10.7ω∕c, and for the expansions of the integrands of the IγjQ ∥ integrals, we used the first n 20 terms. With these parameters, solving the Rayleigh equation for one surface realization required about 12 GB of memory and took approximately 18 min on a machine with two six-core AMD Opteron 2.1 GHz processors. The majority of this time was spent solving the linear set of equations (by the LU-factorization method), which was done by the optimized and highperformance numerical parallel library ScaLAPACK [27] and Intel MKL [28]. Figure 2 presents, for angles of incidence θ 0 ; ϕ 0 0°; 45°, the angular distributions of the mean differential reflection coefficients that result from the incoherent component of the scattered light. Results obtained for 5000 realizations of the surface profile function were averaged to obtain these figures. The incident light was p polarized when produc-  Fig. 2 are observed for the copolarized or cross-polarized scattering configurations where the scattered intensity distributions show "dipole-like" patterns. In particular, in copolarized scattering [Figs. 2(b) and 2(f)], the maximum of the angular intensity distribution is oriented along the plane of incidence [29], while in the cross-polarized scattering configuration [Figs. 2(c) and 2(e)], the maximum intensity is found in a direction that is perpendicular to this plane. The reason for this behavior can be explained by noting that the main contribution to the scattered intensity, for the roughness parameters assumed here, comes from single scattering processes that for copolarization (cross polarization) is zero for the out-of-plane (in-plane) scattering [10,11].
When the polarization of the scattered light is not recorded [Figs. 2(a) and 2(d)], in contrast to when it is, the resulting intensity distributions have much less dependence on the azimuthal angle of scattering ϕ s . However, there is some dependence on ϕ s , and a close inspection of Figs. 2(a) and 2(d) reveals that the equi-intensity lines have elongated shapes (ellipse-like) with the long axis being in the direction of the incident electric field. It should be noted that many of the characteristics of Fig. 2 are similar to what has been recently reported in numerical studies of related scattering geometries using either similar or different numerical approaches for their study. Here we mention light scattering from potentially strongly rough perfectly conducting surfaces [10,19], light scattering from penetrable rough surfaces of in principle any degree of roughness [11,19,24], and, finally, weakly rough surfaces that are part of a clean surface or film geometry [20,21,25].
An example of the angular dependence of the incoherent component of the mean differential reflection coefficients for non-normal incidence and the same scattering system is given in Fig. 3, where it was assumed that θ 0 ; ϕ 0 22°; 45°. In this figure, the meaning of the various subfigures is the same as those of Fig. 2. The first observation to be made from Fig. 3 is that the scattered intensity is now predominantly in forward scattering directions, ϕ s − ϕ 0 ∈ −90°; 90°, and only a small fraction of the incident light is scattered in other directions. Also for non-normal incidence, we find that the maxima of the copolarized intensity distributions are oriented along the plane of incidence. For cross polarization, on the other hand, the "wings" of the "dipole-like" intensity patterns  1) and (2), with δ λ∕10 and a 3λ∕5, where λ is the wavelength of the incident radiation (in vacuum). An ensemble of 5000 surface realizations was used to produce these results. Both edges of the surfaces were L 15λ, and the surface was discretized at 321 × 321 points.
are pushed toward the forward scattering directions, but still show local minima in the plane of incidence. We note that as the local slopes of the surfaces increase, it is expected that the fraction of incident intensity that is scattered into the forward directions (ϕ s − ϕ 0 ∈ −90°; 90°) will decrease, and the scattering into the backscattering directions (ϕ s − ϕ 0 ∉ −90°; 90°) will increase, while, at the same time, the wings of the "dipole-like" angular intensity patterns observed for the crosspolarized scattering will gradually shift toward the backscattering directions [10,16,24]. An additional observation worth making from Figs. 3(a) and 3(d) is that the positions of the maxima of the incoherent component of the mean differential reflection coefficients are not coinciding with the specular directions (indicated by white dots in Figs. 2 and 3). In particular, for the roughness parameters assumed here, it is found that the polar angles of scattering for the these maxima, θ s , are smaller than those of the specular direction; that is, θ s < θ 0 . In neither of the simulation results presented in Figs. 2 and 3 was the enhanced backscattering phenomenon [16,22,30,31] observed; see also Figs. 6 and 7. However, for the roughness parameters assumed in producing the results shown in these figures, this is not to be expected [16,30,31].
It should be noted that since the statistical properties of the surfaces are isotropic, the quantity h∂R αβ ∕∂Ω s i will only have a trivial dependence on the azimuthal angles of incidence ϕ 0 , and, therefore, the value ϕ 0 45°assumed in the simulations is somewhat arbitrary. By trivial, we here mean features that can be compensated for by a simple rotation of the resulting angular intensity distributions.
In Fig. 4 we present the reflectivity, R α θ 0 [see Eq. (30)], as a function of the polar angle of incidence θ 0 of the α-polarized incident light. The roughness parameters assumed in producing the results of Fig. 4 are identical to those assumed in producing the results of Fig. 2. Figure 4 shows that in the case of s-polarized incident light, the reflectivity R s θ 0 is a monotonically increasing function of the polar angle of incidence θ 0 . On the other hand, when the incident light is p polarized, the corresponding reflectivity satisfies R p θ 0 ≤ R s θ 0 ; for the smallest polar angles of incidence (θ 0 ≤ 30°), the two reflectivities are equal, while for larger angles of incidence, R p θ 0 is the smaller of the two. The results in Fig. 4 show that R p θ 0 goes through a maximum around θ 0 ∼ 65°, and a minimum around θ 0 ∼ 88°. The decrease of R p θ 0 with angle of incidence we believe is analogous to the pseudo-Brewster (or quasi-Brewster) effect, which has been observed in the reflectivity of lossy metal surfaces in p polarization [32,33].
Since we are dealing with perfectly conducting surfaces, all electromagnetic power flux incident on a surface will be reflected. As explained in Section 4, the conservation of energy implies U β θ 0 ; ϕ 0 ≡ 1 [see Eq. (32)]. This condition should in principle be satisfied for any angles of incidence and any state of polarization of the incident beam, including the linear polarizations β p; s. By numerically performing the Ω s integrals (or corresponding integrals over q ∥ ) of the mean differential reflection coefficients obtained numerically, U β θ 0 ; ϕ 0 can be calculated [Eq. (32)].
In this way, the numbers given in Table 1 were obtained from the results presented in Figs. 2 and 3. The results of Table 1 demonstrate that energy conservation is satisfied to within a fraction of a percent for all the simulation data reported. This testifies to the correctness, self-consistency, and quality of the simulation results.
The energy conservation condition [Eq. (32)] is also useful for probing the range of roughness parameters for which the Rayleigh equations, or more correctly, our implementation and solutions of them, are adequate for a consistent   4. Reflectivity, R α θ 0 , as defined in Eq. (30), as a function of polar angle of incidence θ 0 and for given polarization α of the incident light. The roughness parameters assumed in obtaining these results are identical to those assumed in producing the results of Fig. 2, but the Rayleigh equation has here been solved for a set of 13 slightly different surface lengths, in order to achieve the desired resolution in the angles close to 90°. description of the scattering processes involved. Since the Rayleigh equations are derived under the assumption of the Rayleigh hypothesis [34], this may indirectly also be a test of the validity of this hypothesis. It is outside the scope of this paper to systematically study these issues, and readers interested in this topic are referred elsewhere [35,36]. However, for the purpose of this paper, in Fig. 5 we present the dependence of Uθ 0 ; ϕ 0 U p θ 0 ; ϕ 0 U s θ 0 ; ϕ 0 ∕2 on the lateral correlation length a, or, equivalently, the rms slope s of the surface, for given angles of incidence θ 0 ; ϕ 0 0°; 45°, a rms roughness of δ λ∕10, and with the remaining parameters being similar to those used to produce the results of Fig. 2.
When both the surface height distribution and the lateral correlation function are Gaussian, as assumed here, the rms slope is given by s 2 p δ∕a [16,23]. The reason that we present Uθ 0 ; ϕ 0 , and not U p and U s separately, is only to improve the statistics. The inset in Fig. 5 indicates that the energy is conserved to within 1%, or better, for a rms slope of s ≈ 2∕3, or smaller (at least for δ λ∕10 assumed here). Moreover, from the same figure energy is conserved to within 10 −5 , or better, for s ≤ 0.28 (a ≥ λ∕5), again assuming δ λ∕10. We note that in [35] it was shown that in the δ; a-parameter region, for which s ≤ 0.28, energy was conserved with an error smaller than 10 −4 , or better, if δ was not too large. However, as δ∕λ became larger than about 0.15, the condition U ≡ 1 seemed to not be satisfied, i.e., had an error greater than 10 −2 , for any correlation length. This lack of energy conservation was interpreted in [35] not as an indication that the Rayleigh equations themselves are not valid, but instead it was believed to be caused by inaccuracies in the evaluation of the IγjQ ∥ integrals by the expansion method used to calculate them (see [35] for details).
Finally, we will present a comparison of mean differential reflection coefficients either obtained by solving the Rayleigh equations by the methods presented in this paper, or obtained from the Green's function surface integral approach of [10] that is based on the Stratton-Chu formula [18]. An attractive practical advantage of the former approach over the latter is that when we apply the former, it requires about 25 times less memory than the latter [25]. While the Rayleigh equation approach is expected to be fairly accurate only for surfaces of small local slopes since it assumes the Rayleigh hypothesis [17,[34][35][36], the Green's function surface integral approach is rigorous and should therefore in principle be valid for any surface roughness. In Figs. 6 and 7, we compare results obtained when these two approaches are applied to the same scattering system (whose properties are defined in the caption of Fig. 2). In particular, what is compared is the polar angular Table 1. Values for U β θ 0 ; ϕ 0 and U incoh β θ 0 ; ϕ 0 as Defined by Eqs. (32) and (33) 5. Unitarity, Uθ 0 ; ϕ 0 U p θ 0 ; ϕ 0 U s θ 0 ; ϕ 0 ∕2, as a function of the correlation length a, or, equivalently, the rms slope s of the surface, for the angles of incidence θ 0 ; ϕ 0 0°; 45°. Since the surfaces used to produce these results had both Gaussian height distributions and Gaussian correlation functions, the rms slope of the surfaces is given by s 2 p δ∕a [16,23]. The inset depicts a semilogarithmic scale jU − 1j as a function of correlation length in order to better identify the deviation of U from unity. In the simulations performed to produce the data presented in this figure, we kept the rms height of the surface fixed at δ λ∕10 while varying its transverse correlation length a. For each set of roughness parameters the results were averaged over an ensemble consisting of 20 surface realizations, which was sufficient to produce converged results. The remaining parameters used in obtaining these results are given in the caption of Fig. 2.  Fig. 6. Comparison of the results for the incoherent component of the mean differential reflection coefficient as obtained by numerically solving the Rayleigh equations (lines) and those resulting from applying the (rigorous) Green's function surface integral method [10] (open symbols) to the same scattering system. Both the co-and crosspolarized components of the mean differential reflection coefficients coming from the light that has been scattered incoherently by the surface are presented, but only their in-plane (left panel) and out-of-plane (right panel) dependencies are shown. The incident light was p polarized, and the remaining roughness parameters were identical to those assumed in producing the results of Fig. 2. The curves corresponding to the Rayleigh equation approach were actually obtained by making appropriate cuts through the corresponding angular intensity distributions of Fig. 2. When applying the Green's function surface integral approach, the simulation parameters were identical to those used in the Rayleigh approach (see caption of Fig. 2) with the following exceptions: the spatial discretization interval was Δx 2λ∕17, the incident beam had a finite size of full width w 4λ, and the number of surface realizations used was N p 10 000.
dependence of the co-and cross-polarized incoherent components of the mean differential reflection coefficients in the plane of incidence (in-plane) or in a plane perpendicular to it (out-of-plane). Figures 6 and 7 show excellent agreement between the results obtained by the two approaches. We have also verified that similar agreement exists for the other parameters of the incident beam used in producing Figs. 2 and 3, but these results are not shown here.

CONCLUSIONS
In conclusion, we have presented the Rayleigh equations for the scattering of light from a two-dimensional, randomly rough, perfectly conducting surface, and solved them nonperturbatively and purely numerically to obtain the scattered light. In particular, for a selection of angles of incidence, we reported the full angular distribution for the co-and cross-polarized incoherent components of the mean differential reflection coefficients of the scattered light for either p or s polarization of the incident plane waves. The fraction of the incident light that was scattered incoherently and coherently by the randomly rough surface was calculated and reported. The simulation results that were obtained by a purely numerical solution of the Rayleigh equations were, for the same surface parameters, found to give very good agreement with the results obtained by the more computationally costly, but fully rigorous, Green's function surface integral method. The quality of the simulation results was quantified by investigating energy conservation (unitarity). For the main results reported, energy conservation was found to be satisfied with an error smaller than 10 −4 , or better. We also investigated the energy conservation for surfaces of increasing rms slopes, and it was found that the error was smaller than 0.1% (1%) for slopes smaller than about 0.47 (0.71) (at least for the roughness parameters that we assumed in the simulations). In addition to indicating a formally correct numerical solution of the Rayleigh equations, such results also testify to the usefulness of purely numerical, and therefore nonperturbative, solutions of such equations.  Fig. 7. Same as Fig. 6, but now the incident light is s polarized and the polar angle of incidence is θ 0 22°.