Second-order perturbation theory for scattering from
heterogeneous rough surfaces.

Charles-Antoine Guérin, Anne Sentenac

► To cite this version:

Charles-Antoine Guérin, Anne Sentenac. Second-order perturbation theory for scattering from heterogeneous rough surfaces. Journal of the Optical Society of America. A Optics, Image Science, and Vision, 2004, 21 (7), pp.1251-1260. 10.1364/JOSAA.21.001251 . hal-00083195

HAL Id: hal-00083195

https://hal.science/hal-00083195v1

Submitted on 10 Jul 2016

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Second-order perturbation theory for scattering from heterogeneous rough surfaces

Charles-Antoine Guérin and Anne Sentenac

Institut Fresnel, Unité Mixte de Recherche 6133, Centre National de la Recherche Scientifique, Faculté des Sciences de Saint-Jérôme, case 162, F-13397 Marseille Cedex 20, France

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 permittivities 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.

## 1. 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 thin-film optical devices. Many rigorous¹ and approximate²,³ methods have been proposed in the case of rough, homogeneous surfaces. Among them, the small perturbation method⁴⁻⁶ (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⁷⁻⁹ is more restricted than that of other perturbative approaches,¹⁰⁻¹³ and even though higher-order, efficient numerical schemes¹⁴⁻¹⁶ 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¹⁷ or contaminated optical surfaces¹⁸ 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,¹⁹ in the study of extinction spectra of metallic nanoparticle monolayers,²⁰ and in the detection of defects in resonant surface-plasmon sensors.²¹ Another difficulty arising in many of these applications is that the substrate is often a layered medium. It can be a planar waveguide²⁰ or a multilayer stack,²² or it may have a continuous permittivity gradient.¹⁷

Scattering from inhomogeneous surfaces with stratified substrates is a complex issue, and most available rigorous methods are restricted to geometries possessing an axis of

invariance¹⁷,²³,²⁴ or periodicity.²⁵ In three dimensions and for large scatterers the problem is still beyond the reach of existing methods.²⁶,²⁷

On the other hand, scattering from homogeneous rough surfaces placed on layered substrates has been tackled by several approximate models.¹⁷,²²,²⁸⁻³² 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.²⁸,²⁹,³¹ The other one consists of building a Born series from a volume-integral equation involving the Green's tensor of the stratified medium.¹⁷,³⁰,³² Except for some models restricted to stratified substrates with two or three layers at most,²⁸,²⁹ the majority of techniques account only for single-scattering processes on the rough boundaries.²²,³⁰⁻³²

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.¹⁸,³³ In the most sophisticated theories, multiple-scattering phenomena are handled as many-body problems by introducing the scattering operators of the particles and the propagator of the surface.³³,³⁴ This 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.¹⁹

In 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.¹⁷,¹⁸ This 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¹³ in the case of homogeneous surfaces, thereby justifying the same “SPM” terminology. Then we compare the SPM with a rigorous method³⁵ 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.

## 2. HOMOGENEOUS ROUGH DEPOSITS ON LAYERED SUBSTRATE

In this section we evaluate the scattering from a homogeneous deposit of relative permittivity $\epsilon_h$ on a layered substrate (lower medium) with relative permittivity $\epsilon_{ml}(z)$ and surrounded by air (upper medium) with relative permittivity 1. We choose the right Cartesian coordinate $(\hat{\mathbf{x}}, \hat{\mathbf{y}}, \hat{\mathbf{z}})$ system with $z$ axis directed upward. Hereafter we use $\mathbf{R} = \mathbf{r} + z\hat{\mathbf{z}}$ as a point vector and $\mathbf{r} = x\hat{\mathbf{x}} + y\hat{\mathbf{y}}$ as its horizontal projection. For an arbitrary vector $\mathbf{a}$, the notation $a$ will refer to its norm and $\hat{\mathbf{a}}$ to its direction. The profile of the deposit is defined by the region $0 \le z \le h(\mathbf{r}) = h(x, y)$, for some nonnegative function $h$ of finite extent (see Fig. 1). The relative permittivity $\epsilon$ as a function of position is thus given by

$$ \epsilon(\mathbf{R}) = \begin{cases} 1 & \text{if } z > h(\mathbf{r}) \\ \epsilon_h & \text{if } 0 \le z \le h(\mathbf{r}), \\ \epsilon_{ml}(z) & \text{if } z < 0 \end{cases} \quad (1) $$

