Second-order Perturbation Theory for Scattering from Heterogeneous Rough Surfaces

We propose a model to calculate scattering from inhomogeneous three-dimensional, rough surfaces on top of a stratified medium. The roughness is made up of an ensemble of deposits with various shapes and permittivi-ties whose heights remain small with respect to the wavelength of the incident light. This geometry is encountered in the remote sensing of soil surfaces, or in optics wherever there are contaminated planar components. Starting from a volume-integral equation involving the Green's tensor of the stratified medium, we derive a height-perturbative expansion up to second order. Our formulation, which depends explicitly on the profiles of each deposit and on the Fresnel coefficients of the layered substrate, accounts for double-scattering events and permits an evaluation of depolarization in the plane of incidence. Comparisons with rigorous calculations in the simplified case of two-dimensional geometries are presented. It is shown that the second-order scattering term can be much more important for heterogeneous surfaces than for their homogeneous counterparts.


INTRODUCTION
The problem of scattering from rough surfaces has been the subject of ongoing research over the past several decades because of its important applications in various domains such as environmental remote sensing and thinfilm optical devices.Many rigorous 1 and approximate 2,3 methods have been proposed in the case of rough, homogeneous surfaces.Among them, the small perturbation method [4][5][6] (SPM), also referred to as the Rayleigh-Rice or Rayleigh-Fano theory, is the oldest approximate theory.It is still employed on a regular basis at the lowest two orders, although its domain of validity [7][8][9] is more restricted than that of other perturbative approaches, [10][11][12][13] and even though higher-order, efficient numerical schemes [14][15][16] improve drastically the accuracy and the range of the SPM.
In contrast, scattering from inhomogeneous surfaces, namely, deposits of different materials on a flat substrate, has been less intensively studied.Yet these geometries are often encountered in practical situations: Soil surfaces with clods and rocks 17 or contaminated optical surfaces 18 are common types of inhomogeneous surfaces.Scattering models for such geometries are useful not only in remote sensing and optical probing but also in the interpretation of near-field optical microscopy data, 19 in the study of extinction spectra of metallic nanoparticle monolayers, 20 and in the detection of defects in resonant surface-plasmon sensors. 21Another difficulty arising in many of these applications is that the substrate is often a layered medium.It can be a planar waveguide 20 or a multilayer stack, 22 or it may have a continuous permittivity gradient. 17cattering from inhomogeneous surfaces with stratified substrates is a complex issue, and most available rigorous methods are restricted to geometries possessing an axis of invariance 17,23,24 or periodicity. 25In three dimensions and for large scatterers the problem is still beyond the reach of existing methods. 26,279][30][31][32] Two main approaches have been adopted.The first one consists of deriving a height-perturbative expansion from the surface boundary conditions written under the Rayleigh hypothesis. 28,29,31The other one consists of building a Born series from a volume-integral equation involving the Green's tensor of the stratified medium. 17,30,32][32] Scattering from inhomogeneous surfaces, in general terms, has been investigated essentially when the roughness, due either to contaminants or nanoparticles, can be considered randomly or uniformly distributed identical scatterers. 18,33n the most sophisticated theories, multiple-scattering phenomena are handled as manybody problems by introducing the scattering operators of the particles and the propagator of the surface. 33,34This approach is well suited to geometries involving identical particles but is less convenient when the roughness is made up of different aggregates with various permittivities and shapes.Furthermore, it does not give a relationship between the scattering and the topography as is needed in near-field optics or remote sensing. 19n this paper we present a systematic perturbative expansion of the scattering amplitude for heterogeneous deposits on stratified media.The method, again, will be referred to as the SPM.By starting from a volume-integral equation, we reconcile the Born series with a perturbative development to obtain an analytical formulation of the scattered amplitude to first (SPM1) and second (SPM2) order in height.Thus our formulation permits evaluation of depolarization in the plane of incidence. 17,18This derivation cannot be obtained by a simple refinement of the existing methods addressing the homogeneous case.It is, however, consistent with the latter inasmuch as our expressions reduce to those of the classical Rayleigh-Rice theory 13 in the case of homogeneous surfaces, thereby justifying the same ''SPM'' terminology.Then we compare the SPM with a rigorous method 35 and a classical model based on the Born approximation.It appears that the contribution of the second-order term can be much more important in heterogeneous surface scattering than in its homogeneous counterpart.

