Beam implementation in a nonorthogonal coordinate system: application to the scattering from random rough surfaces.

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.


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 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][2][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][5][6][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 corru-gated waveguides and diffraction gratings. Recently, it has been successfully extended to diffraction of plane waves from nonperiodic surfaces [12][13][14][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].

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-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͑it͒, where 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 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 ⌿ satisfies the following set of partial differential equations that is valid in the two cases of polarization: where ⌿ = E z in the TE case and ⌿ = ZH z in the TM case, with Z denoting the impedance of vacuum. I is the identity operator, k is the wave vector, ‫ץ‬ u = ‫ץ‬ / ‫ץ‬u , ‫ץ‬ v = ‫ץ‬ / ‫ץ‬v , and ȧ =da / du . The second step consists in seeking plane-wave solutions under the form ⌿͑u,v͒ = exp͑− ikrv͒͑u͒, ‫ץ‬ v ⌿͑u,v͒ = exp͑− ikrv͒ r ͑u͒.
From these equations, it is worth noticing that the profile function comes into play through its derivative ȧ , 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][18][19]. Equation (2b) is then transposed into the spectral domain by Fourier transformation where f͑u͒ stands for ȧ ͑u͒, ͑u͒,or 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 where L 1 ͑␣͒ * ⌫ ͑␣͒ [resp. L 2 ͑␣͒ * ⌫ ͑␣͒] denotes the Fourier transform of the first (resp. the second) hand of Eq. (2b), * is the convolution product, and ⌫ ͑␣͒ = ͑ ͑␣͒ , r ͑␣͒͒ T .A t this stage, Eq. (4) is solved numerically by sampling the spectral domain into a set discrete points ͑␣ n ͒ over which the fields are evaluated: Then, according to Shannon's first theorem, the discretized functions can be interpolated by Numerically, only a finite sequence of sampled values is kept by first choosing a spectral increment ⌬␣ and a number N =2M + 1 of sample points so that the corresponding sequence is given through the relation ⌬␣ = n⌬␣ (with n ͓−M , M͔). Note that when ͉␣ n ͉ Յ 1 it corresponds to a propagating wave, whereas it corresponds to an evanescent one for ͉␣ n ͉ Ͼ 1. For a given ⌬␣, 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 ⌬␣.
After this operation, Eq. (4) becomes where ␣ = diag͑␣ n ͒, ȧ is a toeplitz matrix formed by the Fourier coeffcients of the derivative of the profile function ͑ȧ nm = ȧ n−m ͒, and I, is the identity matrix, (resp. r )isa column vector of ͑␣͒ components [resp. r ͑␣͒]onb 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 ⌿ ͑␣ , v͒ above the PEC surface 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 ͑ nq ± ͒ are rearranged accordingly. The expression of the total field in the spatial domain is obtained by the inverse Fourier transform of Eq. (9): Now it is time to write the boundary conditions that are satisfied by the tangential components of the fields. From the obtained z component ⌿͑u , v͒ we can determine the other tangential component ⌽͑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][11][12][13][14] ⌽͑u,v͒ = ‫ץ͓‬ v − ȧ ͑u͒͑‫ץ‬ u − ȧ ͑u͒‫ץ‬ v ͔͒⌿͑u,v͒,

͑12͒
Here, = i / kZ for TE polarization and =−iZ/ k for TM polarization.
Substituting the expression of ⌿͑u , v͒ from Eq. (11) into Eq. (12) gives with ⌽ +,− = ͓͑I + ȧ ȧ ͔ +,− r +,− ȧ ␣ +,− ͒, where ⌽ +,− (resp. +,− ) is a matrix whose elements are nq +,− (resp. 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 ⌿ or ⌽ 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 ⌬P s and the bistatic scattering coefficient whose expressions are derived in Appendices B and C:

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 q . This is done by observing numerically that the eigenvalues corresponding to propagative waves tend toward values that are exactly cos q when the number of samples is increased. Note that for these waves ␣ q = k sin 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 intrication 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: Numerically, it is found that each barycenter ␣ q tends toward a spectral point ␣ 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 Here, g denotes the waist of the beam at the origin and ␣ 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͑␣͒ = h 2 l c / 2 ͱ exp͑− ␣ 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͒ = ͚ n=1−P/2 P/2 b n exp͑i 2nx / 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 ȧ through the relation ȧ nm = ik⌬␣͑n − m͒b n−m .

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): 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 ⌬P s ͑␣ 0 ͒ in the specular direction as the number of samples N is increased, and we denote ⌬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: where Int͑␦͒ designates the integer part of ␦. Note that ͑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 ȧ ͑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 parameters remaining unchanged. In all the cases presented below, the surfaces are illuminated under the mean incidence i = −30°by a Gaussian beam with the following characteristics: = 1 and g = 0.25d. The first case consists of a surface with low h͑0.05͒ and l c ͑0.35͒. The shape of the surface of such a realization is shown in Fig.  2(a), and the accuracy function ͑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) 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. 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 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͒ and correlation length l c = 0.35. The computed averaged bistatic scattering coefficient ͗͑͒͘ 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  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 and l c = 0.68. 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͒ and low correlation length ͑l c = 0.2͒ corresponding to a rather important roughness. The averaged bistatic scattering coefficient ͗͑͒͘ (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 = −30°. This phenomenon is known as enhanced backscattering and has been studied extensively in the past; see, e.g., [22,23].

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.

APPENDIX A: EXPRESSION OF THE INCIDENT WAVE IN THE TRANSLATION COORDINATE SYSTEM
The incident wave is described through the plane-wave spectrum concept: The distribution of spectral amplitudes I͑␣͒ determines the nature of the beam. Thus if I͑␣͒ is a delta function ͓I͑␣͒ = ␦͑␣ − ␣ 0 ͔͒, one recovers the plane wave under the incidence 0 = arc sin͑␣ 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 ͑␣ n ͒, which gives for Eq. where ␣ q is defined in Eq. (16), r q − = ␤ q , and nq − = 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 I =1 / 2 Re͓E ∧ H*͔ through the surface represented in Fig. 8. In our study, only the v component is of interest: More precisely, this gives in terms of the quantities ⌿ and ⌽ Similarly, we obtain the expression of the power impinging upon the surface: Thus one can define a normalized expression of the spectral power density by The total reflectivity is then given by In the case of a nonpenetrable surface, the energy conservation is simply expressed by ͚ q/␣ q ͓−1;1͔ ⌬P s ͑␣ q ͒ =1. ͑B14͒

APPENDIX C: ANGULAR POWER DENSITY OR BISTATIC SCATTERING COEFFICIENT
One can also define an angular power density ͑͒, which can be related to the spectral power density in considering ␣ as a function of angle : For propagative waves, we have ␣ = sin͑͒, so that d␣ = cos͑͒d. Then ͑͒ = dP s ͑␣͒ d␣ cos͑͒. ͑C2͒ Finally, taking into account the discretization of the spectral domain, we get This coefficient is often known as the bistatic scattering coefficient.