where $\epsilon_{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 $\epsilon^{\text{ref}}$.

$$ \epsilon^{\text{ref}}(\mathbf{R}) = \epsilon^{\text{ref}}(z) = \begin{cases} 1 & \text{if } z > 0 \\ \epsilon_{ml}(z) & \text{if } z < 0 \end{cases} \quad (2) $$

Throughout this paper, a harmonic time dependence $\exp(-i\omega t)$ is implicitly assumed. The equation satisfied by the electric field is

$$ \nabla \times \nabla \times \mathbf{E}(\mathbf{R}) - K^2 \epsilon^{\text{ref}}(z) \mathbf{E}(\mathbf{R}) = \mathbf{J}(\mathbf{R}) + V(\mathbf{R}) \mathbf{E}(\mathbf{R}), \quad (3) $$

where $\mathbf{J}$ is the source that creates the incident field and $V(\mathbf{R}) = K^2[\epsilon(\mathbf{R}) - \epsilon^{\text{ref}}(z)]$ is the “potential” term. The potential $V(\mathbf{R})$ is null everywhere except within the source region $D = \{\mathbf{R}: 0 < z < h(\mathbf{r})\}$, where it assumes the constant value $V := K^2(\epsilon_h - 1)$. Thanks to the introduction of the multilayer, dyadic Green’s tensor $G(\mathbf{R}, \mathbf{R}') = G(\mathbf{r} - \mathbf{r}', z, z')$ that is the solution of

$$ \nabla_{\mathbf{R}} \times \nabla_{\mathbf{R}} \times G(\mathbf{R}, \mathbf{R}') - K^2 \epsilon^{\text{ref}}(z) G(\mathbf{R}, \mathbf{R}') = I \delta(\mathbf{R} - \mathbf{R}'), \quad (4) $$

which satisfies the radiation condition at infinity, we may rewrite Eq. (3) as

$$ \mathbf{E}(\mathbf{R}) = \mathbf{E}^{\text{ref}}(\mathbf{r}, z) + \int d\mathbf{R}' G(\mathbf{R}, \mathbf{R}') V(\mathbf{R}') \mathbf{E}(\mathbf{R}'), \quad (5) $$

where $\mathbf{E}^{\text{ref}}$ is the field that is created by $\mathbf{J}$ in the reference geometry. If the incident field is a plane wave with polarization $\mathbf{p}_{\text{inc}}$, wave number $K = \omega/c$, and incident wave vector $\mathbf{K}_0 = \mathbf{k}_0 - q_0\hat{\mathbf{z}}$, the reference field for $z > 0$ is the sum of the incident field plus the field reflected by the multilayer,

$$ \mathbf{E}^{\text{ref}}(\mathbf{R}) = \exp(i\mathbf{K}_0^{-} \cdot \mathbf{R})\mathbf{p}_{\text{inc}} + \exp(i\mathbf{K}_0^{+} \cdot \mathbf{R})\mathbf{R}(\mathbf{k}_0)\mathbf{p}_{\text{inc}}, \quad (6) $$

where $\mathbf{K}_0^+ = \mathbf{k}_0 + q_0\hat{\mathbf{z}}$ and $\mathbf{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³

$$ G(\mathbf{R}, \mathbf{R}') = -K^{-2} \hat{\mathbf{z}} \hat{\mathbf{z}} \delta(\mathbf{R} - \mathbf{R}') + \tilde{G}(\mathbf{r} - \mathbf{r}', z, z') $$

with

<figure><img src="image_2.png" /><figcaption>Fig. 1. Homogeneous rough deposit on a multilayer substrate.</figcaption></figure>
$$
\begin{equation}
\begin{split}
\tilde{\mathcal{G}}(\mathbf{r}, z, z') = \frac{i}{8\pi^2} \int_{\mathbb{R}^2} \frac{d\mathbf{k}}{q} \exp(i\mathbf{k}\cdot\mathbf{r}) & \{\exp[iq(z+z')]\mathbb{R}(\mathbf{k}) \\
& + \exp(iq|z-z'|)\mathbb{P}^{\text{sign}}(z-z')(\mathbf{k})\},
\end{split}
\tag{7}
\end{equation}
$$

where $q = (K^2 - k^2)^{1/2}$, $\Im q \ge 0$ and $P^\pm(\mathbf{k}) = I - \hat{\mathbf{K}}^\pm \hat{\mathbf{K}}^\pm$
is the projector on the polarization plane for upward- and downward-going wave vectors $\mathbf{K}^\pm = \mathbf{k} \pm q\hat{\mathbf{z}}$. The Dirac component of the Green's function contributes to Eq. (5) in the source region only (where we may have $\mathbf{R}' = \mathbf{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 ($\mathbf{R} \notin \mathcal{D}$) we may write

$$
\mathbf{E}^{\text{ext}}(\mathbf{R}) = \mathbf{E}^{\text{ref}}(\mathbf{R}) + \int d\mathbf{R}' \tilde{\mathcal{G}}(\mathbf{R}, \mathbf{R}') V(\mathbf{R}') \mathbf{E}^{\text{int}}(\mathbf{R}'), \quad (8)
$$

and the internal field ($\mathbf{R} \in \mathcal{D}$) satisfies

$$
\mathbf{E}^{\text{int}}(\mathbf{R}) = A \mathbf{E}^{\text{ref}}(\mathbf{R}) + \int d\mathbf{R}' A \tilde{\mathcal{G}}(\mathbf{R}, \mathbf{R}') V(\mathbf{R}') \mathbf{E}^{\text{int}}(\mathbf{R}'), \quad (9)
$$

where A is defined by

$$
A = I + (1/\epsilon_h - 1)\hat{z}\hat{z}. \quad (10)
$$

3. PERTURBATIVE EXPANSION

We will now seek the electric field in the form of a system-
atic expansion **E**(**R**) = **E**<sub>0</sub>(**R**) + **E**<sub>1</sub>(**R**) + **E**<sub>2</sub>(**R**) + ...,
where the term **E**<sub>n</sub> is of order **h**<sup>n</sup> in height. For this, we
combine the basic integral equations (8) and (9) to get the
following (exact) expression of the electric field:

$$
\begin{equation}
\begin{split}
\mathbf{E}^{\text{ext}}(\mathbf{r}, z) ={}& \mathbf{E}^{\text{ref}}(\mathbf{r}, z) + V \int d\mathbf{r}' \int_0^{h(\mathbf{r}')} dz' \tilde{\mathcal{G}} \\
& \times (\mathbf{r} - \mathbf{r}', z, z') A \mathbf{E}^{\text{ref}}(\mathbf{r}', z') \\
& + V^2 \int d\mathbf{r}' d\mathbf{r}'' \int_0^{h(\mathbf{r}')} dz' \int_0^{h(\mathbf{r}'')} dz'' \tilde{\mathcal{G}} \\
& \times (\mathbf{r} - \mathbf{r}', z, z'') A \tilde{\mathcal{G}}(\mathbf{r}' - \mathbf{r}'', z', z'') \\
& \times \mathbf{E}^{\text{int}}(\mathbf{r}'', z'').
\end{split}
\tag{11}
\end{equation}
$$

Note that the above equation could be iterated further.
However, another iteration would give terms that are of
order O(h³) because of the occurrence of a triple integral
∫₀^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 itera-
tions of the integral equation is complete in order h^n.

Zeroth order. At order zero in height we have clearly

$$
\mathbf{E}_0(\mathbf{R}) = \mathbf{E}^{\text{ref}}(\mathbf{R}), \quad (12)
$$

which is the total field above the multilayer in the ab-
sence 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, $\tilde{\mathcal{G}}(\mathbf{r} - \mathbf{r}', z, z')$ and $\mathbf{E}^{\text{ref}}(\mathbf{r}', z')$

are differentiable functions of z' in the source region that
can be expanded about z' = 0⁺. Thus

$$
\begin{equation}
\begin{split}
\mathbf{E}_1(\mathbf{r}, z) = V & \int d\mathbf{r}' h(\mathbf{r}') \\
& \times \tilde{\mathcal{G}}(\mathbf{r} - \mathbf{r}', z, 0^+) A \mathbf{E}^{\text{ref}}(\mathbf{r}', 0^+), \quad z \notin \mathcal{D}.
\end{split}
\tag{13}
\end{equation}
$$

