The scattering of a scalar beam from isotropic and anisotropic two-dimensional randomly rough Dirichlet or Neumann surfaces: The full angular intensity distributions
Torstein Storflor Hegge, Torstein Nesse, Alexei A. Maradudin, Ingve Simonsen
TThe scattering of a scalar beam from isotropic and anisotropic two-dimensionalrandomly rough Dirichlet or Neumann surfaces: The full angular intensitydistributions
Torstein Storflor Hegge, Torstein Nesse, Alexei A. Maradudin, and Ingve Simonsen
1, 3, 4, ∗ Department of Physics, NTNU — Norwegian University of Science and Technology, NO-7491 Trondheim, Norway Department of Physics and Astronomy, University of California, Irvine, CA 92697, U.S.A. Department of Petroleum Engineering, University of Stavanger, NO-4036 Stavanger, Norway Surface du Verre et Interfaces, UMR 125 CNRS/Saint-Gobain, F-93303 Aubervilliers, France (Dated: September 20, 2018)By the use of Green’s second integral identity we determine the field scattered from a two-dimensional randomly rough isotropic or anisotropic Dirichlet or Neumann surface when it is illu-minated by a scalar Gaussian beam. The integral equations for the scattering amplitudes are solvednonperturbatively by a rigorous computer simulation approach. The results of these calculationsare used to calculate the full angular distribution of the mean differential reflection coefficient. Forisotropic surfaces, the results of the present calculations for in-plane scattering are compared withthose of earlier studies of this problem. The reflectivities of Dirichlet and Neumann surfaces arecalculated as functions of the polar angle of incidence, and the reflectiveties for the two kinds ofsurfaces of similar roughness parameters are found to be different. For an increasing level of surfaceanisotropy, we study how the angular intensity distributions of the scattered waves are affected bythis level. We find that even small to moderate levels of surface anisotropy can significantly alterthe symmetry, shape, and amplitude of the scattered intensity distributions when Gaussian beamsare incident on the anisotropic surfaces from different azimuthal angles of incidence. ∗ [email protected] a r X i v : . [ phy s i c s . op ti c s ] D ec I. INTRODUCTION
The earliest nonperturbative calculations of the scattering of a field from a two-dimensional randomly rough surfacewere the studies of the scattering of a scalar beam, incident from vacuum, on a Dirichlet [1–3] or a Neumann surface [3]carried out by Tran and Maradudin and by Macaskill and Kachoyan. These calculations were based on Green’s secondintegral identity [4]. The integral equations for the source functions, namely the values of the field in the vacuumor its normal derivative, evaluated on the rough surface, were transformed into matrix equations which were thensolved by iterative approaches. The amplitudes of the scattered field are expressed in terms of these source functions,and the differential reflection coefficient is expressed through the scattering amplitudes. The differential reflectioncoefficient (DRC), an experimentally accessible quantity, gives the fraction of the total time-averaged flux incident onthe rough surface that is scattered into an element of solid angle about a specified direction of scattering. In scatteringfrom a randomly rough surface it is the average of the DRC over the ensemble of realizations of the surface profilefunction that is calculated. The result is called the mean differential reflection coefficient (mean DRC). Multiplescattering effects, in particular enhanced backscattering [5], were present in the results for the dependence of themean DRC for in-plane scattering on the polar angle of scattering and a fixed polar angle of incidence.Although in the years following this pioneering work several nonperturbative calculations of the scattering of vectorfields from impenetrable [6–9] and penetrable [10–15] two-dimensional randomly rough surfaces were carried out, littleattention seems to have been directed at rigorous nonperturbative calculations of the scattering of incident beamsfrom Dirichlet and Neumann surfaces perhaps because they are simpler than the scattering problems studied in thesereferences. Nevertheless, the results of these calculations are relevant, for example, in ocean acoustics in the contextof the scattering of a sonic wave from a rough ocean floor [16, 17].In this paper we revisit the problem of the scattering of a scalar beam from a two-dimensional randomly roughsurface, and investigate properties of the scattered field not considered in the earliest studies of this problem [1–3].Thus, in addition to presenting results for scattering from surfaces whose profiles are isotropic Gaussian randomprocesses we also present results for the scattering from surfaces whose profiles are anisotropic Gaussian randomprocesses. In addition to the contribution to the mean differential reflection coefficient from the field scatteredincoherently in plane, we also present results for the reflectiveties of these surfaces and the full angular distribution ofthe intensity of the scattered field. Moreover, these calculations are carried out by means of improved algorithms thatyield accurate solutions of the integral equations arising in the scattering theory without the use of iterative methodsof the Sturm-Liouville type or modifications thereof [1–3].
II. SCATTERING SYSTEM
The system we consider in this work consists of a medium that supports the propagation of scalar waves withoutabsorption, e.g. a liquid, in the region x > ζ ( x (cid:107) ), where x (cid:107) = ( x , x ,
0) is an arbitrary vector in the plane x = 0,and a medium that is impenetrable to scalar waves in the region x < ζ ( x (cid:107) ) [Fig. 1]. The surface profile function ζ ( x (cid:107) ) is assumed to be a single-valued function of x (cid:107) that is differentiable with respect to x and x , and constitutesa stationary, zero-mean, Gaussian random process. It is defined by (cid:10) ζ ( x (cid:107) ) (cid:11) = 0 (1a) (cid:68) ζ ( x (cid:107) ) ζ ( x (cid:48)(cid:107) ) (cid:69) = δ W ( x (cid:107) − x (cid:48)(cid:107) ) , (1b)where the angle brackets here and in all that follows denote an average over the ensemble of realizations of the surfaceprofile function. The quantity δ , the root-mean-square roughness of the surface, is defined by δ = (cid:10) ζ ( x (cid:107) ) (cid:11) . (2)The function W ( x (cid:107) ) introduced in Eq. (1b) is the normalized surface height autocorrelation function , and has theproperty that, W ( ) = 1. In what follows we will also require the power spectrum of the surface roughness, g ( k (cid:107) ),where k (cid:107) is a two-dimensional wave vector k (cid:107) = ( k , k , g ( k (cid:107) ) = (cid:90) d x (cid:107) W ( x (cid:107) ) exp (cid:0) − i k (cid:107) · x (cid:107) (cid:1) . (3)In this work we will assume the following Gaussian form for W ( x (cid:107) ) [8] W ( x (cid:107) ) = exp (cid:18) − x a − x a (cid:19) , (4) x x x qk q k k k φ s φ θ s θ FIG. 1: Schematics of the scattering geometry.where the characteristic lengths a > a > x and x axes, respectively.For the choice of W ( x (cid:107) ) given by Eq. (4), the power spectrum has the Gaussian form g ( k (cid:107) ) = πa a exp (cid:18) − a k − a k (cid:19) . (5)The surface roughness is said to be anisotropic when a (cid:54) = a , and is called isotropic when a = a . In the latter casethe surface height autocorrelation function W ( x (cid:107) ) and the power spectrum g ( k (cid:107) ) depend on x (cid:107) and k (cid:107) , respectively,only through their magnitudes x (cid:107) and k (cid:107) and not on their directions. III. SCATTERING THEORY
The impenetrable surface x = ζ ( x (cid:107) ) is illuminated from the region x > ζ ( x (cid:107) ) by a scalar field of angularfrequency ω , that produce a scattered field of the same frequency. The field ψ ( x ; t ) in the region x > ζ ( x (cid:107) ) has theform ψ ( x ; t ) = ψ ( x | ω ) exp[ − i ωt ] where the amplitude function ψ ( x | ω ) is the solution of the Helmholtz equation (cid:20) ∇ + ω c (cid:21) ψ ( x | ω ) = 0 , (6)with c the speed of the field. This field satisfies either (a) the Dirichlet boundary condition, which corresponds tozero pressure on the wall, ψ ( x | ω ) | x = ζ ( x (cid:107) ) = 0 , (7)or (b) the Neumann boundary condition, which corresponds to zero normal velocity at the wall, ∂∂n ψ ( x | ω ) (cid:12)(cid:12)(cid:12)(cid:12) x = ζ ( x (cid:107) ) = 0 . (8)In Eq. (8) ∂/∂n is the derivative along the normal to the surface x = ζ ( x (cid:107) ) at each point of it, directed into themedium of incidence, ∂∂n = 1 (cid:104) (cid:8) ∇ ζ ( x (cid:107) ) (cid:9) (cid:105) (cid:20) − ζ ( x (cid:107) ) ∂∂x − ζ ( x (cid:107) ) ∂∂x + ∂∂x (cid:21) (9a) ≡ (cid:104) (cid:8) ∇ ζ ( x (cid:107) ) (cid:9) (cid:105) ∂∂N , (9b)where ζ α ( x (cid:107) ) = ∂ζ ( x (cid:107) ) /∂x α ( α = 1 , ψ ( x | ω ) that is convenient to solve numerically we begin by introducing the scalarGreen’s function that satisfies the inhomogeneous Helmholtz equation (cid:20) ∇ + ω c (cid:21) g ( x | x (cid:48) ) = − π δ ( x − x (cid:48) ) . (10)The Green’s function has the representations g ( x | x (cid:48) ) = exp (cid:104) i ωc | x − x (cid:48) | (cid:105) | x − x (cid:48) | (11a)= (cid:90) d q (cid:107) (2 π ) π i α ( q (cid:107) ) exp (cid:104) i q (cid:107) · (cid:16) x (cid:107) − x (cid:48)(cid:107) (cid:17)(cid:105) exp (cid:2) i α ( q (cid:107) ) | x − x (cid:48) | (cid:3) , (11b)where we have introduced the in-plane component of the wave vector q (cid:107) = ( q , q ,
0) and its corresponding normalcomponent α ( q (cid:107) ) = (cid:104) ω c − q (cid:107) (cid:105) q (cid:107) < ωc i (cid:104) q (cid:107) − ω c (cid:105) q (cid:107) > ωc . (12)We next apply Green’s second integral identity [4, p. 152] (cid:90) Ω d x (cid:107) (cid:0) u ∇ v − v ∇ u (cid:1) = (cid:90) Σ d s (cid:18) u ∂v∂ν − v ∂u∂ν (cid:19) , (13)where u ( x ) and v ( x ) are arbitrary scalar functions of x defined in a volume Ω that is bounded by a closed surface Σ.The derivative ∂/∂ν is taken along the normal to the surface Σ at each point of it, directed away from the volume Ω.We assume that the volume Ω is the region x > ζ ( x (cid:107) ), while the surface Σ is the union of the rough surface x = ζ ( x (cid:107) ), which we denote by S , and a hemispherical cap of infinite radius in the upper half space, which we denoteby S ( ∞ ) . Then, on setting u = ψ ( x | ω ) and v = g ( x | x (cid:48) ) in Eq. (13), and taking into account Eqs. (6) and (10), weobtain − π θ (cid:16) x (cid:48) − ζ ( x (cid:48)(cid:107) ) (cid:17) ψ ( x (cid:48) | ω ) = − (cid:90) S d S (cid:20) ψ ( x | ω ) ∂∂n g ( x | x (cid:48) ) − g ( x | x (cid:48) ) ∂∂n ψ ( x | ω ) (cid:21) + (cid:90) S ( ∞ ) d S (cid:20) ψ ( x | ω ) ∂∂ν g ( x | x (cid:48) ) − g ( x | x (cid:48) ) ∂∂ν ψ ( x | ω ) (cid:21) , (14)where θ ( z ) is the Heaviside unit step function. Because the scattered field satisfies a radiation condition at infinity,its contribution to the surface integral over the hemispherical cap of infinite radius S ( ∞ ) on the right-hand side ofEq. (14) vanishes. This integral therefore yields − πψ ( x | ω ) inc where ψ ( x | ω ) inc is the incident field. By using thesymmetry of g ( x | x (cid:48) ), we can now rewrite Eq. (14) as θ (cid:0) x − ζ ( x (cid:107) ) (cid:1) ψ ( x | ω ) = ψ ( x | ω ) inc + 14 π (cid:90) S d S (cid:48) (cid:20)(cid:18) ∂∂n (cid:48) g ( x | x (cid:48) ) (cid:19) ψ ( x (cid:48) | ω ) − g ( x | x (cid:48) ) ∂∂n (cid:48) ψ ( x (cid:48) | ω ) (cid:21) . (15)Since we have assumed that the surface profile function is a single valued function of x (cid:107) , we can replace the integrationover the surface S by integration over the plane x = 0 with the use of the relation d S = (cid:104) (cid:8) ∇ ζ ( x (cid:107) ) (cid:9) (cid:105) d x (cid:107) , θ (cid:0) x − ζ ( x (cid:107) ) (cid:1) ψ ( x | ω ) = ψ ( x | ω ) inc + 14 π (cid:90) d x (cid:48)(cid:107) (cid:40) (cid:20) ∂∂N (cid:48) g ( x | x (cid:48) ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) x (cid:48) = ζ ( x (cid:48)(cid:107) ) F ( x (cid:48)(cid:107) | ω ) − [ g ( x | x (cid:48) )] | x (cid:48) = ζ ( x (cid:48)(cid:107) ) N ( x (cid:48)(cid:107) | ω ) (cid:41) , (16)where F ( x (cid:107) | ω ) = ψ ( x | ω ) | x = ζ ( x (cid:107) ) (17a) N ( x (cid:107) | ω ) = ∂∂N ψ ( x | ω ) (cid:12)(cid:12)(cid:12)(cid:12) x = ζ ( x (cid:107) ) . (17b)With the use of either the Dirichlet or Neumann boundary condition, we can simplify Eq. (16). A. Dirichlet boundary condition
The Dirichlet boundary condition, Eq. (7), can be called an acoustically soft-wall boundary condition. Its assump-tion simplifies Eq. (16) to θ (cid:0) x − ζ ( x (cid:107) ) (cid:1) ψ ( x | ω ) = ψ ( x | ω ) inc − π (cid:90) d x (cid:48)(cid:107) [ g ( x | x (cid:48) )] | x (cid:48) = ζ ( x (cid:48)(cid:107) ) N ( x (cid:48)(cid:107) | ω ) . (18)The scattered field is the second term on the right-hand side of Eq. (18). With the use of the representation for g ( x | x (cid:48) ) given by Eq. (11b), the scattered field can be written as ψ ( x | ω ) sc = (cid:90) d q (cid:107) (2 π ) R D ( q (cid:107) , ω ) exp (cid:2) i q (cid:107) · x (cid:107) + i α ( q (cid:107) ) x (cid:3) , (19)where R D ( q (cid:107) , ω ) = − i2 α ( q (cid:107) ) (cid:90) d x (cid:107) N ( x (cid:107) | ω ) exp (cid:2) − i q (cid:107) · x (cid:107) − i α ( q (cid:107) ) ζ ( x (cid:107) ) (cid:3) . (20) B. Neumann boundary condition
The Neumann boundary condition, Eq. (8), can be called an acoustically hard-wall boundary condition. Its usesimplifies Eq. (16) to θ (cid:0) x − ζ ( x (cid:107) ) (cid:1) ψ ( x | ω ) = ψ ( x | ω ) inc + 14 π (cid:90) d x (cid:48)(cid:107) (cid:20) ∂∂N (cid:48) g ( x | x (cid:48) ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) x (cid:48) = ζ ( x (cid:48)(cid:107) ) F ( x (cid:48)(cid:107) | ω ) . (21)The scattered field is the second term on the right-hand side of this equation. With the use of the representation of g ( x | x (cid:48) ) given by Eq. (11b), it assumes the form ψ ( x | ω ) sc = (cid:90) d q (cid:107) (2 π ) R N ( q (cid:107) , ω ) exp (cid:2) i q (cid:107) · x (cid:107) + i α ( q (cid:107) ) x (cid:3) , (22)where R N ( q (cid:107) , ω ) = − α ( q (cid:107) ) (cid:90) d x (cid:107) F ( x (cid:107) | ω ) (cid:2) q (cid:107) · ∇ ζ ( x (cid:107) ) − α ( q (cid:107) ) (cid:3) exp (cid:2) − i q (cid:107) · x (cid:107) − i α ( q (cid:107) ) ζ ( x (cid:107) ) (cid:3) . (23) IV. THE MEAN DIFFERENTIAL REFLECTION COEFFICIENT
The scattering amplitudes R D ( q (cid:107) , ω ) and R N ( q (cid:107) , ω ) play a central role in the theory of the scattering of a wavefrom a randomly rough surface. This is because the coherent (specular) and incoherent (diffuse) scattering of thewave can be described quantitatively in terms of these amplitudes, through their presence in the expressions for themean differential reflection coefficient.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 direction ofscattering defined by the polar and azimuthal angles of scattering θ s and φ s , respectively.Because in numerical simulations of scattering from a rough surface only a finite region of the plane x = 0 can becovered by a two-dimensional rough surface, we assume for the incident field a Gaussian beam that illuminates onlythat finite region of the x plane, minimizing edge effects in the scattering thereby. The Gaussian incident field isgiven by ψ ( x | ω ) inc = w π (cid:90) q (cid:107) < ωc d q (cid:107) exp (cid:2) i q (cid:107) · x (cid:107) − i α ( q (cid:107) ) x (cid:3) exp (cid:20) − w (cid:0) q (cid:107) − k (cid:107) (cid:1) (cid:21) . (24)The magnitude of the total time-averaged flux incident on the surface is given by P inc = − A Im L (cid:90) − L d x L (cid:90) − L d x ψ ∗ ( x | ω ) inc ∂ψ ( x | ω ) inc ∂x , (25)where L and L are the lengths of the scattering surface along the x and x axes, respectively, while A is a coefficientthat drops out of the expression for the differential reflection coefficient [18]. The minus sign that appears on theright-hand side of Eq. (25) compensates for the fact that the incident flux is negative. On substituting Eq. (24) intoEq. (25), we obtain P inc = A (cid:18) w π (cid:19) Im i (cid:90) q (cid:107) < ωc d q (cid:107) (cid:90) q (cid:48)(cid:107) < ωc d q (cid:48)(cid:107) α ( q (cid:48)(cid:107) ) exp (cid:110) i (cid:104) α ∗ ( q (cid:107) ) − α ( q (cid:48)(cid:107) ) (cid:105) x (cid:111) × exp (cid:20) − w (cid:0) q (cid:107) − k (cid:107) (cid:1) − w (cid:16) q (cid:48)(cid:107) − k (cid:107) (cid:17) (cid:21) L (cid:90) − L d x L (cid:90) − L d x exp (cid:104) − i (cid:16) q (cid:107) − q (cid:48)(cid:107) (cid:17) · x (cid:107) (cid:105) = Aw (cid:90) q (cid:107) < ωc d q (cid:107) α ( q (cid:107) ) exp (cid:104) − w (cid:0) q (cid:107) − k (cid:107) (cid:1) (cid:105) , (26)in the limit as L , → ∞ . By carrying out the angular integration in Eq. (26) we obtain an expression for P inc as aone-dimensional integral P inc = 2 πAw exp (cid:16) − w k (cid:107) (cid:17) ωc (cid:90) d q (cid:107) q (cid:107) α ( q (cid:107) ) I (cid:0) w q (cid:107) k (cid:107) (cid:1) exp (cid:16) − w q (cid:107) (cid:17) , (27)where I ( z ) is the modified Bessel function of the first kind and order zero. An alternative expression for P inc isobtained when we make the change of variable q (cid:107) = ( ω/c ) sin ϑ , namely P inc = 2 πAw (cid:16) ωc (cid:17) exp (cid:16) − w k (cid:107) (cid:17) π (cid:90) d ϑ sin ϑ cos ϑ I (cid:16) w ωc k (cid:107) sin ϑ (cid:17) exp (cid:20) − (cid:16) wωc (cid:17) sin ϑ (cid:21) . (28)The relation between the wave vector k (cid:107) and the polar and azimuthal angles of incidence, θ and φ , respectively, is k (cid:107) = ωc sin θ (cos φ , sin φ , , (29)so that k (cid:107) = ( ω/c ) sin θ . With this result the expression for P inc , Eq. (28), becomes P inc = P inc ( θ ) = A p inc ( θ ) , (30a)where p inc ( θ ) = 2 πw (cid:16) ωc (cid:17) exp (cid:20) − (cid:16) wωc (cid:17) sin θ (cid:21) π (cid:90) d ϑ sin ϑ cos ϑ I (cid:18) (cid:16) wωc (cid:17) sin θ sin ϑ (cid:19) exp (cid:20) − (cid:16) wωc (cid:17) sin ϑ (cid:21) . (30b)The expressions for P inc given by Eqs. (27) and (30) have to be evaluated numerically.The field scattered from either a Dirichlet or Neumann surface can be written as (see Eqs. (19) and (22)) ψ ( x | ω ) sc = (cid:90) d q (cid:107) (2 π ) R ( q (cid:107) , ω ) exp (cid:2) i q (cid:107) · x (cid:107) + i α ( q (cid:107) ) x (cid:3) , (31)where R ( q (cid:107) , ω ) is either R D ( q (cid:107) , ω ) or R N ( q (cid:107) , ω ). The total time-averaged scattered flux is given by P sc = A Im L (cid:90) − L d x L (cid:90) − L d x ψ ∗ ( x | ω ) sc ∂ψ ( x | ω ) sc ∂x . (32)When we substitute Eq. (31) into Eq. (32) the latter becomes P sc = A Im (cid:90) d q (cid:107) (2 π ) (cid:90) d q (cid:48)(cid:107) (2 π ) i α ( q (cid:48)(cid:107) ) R ∗ ( q (cid:107) , ω ) R ( q (cid:48)(cid:107) , ω ) exp (cid:110) − i (cid:2) α ∗ ( q (cid:107) ) − α ( q (cid:48)(cid:107) ) (cid:3) x (cid:111) × L (cid:90) − L d x L (cid:90) − L d x exp (cid:104) − i (cid:0) q (cid:107) − q (cid:48)(cid:107) (cid:1) · x (cid:107) (cid:105) = A Im i (cid:90) d q (cid:107) (2 π ) α ( q (cid:107) ) (cid:12)(cid:12) R ( q (cid:107) , ω ) (cid:12)(cid:12) exp (cid:2) − α ( q (cid:107) ) x (cid:3) , (33)in the limit where L , → ∞ . The function α ( q (cid:107) ) is real for 0 < q (cid:107) < ω/c , and imaginary for q (cid:107) > ω/c . Thus wefinally obtain for P sc the result P sc = A (cid:90) q (cid:107) <ω/c d q (cid:107) (2 π ) α ( q (cid:107) ) (cid:12)(cid:12) R ( q (cid:107) , ω ) (cid:12)(cid:12) . (34)The relation between the wave vector q (cid:107) and the polar and azimuthal angles of scattering θ s and φ s is q (cid:107) = ωc sin θ s (cos φ s , sin φ s , , (35)so that q (cid:107) = ωc sin θ s (36a) α ( q (cid:107) ) = ωc cos θ s , (36b)while d q (cid:107) = (cid:16) ωc (cid:17) sin θ s dΩ s , (36c)where dΩ s = sin θ s d θ s d φ s is the element of solid angle at ( θ s , φ s ). The total time-averaged scattered flux can thenbe written in the form P sc = (cid:90) dΩ s P sc ( θ s , φ s ) , (37)where P sc ( θ s , φ s ) = A (cid:16) ω πc (cid:17) ωc cos θ s (cid:12)(cid:12) R ( q (cid:107) , ω ) (cid:12)(cid:12) . (38)By definition the differential reflection coefficient is given by ∂R∂ Ω s = P sc ( θ s , φ s ) P inc ( θ ) = 14 π (cid:16) ωc (cid:17) cos θ s p inc ( θ ) (cid:12)(cid:12) R ( q (cid:107) , ω ) (cid:12)(cid:12) . (39)Since we are considering scattering from a randomly rough surface, it is not the differential reflection coefficientitself that we need to calculate, but rather its average over the ensemble of realizations of the surface profile function.The resulting mean differential reflection coefficient is given by (cid:28) ∂R∂ Ω s (cid:29) = 14 π (cid:16) ωc (cid:17) cos θ s p inc ( θ ) (cid:68)(cid:12)(cid:12) R ( q (cid:107) , ω ) (cid:12)(cid:12) (cid:69) . (40)If we write the scattering amplitude R ( q (cid:107) | k (cid:107) ) as the sum of its mean value and of its fluctuation away from themean value, R ( q (cid:107) | k (cid:107) ) = (cid:10) R ( q (cid:107) | k (cid:107) ) (cid:11) + (cid:2) R ( q (cid:107) | k (cid:107) ) − (cid:10) R ( q (cid:107) | k (cid:107) ) (cid:11)(cid:3) , (41)we find that each term contributes separately to the mean differential reflection coefficient, which then takes the form (cid:28) ∂R∂ Ω s (cid:29) = (cid:28) ∂R∂ Ω s (cid:29) coh + (cid:28) ∂R∂ Ω s (cid:29) incoh , (42)where (cid:28) ∂R∂ Ω s (cid:29) coh = 14 π (cid:16) ωc (cid:17) cos θ s p inc ( θ ) (cid:12)(cid:12)(cid:10) R ( q (cid:107) , ω ) (cid:11)(cid:12)(cid:12) (43)and (cid:28) ∂R∂ Ω s (cid:29) incoh = 14 π (cid:16) ωc (cid:17) cos θ s p inc ( θ ) (cid:20)(cid:68)(cid:12)(cid:12)(cid:12) R ( q (cid:107) , ω ) (cid:12)(cid:12)(cid:12) (cid:69) − (cid:12)(cid:12)(cid:12)(cid:68) R ( q (cid:107) , ω ) (cid:69)(cid:12)(cid:12)(cid:12) (cid:21) . (44)The former term gives the contribution to the mean differential reflection coefficient from the coherently (specularly)scattered field, while the latter term gives the contribution from the incoherently (diffusely) scattered field.The reflectivity of the randomly rough surface is given by R ( θ , φ ) = π (cid:90) d θ s sin θ s π (cid:90) − π d φ s (cid:28) ∂R∂ Ω s (cid:29) coh . (45)The dependence of the reflectivity on the azimuthal angle of incidence φ arises only when the randomly roughsurface is defined by a surface profile fucntion that is a stationary anisotropic random process. When the surfaceprofile function is a stationary isotropic random process, the reflectivity is independent of φ . V. EQUATIONS SATISFIED BY THE SOURCE FUNCTIONS
We see from Eqs. (18) and (21) that once the source functions F ( x (cid:107) | ω ) and N ( x (cid:107) | ω ) are known these equationsallow the scattered fields to be determined at any point in the region x > ζ ( x (cid:107) ). To obtain the equations satisfiedby these functions we consider the cases of Dirichlet and Neumann surfaces in turn. A. Dirichlet Surfaces
To obtain the equation for the source function N ( x (cid:107) | ω ) we first assume that x > ζ ( x (cid:107) ) and apply the derivativeoperator ∂/∂N to both sides of Eq. (18) to obtain θ (cid:0) x − ζ ( x (cid:107) ) (cid:1) ∂∂N ψ ( x | ω ) = ∂∂N ψ ( x | ω ) inc − π (cid:90) d x (cid:48)(cid:107) (cid:20) ∂∂N g ( x | x (cid:48) ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) x (cid:48) = ζ ( x (cid:48)(cid:107) ) N ( x (cid:48)(cid:107) | ω ) . (46)We next evaluate this equation at x = ζ ( x (cid:107) ) + η and at x = ζ ( x (cid:107) ) − η , where η is a positive infinitesimal, add theresulting equations, and let η →
0. The result is N ( x (cid:107) | ω ) = 2 N ( x (cid:107) | ω ) inc − π (cid:90) d x (cid:48)(cid:107) (cid:20) ∂∂N g ( x | x (cid:48) ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) x = ζ ( x (cid:107) )+ ηx (cid:48) = ζ ( x (cid:48)(cid:107) ) + (cid:20) ∂∂N g ( x | x (cid:48) ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) x = ζ ( x (cid:107) ) − ηx (cid:48) = ζ ( x (cid:48)(cid:107) ) N ( x (cid:48)(cid:107) | ω ) , (47)where N ( x (cid:107) | ω ) inc = [ ∂ψ ( x | ω ) inc /∂N ] (cid:12)(cid:12) x = ζ ( x (cid:107) ) . With the use of the results (cid:20) ∂∂N g ( x | x (cid:48) ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) x = ζ ( x (cid:107) )+ ηx (cid:48) = ζ ( x (cid:48)(cid:107) ) + (cid:20) ∂∂N g ( x | x (cid:48) ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) x = ζ ( x (cid:107) ) − ηx (cid:48) = ζ ( x (cid:48)(cid:107) ) = 2 P (cid:115) ∂∂N g ( x | x (cid:48) ) (cid:123) , (48a)where P denotes the Cauchy principle value, and (cid:113) f ( x | x (cid:48) ) (cid:121) = f ( x | x (cid:48) ) (cid:12)(cid:12)(cid:12) x = ζ ( x (cid:107) ) x (cid:48) = ζ ( x (cid:48)(cid:107) ) , (48b)we obtain finally the equation satisfied by N ( x (cid:107) | ω ) N ( x (cid:107) | ω ) = 2 N ( x (cid:107) | ω ) inc − π P (cid:90) d x (cid:48)(cid:107) K D ( x (cid:107) | x (cid:48)(cid:107) ) N ( x (cid:48)(cid:107) | ω ) , (49a)where we have defined the (Dirichlet) kernel K D ( x (cid:107) | x (cid:48)(cid:107) ) = (cid:115) ∂∂N g ( x | x (cid:48) ) (cid:123) . (49b) B. Neumann Surfaces
To obtain the equation satisfied by the source function F ( x (cid:107) | ω ), we begin by evaluating Eq. (21) at x = ζ ( x (cid:107) ) + η and x = ζ ( x (cid:107) ) − η , adding the resulting equations and then letting η →
0. The result is F ( x (cid:107) | ω ) = 2 F ( x (cid:107) | ω ) inc + 14 π (cid:90) d x (cid:48)(cid:107) (cid:20) ∂∂N (cid:48) g ( x | x (cid:48) ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) x = ζ ( x (cid:107) )+ ηx (cid:48) = ζ ( x (cid:48)(cid:107) ) + (cid:20) ∂∂N (cid:48) g ( x | x (cid:48) ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) x = ζ ( x (cid:107) ) − ηx (cid:48) = ζ ( x (cid:48)(cid:107) ) F ( x (cid:48)(cid:107) | ω ) , (50)where F ( x (cid:107) | ω ) inc = ψ ( x | ω ) inc (cid:12)(cid:12) x = ζ ( x (cid:107) ) . We next use the result (cid:20) ∂∂N (cid:48) g ( x | x (cid:48) ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) x = ζ ( x (cid:107) )+ ηx (cid:48) = ζ ( x (cid:48)(cid:107) ) + (cid:20) ∂∂N (cid:48) g ( x | x (cid:48) ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) x = ζ ( x (cid:107) ) − ηx (cid:48) = ζ ( x (cid:48)(cid:107) ) = 2 P (cid:115) ∂∂N (cid:48) g ( x | x (cid:48) ) (cid:123) , (51)to obtain finally the equation satisfied by the source function F ( x (cid:107) | ω ) F ( x (cid:107) | ω ) = 2 F ( x (cid:107) | ω ) inc + 12 π P (cid:90) d x (cid:48)(cid:107) K N ( x (cid:107) | x (cid:48)(cid:107) ) F ( x (cid:48)(cid:107) | ω ) , (52a)where the (Neumann) kernel is K N ( x (cid:107) | x (cid:48)(cid:107) ) = (cid:115) ∂∂N (cid:48) g ( x | x (cid:48) ) (cid:123) . (52b) VI. NUMERICAL SOLUTION OF THE EQUATIONS FOR THE SOURCE FUNCTIONS
The integral equations (49) and (52) satisfied by the source functions N ( x (cid:107) | ω ) and F ( x (cid:107) | ω ), respectively, have tobe solved numerically. In this section we show how this is done in scattering from a Dirichlet and Neumann surface.To solve Eq. (49) we first replace integration over the entire x (cid:48) x (cid:48) plane by integration over the finite square regiondefined by − L/ < x (cid:48) < L/ − L/ < x (cid:48) < L/
2, that is much larger than the illuminated region of the surface. Arealization of the surface profile function is then generated numerically, by a two-dimensional extension of the methoddescribed in Refs. [19, 20] on a grid of (2 N +1) points within this square region of the x (cid:48) x (cid:48) plane [8]. The coordinatesof the grid points are x (cid:107) ( (cid:96) ) = ( (cid:96) , (cid:96) , x. (53)0In this expression (cid:96) and (cid:96) are integers that each take the values − N , − N + 1, . . . , N − N and which we denotecollectively by (cid:96) , while ∆ x = L/ (2 N + 1).We next rewrite Eq. (49) as N ( x (cid:107) | ω ) = 2 N ( x (cid:107) | ω ) inc − π N (cid:88) (cid:96) (cid:48) = − N N (cid:88) (cid:96) (cid:48) = − N P ( (cid:96) (cid:48) + )∆ x (cid:90) ( (cid:96) (cid:48) − )∆ x d x (cid:48) (cid:96) (cid:48) + )∆ x (cid:90) ( (cid:96) (cid:48) − )∆ x d x (cid:48) K D ( x (cid:107) | x (cid:48)(cid:107) ) N ( x (cid:48)(cid:107) | ω ) , (54)where K D ( x (cid:107) | x (cid:48)(cid:107) ) = exp (cid:104) i ωc (cid:12)(cid:12) x (cid:107) − x (cid:48)(cid:107) (cid:12)(cid:12)(cid:105) i ω/c (cid:12)(cid:12) x (cid:107) − x (cid:48)(cid:107) (cid:12)(cid:12) − (cid:12)(cid:12) x (cid:107) − x (cid:48)(cid:107) (cid:12)(cid:12) (cid:104) − (cid:16) x (cid:107) − x (cid:48)(cid:107) (cid:17) · ∇ ζ ( x (cid:107) ) + ( x − x (cid:48) ) (cid:105)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = ζ ( x (cid:107) ) x (cid:48) = ζ ( x (cid:48)(cid:107) ) . (55)We assume that N ( x (cid:48)(cid:107) | ω ) is a slowly varying function fucntion of x (cid:48) and x (cid:48) in each of the intervals ( (cid:96) (cid:48) − )∆ x Rigorous computer simulations were carried out to obtain the field that is scattered from isotropic or anisotropic randomly rough Dirichlet or Neumann surfaces. These calculations were performed by numerically solving the inho-mogeneous integral equations (49) and (52) by the method outlined in Sec. VI; that is, we solved the linear systemof equation in Eqs. (58) and (62). From their solutions, the reflection amplitudes for Dirichlet or Neumann surfaceswere calculated from Eqs. (20) and (23), respectively, and the results were subsequently used to obtain the meanDRCs defined in Eq. (40). The randomly rough surfaces were assumed to constitute a Gaussian random process thatis characterized by the Gaussian surface height autocorrelation function W ( x (cid:107) ) of the form (4). Realizations of therandomly rough surfaces were generated by the Fourier filtering method as described in Ref. 8. If nothing is said toindicate otherwise, the edges of the square region of the x x -plane covered by the rough surface were L = 21 λ with λ the wavelength of the incident beam, and the half-width of the incident beam, given by the expression in Eq. (24),was w = L/ λ . Moreover, the resulting linear set of equations was solved using the stabilized biconjugatedgradient (BiCGStab) iterative method [21], with the matrix-vector multiplications that it requires performed usingroutines from BLAS (Basic Linear Algebra Subprograms) [22]. A. Isotropic surfaces In the first set of calculations that we performed it was assumed that the rms-roughness of the surface was δ = λ and the correlation length of the isotropic surface was a = 2 λ . For an isotropic surface, one takes a = a ≡ a inEqs. (4) and (5). These roughness parameters are identical to those used by Tran and Maradudin [1] in their initialstudy of the scattering of scalar waves from rough impenetrable surfaces. We will start by assuming a rough Dirichletsurface so it is the integral equation (49), and the corresponding linear set of equations (58), we want to solve. 1. Isotropic Dirichlet surfaces Figure 2(a) shows the in-plane angular dependence of the mean DRC of the scattered wave when a scalar Gaussianbeam is incident normally on the rough surface. This distribution, up to statistical fluctuations, displays a reflectionsymmetry ( θ s ↔ − θ s ) with respect to the angle of scattering θ s = 0 ◦ . Moreover, this distribution is equal tothe corresponding out-of-plane distribution, for a normally incident beam, and the latter distribution is thereforenot shown. When the beam instead is incident on the rough surface at the polar angle of incidence θ = 20 ◦ ,we obtain the in-plane and out-of-plane mean DRCs depicted in Figs. 2(b)–(c), respectively. The vertical dashedand dotted lines in Fig. 2 correspond to the backscattering and specular directions, respectively, and these lines areincluded as a guide to the eye. Well defined peaks in the scattered intensity distributions about the retroreflectiondirections θ s = − θ are observed in Figs. 2(a)–(b). These peaks are enhanced backscattering peaks that result fromthe constructive interference of volume waves that are scattered multiple times in the grooves of the strongly roughsurface [1, 5, 23]. This wave phenomenon was first observed experimentally in the scattering of light from stronglyrough metal surfaces [5]. We therefore stress that the peak located at θ = 0 ◦ in Fig. 2(a), for instance, is not due tocoherent (specular) scattering; in fact, the contribution from coherently scattered waves to the angular integral of themean DRC [see Eq. (63) below] is less than 0 . φ s = φ ± ◦ ] of the mean DRCfor the polar angle of incidence θ = 20 ◦ . The maximum scattered intensity is found around θ s = 0 ◦ , indicated by thevertical dash-dotted line in Fig. 2(c). From this figure we observe that as | θ s | increases, the scattered intensity drops2 R / s (a) (b) (c)0.00.10.20.30.4 R / s (d) (e) (f)-90 -60 -30 0 30 60 90 s [deg]0.00.10.20.30.4 R / s (g) -60 -30 0 30 60 90 s [deg](h) -60 -30 0 30 60 90 s [deg](i) FIG. 2: The in-plane and out-of-plane dependencies of the mean DRCs as functions of the scattering angle θ s forrandomly rough, isotropic, Gaussian correlated Dirichlet surfaces characterized by the surface roughness δ = λ andcorrelation length a = 2 λ [Figs. 2(a–c)], a = 3 λ [Figs. 2(d–f)], and a = 4 λ [Figs. 2(g–i)] when the surfaces areilluminated by Gaussian beams of wavelength λ . These results were obtained by assuming a square region of the x = 0 plane, of edges L = 21 λ , covered by the randomly rough surface that was discretized using the intervals∆ x = λ/ 10 [ a = 2 λ ] and ∆ x = λ/ a > λ ]. The half-width of the Gaussian incident beam was w = L/ λ , thepolar angle of incidence was θ = 0 ◦ or 20 ◦ , and the azimuthal angle of incidence was φ = 0 ◦ . The different panelscorrespond either to the in-plane or out-of-plane dependence of the mean DRCs and to different choices for the polarangle of incidence; the in-plane dependence of the mean DRCs are presented in Figs. 2(a, d, g) and 2(b, e, h) for θ = 0 ◦ and 20 ◦ , respectively, while the out-of plane dependence of the mean DRCs for θ = 20 ◦ depicted inFigs. 2(c, f, i). The results for the mean DRCs that we report were obtained by averaging over an ensembleconsisting of at least N ζ = 4000 realizations of the surface profile function. As a guide to the eye, in each of thein-plane panels we have indicated the backscattering (retroreflection) and the specular directions by vertical dashedand dotted lines, respectively. Moreover, in out-of-plane panels the directions corresponding to θ s = 0 ◦ have beenmarked by vertical dash-dotted lines.off from its maximum value at θ s = 0 ◦ and the distribution is approximately reflection symmetric with respect to thenormal scattering direction; the scattered intensity distribution from an isotropic surface is expected to display suchleft-right symmetry with respect to the plane of incidence.The next set of roughness parameters we consider consists of δ = λ (unchanged) and a = 3 λ . The resulting in-planeand out-of-plane dependencies of the mean DRCs are presented in Figs. 2(d)–(e) and 2(f), respectively, for the polarangles of incidence θ s = 0 ◦ and 20 ◦ . Qualitatively these results show several features that are similar to what isobserved when the correlation length is a = 2 λ [Fig. 2(a)–(c)]. For instance, for both of the polar angles of incidence θ = 0 ◦ and 20 ◦ , enhanced backscattering peaks are observed at θ s = − θ . However, there are also differences betweenthe mean DRCs obtained for the two sets of roughness parameters. For instance, Fig. 2(e) shows that the in-planescattered intensity has a maximum in the forward scattering plane ( θ s > ◦ ); this is not the case when a = 2 λ [Fig. 2(b)] for which the maximum of the in-plane scattered intensity is found in the backscattering direction and3therefore in the backscattering plane ( θ s < ◦ ). For normal incidence, a comparison of the in-plane mean DRCs inFigs. 2(a) and 2(d) reveals that the backscattered intensity is higher when a = 2 λ than what it is when a = 3 λ .A more detailed study of the simulation results in the same figures show that the ratio of the scattered intensityin the retroreflection direction θ s = − θ to the intensity of the background at its position is higher when a = 2 λ then what it is when a = 3 λ . If single-scattering contributions to the scattered intensity can be neglected aroundthe retroreflection direction, the intensity of the backscattering peak is expected to be twice that of the backgroundintensity [19, 23]; this implies that the intensity ratio we defined above should have the value two. From the resultsin Figs. 2(a)–(b) [ a = 2 λ ] and 2(d)–(e) [ a = 3 λ ] it is observed that the intensity ratio is a little smaller than two when a = 2 λ and smaller than this value when a = 3 λ . These findings we take as an indication that multiple scatteringprocesses contribute more significantly to the scattered field for the former set of roughness parameters than for thelatter set of roughness parameters. We find also, for both θ = 0 ◦ and 20 ◦ , that the enhanced backscattering peakis wider in Figs. 2(a)–(b) than in Figs. 2(d)–(e), that is, for the shorter transverse correlation length of the surfaceroughness. This is consistent with what has been observed previously for the scattering of light from one-dimensionalrandomly rough surfaces [19]. Furthermore, the results in Figs. 2(a) and 2(d) show that the entire in-plane scatteredintensity distribution is broader for the case when the correlation length is a = 3 λ than for the case when it is a = 2 λ .Finally, Figs. 2(g)–(i) present computer simulation results for the in-plane or out-of-plane mean DRCs obtainedwhen the correlation length of the rough Dirichlet surface is a = 4 λ with the remaining roughness and numericalparameters having unchanged values; this is twice the value of the correlation length assumed in obtaining the resultspresented in, for instance, Figs. 2(a)–(c). Contrary to what was found when a = 2 λ and a = 3 λ (with δ = λ ), theresults in Figs. 2(g)–(i) for a = 4 λ show no well-defined enhanced backscattering peaks in the scattered intensitydistributions; in this case, the local slopes [23] of the Gaussian surface are simply too small to allow for any significantcontribution from multiple scattering of volume waves. Instead the in-plane and out-of-plane scattered intensitydistributions are found to be featureless around the backscattering and specular directions. Moreover, when θ = 20 ◦ the maximum scattered intensity is located in the forward scattering plane ( θ s > ◦ ), see Fig. 2(h).Based on the angular dependence of the in-plane mean DRCs in Fig. 2 alone, one can strictly speaking not attributethe features observed at θ s = − θ in these distributions to the backscattering phenomenon. To reach this conclusion,one is required to demonstrate that these features are actually peaks in the retroreflection direction and not, forinstance, ridges as were recently observed in the mean DRC for cross-polarized light scattered from a rough perfectelectric conducting surface [7]. The full angular distribution of the mean DRCs for the scattering of scalar waves arepresented in Fig. 3 for polar angles of incidence θ = 0 ◦ , 20 ◦ , 40 ◦ , and for the three sets of roughness parameters forthe Dirichlet surfaces used to produce the results in Fig. 2. Figure 3 shows several examples of well-defined peaksin the mean DRCs about the retroreflection direction [ q (cid:107) = − k (cid:107) ]; in particular, these results demonstrate explicitlythat the peaks at θ s = − θ in Figs. 2(a)–(b) and Figs. 2(d)–(e) are enhanced backscattering peaks. Moreover, thecomments that were made about the results in Fig. 2 regarding the width of the angular distributions are readilyseen to apply to the results presented in Fig. 3; indeed the distributions in Figs. 3(a)–(b) that correspond to thecorrelation length a = 2 λ are less wide than the angular distributions in Figs. 3(a)–(b) [ a = 3 λ ]. However, the mostinteresting observation to be made from the results in Fig. 3 is how the backscattering peaks depend on the polarangle of incidence θ and the correlation length of the surface when the surface roughness is the same. The resultspresented in Figs. 3(a)–(c) correspond to a = 2 λ and display well defined enhanced backscattering peaks for all thepolar angles of incidence θ = 0 ◦ , 20 ◦ , and 40 ◦ . On the other hand, when the correlation length of the surface isincreased from a = 2 λ to a = 3 λ , the mean DRCs in Figs. 3(d)–(e), corresponding to the polar angles of incidence θ = 0 ◦ , 20 ◦ , and 40 ◦ , respectively, show enhanced backscattering peaks that gradually disappear with increasingpolar angles of incidence. For instance, in Fig. 3(f) [ θ = 40 ◦ ] no backscattering peak is observed, while such a peakis observed in the mean DRC for the same polar angle of incidence when a = 2 λ [Fig. 3(c)].At this stage it should be commented that from the results of the full angular dependence of the mean DRC in Fig. 3one can check the satisfaction of the energy conservation; this is often referred to as unitarity . Since the Dirichlet andNeumann surfaces are impenetrable to scalar waves, all energy incident on them has to be reflected away from them.From the definition of the DRCs in Eqs. (39) and (40) it follows that U ( k (cid:107) ) = (cid:90) q (cid:107) <ω/c d Ω s (cid:28) ∂R ( q (cid:107) | k (cid:107) ) ∂ Ω s (cid:29) = 1 . (63)The relation in Eq. (63) is a consequence of energy conservation. Strictly speaking energy conservation should besatisfied individually for each of the DRCs that enters into the calculation of the mean DRC, but this possibilitywe will not explore here. It should be noted, that relation (63) is a necessary but not a sufficient condition forcorrect results. However, we have found that the satisfaction of Eq. (63) is a good criterion for gauging the quality ofsimulation results and, for instance, if the discretization interval used in performing the calculations is small enoughfor the assumed roughness parameters. For the results for normal incidence presented in Fig. 3, and therefore also in4 q / ( / c ) (a) (b) (c)1.00.50.00.51.0 q / ( / c ) (d) (e) (f)1.0 0.5 0.0 0.5 1.0 q /( / c )1.00.50.00.51.0 q / ( / c ) (g) 0.5 0.0 0.5 1.0 q /( / c )(h) 0.5 0.0 0.5 1.0 q /( / c )(i) 0.000.070.140.210.280.350.420.490.560.63 R / s R / s R / s FIG. 3: The full angular dependence of the mean DRCs for randomly rough, isotropic, Gaussian correlated Dirichletsurfaces as functions of the in-plane scattered wave vector q (cid:107) . In all cases the azimuthal angle of incidence was φ = 0 ◦ and the polar angles of incidence were θ = 0 ◦ [Figs. 3(a, d, g)]; 20 ◦ [Figs. 3(b, e, h)]; and 40 ◦ [Figs. 3(c, f, i)].The roughness parameters were those of Fig. 2, that is, for all cases the roughness of the surface was δ = λ , and thecorrelation length was a = 2 λ [Figs. 3(a–c)]; a = 3 λ [Figs. 3(d–f)]; and a = 4 λ [Figs. 3(g–i)]. The remainingscattering and numerical parameters are identical to those presented in the caption of Fig. 2. As a guide to the eye,the positions of the backscattering and the specular directions have been indicated by black crosses and filled circles,respectively.Fig. 2, we found that energy conservation was satisfied to within an error of no more than 2 × − for the numericalparameters [see caption in Fig. 3] assumed in performing these calculations. For non-normal incidence, the error inthe satisfaction of the energy conservation condition was somewhat larger, partly due to the increased footprint ofthe incident beam on the mean surface; for all the calculations reported in Fig. 3 this error was always smaller than10 − which testifies to the quality of the simulation results that we have obtained.At this stage we should mention the computational resources needed to produce the simulation results presentedin Fig. 3. For instance, the results in Fig. 3(a) took 598 s, or almost 10 min, of computer time to produce per surfacerealization when the simulations were performed on a single core of an Intel i7-960 processor (8M Cache, 3 . 20 GHz).This time was spent in the following manner: 385 s for setting up the linear system of equations; 125 s to solve it usingthe BiCGStab iterative solver (for one angle of incidence); 44 s to calculate the reflection amplitudes in Eqs. (20) and(23); and finally another 44 s to obtain the DRC and mean DRC defined by Eqs. (39) and (40). For the simulationsthat we performed, the code was OpenMP parallelized, in which case the wall time for the same simulation wasreduced; using the processor given above and running on all 8 cores resulted in a wall time of 5 . R / s (a) (b) (c)0.00.10.20.30.4 R / s (d) (e) (f)-90 -60 -30 0 30 60 90 s [deg]0.00.10.20.30.4 R / s (g) -60 -30 0 30 60 90 s [deg](h) -60 -30 0 30 60 90 s [deg](i) FIG. 4: Same as Fig. 2 but for Neumann surfaces. 2. Isotropic Neumann surfaces We now address rough Neumann surfaces and the scattering of a scalar beam from them. The parameters of theisotropic randomly rough surfaces we will take to be identical to those assumed for the Dirichlet case [Figs. 2 and 3].The in-plane and out-of-plane angular dependencies of the mean DRCs for rough Neumann surfaces, obtained on thebasis of solving the linear system Eq. (62), are presented in Fig. 4. The corresponding full angular dependencies ofthe mean DRCs are depicted in Fig. 5 where results for θ = 40 ◦ have also been included. At least for the roughnessparameters and polar angles of incidence that we have assumed here, the mean DRCs for both Neumann and Dirichletsurfaces are qualitatively rather similar. For instance, an enhanced backscattering peak is observed around q (cid:107) = − k (cid:107) in the mean DRC for a Neumann surface if it also is observed in the mean DRC for a Dirichlet surface with thesame roughness parameters. Moreover, the behavior we reported above for the mean DRCs of Dirichlet surfaces forincreasing values of a and θ , we also find in the case of scattering from Neumann surfaces. The only minor differencesbetween the mean DRCs for Neumann and Dirichlet surfaces we find worthy of a comment is the difference in shapeof the in-plane and out-of-plane distributions for these two kinds of surfaces. In particular, a detailed comparison ofthe results in Figs. 2(d)–(f) and Figs. 4(d)–(f), corresponding to the scattering from a Dirichlet and Neumann surfaceof the same surface statistics, reveals that the tails and amplitudes of the distributions are somewhat different. Thedistributions in the former case are semi-circular, while in the latter case the distributions are more triangular. 3. Reflectivity For the isotropic surfaces that we have studied until now, the roughness parameters were such that the meanDRCs for both the rough Dirichlet and Neumann surfaces were fully diffuse. In view of Eq. (45) this means that thereflectivity of such surfaces is neglectable and smaller than the error found in the satisfaction of the energy conservationcondition (63). In order to have a non-vanishing reflectivity, we will now consider less rough surfaces; in particular, we6 q / ( / c ) (a) (b) (c)1.00.50.00.51.0 q / ( / c ) (d) (e) (f)1.0 0.5 0.0 0.5 1.0 q /( / c )1.00.50.00.51.0 q / ( / c ) (g) 0.5 0.0 0.5 1.0 q /( / c )(h) 0.5 0.0 0.5 1.0 q /( / c )(i) 0.000.070.140.210.280.350.420.490.560.63 R / s R / s R / s FIG. 5: Same as Fig. 3 but for Neumann surfaces.work with the roughness parameters δ = λ/ 10 and a = λ , still assuming an isotropic Gaussian correlation function.For such roughness parameters the mean DRC was calculated for a sequence of lateral wave vector k (cid:107) = k (cid:107) ˆk (cid:107) ofthe incident Gaussian beam where k (cid:107) = ( ω/c ) sin θ . In particular, the calculation were performed for polar anglesof incidence from θ = 0 ◦ to 80 ◦ in steps of 10 ◦ . For each polar angle of incidence, the reflectivity was calculatedfrom Eq. (45) by using the coherent component of the mean DRC. The mean DRCs were obtained on the basis ofsufficiently many surface realizations so that the calculated reflectivity had converged to at least four decimal places;for the values of the roughness and numerical parameters assumed, 75 surface realizations were sufficient to obtainsuch a convergence. In this way, we obtain the reflectivity curves presented in Fig. 6, where the error bars represent | −U ( k (cid:107) ) | . From the results presented in this figure one observes that the reflectivity of the Dirichlet surface is alwayshigher than the reflectivity of the Neumann surface with the same roughness parameters, and this is even the casefor normal incidence. Moreover, increasing the polar angle of incidence seems to increase the reflectivity; only for thelargest angle of incidence that we consider, θ = 80 ◦ , there may be an exception to this rule in the Neumann case.For such large polar angles of incidence the precision in the simulations (see error bars) is simply not good enough toreach a definite conclusion on this issue. In any case, we remark that when the polar angle of incidence is approaching90 ◦ the reflectiveties of both the Dirichlet and the Neumann surfaces should approach unity. Based on the resultspresented in Fig. 6, it is tempting to speculate that the reflectivity of the Dirichlet surface will smoothly approachunity, while the reflectivity of the Neumann surface will go through a local minimum before again increasing andreaching unity for grazing angles of incidence. Such behavior of the reflectivity was recently observed in the scatteringof s- and p-polarized light from randomly rough perfectly conducting surfaces [9]. The determination of whether ourspeculations are well founded we will leave for a dedicated study that probably will benefit from the use of a planeincident wave.7 θ [deg] R e f l ec ti v it y DirichletNeumann δ = λ/10 a = λ FIG. 6: The reflectivity as a function of the polar angle of incidence calculated on the basis of Eq. (45) for isotropicGaussian correlated Dirichlet and Neumann surfaces of roughness δ = λ/ 10 and correlation length a = λ . Here λ denotes the wavelength of a Gaussian beam of width w = L/ L = 32 λ . Theremaining numerical parameters are identical to those used in obtaining the results of Fig. 3. The error bars, onlyclearly visible for θ = 80 ◦ , indicate the absolute deviation from unity of the unitarity obtained in the numericalsimulations. The reported results were obtain by averaging the results of 50 surface realizations, which was sufficientto achieve convergent results. 4. Comparison to previous results Before leaving the scattering from isotropic surfaces, it should be remarked that more than 20 years ago, Tran andMaradudin published the initial rigorous computational results for the in-plane mean DRC obtained when normallyincident scalar beams are scattered from strongly rough Dirichlet and Neumann surfaces [1, 3]. Their calculations werebased on the integral equations Eqs. (49) and (52), which were transformed into a linear set of equations when eachrealization of the surface profile function ζ ( x (cid:107) ) of edges L = 16 λ was discretized onto a square grid of N = 64 pointsper side (∆ x = λ/ N = 4096) equations was solved iteratively by a methodbased on the Liouville-Neumann series [24]. In the calculation results reported in Refs. [1, 3] an isotropic Gaussiancorrelation function was used, while the values of the roughness parameters, δ = λ and a = 2 λ , assumed were identicalto those assumed in performing the calculations whose results are presented as the first row of sub-figures in Figs. 2–5.The main differences between the results that we report and those reported in Refs. [1, 3] are that we use a largersurface, a wider incident beam, and a smaller discretization interval in performing the calculations, relative to thecorresponding parameters used in the initial calculations. Moreover, the linear equation system is solved in a ratherdifferent manner. A direct comparison of the results from Figs. 2(a) and 4(a) to the results of Figs. 2(b) and 3(b)in Ref. 3 reveals a fair agreement between the former and latter sets of results. For instance, all simulation resultsfor θ = 0 ◦ predict the amplitude of the mean DRC at θ s = 0 ◦ to be somewhat above 0 . 6. However, the results inRef. 3 seem to predict a less broad angular intensity distribution than what we find in the calculations reported here.The observed difference we suspect is caused by the significantly smaller discretization interval used in performingthe calculations reported in Ref. 3. With the use of the numerical parameters assumed in this publication, we wereable to reproduce rather well the results for a Dirichlet surface reported by Tran and Maradudin [3].It ought to be remarked that energy conservation, and thus the quality check on the simulation results that followsfrom it, could not have been performed on the simulation results produced by Tran and Maradudin [1, 3] or byMacaskill and Kachoyan [2]. In none of these publications was the full angular distribution of the scattered intensitycalculated, which is required to obtain U ( k (cid:107) ) defined in Eq. (63). Finally, we remark that in the calculations thatwe performed for θ = 0 ◦ assuming the numerical parameters (with ∆ x = λ/ 4) from Ref. [3], energy conservationwas satisfied within an error of 6 . x = λ/ 10) resulted in an error in the satisfaction of the energyconservation of no more than 0 . ζ ( x (cid:107) ) defined by the roughnessparameters δ = λ/ a = λ and (a) a = λ (isotropic surface); (b) a = 1 . λ ; (c) a = 2 λ and (d) a = 3 λ . Thesurface realizations were generated by the Fourier filtering method [8]. The generation of the surface realizationsassumed the same uncorrelated random numbers so that one can follow how surface features are affected by theincreasing level of anisotropy. B. Anisotropic surfaces Until now we have exclusively been dealing with isotropic surfaces. However, many naturally occurring or manmade surfaces are anisotropic. Therefore, we now turn our attention to the scattering from anisotropic surfaces. Inthis case the correlation function W ( x (cid:107) ) that we assume has the form (4) with a (cid:54) = a . We remark that this form ofanisotropy is not the most general one; for instance, the main axes of the anisotropy do not have to be orthogonal toeach other, but such a more general case will not be addressed here.In the computer simulations that we performed, the polar angle of incidence was θ = 20 ◦ , the surface roughness was δ = λ/ 2, and the correlation length along the x -direction had the constant value a = λ . The correlation length alongthe x -direction was then varied so that a ∈ { λ, . λ, λ, λ } . A realization of the rough surface for each of theseroughness parameters is presented in Fig. 7. To better be able to follow how surface features are affected by increasingthe level of anisotropy, each of these surface realizations was generated on the basis of the same uncorrelated randomnumbers by the method described in Ref. 8. We learn from the surface topographies depicted in Fig. 7 that when a is several times larger than a , the resulting surface topography starts taking the form of a “randomly corrugatediron” structure where the grooves of the structure are along the x -direction. In the limit that a /λ → ∞ , or moreprecisely that a (cid:29) a with a (cid:29) λ , the surface starts approaching a one-dimensional surface for which the surfaceprofile function ζ ( x (cid:107) ) will be independent of the spatial coordinate x . Moreover, the numerical parameters assumedin the simulations and the parameters characterizing the incident beam were identical to those used previously inthe study of the scattering from isotropic surfaces [see captions of Figs. 2 and 3]. The value of the azimuthal angleof incidence, φ , we assumed to be 0 ◦ , 45 ◦ , and 90 ◦ . Notice that when we are dealing with the mean DRC there isno need to consider values of φ outside the interval 0 ◦ to 90 ◦ since the mean DRC for φ / ∈ [0 ◦ , ◦ ] can be relatedby symmetry to the mean DRC for a value of φ ∈ [0 ◦ , ◦ ]. However, note that for the DRC that a single surfacerealization gives rise to, this is no longer true.For later comparison, we start by presenting the full angular dependence of the mean DRC for the isotropic case[ a = a = λ ] for three values of the azimuthal angle of incidence, φ = 0 ◦ , 45 ◦ and 90 ◦ [Figs. 8(a)–(c)]. In eachpanel of Fig. 8 the direction of the plane of incidence has been indicated by a thin black dashed line. Furthermore,in the same figure, the black crosses and filled black circles represent the backscattering and specular directions,respectively. The results presented in Figs. 8(a)–(c) all display well-defined enhanced backscattering peaks, and theseresults are related to each other by azimuthal rotations (up to fluctuations which are caused by the use of a finitenumber of surface realizations); for instance, the mean DRC corresponding to φ = 45 ◦ can be obtained from themean DRC for φ = 0 ◦ (90 ◦ ) by a counter-clockwise rotation through an azimuthal angle of 45 ◦ ( − ◦ ) about the q -axis. This is a consequence of the isotropy of the surface and that the polar angle of incidence is the same for theresults in Figs. 8(a)–(c). It should be mentioned that for the roughness parameters we assume the scattering was fullyincoherent (diffusive); in fact, the reflectivity was of the order of 10 − , which is on a par with the error level in thesesimulations.We now turn to Gaussian correlated anisotropic rough surfaces defined by the correlation lengths a = 1 . λ [and a = λ ]. From the contour plot of the surface realization of such surfaces presented in Fig. 7(b), it is not immediatelyapparent that the surface is anisotropic; to realize this, a close inspection of the topography is needed, and as a result,we will in the following refer to this class of surfaces as being slightly anisotropic. Figures 8(d)–(f) display the angular9FIG. 8: The full angular dependence of the mean DRCs for randomly rough, anisotropic , Gaussian correlatedDirichlet surfaces as functions of the in-plane scattered wave vector q (cid:107) . The surfaces were illuminated by Gaussianincident beams of wavelength λ and angles of incidence ( θ , φ ) were θ = 20 ◦ and φ = 0 ◦ [1st column;Figs. 8(a, d, g, j)], 45 ◦ [2nd column; Figs. 8(b, e, h, k)], and 90 ◦ [3rd column; Figs. 8(c, f, i, l)]. The surface roughness is δ = λ/ a = λ and a = λ [1st row]; a = 1 . λ [2nd row]; a = 2 λ [3rd row]; and finally a = 3 λ [4th row]. The sampling interval usedin performing the numerical simulations was ∆ x = λ/ 8. All remaining numerical parameters and the property of theincident beams are identical to those assumed in obtaining the results in Fig. 3. The reported results were obtainedby averaging over an ensemble consisting of N ζ = 4000 of surface realizations. The thin dashed black lines thatappear in each of the panels indicate the direction of the plane of incidence. The backscattering and speculardirections are marked by black crosses and filled circles. If the projection of the wave vector of incidence is k (cid:107) = ( k , k , 0) then the directions q (cid:107) = ( − k , k , 0) are indicated by black filled triangles in the 2nd column.0dependence of the mean DRCs for the azimuthal angles of incidence φ = 0 ◦ , 45 ◦ and 90 ◦ , respectively. The first thingto notice from these results is that the anisotropy of the surface roughness alters the scattered intensity distributionsrelative to the intensity distributions obtained for the corresponding isotropic surface for which a = λ [Figs.8(a)–(c)].It is also found that the mean DRCs for this anisotropic surface have their highest scattered intensity at, or closeto, the backscattering (retroreflection) direction [black crosses in Fig. 8]. However, the angular dependencies of thescattered intensity around these directions are different for the isotropic and anisotropic surfaces. For instance, bycomparing the mean DRCs in Figs. 8(d) and 8(a) it is observed that what appears as an almost isotropic featurearound the backscattering direction in the isotropic case [Fig. 8(a)], is transformed into an elliptic-like intensitydistribution [Fig. 8(d)] about the same direction that is elongated along the ˆq direction, that is, along the directionfor which the correlation length is the longest. At the same time, the scattered intensity in the backscattering directionis lower in the anisotropic case than in the isotropic case.A similar kind of elongation of the scattered intensity distribution is observed for the anisotropic case when φ = 90 ◦ [Fig. 8(f)]; however, in this case the elongation is predominantly in a positive ˆq -direction away fromthe retroreflection direction and extending all the way to the specular direction. A comparison of the mean DRCs inFigs. 8(f) and 8(c) also reveals that the anisotropy of the surface enhances the scattering into the forward scatteringplane, that is, into directions for which q > 0. It ought to be remarked that in the scattering from isotropic surfaceswe also observed that increasing the correlation length a > λ caused broadening of the mean DRCs as well as enhancedscattering into the forward scattering plane; see Figs. 2 and 3.Figure 8(e) presents the mean DRC for an anisotropic surface for which a = 1 . λ and φ = 45 ◦ . In this case wealso observe an elongation and shift of the scattered intensity distribution around the retroreflection direction intodirection ˆq . The red elliptic-like structures seen in Figs. 8(e)–(f) are rather similar. Moreover, also for the casewhen φ = 45 ◦ [Fig. 8(e)] an enhancement of the scattering into the forward scattering plane is observed; this is quitesimilar to what was found previously when φ = 90 ◦ [Fig. 8(f)].The mean DRCs for the isotropic surface shown in Figs. 8(a)–(c) all have their highest intensity in the backscatteringdirection. For the anisotropic surface of correlation length a = 1 . λ [Figs. 8(d)–(f)] we find that this is also the casewhen φ = 0 ◦ [Fig. 8(d)]. However, for the same roughness parameters and when φ = 45 ◦ and 90 ◦ , a close inspectionof the mean DRCs in Figs. 8(e)–(f) reveals that this is no longer the case; here the directions of maximum intensityare instead shifted away from the retroreflection direction by a small amount in the positive ˆq -direction. Such shiftsof the maximum diffusely scattered intensity are even more apparent when the correlation length a is increasedfurther. Figures 8(g)–(i) and 8(j)–(l) present the mean DRCs corresponding to the correlation lengths a = 2 λ and a = 3 λ , respectively, when the other roughness and numerical parameters are the same. For both these values of a , and when φ = 45 ◦ and 90 ◦ , one finds that the maxima of the scattered intensity are found in the region q > φ = 90 ◦ we find from Figs. 8(i, l) that the directionof maximum scattered intensity is at, or close to, the specular direction. This we speculate is caused by the wavescattering from the grooves of the surface roughness [see Figs. 7(c)–(d)]; such scattering will cause most of the intensitybeing scattered into the forward direction when a > λ , similarly to what we saw previously for the scattering fromisotropic surfaces. Furthermore, when φ = 0 ◦ , we observe from the results in Figs. 8(g, j) that the highest scatteredintensity remains located in the backscattering direction also for these strongly anisotropic surface; it is speculatedthat this will be true when φ = 0 ◦ for any value of a assuming a = λ and δ = λ/ a = 2 λ [Figs. 8(g)–(i)] and a = 3 λ [Figs. 8(j)–(l)], several observations can be make. First, the scattered intensity distributions are found to besignificantly narrower in the ˆq -direction when a = 3 λ than what it is when a = 2 λ . On the other hand, the widthsof the same distributions in the ˆq -direction are only very little affected by the change in the a correlation length.This behavior we attribute to the reduced height difference over a wavelength the surface has in the x directionwhen the correlation length a is several times larger than the wavelength λ . In the extreme limit that a → ∞ ,the mean DRC for φ = 0 ◦ will be proportional to δ ( q ), which has zero width in the ˆq direction, and it is thetransition towards this limit that can be observed in the first column of Fig. 8. Second, the angular dependence of thescattered intensity about the direction of highest scattered intensity is significantly more isotropic for the case when a = 3 λ [Figs. 7(j)–(l)] than what it is for the shorter correlation length a = 2 λ [Figs. 7(g)–(i)]. In this regard, thesituation seen when a = 3 λ resembles more the case of the isotropic surface in Figs. 8(a)–(c); however, the directionsof maximum scattered intensity are rather different in these two cases. Finally, when φ = 45 ◦ it is readily observedfrom the results in Figs. 8(h, k) that the directions where the mean DRCs have their maxima are not in the planeof incidence. This possibility was already alluded to earlier when discussing the result in Fig. 8(e). However, whatdetermines the direction of the highest scattered intensity? A hint towards an explanation is found by observingfrom the results in Fig. 8 that the q coordinate of the directions of the highest scattered intensity is only marginallyaffected, if at all, by the anisotropy of the surface introduced by increasing the correlation length a to values larger1then λ . Moreover, when the correlation length a is significantly larger then a = λ , the incident beam will scatterin a more specular-like fashion from the ridges and grooves of the strongly anisotropic surface. For the roughnessparameters that we assume, the maxima of the mean DRCs for an isotropic surface are in the backscattering directionsdefined by q (cid:107) = − k (cid:107) = ( − k , − k , a > a ≥ λ , predicts that the maximum of the mean DRC will be in the direction q (cid:107) = ( − k , k , a = λ and a = 3 λ that are presented in Figs. 8(j)–(l). A rigorous derivation of the above result for these and other roughnessparameters, we will leave for a dedicated study.One final observation should be made about the angular intensity distributions depicted in Fig. 8. When the planeof incidence contains one of the two main axes of the anisotropy, in our case, ˆq or ˆq , then the mean DRCs areexpected to display a reflection symmetry with respect to the plane of incidence. That this is the case can be observedfrom the mean DRCs presented in the first and last columns of Fig. 8; the error in the satisfaction of the reflectionsymmetry in these results we ascribe to the use of a finite number of surface realization in calculating these meanDRCs. When φ does not equal 0 ◦ or 90 ◦ (and 0 ◦ ≤ φ ≤ ◦ ), this reflection symmetry with respect to the plane ofincidence is not expected to hold, something that can be observed from the panels forming the 2nd column of Fig. 8.Until now not much attention have been given to the amplitudes of the mean DRCs for anisotropic surfaces. Toaddress this and other issues, in Fig. 9 we present the in-plane and out-of-plane angular dependencies of the meanDRCs from Fig. 8. The panels of this figure show the in-plane angular dependencies, except for the 3rd column thatdepicts the out-of-plane dependence. The azimuthal angles of incidence are φ = 0 ◦ [column 1]; φ = 45 ◦ [columns 2and 3]; and φ = 90 ◦ [column 4] while for all cases the polar angle of incidence is θ = 20 ◦ . For the azimuthal angle ofincidence φ = 0 ◦ , the amplitudes of the in-plane mean DRCs do depend on the level of anisotropy. It is observed fromthe results in column 1 of Fig. 8 that the amplitude of the in-plane mean DRCs initially drops with increasing value of a before it starts increasing again when the same parameter is increased further. Such behavior can be understood interms of the observation we did previously regarding the width of the angular dependencies of the mean DRC in theout-of-plane directions when discussing Fig. 8; the widths of these distributions initially increase with increasing a before they gradually decrease with the increase of the same parameters as the full angular distributions of the meanDRCs become centered around the plane of incidence. Note that for φ = 0 ◦ there are well pronounced backscatteringpeaks present in all the in-plane dependencies of the mean DRCs presented in column 1 of Fig. 9.When the azimuthal angle of incidence is φ = 90 ◦ , column 4 of Fig. 9, the dependence of the amplitudes of thein-plane angular distributions on the correlation length a is found to be similar to what was found when φ = 0 ◦ .However, there are also significant differences between the two cases. From column 4 of Fig. 9 it is rather apparent howan increase from unity of the anisotropy ratio a /a causes a shift of the maxima of the in-plane scattered intensitydistributions from the backscattering direction [Fig. 9(d)] to the specular direction [Fig. 9(p)]. For instance, when a /a = 3 the in-plane dependence of the mean DRC [Fig. 9(p)] is already well centered around the specular direction.The reason for this behavior is explained in our discussion of Fig. 8.We now turn to the situation for which φ = 45 ◦ presented in columns 2 and 3 of Fig. 9 for in-plane and out-of-planedistributions, respectively. The amplitudes of the in-plane results for this azimuthal angle of incidence, seem not toincrease with the value of a as we saw previously for the cases φ = 0 ◦ and 90 ◦ . This is partly caused by themaxima of the mean DRC moving out of the plane of incidence, see Fig. 8. However, the most interesting results for φ = 45 ◦ are found for the out-of-plane distributions [Fig. 9, column 3]. It is found that the reflection symmetry ofthis distribution corresponding to an isotropic surface is lost as the anisotropy ratio a /a is increased from unity.For the parameters that we have assumed the maxima of these out-of-plane distributions are found for θ s < ◦ andtheir amplitudes are found to increase with increasing a /a .It should be mentioned that we also did simulations for less rough anisotropic surfaces for which δ = λ/ 10, as inFig. 6, and assuming the same correlation lengths and polar angle of incidence as were used in producing the results inFigs. 8 and 9. The aim was to investigate the dependence of the reflectivity on the azimuthal angle of incidence φ . Itwas found that the variation of the reflectivity with φ was of the order 10 − , which is too small to be significant withthe precision that we have in our simulations, at least, this was the case for the roughness and numerical parametersthat we assumed in performing them. VIII. CONCLUSION In conclusion, we have by numerical simulations studied the scattering of incident Gaussian scalar beams fromisotropic and anisotropic, strongly rough Dirichlet and Neumann surfaces. To this end, we present the rigorous,inhomogenious integral equations for the field (Neumann surface) or the normal derivatives of the field (Dirichletsurface) that are obtained by the use of the Green’s function surface integral method. By a nonperturbative and2 R / s (a) (b) (c) (d)0.00.20.40.6 R / s (e) (f) (g) (h)0.00.20.40.6 R / s (i) (j) (k) (l)-90 -60 -30 0 30 60 90 s [deg]0.00.20.40.6 R / s (m) -60 -30 0 30 60 90 s [deg](n) -60 -30 0 30 60 90 s [deg](o) -60 -30 0 30 60 90 s [deg](p) FIG. 9: Various in-plane and out-of-plane behaviors with the polar scattering angle θ s for the angular dependentmean DRCs presented in Fig. 8. Out-of-plane plots are presented in the 3rd column while the plots in the othercolumns are in-plane plots. The polar angle of incidence assumed is θ = 20 ◦ , while the azimuthal angle of incidenceis φ = 0 ◦ [column 1]; φ = 45 ◦ [columns 2 and 3]; and φ = 45 ◦ [column 4]. The results that appear in a given rowof this figure were obtained from the mean DRCs appearing in the same row of Fig. 8. Therefore the surfaceroughness parameters are a = λ [1st row]; a = 1 . λ [2nd row]; a = 2 λ [3rd row]; a = 3 λ [4th row], and in allcases δ = λ/ a = λ . The remaining numerical parameters are identical to those of Fig. 8. The vertical linesmark the positions of the backscattering direction (dashed lines), specular direction (dotted lines) and the θ s = 0 ◦ direction (dash-dotted lines).purely numerical solution of these integral equations, we obtain the fields scattered from the rough surfaces. For aset of roughness parameters and angles of incidence, we calculate the full angular distribution of the mean DRCs forisotropic, Gaussian correlated Dirichlet and Neumann surfaces. The surface parameters were chosen so that some ofthe results showed enhanced backscattering peaks, which is the hallmark of multiple scattering processes. The qualityof the simulation results was quantified by investigating energy conservation (unitarity), and it was found to be satisfiedwith an error smaller than 2 × − , or better, for the main results presented. We also calculated the dependence ofthe reflectivity on the polar angle of incidence. It was found that for the same parameters characterizing the isotropic3rough surface, the reflectivity of a Dirichlet surface was always higher then the reflectivity of the correspondingNeumann surface independent of the polar angle of incidence.For anisotropic, Gaussian correlated, strongly rough surfaces we calculated the full angular distribution of themean DRCs for both Dirichlet and Neumann surfaces for a given polar angle of incidence [ θ = 20 ◦ ] and threeazimuthal angles of incidence φ = 0 ◦ , 45 ◦ , or 90 ◦ . We found that even for moderate levels of surface anisotropy,as characterized by the ratio of the two correlation lengths, a /a , the full angular distributions of the scatteredintensity (mean DRCs) were affected in a profound manner by the surface anisotropy. For a set of four surfaceroughness parameters, corresponding to an increasing ratio of surface anisotropy, we calculated and presented the fullangular distribution of the mean DRCs. The features of the scattered intensity distributions obtained in this waywere discussed and rationalized. ACKNOWLEDGMENTS