HOMOGENEOUS ROUGH DEPOSITS ON LAYERED SUBSTRATE
In this section we evaluate the scattering from a homogeneous deposit of relative permittivity h on a layered substrate (lower medium) with relative permittivity ml (z) and surrounded by air (upper medium) with relative permittivity 1.We choose the right Cartesian coordinate (x ˆ, y ˆ, z ˆ) system with z axis directed upward.Hereafter we use R ϭ r ϩ zz ˆas a point vector and r ϭ xx ˆϩ yy ˆas its horizontal projection.For an arbitrary vector a, the notation a will refer to its norm and a ˆto its direction.The profile of the deposit is defined by the region 0 р z р h(r) ϭ h(x, y), for some nonnegative function h of finite extent (see Fig. 1).The relative permittivity as a function of position is thus given by where ml (z) is a piecewise constant function that describes the multilayer.The deposit is considered a perturbation of a reference geometry defined by its relative permittivity ref : Throughout this paper, a harmonic time dependence exp(Ϫit) is implicitly assumed.The equation satisfied by the electric field is ) where J is the source that creates the incident field and ͔ is the ''potential'' term.The potential V(R) is null everywhere except within the source region D ϭ ͕R:0 Ͻ z Ͻ h(r)͖, where it assumes the constant value V ª K 2 ( h Ϫ 1).Thanks to the introduction of the multilayer, dyadic Green's tensor G(R, RЈ) ϭ G(r Ϫ rЈ, z, zЈ) that is the solution of which satisfies the radiation condition at infinity, we may rewrite Eq. (3) as (5) where E ref is the field that is created by J in the reference geometry.If the incident field is a plane wave with polarization p inc , wave number K ϭ /c, and incident wave vector K 0 Ϫ ϭ k 0 Ϫ q 0 z ˆ, the reference field for z Ͼ 0 is the sum of the incident field plus the field reflected by the multilayer, where K 0 ϩ ϭ k 0 ϩ q 0 z ˆand R is the reflection operator whose expression, given in Appendix B, depends solely on the Fresnel reflection coefficients of the multilayer.For z Ͼ 0 and zЈ Ͼ 0, the multilayer Green's tensor can be given by 3 G ˜͑r, z, zЈ͒ ϭ i where q ϭ (K 2 Ϫ k 2 ) 1/2 , Iq у 0 and P Ϯ (k) ϭ I Ϫ K ˆϮK ˆϮ is the projector on the polarization plane for upward-and downward-going wave vectors K Ϯ ϭ k Ϯ qz ˆ.The Dirac component of the Green's function contributes to Eq. ( 5) in the source region only (where we may have RЈ ϭ R).This leads to two different integral equations according to whether the electric field is estimated inside or outside the source region.For the external field (R D) we may write and the internal field (R D) satisfies where A is defined by

PERTURBATIVE EXPANSION
We will now seek the electric field in the form of a systematic expansion where the term E n is of order h n in height.For this, we combine the basic integral equations ( 8) and ( 9) to get the following (exact) expression of the electric field: Note that the above equation could be iterated further.However, another iteration would give terms that are of order O(h 3 ) because of the occurrence of a triple integral ͐ 0 h .Thus the first iteration is sufficient to derive the second-order expansion, which is obtained on replacing E int in the last integral by the zeroth-order internal field.Similarly, the approximation obtained after n Ϫ 1 iterations of the integral equation is complete in order h n .
Zeroth order.At order zero in height we have clearly which is the total field above the multilayer in the absence of roughness.

