# Beam implementation in a nonorthogonal coordinate system: Application to the scattering from random rough surfaces

K. Edee, B. Guizal, G. Granet, A. Moreau

► To cite this version:

K. Edee, B. Guizal, G. Granet, A. Moreau. Beam implementation in a nonorthogonal coordinate system: Application to the scattering from random rough surfaces. Journal of the Optical Society of America. A Optics, Image Science, and Vision, 2008, 25 (3), pp.796-804. 10.1364/JOSAA.25.000796. hal-00256005

HAL Id: hal-00256005

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

Submitted on 22 Apr 2021

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
Beam implementation in a nonorthogonal
coordinate system: Application to the scattering
from random rough surfaces

K. Edee,¹ B. Guizal,²,* G. Granet,¹ and A. Moreau¹

¹Laboratoire des Sciences et Matériaux pour l'Electronique et d'Automatique, CNRS/UMR 6602,
Université Blaise Pascal, Les Cézeaux, 63177 Aubière Cedex, France

²Département d'Optique P.M. Duffieux, Institut Franche-Comté Electronique, Mécanique,
Thermique et Optique-Sciences et Technologies, UMR 6174 du CNRS Université de Franche-Comté,
16 Route de Gray 25030 Besançon Cedex, France

*Corresponding author: bguizal@univ-fcomte.fr

The C method is known to be one of the most efficient and versatile tools established for modeling diffraction gratings. Its main advantage is the use of a coordinate system in which the boundary conditions apply naturally and are, *ipso facto*, greatly simplified. In the context of scattering from random rough surfaces, we propose an extension of this method in order to treat the problem of diffraction of an arbitrary incident beam from a perfectly conducting (PEC) rough surface. For that, we were led to revisit some numerical aspects that simplify the implementation and improve the resulting codes.

# 1. INTRODUCTION

Actual surfaces are necessarily rough, and the main difficulty encountered in the modeling of their interaction with electromagnetic waves is related to their level of roughness. This latter is dependent on the geometry and the wavelength of the incident radiation (i.e., how the surface is seen by the wave). Two approximate methods have been developed to study rough surfaces with weak roughness and led to analytic expressions of the diffracted intensity in terms of the statistical parameters of the surfaces such as the rms $\sigma$ and correlation length $l_c$. Among these methods, the most known is the perturbation theory applied to the integral equation of the field [1]. Other approaches use rather heuristic hypotheses concerning the diffraction phenomenon. This is the case for the Kirchhoff approximation [1–3], which can be applied to either shallow surfaces or surfaces with low radii of curvature. In the so-called resonant domain where the geometrical features of the scattering surface are of the order of magnitude of the wavelength, rigorous solutions of the diffraction problem are necessary. This is why exact formalisms have been developed, and one can roughly divide them in two main categories: integral methods [4–7] and differential methods [8]. Integral methods are based on the Green theorem that is used to derive an integral equation that is rigorously solved. As for the differential methods, they rely on solving the partial differential equations directly deduced from Maxwell's equations and transposed into Fourier space. This is the case for Fourier methods, among which one can find the differential method [8], the Fourier modal method [9], and finally the C method [10,11], which will be used in this work. This approach has been introduced, originally, for the study of corrugated waveguides and diffraction gratings. Recently, it has been successfully extended to diffraction of plane waves from nonperiodic surfaces [12–15].

In the present paper we focus, for the sake of conciseness, on the case of a Gaussian rough surface with Gaussian spectrum made of a perfectly conducting (PEC) material and illuminated by a limited beam. This is by no means a limitation of our approach; the case of dielectric or metallic surfaces can be obtained mechanically from the present work. We will use an improved version of the C method to treat this problem. The article is organized as follows. In Section 2, we outline the main steps of the C method as applied to this problem, and in Section 3, we will specify some properties concerning the eigensolutions that will be exploited in Section 4 in the construction of the incident and diffracted beams. Finally, the method is tested numerically against usual criteria of energy conservation and our computations are compared to results taken from the literature and based on both integral methods [15] and differential methods [8].

# 2. STATEMENT OF THE PROBLEM AND INTRODUCTION TO THE C METHOD

We consider (Fig. 1) a rough PEC surface S whose shape is described by a function $y=a(x)$ in the orthonormal coordinate system $(O,x,y,z)$. This surface is invariant along the $z$ direction and is illuminated by an incident wave that can be either TE (the only non-null components of the fields are $E_z$, $H_x$, and $H_y$) or TM (the only non-null components of the fields are $H_z$, $E_x$, and $E_y$) polarized such that we are in the case of the so-called classical incidence [16]. It is well known that in these cases of polar-
<figure><img src="image_2.png" /><figcaption>Fig. 1. PEC rough surface illuminated by an incident beam under the mean incidence $\theta_i$.</figcaption></figure>

ization, all the field components can be expressed in terms of $E_z$ in the TE case or in terms of $H_z$ in the TM case. Let us add that throughout this paper, time dependence of the fields will be held by the term $\exp(i\omega t)$, where $\omega$ denotes the circular frequency of the monochromatic incident radiation.

The starting point of the C method is the choice of a coordinate system in which the diffracting surface (S) coincides exactly with a coordinate surface. This has the big advantage of greatly simplifying the process of writing the boundary conditions. The simplest coordinate system is the so-called translation coordinate system (TCS), which is given by

