Numerical studies of the scattering of light from a two-dimensional randomly rough interface between two dielectric media

The scattering of polarized light incident from one dielectric medium on its two-dimensional randomly rough interface with a second dielectric medium is studied. A reduced Rayleigh equation for the scattering amplitudes is derived for the case where p- or s-polarized light is incident on this interface, with no assumptions being made regarding the dielectric functions of the media. Rigorous, purely numerical, nonperturbative solutions of this equation are obtained. They are used to calculate the reflectivity and reflectance of the interface, the mean differential reflection coefficient, and the full angular distribution of the intensity of the scattered light. These results are obtained for both the case where the medium of incidence is the optically less dense medium, and in the case where it is the optically more dense medium. Optical analogues of the Yoneda peaks observed in the scattering of x-rays from metal surfaces are present in the results obtained in the latter case. Brewster scattering angles for diffuse scattering are investigated, reminiscent of the Brewster angle for flat-interface reflection, but strongly dependent on the angle of incidence. When the contribution from the transmitted field is added to that from the scattered field it is found that the results of these calculations satisfy unitarity with an error smaller than $10^{-4}$.


I. INTRODUCTION
In the great majority of the theoretical studies of the scattering of light from a two-dimensional randomly rough surface of a dielectric medium, the medium of incidence has been vacuum. Recent reviews of such studies can be found in Refs. 1 and 2. As a result of this restriction, effects associated with total internal reflection, which requires that the medium of incidence be optically more dense than the scattering medium, were not considered in these studies. There have been exceptions to this general practice, however.
By the use of the stochastic functional approach [3], Kawanishi et al. [4] studied the coherent and incoherent scattering of an electromagnetic wave from a two-dimensional randomly rough interface separating two different dielectric media. The light could be incident on the interface from either medium. The theoretical approach used in this work [4] is perturbative in nature, and applicable only to weakly rough interfaces. Nevertheless its use yielded interesting results, including the presence of Yoneda peaks in the angular dependence of the intensity of the light scattered back into the medium of incidence when the latter was the optically more dense medium. These are sharp, asymmetric peaks occurring at the critical angle for total internal reflection for a fixed angle of incidence for both p-and s-polarization of the incident light. These peaks were first observed experimentally in the scattering of xrays incident from air on a metal surface [5] and have subsequently been studied theoretically in the context of the scattering of x-rays [6][7][8] and neutrons [7] from rough surfaces.
In Ref. 4, Kawanishi et al. also observed angles of zero scattering intensity, to first order in their approach, in the distributions of the intensity of the incoherently scattered light, when the incident light was p-polarized. Due to their resemblance to the Brewster angle in the reflectivity from a flat interface, they dubbed these angles the "Brewster scattering angles". These were observed, in both reflection and transmission, for light incident from either medium.
Both the Yoneda peaks and the Brewster scattering angles seem to have had their first appearance in optics in the paper by Kawanishi et al. [4]. They have yet to be observed experimentally in this context. It should be mentioned that in an earlier numerical investigation of light scattering from one-dimensional dielectric rough surfaces, Nieto-Vesperinas and Sánchez-Gil [9] observed "sidelobes" in the angular intensity distributions. However, these authors did not associate these features with the Yoneda peak phenomenon, even though we believe doing so would have been correct.
In a subsequent paper Soubret et al. [10] derived a reduced Rayleigh equation for the scattering amplitudes when an electromagnetic wave is incident from one dielectric medium on its two-dimensional randomly rough interface with a second dielectric medium. The solution of this equation was obtained in the form of expansions of the scattering amplitudes in powers of the surface profile function through terms of third order. However, in obtaining the numerical results presented in this paper [10], the medium of incidence was assumed to be vacuum.
In this paper we present a study of this problem free from some of the restrictive assumptions and approximations present in the earlier studies of scattering of polarized light from two-dimensional randomly rough dielectric surfaces. We first derive a reduced Rayleigh equation for the scattering amplitudes when p-or s-polarized light is incident from a dielectric medium whose dielectric constant is ε 1 , on its two-dimensional randomly rough interface with a dielectric medium whose dielectric constant is ε 2 . The dielectric constant ε 1 can be smaller or larger than ε 2 . This equation is then solved by a rigorous, purely numerical, nonperturbative approach. The scattering amplitudes obtained in this way are then used to calculate the reflectivity and reflectance of the interface as a function of the angle of incidence, and also the effect of surface roughness on the contribution to the mean differential reflection coefficient from the light scattered incoherently (diffusely), and the full angular dependence of the intensity of the incoherently scattered light. It is hoped that the presentation of these results will stimulate and motivate experimental studies of such scattering systems.

II. THE SCATTERING SYSTEM
The system we study in this paper consists of a dielectric medium (medium 1), whose dielectric constant is ε 1 , in the region x 3 > ζ(x ), and a dielectric medium (medium 2), whose dielectric constant is ε 2 , in the region x 3 < ζ(x ) [ Fig. 1]. Here x = (x 1 , x 2 , 0) is an arbitrary vector in the plane x 3 = 0, and we assume that both ε 1 and ε 2 are real and positive. The surface profile function ζ(x ) is assumed to be a single-valued function of x that is differentiable with respect to x 1 and x 2 , and constitutes a stationary, zero-mean, isotropic, Gaussian random process defined by where W (x ) is the normalized surface height autocorrelation function, with the property that W (0) = 1. The angle brackets here and in all that follows denote an average over the ensemble of realizations of the surface profile function.
x 1 x 2 1. A sketch of the scattering geometry assumed in this work. The figure also shows the coordinate system used, angles of incidence(θ0, φ0) and scattering (θs, φs), and the corresponding lateral wave vectors k and q , respectively.
The root-mean-square height of the surface is given by The power spectrum of the surface roughness g(k ) is defined by For W (x ) we assume the Gaussian function W (x ) = exp −x 2 /a 2 , where the characteristic length a is the transverse correlation length of the surface roughness. The corresponding power spectrum is given by