Introducing the Fourier transform of the roughness,

$$
\hat{h}(\mathbf{k}) = \frac{1}{4\pi^2} \int \exp(-i\mathbf{k} \cdot \mathbf{r}) h(\mathbf{r}) d\mathbf{r}, \quad (14)
$$

and using the spectral representation of Eqs. (16) and (17)
we obtain, for $z > \max h$,

$$
\begin{equation}
\begin{split}
\mathbf{E}_1(\mathbf{r}, z) = \int d\mathbf{k} \hat{h}(\mathbf{k} - \mathbf{k}_0) & \\
\times \exp(i\mathbf{k} \cdot \mathbf{r} + iqz) B_1(\mathbf{k}, \mathbf{k}_0) \mathbf{p}_{\text{inc}}(\mathbf{k}_0),
\end{split}
\tag{15}
\end{equation}
$$

with

$$
B_1(\mathbf{k}, \mathbf{k}_0) = -\frac{V}{2iq} C^+(\mathbf{k}) A C^-(\mathbf{k}_0). \quad (16)
$$

Here we have introduced the tensors

$$
C^{\pm}(\mathbf{k}) = \mathbb{R}(\mathbf{k}) + \mathbf{P}^{\pm}(\mathbf{k}). \quad (17)
$$

Second order (SPM2). The second-order field can be
written in the form **E**₂(**r**, *z*) = *I* + *J*, with

$$
I = V \left[ \int d\mathbf{r}' \int_0^{h(\mathbf{r}')} dz' \tilde{\mathcal{G}}(\mathbf{r} - \mathbf{r}', z, z') A \mathbf{E}^{\text{ref}}(\mathbf{r}', z') \right]_2, \quad (18a)
$$