First order (SPM1).
The unique contribution to first order in height is given by the first iteration term in Eq. (11).When the observation point z is above the maximal excursion of the surface, G ˜(r Ϫ rЈ, z, zЈ) and E ref (rЈ, zЈ) are differentiable functions of zЈ in the source region that can be expanded about zЈ ϭ 0 ϩ .Thus Introducing the Fourier transform of the roughness, (14) and using the spectral representation of Eqs. ( 16) and ( 17) we obtain, for z Ͼ max h, with Here we have introduced the tensors Second order (SPM2).The second-order field can be written in the form E 2 (r, z) ϭ I ϩ J, with where the notation ͓ • ͔ n stands for the nth order in height of a functional of h.I corresponds to the secondorder contribution of the first iteration, while J stems from the second iteration and describes double-scattering phenomena.The calculation of the second-order field is more complicated than that of the first-order field and is given in Appendix A. In particular, special care must be taken to handle the discontinuity of the Green's tensor on the diagonal zЈ ϭ zЉ.We obtain with Above the maximal excursion of the surface, the scattered field can be written as a superposition of upward-going plane waves: where the dyad S(k, k 0 ) is the scattering operator.The perturbative expansion S ϭ S 0 ϩ S 1 ϩ S 2 ϩ ¯for the scattering operator is an immediate consequence of the previous results: The expanded expressions of the tensors S 1 and S 2 in the canonical polarization basis are given in Appendix B.
The particular case of a rough interface between two homogeneous media is realized with ml (z) ϭ h ϭ , and this can be used as a consistency check.The corresponding formulas actually coincide with those of the classical Rayleigh-Rice perturbation theory to first order. 2,3An analytical identification of the second-order kernel B 2 (k, k 0 , ) is very difficult for it possesses several equivalent expressions.Indeed, any dyad of the form ) leaves the integral of Eq. ( 24) unchanged.However, the numerical comparisons that have been performed (for instance, the ''SPM2'' curves in Figs. 3 and 4) show a perfect agreement with the Rayleigh-Rice, second-order formulas. 13

HETEROGENEOUS ROUGH DEPOSITS
We now consider the case of surface roughness on a substrate made of different materials, as exemplified in Fig. 2. The geometry is defined by a positive function where the domains D i for which h i is non-null are isolated (i.e., not connected) for all i ϭ 1, N. The dielectric constant for z у 0 is given by The zeroth-order scattering operator is unchanged, as it corresponds to reflection by the multilayer, while the first-order scattering operator of the inhomogeneous surface is easily found as the sum of the first-order scattering operators S 1 (k, k 0 , h i , h i ) of each of the homogeneous surfaces h i deposited alone on the substrate: 2. Heterogeneous rough deposit on a multilayer substrate.
Proceeding as before, we obtain the second-order scattering operator as the sum of that of each homogeneous surface taken separately plus a coupling term.On introducing the potential V i associated with h i , V i ϭ ( h i Ϫ 1)K 2 , and the tensor A i where is replaced by h i in Eq. ( 10) we have where It is worth comparing the above result with the Born approximation, which is often resorted to in the heterogeneous case.The Born approximation consists of replacing the unknown internal field in the integral of Eq. ( 8) by the solution of the unperturbed problem, i.e., the reference field Comparing Eq. ( 28) with Eq. ( 16), we see that the Born approximation is consistent with the SPM inasmuch as it yields the same result at first order in height.The second-order term stemming from the Born approximation is, however, incomplete, since it is missing the contribution J of Eq. (18b) of the second iteration that accounts for double-scattering events.Figures 6 and 7 below show how this deficiency can result in a drastic deterioration of accuracy as compared with SPM2.The Born approximation is sometimes alternatively defined 19,36 as as suggested by the form of Eq. ( 8).This has little consequence for small contrast (A Ӎ I), but leads to an important difference for larger contrasts.In particular, the first perturbative order is not recovered in the limit of small heights.