III. THE REDUCED RAYLEIGH EQUATION
The interface x 3 = ζ(x ) is illuminated from the region x 3 > ζ(x ) (medium 1) by an electromagnetic wave of frequency ω. The total electric field in this region is the sum of an incoming incident field and an outgoing scattered field, while the electric field in the region x 3 < ζ(x ) is an outgoing transmitted field, In writing these equations we have introduced the functions At this point it is convenient to introduce the representation On differentiating both sides of Eq. (15) with respect to x j (j = 1, 2) we obtain the result Finally, to be able to evaluate the function I(γ|Q ) we need the inverse of Eq. (15), namely On combining Eqs. (14)- (17) with Eqs. (7) and (12) we obtain the results When the results given by Eq. (19) are substituted into Eq. (13), the latter becomes Thus, the amplitude of the transmitted field B(q ) has been eliminated from the problem and we have obtained an equation for the scattering amplitude A(q ) alone. To transform Eq. (20) into a more useful form we first introduce three mutually perpendicular unit vectors: In terms of these vectors Eq. (20) becomes: We now write the vectors E 0 (k ) and A(q ) in the forms and In these expressions E 0p (k ) and E 0s (k ) are the amplitudes of the p-and s-polarized components of the incident field with respect to the plane of incidence, defined by the vectorsk andx 3 . Similarly, A p (q ) and A s (q ) are the amplitudes of the p-and s-polarized components of the scattered field with respect to the plane of scattering defined by the vectorsq andx 3 . Equation (22) is a vector equation: it is a set of three coupled equations. However, there are only two unknowns, namely A p (q ) and A s (q ). Consequently, one of these equations is redundant. To obtain A p (q ) and A s (q ) in terms of E 0p (k ) and E 0s (k ) we proceed as follows. We take the scalar product of Eq. (22) with each of the three unit vectors given by Eq. (21) in turn. The results are: Equations (25b) and (25c) are the two equations we seek. With the use of Eqs. (21) and (23)-(24), Eqs. (25b) and (25c) can be rewritten in the form (α = p, s, β = p, s) On combining Eqs. (25b) and (25c) with Eq. (26) we find that the scattering amplitudes {R αβ (q |k )} are the solutions of the equation where and Equation (27) is the reduced Rayleigh equation for the scattering amplitudes.

IV. THE MEAN DIFFERENTIAL REFLECTION COEFFICIENT
From the knowledge of the scattering amplitudes the mean differential reflection coefficient, the reflectivity, and the reflectance can be calculated. The differential reflection coefficient ∂R/∂Ω s is defined such that (∂R/∂Ω s )dΩ s is the fraction of the total time-average flux incident on the interface that is scattered into the element of solid angle dΩ s about the scattering direction defined by the polar and azimuthal scattering angles (θ s , φ s ). To obtain the mean differential reflection coefficient we first note that the magnitude of the total time-averaged flux incident on the interface is given by In this result S is the area of the x 1 x 2 plane covered by the randomly rough surface. The minus sign on the right-hand side of the first equation compensates for the fact that the 3-component of the incident flux is negative, and we have used the fact that α 1 (k ) is real, so that Q 0 (k ) is real, and E * 0 (k ) · Q 0 (k ) = 0. In a similar fashion we note that the total time-averaged scattered flux is given by The integral in the second term is purely imaginary. Thus we have The wave vectors k and q can be expressed in terms of the polar and azimuthal angles of incidence (θ 0 , φ 0 ) and scattering (θ s , φ s ), respectively, by From these results it follows that where dΩ s = sin θ s dθ s dφ s . The total time-averaged scattered flux therefore becomes Similarly, the total time averaged incident flux, Eq. (29), becomes Thus by definition, the differential reflection coefficient is given by From this result and Eq. (26) we find that the contribution to the differential reflection coefficient when an incident plane wave of polarization β, the projection of whose wave vector on the mean scattering plane is k , is reflected into a plane wave of polarization α, the projection of whose wave vector on the mean scattering plane is q , is given by As we are dealing with scattering from a randomly rough interface, it is the average of this function over the ensemble of realizations of the surface profile function that we need to calculate. This is the mean differential reflection coefficient, which is defined by If we write the scattering amplitude R αβ (q |k ) as the sum of its mean value and the fluctuation from this mean, then each of these two terms contributes separately to the mean differential reflection coefficient, where and The former contribution describes the coherent (specular) reflection of the incident field from a randomly rough surface, while the latter contribution describes the incoherent (diffuse) component of the scattered light.

A. Reflectivity and Reflectance
Equation (41) is the starting point for obtaining the reflectivity of the two-dimensional randomly rough interface. We begin with the result that The presence of the delta function is due to the stationarity of the randomly rough surface; the Kronecker symbol δ αβ arises from the conservation of angular momentum in the scattering process; and the result that R α (k ) depends on k only through its magnitude is due to the isotropy of the random roughness.
With the result given by Eq. (43), the expression for ∂R αβ (q |k )/∂Ω s coh given by Eq. (41), becomes where we have used the result in obtaining this expression. We next use the relation together with the relations to obtain The reflectivity, R α (θ 0 ), for light of α polarization is defined by The function R α (k ) is obtained from Eq. (43), with the aid of the result that (2π) 2 δ(0) = S, in the form In addition to the reflectivity (49) that depends only on the co-polarized light reflected coherently by the rough interface, it is also of interest to introduce the reflectance for β-polarized light defined as where In short, the reflectance measures the fraction of the power flux incident on the rough surface that was reflected by it, taking both specularly and diffusely reflected light into account: In view of Eq. (40), the reflectance is the sum of a contribution from light that has been reflected coherently and a contribution from light that has been reflected incoherently by the rough interface, R β (θ 0 ) = R β (θ 0 ) coh +R β (θ 0 ) incoh , and both co-and cross-polarized reflected light contribute. Since cross-polarized coherently reflected light is not allowed [see Eq. (43)], the coherent contribution to the reflectance for β-polarized light equals the reflectivity for β-polarized light; R β (θ 0 ) coh = R β (θ 0 ). Equation (51a) can therefore also be written in the form If the incident light is not purely p-or s-polarized, the reflectance and the reflectivity of the rough surface will have to be calculated on the basis of weighted sums of the expressions in Eqs. (49) and (52), where the weights reflect the fractions of the different polarizations contained within the incident light.