$$
J = V^2 \left[ \int d\mathbf{r}' d\mathbf{r}'' \int_0^{h(\mathbf{r}')} dz' \int_0^{h(\mathbf{r}'')} dz'' \right. \\
\times \tilde{\mathcal{G}}(\mathbf{r} - \mathbf{r}', z, z') A \tilde{\mathcal{G}}(\mathbf{r}' - \mathbf{r}'', z', z'') \\
\left. \times A \mathbf{E}^{\text{ref}}(\mathbf{r}', z'') \right]_2, \qquad (18\text{b})
$$

where the notation $[ \cdot ]_n$ stands for the $n$th order in
height of a functional of $h$. $I$ corresponds to the second-
order 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

$$
\begin{equation}
\begin{split}
\mathbf{E}_2(\mathbf{r}, z) = \int d\mathbf{k} d\xi \hat{h}(\mathbf{k} - \xi) \hat{h}(\xi - \mathbf{k}_0) \exp(i\mathbf{k} \cdot \mathbf{r} + iqz) \\
\times B_2(\mathbf{k}, \mathbf{k}_0, \xi) \mathbf{p}_{\text{inc}}(\mathbf{k}_0),
\end{split}
\tag{19}
\end{equation}
$$

with
$$
\begin{equation}
\begin{aligned}
B_2(\mathbf{k}, \mathbf{k}_0, \xi) = & -\frac{V}{4q} \{ q_0 C^+(\mathbf{k}) A[\mathrm{R}(\mathbf{k}_0) - \mathrm{P}^-(\mathbf{k}_0)] \\
& + q[\mathrm{R}(\mathbf{k}) - \mathrm{P}^+(\mathbf{k})] \mathrm{AC}^-(\mathbf{k}_0) \} \\
& - \frac{V^2}{8qK^2} C^+(\mathbf{k}) A[(\mathbf{k} - \mathbf{k}_0)\hat{\mathbf{z}} \\
& + \hat{\mathbf{z}}(\mathbf{k} - \mathbf{k}_0)] \mathrm{AC}^-(\mathbf{k}_0) \\
& - \frac{V^2}{8qq_\xi} C^+(\mathbf{k}) A[C^+(\xi) + C^-(\xi)] \mathrm{AC}^-(\mathbf{k}_0).
\end{aligned}
\tag{20}
\end{equation}
$$

Above the maximal excursion of the surface, the scattered field can be written as a superposition of upward-going plane waves:

$$
\mathbf{E}(\mathbf{r}, z) - \mathbf{E}^{\text{inc}}(\mathbf{r}, z) = \int d\mathbf{k} \exp(i\mathbf{k} \cdot \mathbf{r} + i\mathbf{q}z) \\
\times S(\mathbf{k}, \mathbf{k}_0)\mathbf{p}_{\text{inc}}(\mathbf{k}_0), \quad (21)
$$

where the dyad $S(\mathbf{k}, \mathbf{k}_0)$ is the scattering operator. The perturbative expansion $S = S_0 + S_1 + S_2 + \dots$ for the scattering operator is an immediate consequence of the previous results:

$$
S_0(\mathbf{k}, \mathbf{k}_0) = \mathbb{R}_0(\mathbf{k}_0) \delta(\mathbf{k} - \mathbf{k}_0), \quad (22)
$$

$$
S_1(\mathbf{k}, \mathbf{k}_0) = \hat{h}(\mathbf{k} - \mathbf{k}_0) B_1(\mathbf{k}, \mathbf{k}_0), \quad (23)
$$

$$
S_2(\mathbf{k}, \mathbf{k}_0) = \int d\xi \hat{h}(\mathbf{k} - \xi) \hat{h}(\xi - \mathbf{k}_0) \\
\times B_2(\mathbf{k}, \mathbf{k}_0, \xi). \qquad (24)
$$

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 $\varepsilon_{ml}(z) = \varepsilon_h = \varepsilon$, 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.²,³ An

analytical identification of the second-order kernel $B_2(\mathbf{k}, \mathbf{k}_0, \xi)$ is very difficult for it possesses several equivalent expressions. Indeed, any dyad of the form $\mathbb{L}(\mathbf{k}, \mathbf{k}_0, \mathbf{k} - \xi) - \mathbb{L}(\mathbf{k}, \mathbf{k}_0, \xi - \mathbf{k}_0)$ added to $B_2(\mathbf{k}, \mathbf{k}_0, \xi)$ 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.¹³

4. HETEROGENEOUS ROUGH DEPOSITS

We now consider the case of surface roughness on a sub-
strate made of different materials, as exemplified in Fig.
2. The geometry is defined by a positive function

$$
h(\mathbf{r}) = \sum_{i=1}^{N} h_i(\mathbf{r}), \quad (25)
$$

where the domains $\mathcal{D}_i$ for which $h_i$ is non-null are isolated
(i.e., not connected) for all $i = 1, N$. The dielectric con-
stant for $z \ge 0$ is given by

$$
\epsilon(\mathbf{R})
$$

$$
= \begin{cases}
\epsilon_{h_i} & \text{if } \mathbf{r} \in \mathcal{D}_i, \ 0 < z < h_i(\mathbf{r}) \ (i = 1, \dots, N) \\
\epsilon_{ml}(z) & \text{if } z < 0 \\
1 & \text{elsewhere}
\end{cases}
.
$$

The zeroth-order scattering operator is unchanged, as
it corresponds to reflection by the multilayer, while the
first-order scattering operator of the inhomogeneous sur-
face is easily found as the sum of the first-order scattering
operators S<sub>1</sub>(<b>k</b>, <b>k</b><sub>0</sub>, ε<sub>h<sub>i</sub></sub>, h<sub>i</sub>) of each of the homogeneous
surfaces h<sub>i</sub> deposited alone on the substrate:

$$
S_0(\mathbf{k}, \mathbf{k}_0) = \mathbb{R}_0(\mathbf{k}_0) \delta(\mathbf{k} - \mathbf{k}_0),
$$

<figure><img src="image_3.png" /><figcaption>Fig. 2. Heterogeneous rough deposit on a multilayer substrate.</figcaption></figure>
$$S_1(\mathbf{k}, \mathbf{k}_0) = \sum_{i=1}^{N} S_1(\mathbf{k}, \mathbf{k}_0, \varepsilon_{h_i}, h_i).$$

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 = (\varepsilon_{h_i} - 1)K^2$, and the tensor $A_i$ where $\varepsilon$ is replaced by $\varepsilon_{h_i}$ in Eq. (10) we have

$$S_2(\mathbf{k}, \mathbf{k}_0) = \sum_{i=1}^{N} S_2(\mathbf{k}, \mathbf{k}_0, \varepsilon_{h_i}, h_i) + \sum_{i \neq j} \int d\xi \hat{h}_i \times (\mathbf{k} - \xi) \hat{h}_j(\xi - \mathbf{k}_0) B_0^{ij}(\mathbf{k}, \mathbf{k}_0, \xi), \quad (26)$$

where

$$B_0^{ij}(\mathbf{k}, \mathbf{k}_0, \xi) = -\frac{V_i V_j}{8 q q_{\xi}} C^+(\mathbf{k}) A_i [C^+(\xi) + C^-(\xi)] A_j C^-(\mathbf{k}_0). \quad (27)$$

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

$$\mathbf{E}_{\text{Born}}^{\text{ext}}(\mathbf{R}) = \mathbf{E}^{\text{ref}}(\mathbf{R}) + \int d\mathbf{R}' \tilde{G}(\mathbf{R}, \mathbf{R}') V(\mathbf{R}') A \mathbf{E}^{\text{ref}}(\mathbf{R}). \quad (28)$$

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

$$\mathbf{E}_{\text{Born}}^{\text{ext}}(\mathbf{R}) = \mathbf{E}^{\text{ref}}(\mathbf{R}) + \int d\mathbf{R}' \tilde{G}(\mathbf{R}, \mathbf{R}') V(\mathbf{R}') \mathbf{E}^{\text{ref}}(\mathbf{R}), \quad (29)$$

as suggested by the form of Eq. (8). This has little consequence for small contrast ($A \approx 1$), but leads to an important difference for larger contrasts. In particular, the first perturbative order is not recovered in the limit of small heights.

## 5. NUMERICAL APPLICATIONS

The method has been numerically tested for one-dimensional surfaces (i.e., invariant along the y axis) because this allows a comparison to be made with a volume-integral rigorous method.$^{35}$ A 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.$^{26}$ In 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 time-consuming. The test geometry consists of two trapezoidal rods on a homogeneous half-plane (permittivity $\varepsilon$) as depicted in Fig. 2. The rods have identical dimensions (height = $\lambda/50$, lower basis = $\lambda/2$, upper basis = $0.275\lambda$, separation = $\lambda/2$) but different permittivity $\varepsilon_{h_1}$ and $\varepsilon_{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}(\mathbf{k}, \mathbf{k}_0)|^2$ in $s$ ($i = 2$) and $p$ ($i = 1$) polarization for an incidence angle of $20^\circ$. 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,9}$ Yet 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-

<figure><img src="image_4.png" /><figcaption>Fig. 3. Scattering cross section at $20^\circ$ incidence for the geometry of Fig. 2 with $\varepsilon_{h_1} = \varepsilon_{h_2} = \varepsilon_{ml}(z) = 2.25$ (glass on glass).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Same as Fig. 3 with $ε_{h_1} = ε_{h_2} = ε_{ml}(z) = -3 + 0.8i$ (metal on metal).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Same as Fig. 3 with $ε_{h_1} = ε_{h_2} = -3 + 0.8i$, $ε_{ml}(z) = 2.25$ (metal on glass).</figcaption></figure>

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 $\approx \lambda/50$.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Same as Fig. 3 with $ε_{h_1} = ε_{h_2} = 2.25$, $ε_{ml}(z) = -3 + 0.8i$ (glass on metal).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Same as Fig. 3 with $ε_{h_1} = 4$, $ε_{h_2} = 2.25$, $ε_{ml}(z) = -3 + 0.8i$ (silicon + glass on metal).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 8. Same as Fig. 3 with $\epsilon_{h_1} = -3 + 0.8i$, $\epsilon_{h_2} = 2.25$ and a bilayer $\epsilon_{ml}(z) = 4$ for $0 > z > -\lambda/4$ and $\epsilon_{ml}(z) = -3 + 0.8i$ for $-\infty < z < -\lambda/4$ (metal + glass on a bilayer of silicon-metal).</figcaption></figure>

## 6. 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 $\approx \lambda/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