NUMERICAL APPLICATIONS
The method has been numerically tested for onedimensional surfaces (i.e., invariant along the y axis) because this allows a comparison to be made with a volumeintegral rigorous method. 35A rigorous computation of the scattered field requires a tight sampling of the inhomogeneous region and a matrix inversion for the corresponding number of unknowns.Therefore a rigorous numerical solution of scattering from rough heterogeneous surfaces is still out of reach of current computer capacities. 26In such cases, it is important to have a re-liable and tractable analytical approximate method at hand, at least for surfaces of small height.In addition, an analytical formulation such as SPM allows for statistical expressions of the cross section in the case of randomly rough surfaces.This keeps us from having to resort to Monte Carlo averages, which are extremely timeconsuming.
The test geometry consists of two trapezoidal rods on a homogeneous half-plane (permittivity ) as depicted in Fig. 2. The rods have identical dimensions (height ϭ /50, lower basis ϭ /2, upper basis ϭ 0.275, separation ϭ /2) but different permittivity h 1 and h 2 .We have chosen three typical values for the permittivity: 2.25 (glass), 4 (silicon), and Ϫ3 ϩ 0.8i (metal).We have plotted the scattering bistatic cross section q 2 ͉S ii (k, k 0 )͉ 2 in s (i ϭ 2) and p (i ϭ 1) polarization for an incidence angle of 20°.In all plots we present the results given by SPM1 (dashed curve) SPM2 (solid curve), Born approximation (squares) and the rigorous simulations (dots).We recall that Born approximation accounts for the variations of the reference field inside the defects but involves only single-scattering mechanisms (the Green's tensor appears once in its expression).The difference between the Born approximation and SPM2 gives insight into the presence of multiple-scattering phenomena.
In Figs. 3 and 4 the rough surface is homogeneous and dielectric or metallic, respectively.In both cases, SPM1 is quite accurate whatever the polarization state.This is in agreement with the numerous studies on the validity of SPM1 for homogeneous surfaces. 8,9Yet the scattering mechanisms are quite different in these two cases.When the surface is dielectric, the h 2 term stemming from the single-scattering expression I and that stemming from the double-scattering expression J are negligible compared with the first-order term.As a result, SPM2 and the Born approximation are similar to SPM1.On the contrary, when the surface is metallic, the contribution due to the field variation I is no longer negligible with respect to SPM1 but interferes destructively with that due to multiple scattering J. Hence the Born ap-Fig.3. Scattering cross section at 20°incidence for the geometry of Fig. 2 with h 1 ϭ h 2 ϭ ml (z) ϭ 2.25 (glass on glass).proximation, which does not account for double scattering, differs significantly from SPM1 and SPM2.
In Fig. 5 the rods are metallic and deposited on glass.In s polarization, multiple scattering is important, making SPM2 more accurate than SPM1 and Born.In p polarization, only single scattering occurs and the three approximations (SPM1, SPM2, and Born) are equally accurate.The absence of multiple scattering in p polarization is probably due to the strong anisotropy of the field radiated by the induced dipoles in the rods, which diminishes the coupling between the two defects.
In Fig. 6 the rods are made of glass and deposited on a metal surface.We observe that SPM1 is inaccurate and that the Born approximation significantly improves the results.Hence it appears necessary to account for the field variations inside the rods to evaluate correctly the diffracted amplitudes.This is not surprising since the tangential field of reference is almost null at the interface and increases as one moves away from the metallic surface.Yet it is also important to account for some multiple scattering, as with SPM2, to get a satisfactory evaluation.Note that double scattering is more important in p polarization than in s polarization.This is possibly due to the excitation of surface plasmons.In Fig. 7 the rods are dielectric with different permittivity and are deposited on a metallic surface, while Fig. 8 is an example of a heterogeneous deposit on a multilayer (actually a bilayer) substrate.The same observations as those noted for Fig. 5 apply for these more complex cases.It is worth noting that, in the heterogeneous case, the importance of multiple scattering and of field variations inside the rods make the use of SPM2 mandatory even for very small deposits whose height is Ϸ/50.Fig. 4. Same as Fig. 3 with h 1 ϭ h 2 ϭ ml (z) ϭ Ϫ3 ϩ 0.8i (metal on metal).Fig. 5. Same as Fig. 3 with h 1 ϭ h 2 ϭ Ϫ3 ϩ 0.8i, ml (z) ϭ 2.25 (metal on glass).Fig. 6.Same as Fig. 3 with h 1 ϭ h 2 ϭ 2.25, ml (z) ϭ Ϫ3 ϩ 0.8i (glass on metal).Fig. 7. Same as Fig. 3 with h 1 ϭ 4, h 2 ϭ 2.25, ml (z) ϭ Ϫ3 ϩ 0.8i (silicon ϩ glass on metal).

CONCLUSION
We have presented a perturbative expression up to second order in height of the field scattered by three-dimensional inhomogeneous deposits on a layered media.Our expression, which was obtained from a volume-integral equation, depends explicitly on the Fresnel reflection factors of the multilayer and on the profiles of the deposits.This perturbation theory is original in that it cannot be obtained by an adaptation of the usual Rayleigh-Rice procedure to the heterogeneous case.It is, however, consistent with the latter in the homogeneous case.We have compared our model to rigorous calculations and to the Born approximation for various geometries.We have observed that the applicability domain of first-order perturbation theory is more restricted in the inhomogeneous case than in the homogeneous one.In our examples, the contribution of multiple scattering is not negligible even when the height of the defects is Ϸ/50.Our formulation permits evaluation of the influence of the coupling between the deposits in the far-field amplitude.

APPENDIX A: DERIVATION OF THE SECOND-ORDER FIELD
Let us proceed to the explicit calculation of I and J from Eqs. (18a) and (18b).By the same argument as previously (with z outside the source region) we have The second term J requires more caution since the Green's tensor is discontinuous on the diagonal zЈ ϭ zЉ.However, we can decompose the latter into a symmetric and an antisymmetric part, (A3b) Note that the antisymmetric part of the Green's tensor does not depend on the permittivity and assumes the simple expression as can be seen from the identity We use J s and J as , respectively, to denote the two integrals resulting from this decomposition, with an obvious notation, so that J ϭ J s ϩ J as .The symmetric part of the Green's tensor is unambiguously defined at zЈ ϭ zЉ ϭ 0 ϩ and thus We now employ ⌿(rЈ, zЈ) to denote the innermost integral in J as to first order in h as An elementary calculation then gives ⌿͑rЈ, zЈ͒ ϭ where (rЉ, zЈ) ª min͕zЈ,h(rЉ)͖ Ϫ h(rЉ)/2, yielding the following expression for J as : Resorting to the spectral representation in Eq. ( 7) of the Green's tensor we obtain the following final expression for the scattered field to the second order in height E 2 (r, z) (z D):

APPENDIX B: EXPANDED EXPRESSIONS THE CANONICAL POLARIZATION BASIS
The vectors defining the fundamental polarization basis (V, H) are given by The operators R and P Ϯ assume the following simple expressions in this basis: where the r j (k) are the Fresnel reflection coefficients of the multilayer in different polarizations.We refer to, e.g., Refs.37 and 38 for the expressions of these coefficients for various multilayer configurations, and we state explicitly the relations in two important particular cases, namely, the semi-infinite monolayer and bilayer.For a homogeneous half-space z Ͻ 0 with relative permittivity we have with qЈ ϭ (K 2 Ϫ k 2 ) 1/2 , where the complex square root with nonnegative, imaginary part has to be taken.For a bilayer ml (z) ϭ 1 for 0 Ͼ z Ͼ ϪL, ml (z) ϭ 2 for z Ͻ ϪL we have where r j (01) (k) [resp.r j (12) (k)] is the Fresnel reflection coefficient from the vacuum [resp.medium 1 ] to a homogeneous medium 1 [resp.2 ], and Note that in some textbooks the Fresnel coefficients in V polarization are defined on the magnetic field, which implies an extra 1/ͱ factor.The kernels B 1 and B 2 can also be decomposed in the canonical polarization basis as where the coefficients (B n Ϯ ) ij are obtained by performing the scalar products p j (k) • ͓B n (k, k 0 )p i (k 0 )͔.Long but straightforward calculations lead to the following expressions for the different kernels.Note, however, that the matricial formulations of Eqs. ( 16), (20), and ( 27) can be advantageously used for a numerical implementation.
First-order kernel.
Second-order kernels for heterogeneous deposits.
In these formulas it is understood that the ␤ ij of the homogeneous case have to be taken with appropriate replacement of the h in the order of appearance.