V. NUMERICAL SOLUTION OF THE REDUCED RAYLEIGH EQUATION
The simulation results to be presented in this work were obtained by a nonperturbative numerical solution of the reduced Rayleigh equation (27), which was carried out in the following manner; 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 edges L. This surface profile enters Eq. (27) through the function I(γ|Q ), given by Eq. (17). Utilizing the Taylor expansion detailed in Eq. (17), the Fourier transform of ζ n (x ) was calculated by use of the fast Fourier transform [2], and the Taylor series was truncated at the finite order N T . In evaluating the q integral in Eq. (27), the infinite limits of integration were replaced by finite limits |q | < Q/2, and the integration was carried out by a two-dimensional version of the extended midpoint rule [11, p. 135] applied to a circular subsection of a grid of N q × N q points in the q 1 q 2 -plane, whose size and discretization was determined by the Nyquist sampling theorem [11, p. 494] and the properties of the discrete Fourier transform [2]. In momentum space, these limits lead to discretization intervals of ∆q = 2π/L along the orthogonal axes of the q 1 q 2 -plane, and upper limits on the magnitude of resolved wave vectors are given by Q = ∆q N x /2 , where · denotes the floor function [11, p. 948]. The resulting linear system of equations was solved by LU factorization and back substitution. These calculations were performed simultaneously for incident light of both p-and s-polarization, and they were performed for a large number N p of realizations of the surface profile function ζ(x ). The resulting scattering amplitude R αβ (q |k ) and its squared modulus |R αβ (q |k )| 2 were obtained for each realization. An arithmetic average of the N p results for these quantities yielded the mean values R αβ (q k ) and |R αβ (q |k )| 2 that enter Eqs. (50), (51) and (42) for the reflectivity, reflectance and the mean differential reflection coefficient, respectively. A more detailed description of the numerical method can be found in Ref. 2.

VI. RESULTS AND DISCUSSIONS
The two-dimensional randomly rough dielectric interfaces we study in this work were defined by an isotropic Gaussian height distribution of rms height δ = λ/40, and an isotropic Gaussian correlation function of transverse correlation length a = λ/4 [Eq. (4)]. They covered a square region of the x 1 x 2 -plane of edge L = 25λ, giving an area S = L 2 = 25λ × 25λ. The incident light was assumed to be a p-or s-polarized plane wave of wavelength λ in vacuum. One of the two media in our configuration was assumed to be vacuum with a dielectric constant ε = 1.0, and the other medium was assumed to be a photoresist defined by the dielectric constant ε = 2.64. Since the dielectric constants entering the calculations are independent of the wavelength, all lengths appearing in them can be scaled with respect to λ. The angles of incidence were (θ 0 , φ 0 ), where the azimuthal angle of incidence was set to φ 0 = 0 • , without loss of generality. We remark that this value of φ 0 was chosen since it coincided with one of the two axes of the numerical grid, but that it is, due to the isotropy of the roughness, an arbitrary choice in the sense that results for any other value of φ 0 can be obtained from the results presented through a trivial rotation. The surface profiles were generated by the Fourier filtering method (see Refs. 12 and 13) on a grid of N x × N x = 321 × 321 points. The values used for N x and L correspond to Q = 6.4 ω/c, where Q is the limit in the I(γ|Q )-integrals, Eq. (17), and we used the first N T = 18 terms of the Taylor expansion in the calculation of these integrals.
Investigating the energy conservation of our simulation results can be a useful test of their accuracy. In combining simulation results from the current work with corresponding results obtained for the mean differential transmission coefficient ∂T αβ /∂Ω t through the use of computationally similar methods [14,15], we may add the total reflected and transmitted power for any lossless system. When the reflectance is added to the transmittance for any of the systems investigated in the current work, it is found that the results of these calculations satisfy unitarity [14], a measure of energy conservation, with an error smaller than 10 −4 . This testifies to the accuracy of the approach used, and it is also a good indicator of satisfactory discretization. It should be noted, however, that unitarity is a necessary, but not sufficient, condition for the correctness of the presented results. In a separate investigation [16], unitarity was found to be satisfied to a satisfactory degree for surfaces with a root mean square roughness up to about three times larger than the roughness used in obtaining the results presented in this paper.