$$x=u, \quad y=v+a(u), \quad z=w. \tag{1}$$

Thus the boundary conditions can be written in a simple way at the coordinate $v=0$. In this new coordinate system, it can be shown [13] that $\Psi$ satisfies the following set of partial differential equations that is valid in the two cases of polarization:

$$\begin{pmatrix} \partial_u \partial_u + k^2 I & 0 \\ 0 & I \end{pmatrix} \begin{pmatrix} \Psi \\ \partial_v \Psi \end{pmatrix} = \begin{pmatrix} \partial_u \dot{a} + \dot{a} \partial_u - (1+\dot{a}\dot{a}) \\ I & 0 \end{pmatrix} \partial_v \begin{pmatrix} \Psi \\ \partial_v \Psi \end{pmatrix}, \tag{2a}$$

where $\Psi=E_z$ in the TE case and $\Psi=ZH_z$ in the TM case, with $Z$ denoting the impedance of vacuum. $I$ is the identity operator, $k$ is the wave vector, $\partial_u=\partial/\partial u$, $\partial_v=\partial/\partial v$, and $\dot{a}=da/du$.

The second step consists in seeking plane-wave solutions under the form

$$\Psi(u,v) = \exp(-ikrv)\psi(u), \quad \partial_v\Psi(u,v) = \exp(-ikrv)\psi_r(u).$$

Since $\partial_v = -ikr$, Eq. (2a) can be recast to

$$\begin{pmatrix} \partial_u \partial_u + k^2 I & 0 \\ 0 & I \end{pmatrix} \begin{pmatrix} \psi \\ \psi_r \end{pmatrix} = -ikr \begin{pmatrix} \partial_u \dot{a} + \dot{a} \partial_u - (1+\dot{a}\dot{a}) \\ I & 0 \end{pmatrix} \begin{pmatrix} \psi \\ \psi_r \end{pmatrix}. \tag{2b}$$

From these equations, it is worth noticing that the profile function comes into play through its derivative $\dot{a}$, which must exist. However, profiles with rather sharp edges can be treated efficiently by approximating them by smoother ones or by choosing another coordinate system. These issues have been addressed in previous works and the interested reader can refer to [17–19].

Equation (2b) is then transposed into the spectral domain by Fourier transformation

$$\hat{f}(\alpha) = \int_{-\infty}^{+\infty} f(u) \exp(-ik\alpha u) du, \tag{3}$$

where $f(u)$ stands for $\dot{a}(u)$, $\psi(u)$, or $\psi_r(u)$. We assume that the spectrum of the function $f(u)$ is a bounded stand. This hypothesis can be justified by the limited spatial extension of the incident beam. Under these conditions, Eq. (2b) becomes

$$L_1(\alpha) * \hat{\Gamma}(\alpha) = -ikrL_2(\alpha) * \hat{\Gamma}(\alpha), \tag{4}$$

where $L_1(\alpha) * \hat{\Gamma}(\alpha)$ [resp. $L_2(\alpha) * \hat{\Gamma}(\alpha)$] denotes the Fourier transform of the first (resp. the second) hand of Eq. (2b), * is the convolution product, and $\hat{\Gamma}(\alpha)=(\hat{\psi}(\alpha), \hat{\psi}_r(\alpha))^T$. At this stage, Eq. (4) is solved numerically by sampling the spectral domain into a set discrete points ($\alpha_n$) over which the fields are evaluated:

$$\alpha_n = n\Delta\alpha, \quad \text{where } n \text{ is a relative integer.} \tag{5}$$

Then, according to Shannon's first theorem, the discretized functions can be interpolated by

$$\hat{f}(\alpha) = \sum_{n=-\infty}^{+\infty} f_n \sin\left(\frac{\pi(\alpha - n\Delta\alpha)}{\Delta\alpha}\right) = \sum_{n=-\infty}^{+\infty} f_n \hat{b}_n(\alpha). \tag{6}$$

Numerically, only a finite sequence of sampled values is kept by first choosing a spectral increment $\Delta\alpha$ and a number $N=2M+1$ of sample points so that the corresponding sequence is given through the relation $\Delta\alpha=n\Delta\alpha$ (with $n \in [-M, M]$). Note that when $|\alpha_n| \le 1$ it corresponds to a propagating wave, whereas it corresponds to an evanescent one for $|\alpha_n| > 1$.

For a given $\Delta\alpha$, increasing $N$ (or $M$) results in adding more propagative and evanescent waves. In addition, in order to refine the density of samples it suffices to take smaller values of $\Delta\alpha$.

After this operation, Eq. (4) becomes

$$\begin{pmatrix} \alpha\alpha - \mathbf{I} & 0 \\ 0 & \frac{i}{k}\mathbf{I} \end{pmatrix} \begin{pmatrix} \psi \\ \psi_r \end{pmatrix} = r \begin{pmatrix} \alpha\dot{a} + \dot{a}\alpha & -\frac{1}{k}(\mathbf{I} + \dot{a}\dot{a}) \\ \mathbf{I} & 0 \end{pmatrix} \begin{pmatrix} \psi \\ \psi_r \end{pmatrix}, \tag{7}$$

