In the diffraction shadow: Norton waves versus surface plasmon-polaritons in the optical region

Surface electromagnetic modes supported by metal surfaces have a great potential for uses in miniaturised detectors and optical circuits. For many applications these modes are excited locally. In the optical regime, Surface Plasmon Polaritons (SPPs) have been thought to dominate the fields at the surface, beyond a transition region comprising 3-4 wavelengths from the source. In this work we demonstrate that at sufficiently long distances SPPs are not the main contribution to the field. Instead, for all metals, a different type of wave prevails, which we term Norton waves for their reminiscence to those found in the radio-wave regime at the surface of the Earth. Our results show that Norton Waves are stronger at the surface than SPPs at distances larger than 6-9 SPP's absorption lengths, the precise value depending on wavelength and metal. Moreover, Norton waves decay more slowly than SPPs in the direction normal to the surface.

The confinement of the electromagnetic field associated to Surface Plasmon Polaritons (SPPs), and their intrinsic speed, make them very interesting candidates for their use in photonics [1,2]. Due to this, the study of the electromagnetic (EM) fields radiated by localized sources (like defects [3], nano-gratings or apertures [4]) placed on a surface has received a renewed interest in the last decade. This is an old problem, which was put at the forefront in the early 1900's by its possible relevance to the transmission of radio signals. The seminal works of Zenneck [5] and Sommerfeld [6] unveiled the existence of surface waves running along the Earth, which can be considered as a lossy dielectric. The interest in these works fainted after the realization that radio transmission does not occur via the exponentially damped surface modes, but through reflection at the ionosphere. Nevertheless, Norton subsequently showed that, in the long distance limit, radio waves decay algebraically at the surface [7]. This result triggered a debate on the range of validity of Zenneck-Sommerfeld and Norton waves in the radio regime that has propagated to our days (see Ref. [8] for more details and an historical account). Recently, advances in nanofabrication have allowed the scaling down of old radio devices into the optical regime [9]. Metallic surfaces are specially interesting because they support SPPs, which are surface EM modes strongly confined to the plane. The analysis of the surface EM fields created by a localized source in a metal surface has revealed the existence of a near-field region, extending for 3-4 wavelengths, where the field presents a complex dependence [10,11]. SPPs have been thought to dominate the EM field beyond this region. In this work we show that, irrespectively of the metal considered, the long-distance asymptotic limit of the EM field at metal surface is not the SPP but a different type of wave, which we denote as Norton waves (NWs) due to their reminiscence to those found in dielectric * Electronic address: lmm@unizar.es surfaces. We show the range of validity of SPPs and NWs and the distance and field amplitude after which the latter dominate.
Although we will show later how the obtained results apply to dipole sources, let us concentrate first on the EM fields emerging from a subwavelength slit, placed in an optically thick metal film. The film is back-illuminated by normal-incident p-polarized light with wavelength λ (i.e. the wavevector in vacuum is g = 2π/λ ). The frequencydependent dielectric constant of the metal is ε m . This system has been chosen for analytical simplicity (the full EM field can be expressed in terms of the magnetic field along the slit axes, H y (X, Z)) and, also because it is a configuration that has been amply studied both theoretically [12,13,14,15,16,17] and experimentally [18,19,20,21]. Figure 1 renders a snapshot of the radiated H y (X, Z) (computed with the FDTD method) for slit width of A = 100nm and λ = 540nm, for both a Perfect Electrical Conductor (PEC, characterized by |ε m | = ∞) and Au [22]. The choice of metal and wavelength is motivated for proof-of-principle purposes on the existence of NWs, but we will show later on that our results are applicable to other metals and frequency ranges. Our treatment fully takes into account the vectorial nature of the EM fields and, therefore, goes beyond the scalar approximations considered in other works [17,20].
It is apparent from Fig. 1 that the effect of a finite ε m is to strongly modify the radiation pattern close to the metal surface. Although we will provide expressions for the field everywhere, our main focus will be to characterize the fields within the diffraction shadow which, loosely speaking is the region where radiation from a slit in a real metal is strongly reduced with respect to the PEC case.
The Green's dyadic method is more suitable for an analytical study of this problem. Within this method, the field radiated at the point R = (X, Z) by a slit of width A is (see Supplementary Information S1 for the justification of this expression arXiv:0909.4328v1 [cond-mat.mes-hall] 23 Sep 2009 1: Snapshot of the magnetic field radiated by a back-illuminated subwavelength slit in an optically thick metal film. In the top panel the metal is treated as a Perfect Electrical Conductor, while the metal in the lower panel is Au. The wavelength is λ = 540nm. and its validation with numerical calculations) where the "Green's function" G(x, z) is the magnetic field generated by a dipolar source with the electric field pointing along the x-direction, placed at the metal interface, and δ is the skin depth for the metal. In this expression and throughout the paper all distances denoted by lower case letters are expressed in dimensionless units as x = gX, z = gZ and a = gA. Alternatively, given that the fundamental waveguide mode inside the slit is constant in the x-direction, G(x, z) can be seen as the magnetic field radiated by an infinitesimally thin slit. The angular spectrum representation of this function is: where q is the x-component of the wavevector (in units of g), D(q) = q zm / [2π (ε m q z + q zm )], q z = 1 − q 2 and q zm = ε m − q 2 [23]. The solution to this integral is not known in closed form. Fortunately, there are mathematical methods [28] for extracting its long-distance asymptotic expression, G asymp (x, z). The rigorous calculation for G asymp (x, z) is provided in the Supplementary Material S2 and, additionally, a simplified derivation will be given later on. But, before going into the mathematical details, let us now concentrate on the fields at the metal surface and give the result obtained: In this expression G SPP (x, 0) is the SPP contribution where q p = ε m /(1 + ε m ) is the SPP momentum, and C p = q 3 p / [2π(ε m − 1)] is the residue of D(q) at q p . The second term is As will be shown later, this term is the 2D optical analog in metal surfaces of the Norton wave [7] found in the study of the radio-wave radiation of point dipoles on lossy dielectric interfaces. Dimensionality accounts for the difference between the decay laws: x −3/2 (2D dipoles) and the x −2 (3D dipoles). The validity of Eq. (3) and the competition between SPPs and NWs is illustrated in Fig. 2, which shows the magnetic field at the surface radiated by an infinitesimally thin subwavelength slit, for Au at two different wavelengths. In each case, this figure renders the exact result (computed numerically from Eq. (1)) and the the SPP and NW contributions. For the cases considered in this figure, the asymptotic result given by Eq. (3) is virtually indistinguishable from the exact result even for X ≈ 3µm and it is not represented. The field is distance from the source along the metal X (μm) mainly SPP-like at the shorter distances, while NW dominates at sufficiently long distances from the source. Notice that the relative phase of the NW and the SPP contributions at the distance where their modulus are equal changes with wavelength. So, their destructive interference may lead to the cancellation of the field (as in Au, at λ = 540nm at X ≈ 12µm) or, if the cancellation is not complete, to the appearance of small oscillations in the total field amplitude of the field (as in the presented case of Au at λ = 567nm). It is worth noticing that similar oscillations were found in Scanning Near Field Optical Microscope experiments in Au [20], but their origen was unknown.
Beyond the particular examples presented in Fig. 2, the expression given by Eq. (3) is a good approximation for the field at the surface, at sufficiently long distances from the source. In order to quantify this statement, we define x a (β ) as the minimum distance such that |(G(x a , 0) − G asymp (x a , 0))/G(x a , 0)| < β . The discontinuous lines in Fig. 3 show the spectral dependence in the optical and telecom regimes of x a (0.1) for different metals, in units of the corresponding SPP absorption length. Given that the NW decays algebraically and the SPP exponentially with distance, at sufficiently large distances the NW is the main contribution to the field at the surface, for all metals and all wavelengths. The crossover from NW to SPP is represented in Figure. 3, which renders the spectral dependence of the distance at which the NW contribution is larger than the SPP one, x NW , for different metals. This distance strongly depends upon the dielectric permittivity of the metal, being smaller for very lossy metals, as Cu and Au in the region of inter-band transition (close to λ = 500 nm for both metals).
Undoubtedly, the existence of NWs in metal surfaces has passed unnoticed up to now due to their small amplitude. In order to characterize how much has the field decayed when the NW takes over, we consider the ratio |G(X NW , 0)/G(X = x NW , defined as the distance at which the amplitude of the NW is larger than the SPP one. The discontinuous lines render x a (0.1) the minimum distance at which G asymp gives a relative error of 10% with respect to the exact result. The inset renders the spectral dependence SPP absorption length for the considered metals.
λ /4, 0)| (the distance X = λ /4 has been arbitrarily chosen to give a representative reference in the near-field). In the optical regime, when the NW takes over the field has decayed by a factor ranging from 10 −2 for lossy metals (like Ni) to 10 −3 − 10 −4 for Ag. Therefore, the NW is not a good channel for sending information along the surface. Nevertheless, and given that estimations of the field at the surface far away from the source based on the decay of SPPs may be orders of magnitude wrong, NWs may have to be taken into account for precise analysis or design of experiments.
In order to show the origin of NWs, and its relation to other waves discussed in the literature, as creeping waves and SPPs, let us concentrate on the physical interpretation of the field radiated by the slit. Additionally, this will lead to a "poor man's" (yet correct) derivation of some of the main results. It is clear from Eq. (2) that a slit excites the whole range of diffraction modes (both radiative and evanescent) with an amplitude given by D(q), which can loosely speaking be understood as the density of EM modes with a given wavevector q at the slit position [25]. The standard treatment of G(x, z) in the far-field relies on the observation that, although all modes are always present, their contributions cancel out due to destructive interference whenever the phase Φ = qr changes rapidly. Thus, only the region in q-space where the phase presents an extremum contributes to the far field. For a given point (x, z) (or (r, θ ) in polar coordinates with θ defined as the angle from the normal to the surface), the extremum occurs at the condition (q/q z ) min = x/z, i.e. q min = sinθ . Expanding the integrand around this extremum leads to the "ray-optics" (RO) contribution This analytical result reproduces what was observed in Fig. 1: the magnetic field radiated by an infinitesimally thin slit in a PEC is isotropic, but the pattern in a real metal is strongly modified close to the surface, for angles such that cosθ 1/ |ε m |.
However, right at the surface the derivative of the phase Φ = qx never cancels and the saddle point approximation outlined above can not be directly applied. As the integral of the product of a smooth and rapidly oscillating function is very small, only the parts of the angular spectrum where D(q) changes rapidly in the scale of 2π/x will give a net contribution to the integral. For very small x all the "density of states" contribute. As x increases, the smooth long-q region of D(q) is progressively canceled out in the integral, which is eventually dominated by the strong (and rapid) contribution from the pole in D(q). The contribution of this pole gives the SPP field. Notice that, in a lossy metal, the density of states associated to the plasmon pole has a finite width, which causes the exponential decrease of the SPP amplitude with distance (characterized by the SPP propagation length l SPP = Im(q p ) −1 ).
The previous argument explains why the field at the surface is not the SPP for all distances and is expected to have a complex dependence with x. Recently, Lalanne and coworkers have termed "creeping wave" (CW) to the difference between the exact field and the approximation given by the SPP pole [11]. The numerical study of the CW has shown that it is a damped wave which, along the surface, oscillates with the free-space wavevector and decays after a few wavelengths. A point to notice is that despite the e ix dependence, the CW arises from the whole angular spectrum, not only from regions close to q = 1.
However, the SPP pole is not the sharpest feature of D(q): the derivative of D(q) diverges at the branch point q z = 0. This is illustrated in Fig. 4, which shows that |D(q)| has a kink at q = 1. The contribution to the integral from this kink is expected to be small but, as the kink can not be characterized by a typical width in q-space, it is not as strongly suppressed as the SPP contribution when integrated with an oscillatory function. In order to show that the kink originates the NW, it is convenient to integrate by parts G(x, 0). Then, from Eq.(2) we obtain G(x, 0) = (i/x) ∞ −∞ ∆(q) e iqx dq, with ∆(q) = q G (q). This representation has the advantage that the kink in D(q) transforms into a square root singularity (see inset to Fig. 4). The contribution close to q z = 0 can be retrieved by keeping the singularity but setting q z = 0 everywhere else, this is, by defining ∆ NW (q) = (ε m / 2π √ ε m − 1 ) (1/q z ). The inset to Fig. 4 renders the comparison between ∆(q) and ∆ NW (q), for a representative case. Of course, ∆ NW (q) is only a good approximation to ∆(q) close to q = 1, so its use for integration over the whole angular spectrum could seem unjustified. However for very large x this is valid, as only the 0 (x)/x. Recalling that this expression is only valid for large x, we substitute H (1) 0 (x) for its asymptotic value and obtain the result in Eq. (5). The motivation for our terminology on this type of wave is that the asymptotic term found by Norton, for the case of radio waves emitted by a dipole in a dielectric [7], also decays algebraically and originates from the angular spectrum close to q = 1. The SPP contribution can also be extracted from the previous representation by expanding ∆ q close to q = q p . As the NW and the SPP arise from different parts of the angular spectrum, their fields can be directly added up, leading to Eq. (3).
The relevance of NWs with respect to SPPs increases when we move away from the surface. Since for z = 0 the NW originates from q-values close to the light-line, its decay with distance to the surface is expected to be slower than the exponential decay of SPPs. In order to obtain the dependence of the Norton wave on both x and z, we have carried out the asymptotic analysis of the Green's function given by Eq. (1), using the general method described in Ref. [28]. The derivation and full asymptotic form can be found in the Supplementary Material S2. The result up to terms of the order r −1/2 was already presented in Ref. [16], where it was shown that the long-distance asymptotic expressions are excellent approximations even for distances as small as x = 1 (i.e., X = λ /(2π)). However, the asymptotic expression given by Ref. [16] misses some contributions of the order r −3/2 , so it can not be used for the problem discussed in this paper. We find that the expression for G(x, z) in the far-field (r 1) is (see Supplementary Material S2): where G SPP (x, z) = G SPP (x, 0) e iq pz z , with q pz = −1/ (1 + ε m ), G RO (x, z) is given by Eq. (6), and G NW (x, z), defined as the term that goes as r −3/2 , is given by Figure. 5 presents the comparison between the zdependence of the exact G(x, z) and G asymp (x, z), at X = 10µm. These results show that the asymptotic expression is very accurate. Also that, as expected, the NW contribution decays with distance to the surface much more slowly than the SPP one. In both chosen examples the SPP dominates right at the surface, but the NW takes over at a finite distance from it. However, at sufficiently large z the RO contribution always dominates. The inset to Fig. 5 shows, for Au at λ = 540nm, which of the three "asymptotic" terms (SPP, RO and NW) dominates in the X − Z plane. This inset shows that the NW is the largest contribution over a "stripe" close to the surface. As the region where the NW is dominant satisfies z x, the full expression for the NW given by Eq. (8) can be approximated by (for |ε m | 1) The comparison of this expression with that of the RO contribution allows for an estimation of the distance to the surface at which the crossover between NW and RO occurs, z NW . We obtain Notice that the algebraic decay of the NW with z, reflects that this wave arises from the interference of its constituent components.
To summarize, we have shown that, in the asymptotic limit of long distances to the source, the SPPs are not the main channel for EM fields at the metal surface. Instead, after a few SPP absorption lengths, Norton waves take over. This occurs for any metal and any frequency range. Norton waves decay much more slowly than SPPs both along the surface (as x −3/2 for 2D dipoles) and along the perpendicular direction.
Consider a p-polarized electromagnetic wave (with wavelength λ and wavevector g = 2π/λ ) incident onto a metal film with a subwavelength slit. The metal film is optically thick and extends from z = −W (where the EM field impinges) to z = 0 (the exit side). The dielectric constant of the metal is ε m . The slit has width A and we set the origin of the x-axis at the center of the slit.
According to Lippmann-Shwinger integral equation [26], the electric field at any point at exit side of the film (z > 0) is given by the following integral relation where E 0 (r) is the solution without the slit, ∆ε(r) = 1 − ε m in the volume occupied by the slit, V , and zero everywhere else.
In the case of an optically thick film the field E 0 (R) can be neglected at the exit side and the dyadicĜ E (R, R ) can be approximated by the one corresponding to a single metalvacuum interface. In order to obtain the magnetic field from Eq. (A1), we use the Maxwell equation H = (−i/g)∇ R × E and arrive at where we have passed to dimensionless distances r = gR.
The dyadicĜ H connects the magnetic field outside of the slit with the electric field inside the slit. For the considered twodimensional geometry, where only p-polarized waves are involved, the magnetic field H points along the y-direction. Assuming that the electric field inside the slit mainly points along the x-direction only the yz element of the dyadicĜ H needs be computed. We denote this element byĜ H yx (x, z; x , z ). Following Ref. 27, we find that ε m q z + q zm e iq(x−x )+iq z z−iq zm z (A3) where q is the dimensionless x-component of the wavevector, q = k x /g, q z = 1 − q 2 and q z = ε m − q 2 .
The integrant contains the exponential factor e −iq zm z , which decays in the distance of a skin depth δ = 1/Im(q zm ), which is of the order of a few tens of nm in the optical regime. Therefore, the integration limits in z can be extended to [−∞, 0]. Moreover, the variation of the Green's dyadic is much faster than that of the electric field inside the slit, so the electric field inside the slit can be approximated by its value at the distance z = −δ (this is obtained as the average distance to the surface, weighted by the exponential decay of the field). An additional advantage of using the field at a short distance inside the slit is that the numerical problems related to the treatment of corners are eliminated.
The integration over z can be performed in the following way .
We then obtain Since we are interested at the field close to the surface and away from the source (which, as we will show, arises from the region in the angular spectrum close to q = 1) and for optical frequencies (where ε m is large), this result can be related to the field radiated by a dipolar source G(x, y) at the metal surface: where G(x, z) = 1 2π dq q zm ε m q z + q zm e iqx+iq z z .
The relation between Eq. (A5) and Eq. (A7) can be easily seen by noticing that q zm ≈ √ ε m − 1 in the region of interest. While Eq. (A7) involves an additional approximation, we have preferred to work with the Green's dyadic for a dipole source (rather than with G slit ) due to its wider applicability to other problems. In any case, the differences for the slit case are minimal and the methods described in this paper could be straightforwardly applied to G slit .
In order to validate the expression (A7), we first performed full-vectorial computations using the Finite Element Method (FEM) of the magnetic field emerging from slits with the thicknesses 50 nm and 600 nm. The wavelength was chosen to be 576 nm and the metal of the film is gold. Then we computed the integral given by Eq. (A7) extracting the electric field inside the slit at z = −δ from the FEM calculations. The function G(x, z) was substituted by its asymptotic value, see the next section S2. The comparison is shown in Fig. 6. These results clearly show that Eq. (A7) is very accurate for subwavelegth slits (A << λ ) and even provide a good approximation for A ∼ λ . Consider the Green's function for a dipole placed at the metal surface given by Eq. (A8). In this section we sketch the asymptotic analysis performed, which has been done following the general method described in [28] for treating Sommerfeld integrals.
In this method, the integrant is first prolonged into the complex q-plane. Subsequently, the following changes of variable are performed q = sin φ and s = √ 2e i π 4 sin φ −θ 2 . In polar coordinates (x = r sin θ and y = r cos θ ) the integral takes the following form where the contour C in the s-plane corresponds to the real axis in the initial plane q.
Then the integrant is separated into singular and non-  When the pole is crossed during the transformation of the initial path to the steepest-decent one, the residue must be taken into account. singular parts where C p is the residue given by C p = ε √ ε 2π(ε 2 −1) √ 1+ε , and the position of the pole in the complex plane s is s p = √ 2e ıπ/4 sin( φ p −θ 2 ), where φ p = arccos −1 √ ε+1 . Then, after deforming the integral path to the steepest descent one (real axis in the plane s), the singular part yields the complementary error function iπC p e r(i−s 2 p ) erfc(−is p √ r) with the argument being the square root of the "numerical distance" introduced by Sommerfeld. The non-singular part of the integral can be expanded in Taylor series close to the saddle point s = 0 providing the infinite sum of the integrals of the Gaussian type. The final result reads G = iπC p e irq p erfc(−is p √ r) + e ir ∑ n∈even Γ( n+1 2 ) n!r n+1 2 d n Φ 0 ds n | s=0 . (B5) The function Φ 0 is composed of two parts: Φ and C p /(s p − s). Note that the part of the sum in (B5) coming from the term C p /(s p −s) coincides, up to a sign, with the asymptotic expansion of the complementary error function for large arguments without the first term (while the last term coincides with the residue contribution into the integral). This expansion reads where Θ is a Heaviside step function, and θ p is the angle defining the diffraction shadow, θ p = Re(φ p ) − arccos (1/ cosh[Im(φ p )]). This is the critical angle such that for θ > θ p we have Im(s p ) < 0 and the initial transformation of the integration path into the steepest-decent one leads to the crossing of the pole, so that the residue must be taken into account, see Fig. 7. An example of the critical angle limiting the diffraction shadow is represented in Fig. 1 of the article by dashed lines.
Retaining the residue contribution, G SPP , and the first two terms G RO (proportional to r −1/2 ) and G NW (proportional to r −3/2 ) in the sum coming from Φ we obtain an expression that is correct in the far-field, up to terms ∆G = O(r −5/2 ). The result can be rewritten as G = G SPP + G RO + G NW + ∆G. (B7) The plasmonic term is the "ray optics" term is G RO = e ir−iπ/4 √ 2πr · cos θ ε m − sin 2 θ ε m cos θ + ε m − sin 2 θ , the Norton wave is given by