A. Normal incidence
In Fig. 2 we display the contribution to the in-plane (q k ) incoherent components of the mean differential reflection coefficient (DRC) as a function of the polar angle of scattering when the random surface is illuminated from the vacuum side at normal incidence by p-and s-polarized light, Fig. 2(a), and when it is illuminated from the dielectric medium side, Fig. 2(b). Notice that the unit vectorsq = q /q andk = k /k are well defined  also for θ s = 0 • and θ 0 = 0 • , respectively, as follows from Eqs. (32) and (47). Only results for in-plane co-polarized scattering are presented, since in-plane cross-polarized scattering is suppressed due to the absence of contribution from single-scattering processes. An ensemble of 4500 realizations of the surface profile function was used to produce the numerical results that this figure is based on. This ensemble size is more than adequate in terms of the interpretation of the results and their features, but we note that a larger ensemble size would have reduced the jaggedness that can be observed in all the (solid line) results presented in this work. From Fig. 2(a) it is observed that the curves corresponding to the two polarizations are featureless, and are nearly identical. In contrast, the curves presented in Fig. 2(b) are rather different for the two polarizations; they display both peaks and dips in p → p scattering, and peaks in s → s scattering. The origins of these features can be understood through small amplitude perturbation theory (SAPT). The contribution to the mean differential reflection coefficients from light scattered incoherently can to the lowest nonzero order in the surface profile function ζ(x ) be expressed as (see the appendix for details):  53), reproduces all the important features found in the mean differential reflection coefficients fairly well, but with a discrepancy in the amplitudes. This discrepancy decreases with decreasing surface roughness (results not shown). For example, similar comparisons for surfaces with an rms-roughness of δ = λ/80, but with the same correlation length a = λ/4, show that the ability of Eq. (53) to reproduce the results based on the RRE is excellent for such weakly rough surfaces. However, for surfaces of rms-roughness δ = λ/20 and the same correlation length, significant discrepancies are observed both in intensity (or amplitude) and angular dependence between the curves obtained on the basis of SAPT and the corresponding curves resulting from a numerical solution of the RRE [Fig. 3]. For instance, from Fig. 3(b) it is observed that the angular dependence of the s → s scattered intensity around the normal scattering direction is not correctly reproduced by SAPT; in this angular interval the numerical simulation results are almost constant and therefore essentially independent of θ s . These results illustrate the importance and necessity to go beyond lowest order SAPT or to do numerical simulations. We therefore stress the point that even if we in the following often turn to SAPT for interpretation of the nonperturbative solutions to the RRE, any conclusion drawn on the basis of Eq. (53) is correct only to the lowest non-zero order in the surface profile function. Results similar to those presented in Figs. 2 and 3 but for scattering out-of-plane [q ·k = 0] are not presented, since, for normal incidence, the results for co-polarized in-plane scattering are the same as the results for cross-polarized out-of-plane scattering when the polarization of the scattered light is the same in the two cases. This symmetry is expected for isotropic surfaces like the ones we are investigating when the lateral momentum of the incident light is zero, supported by the observation that Eq. (53a) evaluated in-plane equals Eq. (53c) evaluated out-of-plane, and correspondingly for Eqs. (53d) and (53b), when k = 0 and θ 0 = 0.
In order to simplify the subsequent discussion we here express d α (q ) and d α (k ) in terms of the polar angles of incidence and scattering using the relations in Eq. (47): We see from Eq. (54) that when ε 2 is greater than ε 1 , both d p (q ) and d s (q ) are real continuous functions of θ s in the interval 0 < |θ s | < π/2, and so therefore are |d p (q )| 2 and |d s (q )| 2 . Hence, no features are introduced into the corresponding mean differential reflection coefficients by these functions. However, when ε 1 is greater than ε 2 , the function (ε 2 /ε 1 ) − sin 2 θ s 1 2 present in d p (q ) and d s (q ) vanishes when |θ s | equals the critical angle θ c = sin −1 ε 2 /ε 1 for total internal reflection for the corresponding flat-surface system, and becomes purely imaginary as |θ s | increases beyond this angle. The functions |d p (q )| 2 and |d s (q )| 2 therefore both have minima at |θ s | = θ c . In the case of s → s in-plane scattering at normal incidence, the minima in the function |d s (q )| 2 for ε 1 > ε 2 lead to sharp peaks at |θ s | = θ c in ∂R ss (q |k )/∂Ω s incoh , as displayed in Fig. 2(b). These same peaks will then also be present for p → s out-of-plane scattering at normal incidence. However, for p → p in-plane scattering, while there are still minima in the function |d p (q )| 2 , we see from Eq. (53a) that ∂R pp (q |k )/∂Ω s incoh , to lowest order in the surface profile function, vanishes when the function vanishes. For normal incidence (k = 0) and in-plane scattering (q k ), we see from this expression that ∂R pp (q |k )/∂Ω s incoh will vanish when α 2 (q ) = 0, which is when q = √ ε 2 ω/c [see Eq. (47)]. This will be the case for θ s < 90 • only when ε 1 > ε 2 , and the expression for ∂R pp (q |k )/∂Ω s incoh in Eq. (53a) will in this case be zero for |θ s | = θ c , explaining the dips shown in Fig. 2(b) for p → p scattering. We note in passing that the out-of-plane distribution of ∂R ps (q |k )/∂Ω s incoh also shows dips at the same polar angles, due to the factor α 2 (q ) in Eq. (53c), but these dips will be present regardless of the angle of incidence. The peaks observed for θ s ≥ θ c in Figs. 2 and 3 for ε 1 > ε 2 are optical analogues of the Yoneda peaks observed by Y. Yoneda in the scattering of x-rays from a metal surface [5] and described as "quasi-anomalous scattering peaks" in the two-dimensional work by Kawanishi et al. [4]. The Yoneda peaks were originally observed as sharp peaks for incidence close to the grazing angle, as the difference in the dielectric constants of the two scattering media is very small at x-ray frequencies. In the following, by Yoneda peaks we will mean well-defined maxima in the angular distribution of the intensity of the scattered light at, or slightly above, the critical polar angle for total internal reflection, when ε 1 > ε 2 .
Although the mathematical origin of the Yoneda peaks is clear from Eqs. (53) and (54), namely, they are associated with the minimum of the functions |d p,s (q )| 2 and |d p,s (k )| 2 , a physical interpretation of them is still under discussion. Thus, Warren and Clarke [17] in a study of the reflection of x-rays from a polished surface (mirror), proposed that these peaks arise when the incident beam at a grazing angle of incidence θ that is greater then the grazing critical angle for total internal reflection θ c is scattered through a small angle β by something just above the mirror surface. The scattered field falls upon the mirror at a grazing angle α, and strong reflection occurs when α < θ c . This reflection is cutoff sharply for α > θ c and less sharply for α < θ c by the rapidly decreasing intensity of small-angle scattering. This produces an asymmetric peak in the intensity of the scattered field, whose maximum occurs at the critical angle. It was suggested that the scatterer could be a projection on an irregular surface.
In a subsequent study of the grazing-angle reflection of x-rays from rough metal surfaces with the use of the distorted-wave Born approximation [18], Vineyard [6] noted from an examination of the angular dependence of the magnitude of the Fresnel coefficient for transmission through a planar vacuum-metal interface, that it produces a transmitted field on the surface whose angular dependence has the form of an asymmetric peak. The peak maximum occurs at the critical angle, and has a magnitude that is twice that of the incident electric field on the surface, leading to enhanced diffuse scattering at this angle. This "Vineyard effect" was invoked by Sinha et al. [7] as the origin of the Yoneda peaks. This result is mathematically similar to the explanation provided by Eqs. (53) and (54), but it is not a physical explanation for these peaks.
Such an interpretation was offered by Kawanishi et al. [4], who suggested that the Yoneda peaks may be associated with the presence of lateral waves [19] propagating along the interface in the optically less dense medium. This wave satisfies the condition for refraction back into the optically more dense medium, and it therefore leaks energy at every position along the interface, along rays whose scattering angle θ s equals θ c [20]. This radiation is restricted to the range θ c < θ s < π/2. This is an attractive explanation, but it needs to be explored more through additional calculations.
We have also calculated the full angular intensity distributions of the reflected light. Figures 4 and 5 present such simulation results for the contribution to the mean differential reflection coefficient from the light that has been scattered incoherently by the randomly rough interface. The angles of incidence were set to (θ 0 , φ 0 ) = (0 • , 0 • ); it was cross-sectional cuts along the plane of incidence of these angular intensity distributions that resulted in the solid curves presented in Fig. 2 Fig. 2(a). The star notation, e.g. p → , indicates that the polarization of the scattered light was not recorded, hence Fig. 4(a) Fig. 2(b). Notice the rapid changes in intensity around the polar angle θs = θc = sin −1 ε1/ε2 [or q = √ ε2ω/c].
normal incidence, the scattering distributions are independent of the azimuthal angle of scattering φ s . This rotational symmetry is expected for isotropic surfaces like the ones we are investigating when the lateral momentum of the unpolarized incident light is zero. However, for ε 1 < ε 2 and when the incident light is linearly polarized but the polarization of the reflected light is not recorded, Figs. 4(a) and (d), we observe a slight skew in the distributions. This is similar to results presented in other, similar work [2,21], and is due to the subtle differences between the distributions of p → p and s → s scattered light, as presented for in-plane scattering in Fig. 2. When ε 1 > ε 2 the Yoneda peaks form a circle of equal intensity at the polar angle θ s = θ c [or q = √ ε 2 ω/c] in  Fig. 5(h). The position and circular symmetry of this groove can be understood through the previously mentioned factors of α 2 (q ) and k present in Eq. (53a) for p → p polarization and the factor α 2 (q ) present in Eq. (53c) for s → p polarization, since α 2 (q ) becomes zero when q = √ ε 2 ω/c and k is zero for normal incidence. It can be of interest to note that we also, as a consequence, observe a φ s -independent peak in Fig. 5(h), at a polar scattering angle significantly larger than θ c : the same peak as seen for p → p scattering in Fig. 2(b). However, this peak is not as sharp as the peak found at θ c in Fig. 5(i), and according to our definition it is not a Yoneda peak. Equations (53) demonstrate that the angular intensity distributions we are investigating can, to lowest order in the surface profile function, be explained through different factors in these equations with good approximation. As an aid in the interpretation of the results presented here and in the following, we notice that the power spectrum of the surface, g(|q − k |) is common for all equations in Eq. (53). As such, the mean DRC in SAPT to lowest order is essentially a distorted Gaussian on which critical angle effects are superposed.  As a starting point for our discussion of results for non-normal incidence, we in Fig. 6(a) present the angular dependence of the light scattered incoherently for a grazing angle of incidence from vacuum: θ 0 = 66.9 • . The scattering distribution for s → s scattered light can be seen to have retained its general shape from Fig. 2(a), but for p → p scattering we now observe a new feature: a local minimum at θ s ≈ 50 • . In the case of small amplitude perturbation theory, represented in Fig. 6(a) by dashed curves, ∂R pp (q |k )/∂Ω s incoh goes to zero at the position of this minimum.
The scattering angles defined by Θ B were first mentioned in the literature by Kawanishi et al. [4], where the angular values of Θ B were explored through a stochastic functional approach for two-dimensional surfaces. They chose to call the angles at which the first order contribution (according to their approach) to ∂R αp (q |k )/∂Ω s incoh vanishes the Brewster scattering angles, as a generalization of the Brewster angle for a flat surface. In what follows, following Kawanishi et al., we will call the polar angles of scattering in the plane of incidence at which p-and s-polarized light is scattered diffusely (incoherently) into light of any polarization with zero, or nearly zero, intensity, the Brewster scattering angles.
The Brewster angle θ B is defined by the zero in the reflectivity from a flat surface (coherent reflection in the specular direction) for p-polarization at the angle of incidence given by θ 0 = θ B = tan −1 ( ε 2 /ε 1 ). For one set of {ε 1 , ε 2 }, there is hence only one Brewster angle for incidence in a given medium. However, in contrast, we would like to stress the fact that the Brewster scattering angles for p → p scattering are present for a wide range of angles of incidence, given by Eq. (56) for in-plane scattering. From Eq. (56) it is also of interest to note that for light incident at the Brewster angle (for the corresponding flat-surface system), θ 0 = θ B , we find that Θ B (θ B ) = θ B ; the scattering intensity for light scattered incoherently vanishes for a scattering angle equal to the Brewster angle. This attests to the close relation between the Brewster angle for coherent reflection and the Brewster scattering angle Θ B for diffuse reflection, and is consistent with the findings of Kawanishi et al. [4]. Figure 6(b) presents simulation results for the same configuration as in Fig. 6(a), but for light scattered out-of-plane [q ·k = 0]. The dot productq ·k in Eq. (53d) indicates that, to lowest non-zero order in SAPT, we should not expect any contribution to the mean DRC from s → s out-of-plane incoherently scattered light. However, this is not the case for p → p scattered light, where, even forq ·k = 0, a closer look at Eq. (53a) indicates that the out-of-plane scattered intensity is zero only for θ s = 0 [q = 0]. This is precisely what we observe for ∂R pp (q |k )/∂Ω s incoh in Fig. 6(b). Figure 7 depicts results similar to those presented in Fig. 6 but for an increased surface rms-roughness of δ = λ/20 with the remaining parameters unchanged. As for normal incidence, it is found, not surprisingly, that small amplitude perturbation theory is most accurate for the smallest surface roughness. However, the most interesting feature to notice from Fig. 7(a), as compared to Fig. 6(a), is the angular position and amplitude of the local minimum of the in-plane p → p intensity distribution. In the former figure [Fig. 7(a)], the intensity at the position of the minimum is non-zero and it is located at an angle that is smaller than the Brewster scattering angle Θ B (θ 0 ) predicted by Eq. (56). We speculate that this shift in the Brewster scattering angle is roughness induced in a way similar to how the "normal" Brewster angle is shifted by the introduction of surface roughness.
The results presented in Fig. 6 were in-plane and out-of-plane cuts from Fig. 8, which presents the full angular distribution of the contributions to the mean DRC from incoherently scattered light for the angles of incidence (θ 0 , φ 0 ) = (66.9 • , 0 • ). Here the white dots indicate the lateral wavevector of the specular reflection, k . Compared to the results presented in Fig. 4, Fig. 8 displays many interesting features that are strongly dependent on both incoming and outgoing polarization, and we are in Fig. 8 left with symmetry in the distributions only about the plane of incidence. For p → p polarized reflection, Fig. 8(b), we observe that a significant fraction of the incoherently scattered light has shifted into the backscattering portion of the q -plane as the angle of incidence has increased. The opposite is true for s → s polarized reflection, Fig. 8(f), where the majority of the incoherently scattered light is scattered into the forward portion of the q -plane. This can be understood through small amplitude perturbation theory: In Eq. (53), the function F (q |k ) [Eq. (55)] constitutes the main difference between s → s and p → p polarized scattering, and it is easy to see that this term will enhance the backward scattering and reduce the forward scattering for p → p polarization. Additionally, the Brewster scattering angle, which for θ 0 = 66.9 • was given by Eq. (56) and found to be at θ s = 51.7 • for the parameters assumed, can now be seen to be part of a more general but still localized minimum in both Fig. 8(a) and Fig. 8(h), i.e., for p → and • → p scattering respectively. Figure 8(b) shows that the Brewster scattering angle for p → p polarized scattering can be found to be part of an interestingly shaped minimum in the q -plane. The shape of this minimum can, however, be extracted in a straightforward manner from Eq. (55). More interesting still is scattering in the inverse configuration, where light is incident from the dielectric side of the rough interface [ε 1 = 2.64, ε 2 = 1.0]: solutions of the RRE for this configuration and angles of incidence (θ 0 , φ 0 ) = (34.5 • , 0 • ), but for otherwise identical parameters as in Figs. 6 and 8, are presented in Fig. 9. Analogous with Fig. 6, Fig. 9(a) shows the incoherent component of the mean DRC for in-plane scattering, and Fig. 9(b) shows the corresponding curves for out-of-plane scattering.
In Fig. 9(a), we now observe that the two dips in ∂R pp (q |k )/∂Ω s incoh at |θ s | = θ c observed in Fig. 2(b) have both turned into Yoneda peaks, albeit with different peak intensities, and that the sharp dip at the same angle for forward scattering have turned into a less sharp local minimum at θ s ≈ 27 • . In order to understand these features, we see from Eq. (56) that, for ε 1 > ε 2 , ∂R pp (q |k )/∂Ω s incoh vanishes for θ s = Θ B (θ 0 ) when θ 0 is in the interval 0 < θ 0 < sin −1 ε 2 /ε 1 . This minimum in ∂R pp (q |k )/∂Ω s incoh will shift its polar position towards θ s = 0 • for increasing θ 0 , eventually "releasing" the Yoneda peaks in the forward scattering plane originating in the |d α (q )| 2 functions also for p → p scattering. In the backscattering plane, we observe through the function given in (55) that the negative sign of (q ·k ) will lead to a monotonic increase in the contribution from Eq. (55) to Eq. (53a) as θ 0 increases, eventually producing a Yoneda peak also forq ·k < 0. The overall distribution of s → s incoherent scattering in Fig. 9(a) also show a strong forward shift in its scattering intensities, which, as we look to Eq. (53d), can be attributed solely to the shifted power spectrum g(|q − k |).
Looking at Fig. 9(b), we observe several features for out-of-plane scattering that warrant a comment. Overall, we observe that the scattering distributions are again symmetric about θ s = 0, as is expected for out-of-plane scattering when the surface roughness is isotropic. Moreover, the distribution of ∂R ps (q |k )/∂Ω s incoh appears similar in shape to the distribution of ∂R pp (q |k )/∂Ω s incoh in Fig. 2(b). Their similarity can, to lowest non-zero order in SAPT, be attributed to their shared factor of α 2 (q ) in Eqs. (53a) and (53c), which in both cases vanishes for q = √ ε 2 ω/c, thereby suppressing the Yoneda peaks at this polar angle. There are no such suppressing factors present in Eq. (53b), and the distribution of ∂R sp (q |k )/∂Ω s incoh therefore displays Yoneda peaks at |θ s | = θ c . Similar to what we observed in Fig. 6(b), we see that the distribution of ∂R pp (q |k )/∂Ω s incoh has a local minimum at θ s = 0; both this minimum and the Yoneda peaks found at |θ s | = θ c are readily understood through the function in Eq. (55) and the factor |d α (q )| 2 , respectively. We now turn to a scattering system for which the rms-roughness of the surface is increased to δ = λ/20, i.e., twice of the roughness assumed in obtaining the results of Fig.9. Results for the in-plane and out-of-plane scattered intensity distributions for different combinations of the polarizations of the incident and scattered light are presented in Fig. 10. Overall the results in Fig. 10 are in qualitative agreement with those of Fig. 9 for the equivalent but less rough scattering system. In general, the increase in surface roughness is again found to result in a poorer agreement between the results obtained on the basis of SAPT and those obtained by a direct numerical solution of the RRE. However, it is interesting to observe that for the case of in-plane as well as out-of-plane p → p scattering, SAPT seems to give a fair representation of the simulated scattered intensity distributions for both levels of roughness considered in Figs. 9 and 10. For other combinations of the polarizations of the incident and scattered light this is not the case. The results presented in Fig. 9 were, as for Fig. 6, in-plane and out-of-plane cuts from Fig. 11 which displays the full angular distribution of the contribution to the mean DRC from the incoherently scattered light for the angles of incidence (θ 0 , φ 0 ) = (34.5 • , 0 • ). In contrast to what was observed in Fig. 8, all four of the lower left 2 × 2 subfigures in Fig. 11 now have significantly differing appearances. Similar to our observations in the case of incidence from vacuum, we observe that the Brewster scattering angle described by Eq. (56) can be seen to be part of a more general but still localized minimum in both Fig. 11(a) and Fig. 11(h), for p → and • → p scattering, respectively. Further, we still, as in Fig. 5, observe Yoneda peaks for all azimuthal angles of scattering. The intensities of these peaks are now, however, significantly stronger in the forward scattering plane, closer to the direction of specular reflection.
We now turn to Fig. 12, which is identical to Fig. 11 but for the angles of incidence (θ 0 , φ 0 ) = (45.5 • , 0 • ). For these angles of incidence, the light incident on a flat surface would exhibit total internal reflection. Incoherent scattering is, as before, greatly enhanced for q ≥ √ ε 2 ω/c, the part of wavevector-space that is evanescent in the medium of transmission. The intensity of the light scattered diffusely into this region is now comparable for s-and p-polarized light, and we see Yoneda peaks in both forward and backward scattering, for a fairly wide range of azimuthal angles. This can, as before, be understood to lowest non-zero order in SAPT through Eqs. (53) and (54). The factors |d p (k )| −2 and |d s (k )| −2 will both have their maxima at θ 0 = sin −1 ( ε 2 /ε 1 ), maxima that coincide with the corresponding maxima for the previously mentioned factors |d p (q )| −2 and |d s (q )| −2 . The contribution from these factors will be the same for all φ s , but common for all combinations of polarized scattering in Eq. (53) is that the multiplicative factor of the power spectrum will have its principal weight at |q − k | = 0; explaining the asymmetry about q 1 = 0 and the consequent shift of scattering to the forward scattering portion of the q -plane. While there is no Brewster scattering angle for the angle of incidence in Fig. 12, we still observe a local minimum in the backward scattering direction close to the critical angle for p-polarized incident light, Fig. 12(a).