where $\alpha=\text{diag}(\alpha_n)$, $\dot{a}$ is a toeplitz matrix formed by the Fourier coefficients of the derivative of the profile function ($\dot{a}_{nm}=\dot{a}_{n-m}$), and $\mathbf{I}$, is the identity matrix, $\psi$ (resp. $\psi_r$) is a column vector of $\hat{\psi}(\alpha)$ components [resp. $\hat{\psi}_r(\alpha)$] on $\hat{b}_n$-basis series.

In the spectral domain, the last equations correspond to an eigenvalue problem that can be solved by standard numerical routines. It is of fundamental importance to understand that solving Eq. (4) gives the elementary solutions “living” in the new coordinate system.

We are now going to use these elementary solutions to build the fields that correspond to our diffraction problem. Above the surface S, the field can be seen as the sum of the incident beam and the diffracted one. Each of these beams is constructed as the sum of elementary waves. This gives for $\hat{\Psi}(\alpha, v)$ above the PEC surface

$$\hat{\Psi}(\alpha, v) = \hat{\Psi}_{\text{inc}}(\alpha, v) + \hat{\Psi}_{\text{diff}}(\alpha, v), \tag{8}$$
$$ \hat{\Psi}(\alpha, v) = \sum_{q=-M}^{M} I_q \hat{\Psi}_q^-(\alpha, v) + \sum_{q=-M}^{M} R_q \hat{\Psi}_q^+(\alpha, v), \quad (9) $$

with

$$ \hat{\Psi}_{q}^{\pm}(\alpha, v) = \exp(-ikr_{q}^{\pm}v) \sum_{n=-M}^{M} \psi_{nq}^{\pm} \hat{b}_{n}(\alpha). \quad (10) $$

$I_q$ and $R_q$ are the spectral amplitudes of the incident and the diffracted beams, respectively. According to the $\exp(-ikr_q v)$ dependence and in order to satisfy the outgoing Summerfeld condition, the eigenvalues $r_q$ are sorted such that their imaginary parts are negative. Positive real parts correspond to diffracted waves whereas negative real parts correspond to incident waves on the diffracted surface. Naturally, the corresponding eigenvectors ($\psi_{nq}^\pm$) are rearranged accordingly.

The expression of the total field in the spatial domain is obtained by the inverse Fourier transform of Eq. (9):

$$ \begin{aligned} \Psi(u, v) = \Delta\alpha\Pi_{\Delta a}(u) \sum_q I_q \exp(-ikr_q^-v) \sum_{n=-M}^{M} \psi_{nq}^- \exp(-ik\alpha_n u) \\ + \Delta\alpha\Pi_{\Delta a}(u) \sum_q R_q \exp(-ikr_q^+v) \\ \times \sum_{n=-M}^{M} \psi_{nq}^+ \exp(-ik\alpha_n u), \end{aligned} \quad (11) $$

with

$$ \Pi_{\Delta a}(u) = \begin{cases} 1 & \text{if } u \in [-1/2\Delta\alpha \quad 1/2\Delta\alpha] \\ 0 & \text{otherwise} \end{cases}. $$

Now it is time to write the boundary conditions that are satisfied by the tangential components of the fields. From the obtained $z$ component $\Psi(u, v)$ we can determine the other tangential component $\Phi(u, v)$ ($E_u$ for TM polarization and $H_u$ for TE polarization), which is necessary to write the boundary conditions, through the expression [10–14]

$$ \Phi(u, v) = \eta[\partial_v - \dot{a}(u)(\partial_u - \dot{a}(u)\partial_v)]\Psi(u, v), \quad (12) $$

Here, $\eta=i/kZ$ for TE polarization and $\eta=-iZ/k$ for TM polarization.

Substituting the expression of $\Psi(u, v)$ from Eq. (11) into Eq. (12) gives

$$ \begin{aligned} \Phi(u, v) = \Delta\alpha\Pi_{\Delta a}(u) \sum_q I_q \exp(-ikr_q^-v) \sum_{n=-M}^{M} \phi_{nq}^- \exp(-ik\alpha_n u) \\ + \Delta\alpha\Pi_{\Delta a}(u) \sum_q R_q \exp(-ikr_q^+v) \\ \times \sum_{n=-M}^{M} \phi_{nq}^+ \exp(-ik\alpha_n u), \end{aligned} \quad (13) $$

with $\Phi^{+,-} = \eta([\mathbf{I}+\hat{\mathbf{a}}\hat{\mathbf{a}}]\psi^{+,-}\mathbf{r}^{+,-}\hat{\mathbf{a}}\alpha\psi^{+,-})$, where $\Phi^{+,-}$ (resp. $\psi^{+,-}$) is a matrix whose elements are $\phi_{nq}^{+,-}$ (resp. $\psi_{nq}^{+,-}$).

Finally, as the lower medium (under the surface S) is perfectly conducting, the electric fields are null inside it, and writing the boundary conditions consists simply in nullifying the tangential components $\Psi$ or $\Phi$ according to

the polarization, on the surface S given by $v=0$. This leads to a set of algebraic equations linking the unknown spectral amplitudes $R_q$ of the diffracted wave to the known spectral amplitudes $I_q$ of the incident one. Once this system is solved, one can compute the spectral power density $\Delta\tilde{P}^s$ and the bistatic scattering coefficient $\sigma$ whose expressions are derived in Appendices B and C:

$$ \Delta\tilde{P}^s(\alpha_q) = \frac{|R_q|^{2\Re} \left( \sum_{n=-M}^{M} \psi_{nq} \Phi_{nq}^* \right)}{\sum_{q/\tilde{\alpha}_q \in [-1;1]} |I(\tilde{\alpha}_q)|^{2\Re} \left( \sum_{n=-M}^{M} \psi_{nq} \Phi_{nq}^* \right)}, \quad (14) $$

$$ \sigma(\theta_q) = \frac{\Delta \tilde{P}^s(\alpha_q)}{\Delta \alpha} \cos(\theta_q). \quad (15) $$

### 3. NUMERICAL IMPLEMENTATION

One of the important keys of the method presented above is the operation of sorting the eigenvalues and assigning a given eigenvalue $r_q$ to a direction of propagation given by its angle $\theta_q$. This is done by observing numerically that the eigenvalues corresponding to propagative waves tend toward values that are exactly cos $\theta_q$ when the number of samples is increased. Note that for these waves $\alpha_q =k \sin \theta_q$. This provides a simple way of sorting the eigenvalues only if there is no degeneracy, which is not the case if the spectral domain is discretized symmetrically. The eigenvalues are then degenerate two by two, and this indication ruins any hope of succeeding to sort the set of the computed eigenvalues. To overcome this difficulty one can, in a naïve approach, discretize the spectral range asymmetrically, but this leads to some difficulties related to the sorting precision. Another way, used in all the previous implementations of the C method, is to replace all the eigensolutions corresponding to real eigenvalues by plane waves (for which the propagation angle is known) expressed in the curvilinear coordinate system. This approach is known as the modified C method.

One elegant way, to our knowledge never reported before for the C method, is to sort the eigenvalues by using their corresponding eigenvectors. Indeed the direction of propagation corresponding to a fixed eigenvalue $r_q$ is given exactly by the barycenter of its associate eigenvector:

$$ \tilde{\alpha}_q = \frac{\int \alpha |\hat{\Psi}_q(\alpha)|^2 d\alpha}{\int |\hat{\Psi}_q(\alpha)|^2 d\alpha}. \quad (16) $$

Numerically, it is found that each barycenter $\tilde{\alpha}_q$ tends toward a spectral point $\alpha_q$ as the number of samples $N$ is increased. It is important to emphasize that all barycenters are found to be different, which makes it possible to associate, without any ambiguity, a given eigenvalue to a propagation direction.
We turn now to describing the incident field. In what follows, we will assume that the incident radiation is a Gaussian beam with spectral amplitudes given by the distribution

$$I(\alpha) = \frac{g}{2\sqrt{\pi}} \exp\left[-\frac{k^2(\alpha - \alpha_0)^2 g^2}{4}\right]. \quad (17)$$

Here, g denotes the waist of the beam at the origin and $\alpha_0$ the mean x component of the incident wave vector.

The method has been implemented for random rough surfaces with Gaussian correlation functions and with power spectrums given by $W(\alpha) = h^2 l_c / (2\sqrt{\pi} \exp(-a^2 l_c^2 / 4))$. The surface is then assumed to possess Gaussian height probability distribution with a rms h and a correlation length $l_c$. The numerical procedure used to generate the surface is a spectral method previously reported by Tsang et al. [20]. The profile is expressed as a Fourier expansion $a(x) = \sum_{n=1-P/2}^{P/2} b_n \exp(i 2 \pi n x / d)$ whose Fourier coefficients $b_n$ are evaluated from a distribution of a set of P Gaussian uncorrelated random numbers with zero mean and unit variance passed through a Gaussian filter. As the C method is based on Fourier expansions of the fields and of the derivative of the profile function, the $b_n$ coefficients allow a direct determination of matrix $\hat{a}$ through the relation $\hat{a}_{nm} = ik\Delta\alpha(n-m)b_{n-m}$.

## 4. RESULTS AND DISCUSSION

In this section, we first check our numerical code against usual criteria of convergence and energy conservation. This is done through a discussion on some key parameters that influence the accuracy of the method, namely, the rms height of the surface and its correlation length. After that, a comparison of our results with some published data will be performed.

### A. Convergence and Accuracy

First of all, it is necessary to get an idea about the domain of applicability of the C method when studying rough random surfaces. One key indicator of the accuracy is known to be the energy conservation test expressed through the following relation (see Appendix B):

$$\sum_{q/\alpha_q \in [-1;1]} \Delta \tilde{P}^s(\alpha_q) = 1. \quad (18)$$

Nevertheless, it is well established that in the case of the C method, the last equality is valid whatever the truncation order. This means that the energy conservation is valid intrinsically and then cannot be used to estimate the convergence rate. In order to measure the convergence speed of our implementation, we compute the power density $\Delta \tilde{P}^s(\alpha_0)$ in the specular direction as the number of samples N is increased, and we denote $\Delta \tilde{P}^{s*}$ the value corresponding to maximum N. If we assume that this last value is the “exact” one, then we can estimate the convergence speed by introducing the following error function:

$$\xi(N) = - \mathrm{Int}(\log_{10} |\Delta \tilde{P}^s(\alpha_0) - \Delta \tilde{P}^{s*}|), \quad (19)$$

where $\mathrm{Int}(\delta)$ designates the integer part of $\delta$. Note that $\xi(N)$ directly gives an idea about the number of correct digits in the computed quantity.

As pointed out above, there are two important parameters that may influence the accuracy of the C method in the present implementation: the rms height and the correlation length of the diffracting rough surface. These two parameters are related to the fact that the profile function $a(x)$ appears in the field equations through its derivative $\hat{a}(x)$ and more precisely through the Fourier coefficients of this latter. This means that the computational load will follow the sharpness of the derivative, which appears to be important for large rms heights and short correlation lengths. To support these assertions, we present three cases with realizations chosen upon a set of surfaces having different statistical parameters h and $l_c$, all the other

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

<figure><img src="image_4.png" /><figcaption>Fig. 2. (Color online) (a) Single realization of a Gaussian rough surface profile and its derivative. $h=0.05\lambda$, $l_c=0.35\lambda$, $d=25.6\lambda$. (b) Accuracy in reflectivity for one profile realization [shown in (a)]. $\theta_i=-30^\circ$, $\lambda=1$, $g=0.25d$, TE and TM polarizations.</figcaption></figure>
parameters remaining unchanged. In all the cases presented below, the surfaces are illuminated under the mean incidence $\theta_i=-30^\circ$ by a Gaussian beam with the following characteristics: $\lambda=1$ and $g=0.25d$. The first case consists of a surface with low $h(0.05\lambda)$ and $l_c(0.35\lambda)$. The shape of the surface of such a realization is shown in Fig. 2(a), and the accuracy function $\xi(N)$ is plotted in Fig. 2(b) for both TE and TM polarizations. It can be seen that it is necessary to take $N$ larger than 150 if one wants to stabilize the fourth digit. In the second example, presented in Figs. 3(a) and 3(b), we increase $h$ (to $0.2\lambda$) and keep the same correlation length. As expected, the number of samples necessary to reach the same accuracy level as in the first case increases to roughly 300. Next we keep the new value of $h$ and increase the correlation length to $l_c = 1.5\lambda$. This time, the needed number of samples is about 50 [Fig. 4(b)], and this low value is in complete concordance with the fact that the profile is smoother (larger radii of curvature) as can be seen from Fig. 4(a). Notice here

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

<figure><img src="image_6.png" /><figcaption>Fig. 3. (Color online) (a) Single realization of a Gaussian rough surface profile and its derivative. $h=0.2\lambda$, $l_c=0.35\lambda$, $d=25.6\lambda$. (b) Accuracy in reflectivity for one profile realization [shown in (a)] $\theta_i=-30^\circ$, $\lambda=1$, $g=0.25d$, TE and TM polarizations.</figcaption></figure>

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

<figure><img src="image_8.png" /><figcaption>Fig. 4. (Color online) (a) Single realization of a Gaussian rough surface profile and its derivative. $h=0.2\lambda$, $l_c=1.5\lambda$, $d=25.6\lambda$. (b) Accuracy in reflectivity for one profile realization [shown in (a)]. $\theta_i=-30^\circ$, $\lambda=1$, $g=0.25d$, TE and TM polarizations.</figcaption></figure>

that, as it is known for the C method, there is no fundamental difference between the convergence behavior for the two cases of polarization.

## B. Comparison with Results Obtained by an Integral Method

Now that we have an idea about the control of the accuracy of the method, we are ready to compare it with some other approaches. We will primarily confront our results with those obtained from the integral method used by Tsang *et al.* in [20] and from the finite-difference time domain (FDTD) method used by Hastings *et al.* in [21]. We begin with an example taken from [20] with low rms height ($h=0.05\lambda$) and correlation length $l_c=0.35\lambda$. The computed averaged bistatic scattering coefficient $\langle\sigma(\theta)\rangle$ over 100 realizations is shown in Fig. 5(a) for the TE case and in Fig. 5(b) for the TM case. It can be seen that there
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 5. (a) Angular power density averaged through 100 realizations, TE polarization. $\theta_i=-30°$, $\lambda=1$, $d=25.6\lambda$, $g=0.25d$, $h=0.05\lambda$, $lc=0.35\lambda$, $M=30$. (b) Angular power density averaged through 100 realizations, TM polarization. $\theta_i=-30°$, $\lambda=1$, $d=25.6\lambda$, $g=0.25d$, $h=0.05\lambda$, $lc=0.35\lambda$, $M=30$.</figcaption></figure>

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

<figure><img src="image_12.png" /><figcaption>Fig. 6. (a) Angular power density averaged through 50 realization, TE polarization. $\theta_i=-45°$, $\lambda=1$, $d=50.6\lambda$, $g=0.25d$, $h=0.16\lambda$, $lc=0.68\lambda$, $M=120$. (b) Angular power density averaged through 50 realizations, TE polarization. $\theta_i=-70°$, $\lambda=1$, $d=50.6\lambda$, $g=0.25d$, $h=0.16\lambda$, $lc=0.68\lambda$, $M=150$.</figcaption></figure>

