Surface waves on two-dimensional rough Neumann surfaces

dispersion relation is obtained and solved for surface waves on a doubly periodic Neumann surface. A dispersion relation is obtained and solved for surface waves on a two-dimensional randomly rough Neumann surface. These results are applicable to a system of a liquid in contact with a rough hard substrate. • These are the first studies of surface waves in such systems. A planar surface on which a scalar wave satisfies the Neumann boundary conditions does not support a surface wave. However, a structured Neumann surface can support such a wave. By means of a Rayleigh equation for a scalar field in the region above the two- dimensional rough surface of a semi-infinite medium on which the Neumann boundary condition is satisfied, we derive the dispersion relation for surface waves on both doubly periodic and randomly rough surfaces. Dispersion curves for these waves on doubly periodic surfaces with three forms of the surface profile function are presented together with dispersion curves for surface waves on a two-dimensional randomly rough surface.


Introduction
In recent years interest has arisen in surface waves that propagate on impenetrable surfaces, in particular on perfectly conducting surfaces [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. These surface waves cannot exist if the perfectly conducting surface on which they propagate is planar: the surface must be periodically [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18] or randomly [1,15] rough. The interest in such waves on periodically rough perfectly conducting surfaces is due to the fact that their dispersion curves mimic those of surface plasmon polaritons on planar metallic surfaces, even though a perfect conductor has no electronic plasma that can display collective oscillations. In fact these surface waves are bound to the surface by its departure from planarity. Moreover, by varying the periodic structuring of a perfectly conducting surface, for example, by changing its period, or the heights and depths of its protuberances or indentations, or by changing the forms of these features, one can design these surface waves to possess specified properties, such as the frequency range in which they exist. Such ''tuning'' of the properties of these surface waves is not possible for surface plasmon polaritons on planar metallic surfaces. The existence of these waves on structured perfectly conducting surfaces has created the possibility of surface plasmon polaritons on periodically structured metallic surfaces in the gigahertz and terahertz frequency ranges, in which a metal is well represented by a perfect conductor. In these frequency ranges a surface plasmon polariton is so weakly bound to a planar metallic surface that it is in fact a surface current, not a surface wave. The existence of surface plasmon polaritons at gigahertz and terahertz frequencies is important in technological applications [19][20][21][22].
In this paper we study surface waves on a different type of impenetrable surface, namely a surface on which a scalar wave satisfies the Neumann boundary condition, i.e. a surface on which the normal derivative of the wave function of the wave vanishes. A physical situation in which this kind of surface wave can occur is a liquid on a hard substrate. It is straightforward to show that such a wave cannot exist if the surface of the substrate is planar. To see this let us denote the scalar field Ψ (x; t) in the region x 3 ≥ 0 by Ψ (x; t) = ψ(x|ω) exp(−iωt). Then the amplitude function ψ(x|ω) satisfies the Helmholtz equation where c is the wave speed in this region. The solution of this equation that propagates in a wavelike fashion, in directions parallel to the plane x 3 = 0, and vanishes at x 3 = +∞ can be written as where x ∥ = (x 1 , x 2 , 0), k ∥ = (k 1 , k 2 , 0), and (1. 3) The dispersion relation that connects the frequency of the wave ω with the two-dimensional wave vector k ∥ is obtained by using the solution (1.2) in the Neumann boundary condition on the surface x 3 = 0, namely In this way we obtain the condition (1.5) In order to obtain a nontrivial solution, i.e. one for which the amplitude A is nonzero, we have to require that β 0 (k ∥ ) = 0, which yields the dispersion relation ω = ck ∥ . The field ψ(x|ω) in the region x 3 ≥ 0 then becomes This function describes a surface skimming bulk wave that is not bound to the surface (β 0 (k ∥ ) ≡ 0).
The wave described by Eq. (1.6) is ''unstable'' [23], in that even a slight change of the boundary condition on the Neumann surface converts it into a surface wave or a surface shape resonance, both of which are bound to the surface.
In this paper we demonstrate that this is the case by studying theoretically the propagation of surface waves on a structured two-dimensional Neumann surface that is doubly periodically rough or randomly rough. The corresponding problem for a one-dimensional periodically or randomly corrugated Neumann surface was studied several years ago in Refs. [24] and [25], respectively, and will not be considered here.
In the work presented here we first derive the Rayleigh equation for the scalar field in the region above a two-dimensional rough Neumann surface. We then consider a doubly periodic surface and use the Rayleigh equation to obtain the dispersion relation of the surface waves it supports, which is then solved numerically. We also obtain the dispersion relation for surface waves on a two-dimensional randomly rough Neumann surface on the basis of the small roughness approximation to the Rayleigh equation. Thus, for both types of two-dimensional surface roughness we find that a Neumann surface supports surface waves.

The Rayleigh equation
The physical system we consider in this work consists of an isotropic and homogeneous medium characterized by a scalar wave speed c in the region x 3 > ζ (x ∥ ), and an impenetrable medium in the region x 3 < ζ (x ∥ ). The vector x ∥ defined by x = (x 1 , x 2 , 0) is an arbitrary vector in the plane x 3 = 0. The surface profile function ζ (x ∥ ) is assumed to be a single-valued function of x ∥ , that is differentiable with respect to x 1 and x 2 .
We seek a solution of the scalar wave equation where ω is the frequency of the wave in this medium. A time dependence of the form exp (−iωt) has been assumed for this field, but has not been indicated explicitly. We require that the solution of Eq. (2.1) satisfy the Neumann boundary condition on the surface x 3 = ζ (x ∥ ). Finally, we require that ψ(x|ω) vanish as x 3 → +∞.
The solution of Eq. (2.1) that satisfies the boundary condition at infinity can be written in the form We now introduce the representation so that (2.7) By taking the derivative of both sides of Eq. (2.6) with respect to x α (α = 1, 2), we obtain (2.8) On substituting Eqs. (2.6) and (2.8) with γ = α 0 (q ∥ ) into Eq. (2.5) the latter becomes The change of variable q ∥ + Q ∥ = k ∥ transforms Eq. (2.9) into (2.10) On equating to zero the k ∥ Fourier coefficient in this equation we obtain as the equation satisfied by the amplitude A(q ∥ ) (2.11) The solvability condition for this homogeneous linear integral equation is the dispersion relation that gives the frequencies ω for each value of the wave vector k ∥ . We see from Eq. (2.3) that only solutions of the dispersion relation for which the imaginary part of α 0 (q ∥ ) is positive correspond to waves localized to the surface x 3 = ζ (x ∥ ).
We now apply Eq. (2.11) to a study of surface waves on two-dimensional periodically and randomly rough surfaces.

A bigrating
We assume here that the surface profile function ζ (x ∥ ) is a doubly periodic function of ζ (x ∥ ), In Eq. (3.2) a 1 and a 2 are the two noncollinear primitive translation vectors of a two-dimensional Bravais lattice, while ℓ 1 and ℓ 2 are any positive or negative integers, or zero, which we denote collectively by ℓ. The area of the primitive unit cell of the two-dimensional lattice defined by Eq. (3.2) is a c = |a 1 × a 2 |.
We also introduce the lattice reciprocal to the one defined by Eq. (3.2), whose sites are given by where the primitive translation vectors b 1 and b 2 are defined by the conditions and h 1 and h 2 are any positive or negative integers or zero that we denote collectively by h.
For a surface possessing the periodicity property described by Eq. (3.1) the function I(γ |Q ∥ ) becomes . (3.6) In obtaining this expression we have used the result In order that the field in the region where k ∥ is the two-dimensional wave vector of the surface wave. When the results given by Eqs. (3.5) and (3.9) are substituted into Eq. (2.11) we obtain as the equation for the coefficients (3.10) The dispersion relation for surface waves on a bigrating on which the Neumann boundary condition is imposed is then obtained by equating to zero the determinant of the matrix of coefficients in Eq. (3.10).
The solutions of this equation are even functions of k ∥ , ω s (−k ∥ ) = ω s (k ∥ ), where s labels the solutions for a given value of k ∥ in the order of increasing magnitude. They are also periodic functions of k ∥ with the periodicity of the reciprocal lattice, The solutions can therefore be sought for values of k ∥ inside the first Brillouin zone of the bigrating, and inside the nonradiative region of ω and k ∥ values, defined by |k ∥ | > (ω/c).
In the present work we have calculated the solutions of the dispersion relation obtained from Eq. (3.10) for three forms of the surface profile function ζ (x ∥ ). We consider them in turn.

A square lattice of hemiellipsoids
In our first example the surface profile function is represented by a square array of hemiellipsoids on an otherwise planar Neumann surface. The primitive translation vectors of the square lattice are The primitive translation vectors of the corresponding reciprocal lattice are The surface profile function describing a square lattice of hemiellipsoids or radius R (<a) and height HR can be written as where (3.14) Because each hemiellipsoid has circular symmetry the dispersion relation satisfies the relation where S is a 2 × 2 real orthogonal matrix representative of any of the point group operations that leave the square lattice invariant. In the present case they constitute the point group C 4v . The property expressed by Eq. (3.15) has the consequence that all of the independent solutions of the dispersion relation are obtained if the wave vector k ∥ is restricted to the irreducible element of the two-dimensional first Brillouin zone. This is the region of the Brillouin zone that generates the entire zone when transformed by the application of the operations of the point group C 4v to it.
The functionÎ(γ |G ∥ ) for the surface profile function defined by Eq. (3.13)-(3.14) is given by [26] where J ν (x) is a Bessel function of the first kind of order ν, and Γ (x) is the gamma function. For the values of the parameters used in the numerical calculations based on the results obtained in this section it was necessary to include only the first few (≈15) terms in these rapidly convergent series.
To solve the dispersion relation obtained from Eq. (3.10) the infinite sum over G ′ ∥ has to be truncated. This was done by that increasing |H| depresses the dispersion curves for lattices of protuberances and indentations, but more strongly for the latter type of surface than for the former.
It is an open question as to whether increasing the magnitude of H further will lead to a second, higher frequency, branch of the dispersion curve entering the nonradiative region. Numerical instabilities in the present calculations for values of |H| larger than unity have prevented us from answering this question.

An orthogonal superposition of two sinusoidal gratings
The second form of the surface profile function we consider is given by

An orthogonal superposition of two symmetric sawtooth gratings
In our third example the doubly periodic surface is an orthogonal superposition of two identical symmetric sawtooth gratings. The surface profile function ζ (x ∥ ) has the form given by Eq. (3.13), where the function s(x ∥ ) is now defined by The functionÎ(γ |G ∥ ) for this surface profile function iŝ We believe that the small value of N max = 4 used in obtaining the dispersion curve when H/a = 4 is due to the nonanalytic nature of the surface profile function (3.20). It is known that the Rayleigh hypothesis is invalid for a onedimensional periodic surface defined by a nonanalytic profile function [27]. It is assumed, without proof, that the same result holds for a doubly periodic surface defined by a nonanalytic surface profile function. When a one-or two-dimensional nonanalytic periodic surface profile function is expanded in a finite number of plane waves, it is an analytic function of the coordinates, for which the Rayleigh hypothesis is valid when its amplitude to period ratio is sufficiently small. However, as more and more plane waves are included in the expansion of the surface profile function, i.e. when N max is increased, a point is reached where the nonanalytic nature of the surface profile function begins to manifest itself, and the Rayleigh hypothesis breaks down. A further increase in N max simply worsens the result. Thus the convergence of the expansion of the scalar field in plane waves has an asymptotic character. In the calculations based on the profile function (3.20) the optimal value of N max is 4 for the rough surface with H/a = 0.4, while it is N max = 12 for the smoother surface with H/a = 0.2.
Similar considerations apply to the results presented in Fig. 1 for another surface defined by a nonanalytic profile function.

A randomly rough surface
We now assume that the surface profile function ζ (x ∥ ), in addition to being a single-valued function of x ∥ and differentiable with respect to x 1 and x 2 , constitutes a stationary, zero-mean, isotropic, Gaussian random process defined by In what follows we will need the Fourier integral representation of ζ (x ∥ ), The Fourier coefficientζ (Q ∥ ) is also a zero-mean Gaussian random process that is defined by The function g(|Q ∥ |) entering this equation is the power spectrum of the surface roughness, and is defined by To obtain the dispersion relation for surface waves on a randomly rough Neumann surface we will use the small roughness approximation. This consists of expanding the function I(α 0 (q ∥ )|k ∥ − q ∥ ) in Eq. (2.11) in powers of the surface profile function and retaining only the first two terms: Eq. (2.11) can then be rewritten in the form (4.7) The surface profile function entering Eq. (4.7) is a random process. Consequently the solution A(k ∥ ) of Eq. (4.6) is a random process. Just as ζ (x ∥ ) is defined by the moments of its probability density function, so can A(k ∥ ) be defined by the moments of its probability density function. Of these a particularly important one is the first moment ⟨A(k ∥ )⟩, which describes the propagation of the mean wave across the randomly rough surface.
To obtain the equation satisfied by ⟨A(k ∥ )⟩ we introduce the smoothing operator P that averages every function to which it is applied over the ensemble of realizations of the surface profile function [28]. We also introduce the complementary operator Q = 1 − P, that produces the fluctuating part of every function on which it acts. On applying these operators to Eq. (4.6) we obtain the pair of equations (4.9) When we use the results that PV (k ∥ |q ∥ ) = 0, because ζ (x ∥ ) is a zero-mean random process, and that QA(r ∥ ) is of order ζ (x ∥ ), Eqs. (4.8)-(4.9) simplify to On combining Eqs. (4.10) and (4.11) we obtain as the equation satisfied by the mean amplitude ⟨A(k ∥ )⟩ The result that ⟨A(k ∥ )⟩. (4.14) The dispersion relation for surface waves on a randomly rough two-dimensional Neumann surface is therefore .
The power spectrum g(|k ∥ − q ∥ |) can be written in the form where where If we now take Eq. (2.4) into account, we can rewrite Eq. (4.18) as (4.20) We now make the Ansatz [Ag 0 (k ∥ |q ∥ ) + Bg 0 (k ∥ |q ∥ ) + Cg 2 (k ∥ |q ∥ )] (4.23) Since according to Eq. (4.21) the deviation of (ω/c) from k ∥ is of O(δ 4 ) we can replace (ω/c) in the integrals in Eq. (4.23) by k ∥ in obtaining ∆(k ∥ , ω) to the lowest order in δ. The result is  . (4.26) Since from Eq. (2.3) we see that the imaginary part of α 0 (q ∥ ) must be positive for the wave to be bound to the surface, it follows from Eq. (4.26) that ∆ 1 (k ∥ ) must be positive.
In obtaining numerical results for ∆ 1 (k ∥ ) and ∆ 2 (k ∥ ) we assume that the normalized surface height autocorrelation function W (|x ∥ |) has the Gaussian form (4.27) The coefficient function g m (k ∥ |q ∥ ), Eq. (4.17), in this case is given by (4.28) where I m (x) is a modified Bessel function of the first kind of order m. The expressions (4.25a) and (4.25b) for ∆ 1 (k ∥ ) and If we introduce the definition x = k ∥ a, and make the changes of variables q ∥ = (x/a) cosh θ in Eq. (4.29a) and q ∥ = (x/a) sin θ in Eq. (4.29b), we find that (4.31b) In the long wavelength limit x → 0 the universal functions d 1 (x) and d 2 (x) take the forms [4]   If we combine the results given by Eqs. (4.21), (4.24), and (4.30), we find that the dispersion relation for surface waves on a two-dimensional, randomly rough, Neumann surface can be written in the form where Thus, in the long wavelength limit the surface wave dispersion relation is The function α 0 (k ∥ ) obtained from Eqs. (4.26), (4.30), and (4.32) is given by (4.36) In Fig. 4 we have plotted d 1 (x) and d 2 (x) as functions of x, while the functions ω 1 (x) and ω 2 (x) are plotted in Fig. 5. The functions d 1 (x) and d 2 (x) are seen to be positive for all values of x. One of the consequences of these results is that the inverse decay length of the wave into the vacuum, Imα 0 (k ∥ ), is always positive. Thus the wave is bound to the surface for all values of x = k ∥ a. In the long wavelength limit Imα 0 (k ∥ ) is proportional to k 2 ∥ (Eq. (4.36)). This wave is therefore weakly bound to the surface in this limit, but it is bound nonetheless.
A second consequence of the positivity of d 1 (x) and d 2 (x) is that Imω(k ∥ ) = −ck ∥ (δ/a) 4 (x 8 /4)d 1 (x)d 2 (x) is negative for all x. This means that the surface wave is attenuated as it propagates along the randomly rough surface for all values of x.
In the long wavelength limit Imω(k ∥ ) is proportional to k 6 ∥ (Eqs. (4.33), (4.34b)) which, in view of Eq. (4.33) means that it is proportional to the sixth power of its frequency. The explanation for this dependence begins with the observation that the frequency dependence of Rayleigh scattering is ω d+1 , where d is the dimensionality of the scatterer. The protuberances and indentations of the surface responsible for the scattering of the wave in the present case are three dimensional, since they are defined by the equation x 3 = ζ (x ∥ ). The Rayleigh scattering law therefore gives us a ω 4 frequency dependence of the scattering rate in the low frequency long wavelength limit. The remaining factor of ω 2 arises because the amplitude of the surface wave decays exponentially with distance into the region x 3 > ζ (x ∥ ), increasing thereby the volume within which its interaction with the surface roughness occurs. We have seen that the decay length of the wave into the region x 3 > ζ (x ∥ ) is proportional to the square of its wavelength parallel to the surface (Imα 0 (k ∥ ) is proportional to (k ∥ a) 2 in this case: see Eq. (4.36)), and this dependence translates into a factor of ω 2 multiplying the Rayleigh scattering rate .
These surface waves display wave slowing (ω 1 (x) > 0) for small x. However, ω 1 (x) becomes negative for values of x greater than a critical value x c ∼ 1.8. Thus, although in the range x > x c the wave is still bound to the surface, Imα 0 (k ∥ ) > 0, its phase velocity is greater than the speed c of scalar waves in the region x 3 > ζ (x ∥ ), so that the surface wave turns out to be in the radiative region of k ∥ and ω values.
We also note that despite the difference between the physical systems studied, and the difference between the equations to be solved, the dispersion relation for a surface wave on a two-dimensional randomly rough Neumann surface turns out to be the same as the dispersion relation for surface electromagnetic waves on a two-dimensional randomly rough perfectly conducting surface in contact with vacuum [15]. We have included the resulting dispersion curve for completeness.

Conclusions
The results presented here have been obtained for a simple system that demonstrates how the structuring of a surface can make it support a surface wave that does not exist in the absence of the structuring. The structuring can be periodic or it can be random. The resulting modes are dispersive because of the presence of a characteristic length in the problem: the lattice parameter in the case of the bigrating, and the transverse correlation length in the case of the randomly rough surface. They display wave slowing. As the height of the protuberances, or depths of indentations of doubly periodic surfaces and the rms height of a randomly rough surface is increased, the dispersion curve of the corresponding surface waves bends away from the dispersion curve ω = ck ∥ into the nonradiative region of ω and k ∥ values at smaller and smaller values of k ∥ . The frequency range within which these waves exist shrinks correspondingly. Increasing the period of the bigrating for fixed values of the height of the protuberances or the depth of the indentations, or decreasing the value of δ/a for the randomly rough surface, has the opposite effect: the surface wave dispersion curve moves closer to the dispersion curve ω = ck ∥ . In general, the dispersion curves of these surface waves can be modified in desirable ways by suitable choices of the parameters that define the roughness of the doubly periodic or randomly rough surface.