C. Reflectivity and reflectance
The reflectivities for the two configurations of media are presented in Fig. 13. Both Figs. 13(a) and 13(b) show only small deviations from the Fresnel reflection coefficients for a corresponding flat-surface system, the only notable difference being in Fig. 13(b) where the surface roughness prevents total internal reflection for incoming light with θ 0 larger than θ c = sin −1 ε 2 /ε 1 ≈ 38.0 • , the critical angle corresponding to the values of the dielectric constants assumed in these simulations. The overall reflectivities for both systems are slightly smaller in all cases than the corresponding Fresnel coefficients, which is expected for a rough surface system since some light is scattered diffusely away from the specular direction. The rough-surface analogues of the Brewster angles for corresponding flat-interface systems, called analogues because the reflectivity does not reach strict zero in the case of surface roughness, are clearly  13. (a) The reflectivities Rα(θ0) of a two-dimensional randomly rough vacuum-dielectric interface [ε1 = 1.0, ε2 = 2.64] for p-and s-polarized light as functions of the polar angle of incidence. (b) The same as in Fig. 13(a), but for a dielectric-vacuum interface [ε1 = 2.64, ε2 = 1.0]. The quantity R F α (θ0) indicates the Fresnel reflection coefficient (flat surface reflectivity). The critical angle θ0 = θc = sin −1 ε2/ε1 for total internal reflection for an equivalent flat-interface system is indicated by a vertical dashed line in Fig. 13(b). Several simulations were run with small perturbations in the surface length L in order to obtain reflectivity data with higher angular resolution. The roughness parameters assumed in obtaining these results were δ = λ/40 and a = λ/4.
The θ0-dependence of the contribution to the reflectance from p-and s-polarized incident light that has been scattered incoherently from a two-dimensional randomly rough surface. This quantity is for β-polarized incident light defined as R β (θ0) incoh = R β (θ0) − R β (θ0). (a) The reflectances for a vacuum-dielectric interface [ε1 = 1.0, ε2 = 2.64] for p-and s-polarized light as functions of the polar angle of incidence. (b) Same as Fig. 14(a), but for a dielectric-vacuum interface [ε1 = 2.64, ε2 = 1.0]. As in Fig. 13, the critical angle for total internal reflection in a corresponding flat-interface system, θc, is indicated by a vertical dashed line in Fig. 14(b). Several simulations were run with small perturbations in the surface length L in order to obtain reflectance data with higher angular resolution. The roughness parameters assumed in obtaining these results were δ = λ/40 and a = λ/4. seen for p-polarized light in both figures.
The differences between the presented results for the reflectivity and the corresponding Fresnel coefficients can be better understood through Fig. 14, which presents the contribution to the reflectance from the light that has been reflected incoherently by the interface: R β (θ 0 ) incoh [see Eq. (52)]. In both subfigures in Fig. 14 we see that the amount of diffusely scattered light in general decreases with an increasing angle of incidence, if we ignore the effects of total internal reflection. This is consistent with the general notion that a rough surface is perceived as less rough for large angles of incidence [14]. Figure 14(a) shows that the incoherent part of the reflectance for the vacuum-dielectric configuration is a monotonically decreasing function of θ 0 for both polarizations, as expected by inspection of Eqs. (53) and (54), for ε 1 < ε 2 .
The functions |d p (k )| −2 and |d s (k )| −2 , and the factor 1/cos(θ 0 ) are all monotonically increasing functions of k (or θ 0 ), but they do not increase rapidly enough to compensate for the monotonously decreasing factor of α 2 1 (k ). A closer inspection of the numerical results and a more careful evaluation of the different factors in Eq. (53) have shown that the more rapid initial decrease of R p (θ 0 ) incoh is due to the contribution from its cross-polarized term, while its co-polarized term is responsible for the eventual less rapid decrease compared to R s (θ 0 ) incoh .
The incoherent part of the reflectance for the dielectric-vacuum configuration is displayed in Fig. 14(b). We can find here the explanation for why the curve for p-polarization in Fig. 13(b) showed a stronger peak for θ 0 just beyond the critical angle θ c than the curve for s-polarization: less light is scattered incoherently for these angles when the incident light is p-polarized than when it is s-polarized. We can also see that the contribution to the reflectance from the light scattered incoherently increases more than two-fold at the critical angle relative to the contribution at normal incidence. This behaviour can, again, be understood in terms of small amplitude perturbation theory to lowest order in the surface profile function, Eqs. (53) and (54). The functions |d p (k )| −2 and |d s (k )| −2 will both have their maximum at the critical angle θ c , but while R s (θ 0 ) incoh will get monotonously increasing contributions from both its co-and cross-polarized components for 0 < θ 0 < θ c , for R p (θ 0 ) incoh the cross-polarized component will go to zero due to the α 2 (k )-factor present in Eq. (53b). This dip in ∂R sp (q |k )/∂Ω s incoh is hence the main reason for the differences in the incoherent component of the reflectance for this configuration of media.