is a remarkable agreement between our TE results and those of [20] (the TM case was not treated in this reference). In a second example, we retrieve the results obtained by the FDTD as published in [21]. In this case, the rms height and the correlation length are fixed to $h=0.16\lambda$ and $l_c=0.68\lambda$. Figures 6(a) and 6(b) show our computed averaged (over 50 realizations) bistatic scattering coefficient, in two cases of incidence: -45° and -70°, which correspond, respectively, to Figs. 7 and 8 of [21]. Here again one can observe an excellent agreement between results obtained from the C method and those obtained from the FDTD method.

Finally, let us take an example with high rms ($h=0.2\lambda$) and low correlation length ($l_c=0.2\lambda$) corresponding to a rather important roughness. The averaged bistatic scattering coefficient $\langle\sigma(\theta)\rangle$ (over 100 realizations) is presented in Figs. 7(a) and 7(b) for the TE and TM polarization cases, respectively. Unlike what is observed in the previous cases, where there is an apparent peak in the

specular direction, here an important amount of energy is backscattered around the angle of incidence $\theta=-30°$. This phenomenon is known as enhanced backscattering and has been studied extensively in the past; see, e.g., [22,23].

## 5. CONCLUSION

An extension of the C method to the case of diffraction of an arbitrary incident radiation by a PEC rough surface has been presented. The accuracy of the method has been investigated and comparisons with another approach (namely, the integral method) showed a very good agreement. Furthermore, an improvement concerning the problem of eigenvalue sorting was proposed through the computation of the eigenvectors' barycenters. Finally, it is important to emphasize that extending the present work to multilayered dielectric or metallic surfaces can be done mechanically and without any difficulty.
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 7. (a) Angular power density average through 200 realizations, TE polarization. $\theta_i=-30^\circ$, $\lambda=1$, $d=25.6\lambda$, $g=0.25d$, $h=0.2\lambda$, $lc=0.2\lambda$, $M=300$. (b) Angular power density average through 200 realizations, TM polarization. $\theta_i=-30^\circ$, $\lambda=1$, $d=25.6\lambda$, $g=0.25d$, $h=0.2\lambda$, $lc=0.2\lambda$, $M=300$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 8. Surface for Poynting's flux computation.</figcaption></figure>

## APPENDIX A: EXPRESSION OF THE INCIDENT WAVE IN THE TRANSLATION COORDINATE SYSTEM

The incident wave is described through the plane-wave spectrum concept:

$$ \Psi^i(x,y) = \int_{-\infty}^{+\infty} I(\alpha) \exp(-ik\alpha x) \exp(-ik\beta(\alpha)y) d\alpha, $$

$$ \alpha^2 + \beta^2 = 1. \qquad (A1) $$

The distribution of spectral amplitudes $I(\alpha)$ determines the nature of the beam. Thus if $I(\alpha)$ is a delta function $[I(\alpha)=\delta(\alpha-\alpha_0)]$, one recovers the plane wave under the incidence $\theta_0=\arcsin(\alpha_0)$. Our point here is to express the incident beam in the TCS from its known expansion (A1) in the Cartesian coordinate system. This will be done over the same discrete spectral range $(\alpha_n)$, which gives for Eq. (A1)

$$ \Psi^{\text{inc}}(x,y) = \Delta\alpha\Pi_{\Delta\alpha}(x) \sum_{n} I(\alpha_n)\exp(-ik\alpha_n x)\exp(-ik\beta_n y), $$

$$ \alpha_n^2 + \beta_n^2 = 1, \qquad (A2) $$

or in a slightly different form

$$ \Psi^{\text{inc}}(x,y) = \Delta\alpha\Pi_{\Delta\alpha}(x) \sum_{q} I(\alpha_q)\exp(-ik\beta_q y) \\ \times \sum_{n} \delta_{nq} \exp(-ik\alpha_n x), \qquad (A3) $$

where $\delta_{nq}$ is the Kronecker symbol. This last expression recalls that $(\delta_{nq})$ is the matrix of eigenvectors in the Cartesian coordinate system. Then, introducing the expression $y=v+a(x)$ in Eq. (A3) leads to

$$ \Psi^{\text{inc}}(u,v) = \Delta\alpha\Pi_{\Delta\alpha}(u) \sum_{q} I(\alpha_q)\exp(-ik\beta_q v) \\ \times \exp(-ik\beta_q a(u)) \sum_{n} \delta_{nq} \exp(-ik\alpha_n u). \qquad (A4) $$

That finally can be rearranged under the form

$$ \Psi^{\text{inc}}(u,v) = \Delta\alpha\Pi_{\Delta\alpha}(u) \sum_{q} I(\alpha_q)\exp(-ik\beta_q v) \\ \times \sum_{n} L_{n-pq} \exp(-ik\alpha_n u), \qquad (A5) $$

where

$$ \exp(-ik\beta_q a(u)) = \sum_{p} L_{pq} \exp(-ik\alpha_p u), \qquad (A6) $$

$$ L_{pq} = L_q(\alpha_p) = \text{TF}[\exp(-ik\beta_q a(u))](\alpha_p). \qquad (A7) $$

$L_{pq}$ represent the $p$ components on the Fourier basis of plane waves associated to $\beta_q$ and expressed in the curvilinear coordinate system. The incident field of Eq. (A5) may be rewritten to match the form of Eq. (11):