$$I = \frac{V}{2} \int d\mathbf{r}' h^2(\mathbf{r}') \left[ \frac{d\tilde{G}}{dz'}(\mathbf{r} - \mathbf{r}', z, 0^+) A \mathbf{E}^{\text{ref}}(\mathbf{r}', 0^+) + \tilde{G}(\mathbf{r} - \mathbf{r}', z, 0^+) A \frac{d\mathbf{E}^{\text{ref}}}{dz'}(\mathbf{r}', 0^+) \right]. \quad (A1)$$

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,

$$\tilde{G}(\mathbf{r}' - \mathbf{r}'', z', z'') = \tilde{G}_s(\mathbf{r}' - \mathbf{r}'', z', z'') + \tilde{G}_{as}(\mathbf{r}' - \mathbf{r}'', z', z''), \quad (A2)$$

with

$$\tilde{G}_s(\mathbf{r}' - \mathbf{r}'', z', z'') = \frac{1}{2} [ \tilde{G}(\mathbf{r}' - \mathbf{r}'', z', z'') + \tilde{G}(\mathbf{r}' - \mathbf{r}'', z'', z') ], \quad (A3a)$$

$$\tilde{G}_{as}(\mathbf{r}' - \mathbf{r}'', z', z'') = \frac{1}{2} [ \tilde{G}(\mathbf{r}' - \mathbf{r}'', z', z'') - \tilde{G}(\mathbf{r}' - \mathbf{r}'', z'', z') ]. \quad (A3b)$$

Note that the antisymmetric part of the Green's tensor does not depend on the permittivity and assumes the simple expression

$$\tilde{G}_{as}(\mathbf{r}' - \mathbf{r}'', z', z'') = -\mathrm{sign}(z' - z'') \frac{i}{8\pi^2 K^2} \int_{\mathbb{R}^2} d\xi \exp[i\xi \cdot (\mathbf{r}' - \mathbf{r}'')] \times \exp(iq_{\xi}|z' - z'|)(\xi\hat{z} + \hat{z}\xi), \quad (A4)$$

as can be seen from the identity

$$\frac{1}{2q} [C^{-}(\mathbf{k}) - C^{+}(\mathbf{k})] = \frac{1}{2q} [P^{-}(\mathbf{k}) - P^{+}(\mathbf{k})] = \frac{1}{K^2} (\mathbf{k}\hat{\mathbf{z}} + \hat{\mathbf{k}}\mathbf{k}). \quad (A5)$$

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