VII. CONCLUSIONS
We have presented a derivation of the reduced Rayleigh equation (RRE) for the reflection amplitudes of light scattered from a two-dimensional, randomly rough, surface. These equations enable a non-perturbative solution of the scattering problem based on the Rayleigh hypothesis. As an example of its solution by purely numerical means, the full angular distributions for both co-and cross-polarized incoherent components of the mean differential reflection coefficients were reported, for configurations of vacuum and an absorptionless dielectric with a Gaussian surface power spectrum and correlation function.
It was shown that a configuration of reflection within the optically denser medium leads to Yoneda peaks in the angular distributions of the diffusely scattered light, namely peaks at the critical angle for total internal reflection in the denser medium. The behaviour and development of these peaks for a wide range of angles of incidence and scattering were investigated, and the lack of such peaks for light scattered into p-polarization for polar angles of incidence smaller than the critical angle were explained through small amplitude perturbation theory (SAPT).
Brewster scattering angles, angles where scattering into p-polarization is suppressed to strict zero in SAPT to lowest non-zero order in the surface profile function, were found to explain many of the differences in scattering into s-and p-polarization for the scattering systems investigated in the current work. These angles were first mentioned in the literature by Kawanishi et al. in Ref. 4. Our results are in good agreement with their findings.
Small amplitude perturbation theory, to lowest non-zero order in the surface profile function, was overall shown to reproduce our numerical results qualitatively to a fairly high degree of accuracy, both through analytical arguments and a numerical implementation of that theory. This leads us to believe that the features presented in the results are single-scattering effects.
The scattering of light from a transparent dielectric is well described by solutions obtained by means of small amplitude perturbation theory, including the full angular distribution of the mean DRC for all combinations of the polarizations of the incident and scattered light. The reduced Rayleigh equation is a powerful starting point for studies of higher-order scattering features, such as enhanced backscattering, for example. The results presented here show that for the degree of surface roughness and the values of the dielectric constants assumed in this work no higher-order features are observed. Nevertheless, the RRE still gives more accurate numerical results for the mean DRC than does SAPT to lowest nonzero order in the surface profile function when the surface roughness is increased.
As an investigation of the quality of the results, energy conservation (unitarity) was found to be satisfied within 10 −4 when the total scattered energy from both reflection and transmission was added together, for the roughness parameters and configurations used in this paper. An investigation similar to the present one but for light transmitted through the dielectric rough interface will be presented in a separate publication [15].