$$ \Psi^{\text{inc}}(u,v) = \Delta\alpha\Pi_{\Delta\alpha}(u) \sum_{q} I(\tilde{\alpha}_q)\exp(-ikr_q^\alpha v) \\ \times \sum_{n} \psi_{nq}^\alpha \exp(-ik\alpha_n u), \qquad (A8) $$

where $\tilde{\alpha}_q$ is defined in Eq. (16), $r_q^\alpha = \beta_q$, and $\psi_{nq}^\alpha = L_{n-pq}$.
## APPENDIX B: SPECTRAL POWER DENSITY

The total power of the scattered field is computed as the flux of the complex Poynting vector $J = 1/2 \Re[e[\mathbf{E} \wedge \mathbf{H}^*]$ through the surface represented in Fig. 8. In our study, only the $v$ component is of interest:

$$J_v = \frac{1}{2} \Re[e[E_u H_z^* - E_z H_u^*]]. \quad (B1)$$

More precisely, this gives in terms of the quantities $\Psi$ and $\Phi$

$$J_v = \begin{cases} \Re[e[\Psi(u,v)\Phi^*(u,v)]] & \text{for TE polarization} \\ -\Re[e[\Psi^*(u,v)\Phi(u,v)]] & \text{for TM polarization.} \end{cases} \quad (B2)$$

For the sake of conciseness, we derive the spectral power density expression for the case of TE polarization. The TM polarization case follows immediately from this latter. The flux of the complex Poynting vector is

$$P^s = \Re \left[ \int_z dz \int_R \Psi(u,v) \Phi^*(u,v) du \right]. \quad (B3)$$

Since the scattering surface and the incident field are invariant in the $z$-direction, $\int_R dz$ is constant and can be taken equal to 1. Then

$$P^s = \Re \left[ \int_R \Psi(u,v) \Phi^*(u,v) du \right], \quad (B4)$$

and by using Parseval's theorem

$$\int_R \Psi(u,v) \Phi^*(u,v) du = \int_R \hat{\Psi}(\alpha, v) \hat{\Phi}^*(\alpha, v) d\alpha, \quad (B5)$$

this gives

$$P^s = \eta \Re \left( \left[ \sum_{q/r_q \in R} R_q \exp(-ikr_q v) \sum_{p=-M}^{M} \psi_{pq} \right] \times \left[ \sum_{m/r_m \in R} R_m^* \exp(ikr_m^* v) \sum_{n=-M}^{M} \Phi_{nm}^* \right] \right) \delta_{np}, \quad (B6)$$

$$P^s = \eta \Re \left( \left[ \sum_{q/r_q \in R} \sum_{m/r_m \in R} R_m^* R_q \exp(-ik(r_q - r_m)v) \times \sum_{n,p=-M}^{M} \psi_{pq} \Phi_{nm}^* \right] \right) \delta_{np}. \quad (B7)$$

Since the eigensolutions are orthogonal, we have

$$P^s = \eta \Re \left( \sum_{q/r_q \in R} R_q R_q^* \sum_{n,p=-M}^{M} \psi_{pq} \Phi_{nq}^* \right),$$

$$P^s = \eta \sum_{q/r_q \in R} |R_q|^2 \Re \left( \sum_{n=-M}^{M} \psi_{nq} \Phi_{nq}^* \right). \quad (B8)$$

The spectral power density in $\alpha_q$ direction is defined by

$$\Delta P^s(\alpha_q) = |R_q|^2 \Re \left( \sum_{n=-M}^{M} \psi_{nq} \Phi_{nq}^* \right). \quad (B9)$$

Then

$$P^s = \eta \sum_{q/\alpha_q \in [-1;1]} \Delta P^s(\alpha_q). \quad (B10)$$

Similarly, we obtain the expression of the power impinging upon the surface:

$$P^{\text{inc}} = \eta \sum_{q/\tilde{\alpha}_q \in [-1;1]} |I(\tilde{\alpha}_q)|^2 \Re \left( \sum_{n=-M}^{M} \psi_{nq} \Phi_{nq}^* \right). \quad (B11)$$

Thus one can define a normalized expression of the spectral power density by

$$\Delta \tilde{P}^s(\alpha_q) = \frac{\Delta P^s(\alpha_q)}{P^{\text{inc}}}. \quad (B12)$$

The total reflectivity is then given by

$$\tilde{P}^s = \sum_{q/\alpha_q \in [-1;1]} \Delta \tilde{P}^s(\alpha_q). \quad (B13)$$

In the case of a nonpenetrable surface, the energy conservation is simply expressed by

$$\sum_{q/\alpha_q \in [-1;1]} \Delta \tilde{P}^s(\alpha_q) = 1. \quad (B14)$$

## APPENDIX C: ANGULAR POWER DENSITY OR BISTATIC SCATTERING COEFFICIENT

One can also define an angular power density $\sigma(\theta)$, which can be related to the spectral power density in considering $\alpha$ as a function of angle $\theta$:

$$\sigma(\theta) = \frac{d\tilde{P}^s(\alpha(\theta))}{d\theta} = \frac{d\tilde{P}^s(\alpha)}{d\alpha} \frac{d\alpha(\theta)}{d\theta}. \quad (C1)$$

For propagative waves, we have $\alpha = \sin(\theta)$, so that $d\alpha = \cos(\theta)d\theta$. Then

$$\sigma(\theta) = \frac{d\tilde{P}^s(\alpha)}{d\alpha} \cos(\theta). \quad (C2)$$

Finally, taking into account the discretization of the spectral domain, we get

$$\sigma(\theta_q) = \frac{\Delta \tilde{P}^s(\alpha_q)}{\Delta \alpha} \cos(\theta_q). \quad (C3)$$

This coefficient is often known as the bistatic scattering coefficient.

## REFERENCES

1. J. A. DeSanto and G. S. Brown, “Analytical techniques for multiple scattering from rough surfaces,” in *Progress in Optics*, E. Wolf, ed. (North-Holland, 1986), Vol. 23, pp. 3–62.

2. E. I. Thoros, “The validity of the Kirchoff approximation for rough surfaces scattering using a Gaussian roughness spectrum,” J. Acoust. Soc. Am. **83**, 78–92 (1987).
3. H. Faure-Geors and D. Maystre, "Improvement of the Kirchoff approximation for scattering from rough surfaces," J. Opt. Soc. Am. A **6**, 532–542 (1989).

4. J. A. Desanto, "Exact spectral formalism for rough surfaces scattering," J. Opt. Soc. Am. A **2**, 2202–2207 (1985).

5. J. M. Soto-Crespo and M. Nieto-Vesperinas, "Electromagnetic scattering from very rough random surfaces and deep reflection gratings," J. Opt. Soc. Am. A **6**, 367–384 (1989).

6. D. Maystre, "Electromagnetic scattering from perfectly conducting rough surfaces in the resonance region," IEEE Trans. Antennas Propag. **31**, 885–895 (1983).

7. M. Sailllard and D. Maystre, "Scattering from metallic and dielectric rough surfaces," J. Opt. Soc. Am. A **3**, 982–990 (1990).

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

9. B. Guizal, D. Barchiesi, and D. Felbacq, "Electromagnetic beam diffraction by a finite lamellar structure," J. Opt. Soc. Am. A **20**, 2277–2280 (2003).

10. J. Chandezon, D. Maystre, and G. Raoult, "A new theoretical method for diffraction gratings and its numerical application," J. Opt. **11**, 235–241 (1980).

11. L. Li and J. Chandezon, "Improvement of the coordinate transformation method for surface-relief gratings with sharp edges," J. Opt. Soc. Am. A **13**, 2247–2255 (1995).

12. R. Dusséaux and R. Oliveira, "Scattering of plane wave by 1-dimensional rough surface study in a nonorthogonal coordinate system," Prog. Electromagn. Res. **34**, 63–88 (2001).

13. G. Granet, K. Edee, and D. Felbacq, "A new curvilinear coordinate system based approach," Prog. Electromagn. Res. **41**, 235–250 (2003).

14. K. Edee, G. Granet, R. Dusséaux, and C. Baudier, "A curvilinear coordinate based hybrid method for scattering of plane wave from rough surfaces," J. Electromagn. Waves Appl. **15**, 1001–1015 (2004).

15. C. Baudier, R. Dusséaux, K. Edee, and G. Granet, "Scattering of plane wave by one-dimensional dielectric random rough surfaces: study with the curvilinear coordinate method," Waves Random Media **14**, 61–74, (2004).

16. R. Petit, *Electromagnetic Theory of Gratings* (Springer-Verlag, 1980).

17. L. Li and J. Chandezon, "Improvement of the coordinate transformation method for surface-relief gratings with sharp edges," J. Opt. Soc. Am. A **13**, 2247–2255 (1996).

18. J. P. Plumey, B. Guizal, and J. Chandezon, "Coordinate transformation method as applied to asymmetric gratings with vertical facets," J. Opt. Soc. Am. A **14**, 610–617 (1997).

19. G. Granet, J. Chandezon, J.-P. Plumey, and K. Raniriharinosy, "Reformulation of the coordinate transformation method through the concept of adaptive spatial resolution. Application to trapezoidal gratings," J. Opt. Soc. Am. A **18**, 2102–2108 (2001).

20. L. Tsang, J. A. Kong, K. Ding, and C. O. Ao, *Scattering of Electromagnetic Waves: Numerical Simulations*, Wiley Series in Remote Sensing, J. A. Kong, ed. (Wiley, 2001), pp. 113–161.

21. F. D. Hastings, J. B. Schneider, and S. L. Broschat, "A Monte-Carlo FDTD technique for rough surface scattering," IEEE Trans. Antennas Propag. **43**, 1183–1191 (1995).

22. A. Ishimaru, "Experimental and theoretical studies on enhanced backscattering from scatterers and rough surfaces," in *Scattering in Volumes and Surfaces*, M. Nieto-Vesperinas and J. C. Dainty, eds. (North-Holland, 1990).

23. A. A. Maradudin, J. Q. Lu, P. Tran, R. F. Wallis, V. Celli, Z. H. Gu, A. R. McGurn, E. R. Méndez, T. Michel, M. Nieto-Vesperinas, J. C. Dainty, and A. J. Sant, "Enhanced backscattering from one- and two-dimensional random surfaces," Rev. Mex. Fis. **38**, 343–397 (1992).