$$J_s = \int d\mathbf{r}' d\mathbf{r}'' h(\mathbf{r}') h(\mathbf{r}'') \tilde{G}(\mathbf{r} - \mathbf{r}', z, 0^+) \times A \tilde{G}_s(\mathbf{r}' - \mathbf{r}'', 0^+, 0^+) A \mathbf{E}^{\text{ref}}(\mathbf{r}', 0^+). \quad (A6)$$

We now employ $\Psi(r', z')$ to denote the innermost integral in $J_{as}$ to first order in $h$ as

$$\Psi(r', z') := \left[ \int d\mathbf{r}'' \int_0^{h(\mathbf{r}'')} dz'' \tilde{G}_{as}(\mathbf{r}' - \mathbf{r}'', z', z'') \times A \mathbf{E}^{\text{ref}}(\mathbf{r}', z'') \right]_1. \quad (A7)$$

An elementary calculation then gives

$$\begin{aligned} \Psi(r', z') &= \frac{1}{4\pi^2 K^2} \int dr'' d\xi \exp[i\xi \cdot (\mathbf{r}' - \mathbf{r}'')] \exp(i\mathbf{k}_0 \cdot \mathbf{r}'') i(\xi\hat{\mathbf{z}} + \hat{\mathbf{z}}\xi) \phi(r'', z') A C^{-}(\mathbf{k}_0) \mathbf{E}^{\text{inc}} \\ &= (\nabla_{r'} \hat{\mathbf{z}} + \hat{\mathbf{z}} \nabla_{r'}) [ \phi(r', z') \times \exp(i\mathbf{k}_0 \cdot \mathbf{r}') ] A C^{-}(\mathbf{k}_0) \mathbf{E}^{\text{inc}}, \end{aligned} \quad (A8)$$
where $\phi(r'', z') := \min\{z', h(r'')\} - h(r'')/2$, yielding the following expression for $J_{as}$:

$$ 
\begin{aligned}
& \frac{V^2}{K^2} \left\{ \int d\mathbf{r}' \int_0^{h(\mathbf{r}')} dz' \tilde{G}(\mathbf{r} - \mathbf{r}', z, z') A(\nabla_{\mathbf{r}'} \hat{\mathbf{z}} + \hat{\mathbf{z}} \nabla_{\mathbf{r}'}) \right. \\
& \qquad \left. \times [(z' - h r'/2) \exp(i\mathbf{k}_0 \cdot \mathbf{r}')] A\text{C}^-(\mathbf{k}_0) \mathbf{E}^{\text{inc}}_2 \right. \\
& \qquad \left. = \frac{-V^2}{2K^2} \left\{ \int d\mathbf{r}' h(\mathbf{r}') \tilde{G}(\mathbf{r} - \mathbf{r}', z, 0^+) A[\nabla h(\mathbf{r}') \hat{\mathbf{z}} \right. \right. \\
& \qquad \left. \left. + \hat{\mathbf{z}} \nabla h(\mathbf{r}')] \exp(i\mathbf{k}_0 \cdot \mathbf{r}') A\text{C}^-(\mathbf{k}_0) \mathbf{E}^{\text{inc}} \right. \right\}
\end{aligned}
$$

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(\mathbf{r}, z)$ ($z \notin \mathcal{D}$):

$$ 
\begin{aligned}
& \left( -\frac{V}{4} \int \frac{d\mathbf{k}}{q} \exp(i\mathbf{k}\cdot\mathbf{r} + iqz) \hat{h}^2(\mathbf{k} - \mathbf{k}_0) \{q_0 C^+(\mathbf{k}) A[R(\mathbf{k}_0) - P^-(\mathbf{k}_0)] + q[R(\mathbf{k}) - P^+(\mathbf{k})]AC^-(\mathbf{k}_0) \} \right. \\
& \left. - \frac{V^2}{8K^2} \int \frac{d\mathbf{k}}{q} \exp(i\mathbf{k}\cdot\mathbf{r} + iqz) \hat{h}^2(\mathbf{k} - \mathbf{k}_0) C^+(\mathbf{k}) \right. \\
& \left. \times A[(\mathbf{k} - \mathbf{k}_0)\hat{\mathbf{z}} + \hat{\mathbf{z}}(\mathbf{k} - \mathbf{k}_0)]AC^{-a}(\mathbf{k}_0) - \frac{V^2}{8} \int \frac{d\mathbf{k}}{q} \frac{d\xi}{q_\xi} \right. \\
& \left. \times \exp(i\mathbf{k}\cdot\mathbf{r} + iqz) \hat{h}(\mathbf{k} - \xi) \hat{h}(\xi - \mathbf{k}_0) C^+(\mathbf{k}) A[C^+(\xi) + C^-(\xi)]AC^-(\mathbf{k}_0) \right) \mathbf{p}_{\text{inc}}.
\end{aligned}
$$

## APPENDIX B: EXPANDED EXPRESSIONS IN THE CANONICAL POLARIZATION BASIS

The vectors defining the fundamental polarization basis ($V, H$) are given by

$$ \mathbf{p}_1^{\pm}(\mathbf{k}_0) = \frac{k_0 \hat{\mathbf{z}} \mp q_0 \hat{\mathbf{k}}_0}{K}, \quad \mathbf{p}_2^{\pm}(\mathbf{k}_0) = \hat{\mathbf{z}} \times \hat{\mathbf{k}}_0. \quad (B1) $$

The operators $R$ and $P^\pm$ assume the following simple expressions in this basis:

$$ P^{\pm}(\mathbf{k}) = \mathbf{p}_{1}^{\pm}(\mathbf{k})\mathbf{p}_{1}^{\pm}(\mathbf{k}) + \mathbf{p}_{2}^{\pm}(\mathbf{k})\mathbf{p}_{2}^{\pm}(\mathbf{k}), \quad (B2a) $$

$$ R(\mathbf{k}) := r_1(k)\mathbf{p}_1^+(\mathbf{k})\mathbf{p}_1^-(\mathbf{k}) + r_2(k)\mathbf{p}_2^+(\mathbf{k})\mathbf{p}_2^-(\mathbf{k}). \quad (B2b) $$

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 $\epsilon$ we have

$$ r_1(k) = \frac{\epsilon q - q'}{q' + \epsilon q'}, \quad r_2(k) = \frac{q - q'}{q + q'}, \quad (B3) $$

with $q' = (\epsilon K^2 - k^2)^{1/2}$, where the complex square root with nonnegative, imaginary part has to be taken. For a bilayer $\epsilon_{ml}(z) = \epsilon_1$ for $0 > z > -L$, $\epsilon_{ml}(z) = \epsilon_2$ for $z < -L$ we have

$$ r_j(k) = \frac{r_j^{(01)}(k) + r_j^{(12)(k)}}{1 + r_j^{(01)}(k) r_j^{(12)}(k) \exp(2iq'_1L)}, \quad (B4) $$

where $r_j^{(01)}(k)$ [resp. $r_j^{(12)}(k)$] is the Fresnel reflection coefficient from the vacuum [resp. medium $\epsilon_1$] to a homogeneous medium $\epsilon_1$ [resp. $\epsilon_2$], and $q'_1 = (\epsilon_1 K^2 - k^2)^{1/2}$. Note that in some textbooks the Fresnel coefficients in $V$ polarization are defined on the magnetic field, which implies an extra $1/\sqrt{\epsilon}$ factor. The kernels $B_1$ and $B_2$ can also be decomposed in the canonical polarization basis as

$$ B_n(\mathbf{k}, \mathbf{k}_0) = \sum_{i,j=1}^{2} (B_n)_{ji}(\mathbf{k}, \mathbf{k}_0) \mathbf{p}_j(\mathbf{k}) \mathbf{p}_i(\mathbf{k}_0), \quad (B5) $$

where the coefficients $(B_n^\pm)_{ij}$ are obtained by performing the scalar products $\mathbf{p}_j(\mathbf{k}) \cdot [B_n(\mathbf{k}, \mathbf{k}_0)\mathbf{p}_i(\mathbf{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.

$$ 
\begin{aligned}
(B_1)_{11}(\mathbf{k}, \mathbf{k}_0) = & \frac{(\epsilon_h - 1)}{2iq} \left\{ \begin{aligned}[t] & [1 - r_1(k)] \\ & \times [1 - r_1(k_0)] q q_0 \hat{\mathbf{k}} \cdot \hat{\mathbf{k}}_0 - [1 + r_1(k)] \\ & \times [1 + r_1(k_0)] \frac{kk_0}{\epsilon_h} \end{aligned} \right\},
\end{aligned}
$$

$$ 
\begin{aligned}
(B_1)_{21}(\mathbf{k}, \mathbf{k}_0) = & -(B_1)_{12}(-\mathbf{k}_0, -\mathbf{k}) \\ 
& = -\frac{(\epsilon_h - 1)q_0 K}{2iq} \\ 
& \times [1 - r_1(k_0)][1 + r_2(k)]\hat{\mathbf{z}}[\hat{\mathbf{k}}, \hat{\mathbf{k}}_0],
\end{aligned}
$$

$$ 
\begin{aligned}
(B_1)_{22}(\mathbf{k}, \mathbf{k}_0) = & -\frac{(\epsilon_h - 1)K^2}{2iq} [1 + r_2(k_0)] \\ 
& \times [1 + r_2(k)]\hat{\mathbf{k}} \cdot \hat{\mathbf{k}}_0.
\end{aligned}
$$

Second-order kernels for homogeneous deposits. To shorten the formulas we will use the abbreviated notation $r_j^0 = r_j(k_0)$, $r_j = r_j(k)$, and $r'_j = r_j(\xi)$ and decompose the kernel in the form

$$ (B_2)_{ji}(\mathbf{k}, \mathbf{k}_0, \xi) = \frac{\epsilon_h - 1}{4} \alpha_{ji} + \frac{(\epsilon_h - 1)^2}{4} \beta_{ji}. $$

Then we have
$$
\begin{aligned}
\alpha_{11} = & -\frac{(\varepsilon_h - 1)}{2\varepsilon_h q} [kq_0(1 + r_1)(1 - r_1^0)(\hat{\mathbf{k}} \cdot \hat{\mathbf{k}}_0 k - k_0) \\
& + k_0 q(1 - r_1)(1 + r_1^0)(k - k_0 \hat{\mathbf{k}} \cdot \hat{\mathbf{k}}_0)] \\
& - \frac{1}{\varepsilon_h q} [(kk_0 + \varepsilon_h q q_0 \hat{\mathbf{k}} \cdot \hat{\mathbf{k}}_0)(q - q_0)(r_1 - r_1^0) \\
& + (kk_0 - \varepsilon_h q q_0 \hat{\mathbf{k}} \cdot \hat{\mathbf{k}}_0)(q + q_0)(r_1 r_1^0 - 1)],
\end{aligned}
$$

$$
\begin{aligned}
\beta_{11} = & (1 - r_1)(1 - r_1^0) \left[ \frac{K^2 q_0}{q_\xi} \hat{\mathbf{z}}[\hat{\mathbf{k}}, \hat{\xi}] \hat{\mathbf{z}}[\hat{\mathbf{k}}_0, \hat{\xi}] (1 + r_2') \right. \\
& \left. - q_0 q_\xi (1 - r_1') \hat{\xi} \cdot \hat{\mathbf{k}} \hat{\xi} \cdot \hat{\mathbf{k}}_0 \right] - \frac{1}{\varepsilon_h^2 q} (1 + r_1) \\
& \times \left[ \frac{kk_0 \xi^2}{q_\xi} (1 + r_1')(1 + r_1^0) + q_0 k \xi \varepsilon_h \hat{\xi} \cdot \hat{\mathbf{k}}_0 r_1' \right. \\
& \left. \times (1 - r_1^0) \right] + \frac{k_0 \xi}{\varepsilon_h} (1 - r_1)(1 + r_1^0) r_1' \hat{\mathbf{k}} \cdot \hat{\xi};
\end{aligned}
$$

$$
\alpha_{12} = K \hat{\mathbf{z}}[\hat{\mathbf{k}}, \hat{\mathbf{k}}_0] \left[ (q + q_0)(r_1 r_2^0 + 1) + (q - q_0) \right. \\
\times (r_1 + r_2^0) + \frac{k^2}{2q\varepsilon_h} (1 + r_1)(1 + r_2^0) \left. \right],
$$

$$
\beta_{12} = \hat{\mathbf{z}}[\hat{\xi}, \hat{\mathbf{k}}_0](1 + r_2^0) \left\{ \frac{k \xi}{q \varepsilon_h} r_1'(1 + r_1) - (1 - r_1) \times [q_{\hat{\xi}}(1 + r_1')\hat{\mathbf{k}} \cdot \hat{\xi} - K(1 + r_2')\hat{\mathbf{z}}(\hat{\xi}, \hat{\mathbf{k}})\hat{\xi} \cdot \hat{\mathbf{k}}_0] \right\};
$$

$$
\alpha_{22} = \frac{K^2}{q} \hat{\mathbf{k}} \cdot \hat{\mathbf{k}}_0 [q_0(1 + r_2)(1 - r_2^0) + q(1 - r_2) \times (1 + r_2^0)],
$$

$$
\beta_{22} = -(1 + r_2^0)(1 + r_2) \frac{K^2}{q} \left[ \frac{K^2}{q_\xi} (\hat{\xi} \cdot \hat{\mathbf{k}})(\hat{\xi} \cdot \hat{\mathbf{k}}_0)(1 + r_2') + q_\xi (1 - r_1') \hat{\mathbf{z}}[\hat{\mathbf{k}}, \hat{\xi}] \hat{\mathbf{z}}[\hat{\xi}, \hat{\mathbf{k}}_0] \right].
$$

We recall that **a**[**b**, **c**] is a short notation for **a**·(**b** × **c**). The missing cross term is given by

$$
(B_2)_{21}(\mathbf{k}, \mathbf{k}_0, \xi) = -(B_2)_{12}(-\mathbf{k}_0, -\mathbf{k}, -\xi), \quad (B6)
$$

which is the expression of reciprocity for the second-order kernel.

Second-order kernels for heterogeneous deposits.

$$
(B_c^{ij})_{11}(\mathbf{k}, \mathbf{k}_0, \xi) = \frac{(\varepsilon_{h_i} - 1)(\varepsilon_{h_j} - 1)}{4} \Big|_{(\varepsilon_{h_i}^2, \varepsilon_{h_j}, \varepsilon_{h_k}) \to (\varepsilon_{h_i} \varepsilon_{h_j}, \varepsilon_{h_k}, \varepsilon_{h_l})} \beta_{11},
$$

$$
(B_c^{ij})_{12}(\mathbf{k}, \mathbf{k}_0, \xi) = \left. \frac{(\varepsilon_{h_i} - 1)(\varepsilon_{h_j} - 1)}{4} \beta_{12} \right|_{\varepsilon_h \to \varepsilon_{h_j}},
$$

$$
(B_c^{ij})_{22}(\mathbf{k}, \mathbf{k}_0, \xi) = \frac{(\varepsilon_{h_i} - 1)(\varepsilon_{h_j} - 1)}{4} \beta_{22}.
$$

In these formulas it is understood that the $\beta_{ij}$ of the homogeneous case have to be taken with appropriate replacement of the $\varepsilon_h$ in the order of appearance.

ACKNOWLEDGMENTS

Many thanks go to Tanos El Fouhaily for useful comments on the manuscript.

Corresponding author C.-A. Guerin's e-mail address is charles-antoine.guerin@fresnel.fr.

REFERENCES

1. M. Saillard and A. Sentenac, “Rigorous solutions for electromagnetic scattering from rough surfaces,” Waves Random Media **11**, R103–R137 (2001).

2. A. G. Voronovich, *Wave Scattering from Rough Surfaces*, Springer Series on Wave Phenomena (Springer, New York, 1994).

3. T. Tsang, J. A. Kong, and R. Shin, *Theory of Microwave Remote Sensing* (Wiley, New York, 1985).

4. Lord Rayleigh, *The Theory of Sound*, 3rd ed. (MacMillan, London, 1896).

5. U. Fano, “The theory of anomalous diffraction gratings and of quasi-stationary waves on metallic surfaces (Sommerfeld’s waves),” J. Opt. Soc. Am. **31**, 213–222 (1941).

6. S. O. Rice, “Reflection of electromagnetic waves from slightly rough surfaces.” Commun. Pure Appl. Math. **4**, 351–378 (1951).

7. M. F. Chen and A. K. Fung, “A numerical study of the regions of validity of the Kirchhoff and small-perturbation rough surface scattering models,” Radio Sci. **23**, 163–170 (1988).

8. J. M. Soto-Crespo, M. Nieto-Vesperinas, and A. T. Friberg, “Scattering from slightly rough random surfaces: a detailed study on the validity of the small perturbation method,” J. Opt. Soc. Am. **7**, 1185–1201 (1990).

9. J. A. Sanchez, A. A. Maradudin, and E. R. Mendez, “Limits of validity of three perturbation theories of the specular scattering of light from one-dimensional, randomly rough, dielectric surfaces.” J. Opt. Soc. Am. **12**, 1547–1557 (1995).

10. D. P. Winebrenner and A. Ishimaru, “Application of the phase perturbation technique to randomly rough surface,” J. Opt. Soc. Am. A **2**, 2285–2294 (1985).

11. D. M. Milder, “An improved formalism for wave scattering from rough surfaces,” J. Acoust. Soc. Am. **89**, 529–541 (1991).

12. S. Smith, “The operator expansion formalism for electromagnetic scattering from rough dielectric surfaces,” Radio Sci. **31**, 1377–1385 (1996).

13. A. G. Voronovich, “Small-slope approximation for electromagnetic wave scattering at a rough interface of two dielectric half-spaces,” Waves Random Media **4**, 337–367 (1994).

14. O. P. Bruno and F. Reitich, “Numerical solution of diffrac-
tion problems: a method of variation of boundaries,” J.
Opt. Soc. Am. A **10**, 1168–1175 (1993).

15. O. P. Bruno and F. Reitich, “Numerical solution of diffrac-
tion problems: a method of variation of boundaries. II.
Finitely conducting gratings, Padé approximants, and sin-
gularities,” J. Opt. Soc. Am. A **10**, 2307–2316 (1993).
16. O. P. Bruno and F. Reitich, "Numerical solution of diffraction problems: a method of variation of boundaries. III. Doubly periodic gratings," J. Opt. Soc. Am. A **10**, 2551–2562 (1993).

17. K. Sarabandi, O. Yisok, and F. Ulaby, "A numerical simulation of scattering from one-dimensional, inhomogeneous, dielectric, random surfaces," IEEE Trans. Geosci. Remote Sens. **34**, 425–432 (1996).

18. L. Sung, G. Mulholland, and T. Germer, "Polarized light-scattering measurements of dielectric spheres upon a silicon surface," Opt. Lett. **24**, 866–868 (1999).

19. R. Carminati and J.-J. Greffet, "Influence of dielectric contrast and topography on the near-field scattered by an inhomogeneous surface," J. Opt. Soc. Am. A **12**, 2716–2725 (1995).

20. S. Linden, J. Kuhl, and H. Giessen, "Controlling the interaction between light and gold nanoparticles with selective suppression of extinction," Phys. Rev. B **86**, 4688–4691 (2001).

21. L. J. Lévesque and B. E. Paton, "Detection of defects in multiple layer structures by using surface plasmon resonance," Atmos. Ocean. **36**, 7199–7203 (1997).

22. I. Ohlidal and K. Navratil, "Scattering of light from multilayer systems with rough boundary," Prog. Opt. **34**, 251–334 (1995).

23. H. Giovannini, M. Saillard, and A. Sentenac, "Numerical study of scattering from rough inhomogeneous films," J. Opt. Soc. Am. A **15**, 1182–1191 (1998).

24. A. Sentenac and J.-J. Greffet, "Scattering by 2D particles deposited on a dielectric planar waveguide, a near-field and far-field study," Waves Random Media **5**, 145–155 (1995).

25. L. Li, "New formulation of the Fourier modal method for crossed surface-relief gratings," J. Opt. Soc. Am. A **14**, 2758–2767 (1997).

26. P. Chaumet, A. Rahmani, F. de Fornel, and J. P. Dufour, "Evanescent light scattering: the validity of the dipole approximation," Phys. Rev. B **58**, 2310 (1998).

27. "Computational wave issues in remote sensing, imaging and target identification, propagation, and inverse scattering," special issue, IEEE Trans. Geosci. Remote Sens. **38**, (2000).

28. A. Soubret, G. Berginc, and C. Bourrely, "A new application of reduced Rayleigh equations to electromagnetic wave scattering by two-dimensional, randomly rough surfaces," Phys. Rev. B **63**, 245411 (2001).

29. P. Dinesen and J. Hesthaven, "Fast and accurate modeling of waveguide grating couplers, three-dimensional vectorial case," J. Opt. Soc. Am. A **18**, 2876–2885 (2001).

30. N. Zhuk, "Scattering of em waves from a slightly rough surface of a generally anisotropic plane-layered half space," IEEE Trans. Antennas Propag. **45**, 1774–1782 (1997).

31. A. Fuks, "Wave diffraction by a rough boundary of an arbitrary plane-layered medium," IEEE Trans. Antennas Propag. **49**, 630–639 (2001).

32. S. Dietrich and A. Haase, "Scattering of X-rays and neutrons at interfaces," Phys. Rep. **260**, (1995).

33. P. Johansson, "Light scattering from disordered overlayers of metallic nanoparticles," Phys. Rev. B **64**, 165405 (2001).

34. T. Sondergaard and S. Boshelvonyi, "Vectorial model for multiple scattering by surface nanoparticles via surface polariton–polariton interactions," Phys. Rev. B **67**, 165405 (2003).

35. F. Pincemin, A. Sentenac, and J. J. Greffet, "Near field scattered by a dielectric rod below a metallic surface," J. Acoust. Soc. Am. **11**, 1117–1127 (1994).

36. *X-Ray and Neutron Reflectivity: Principles and Applications*, J. Daillant and A. Gibaud, eds. (Springer, New York, 1999).

37. J. Lekner, *Theory of Reflection of Electromagnetic and Particle Waves* (Martinus Nijhoff, Dordrecht, The Netherlands, 1987).

38. L. Tsang, J. A. Kong, and K. H. Ding, *Scattering of electromagnetic waves*, first volume of three in Wiley Series in Remote Sensing (Wiley-Interscience, New York, 2001).