Optical surface waves over metallo-dielectric nanostructures: Sommerfeld integrals revisited

The asymptotic closed-form solution to the fundamental diffraction problem of a linear horizontal Hertzian dipole radiating over the metallo-dielectric interface is provided. For observation points just above the interface, we confirm that the total surface near-field is the sum of two components: a long-range surface plasmon polariton and a short-range radiative cylindrical wave. The relative phases, amplitudes and damping functions of each component are quantitatively elucidated through simple analytic expressions for the entire range of propagation: near and asymptotic. Validation of the analytic solution is performed by comparing the predictions of a dipolar model with recently published data


Introduction
Sparked by the discovery of enhanced optical light transmission through subwavelength holes pierced in a metal film [1], considerable interest has grown in the last decade for the guiding and manipulation of light in metallic nanostructures. A central role is attributed to electromagnetic surface waves (SW) that are launched by light diffraction on a structural nano-defect and subsequently propagated along the metallo-dielectric interface. While novel nano-devices operating on these surface-bound vectors are envisioned in many areas of science and engineering, fundamental-level research is still being performed in order to completely assess and control their intrinsic properties. In this respect, surface plasmon polaritons (SPP) have been identified early on to occupy a predominant role, especially in the asymptotic propagation regime. The SPP, which is a guided electromagnetic wave with fields evanescently extending on either side of the interface, represents a solution stemming from Maxwell's equations provided the metal substrate has a finite conductivity and the magnetic field is polarized transverse (TM) to the plane of incidence. An alternate type of SW denominated "composite diffracted evanescent wave" (CDEW), and characterized by a field damping scaling as 1/x with propagation distance x, has also been proposed [2]. The CDEW model is based on a number of approximations. First of all, it derives from a scalar-wave theory which ignores field polarization, and second, it assumes an opaque and infinitely thin screen as the boundary condition, thus precluding any likely involvement of the SPP mode. A key conceptual difference between the CDEW and the classical SPP model is that the former allows several propagating surface modes to be generated by diffraction while the latter only admits the SPP. Both models were applied to interpret recent experimental investigations of the near field with relative success [3][4][5]: the CDEW is ostensibly better suited in the immediate vicinity of the source whereas the SPP is most accurate further away from it. The indication that each theory presents complementing pictures of the same phenomenon was made apparent by Lalanne and Hugonin [6] when they suggested that the composite diffracted surface wave essentially consists of two components: a SPP and a radiative "creeping wave" characterized by a damping scaling with 1/x 1/2 . In their demonstration, the authors used a rigorous Green's function formalism to describe the field radiated by a line source over a metallo-dielectric half-space. However, no closed-form solution for the creeping wave was provided. New experimental results [7] taken from direct measurements of the surface near field have lent further credence to the thesis of a SW dual composition.
In this theoretical contribution, we demonstrate that the scattering process of incident TM-polarized light by a one-dimensional subwavelength nano-defect on an otherwise flat metal surface is generically modeled through the fundamental diffraction problem of a horizontal Hertzian dipole radiating over the metallo-dielectric interface. Starting directly from Maxwell's equations and enforcing the proper boundary conditions, we write down the exact Sommerfeld-type integral for the electric field perpendicular to the interface. The integral is then rigorously solved via the modified method of steepest descents to obtain the asymptotic closed-form solution. We explicitly show that the total diffracted field along the surface is composed of two distinct contributions as predicted in [6]: a short-range radiative wave and a long-range SPP. The amplitudes, phases and damping functions of both components are quantitatively revealed through simple analytic expressions. Finally, we propose a double-dipole model to characterize the near-field interactions between two nanoobjects and show the excellent agreement between our model's predictions and the experiments, thus validating the analytical solution in the same stretch.

Derivation of the closed-form solution
It was demonstrated that the scattering of incident TM-polarized light on a subwavelength metallic slit induces oscillating electric charges around the sharp slit edges that emulate a horizontal electric dipole [8,9]. It can be argued that this distinctive feature equally applies to the subwavelength groove if it is considered as a partially filled slit. Hence, the physical problem of light scattering on a nano-slit, or nano-groove, can be represented by the basic diffraction problem of an infinitely small horizontal electric dipole radiating over the interface. That assumption is corroborated by other analytical and experimental investigations [10,11]. From our calculations, this approximation is fairly accurate if the width w of the given nano-object is smaller than a half-wavelength, w< /2, and if observation points are located at a distance x from the source larger than x> / . The remaining discussion in this Section follows the contemporary mathematical treatment derived by R. E. Collin for the distinct case of a 3D vertical point dipole radiating over the earth's surface [12]; which was the original diffraction problem famously addressed by Arnold Sommerfeld one century ago. We consider the 2D linear electric dipole oriented along the x-axis and located in the halfspace z 0 at a height h from the interface separating two semi-infinite nonmagnetic and (Fig. 1). The infinitesimal line dipole is described by its surface current density J x = (x) (z-h) normalized to unit electric moment. We will not cover here the details of the straightforward though lengthy derivation of Eq. (1), which can be found in textbooks [13], but rather briefly summarize the procedures therein. After enforcing the boundary conditions at the interface for the E x , E z and H y fields (TM-polarization) in Maxwell's curl equations then taking the inverse Fourier transform of the fields and choosing the appropriate Green's function, one derives the Sommerfeld-type integral for the normal electric field in region 1 (z 0):

Expression of the Sommerfeld integral for a linear horizontal Hertzian dipole
where k is the lateral component of the total wavevector which is related to the normal component in each region "i" via i = k i 2 k 2 with k i = k 0 i and k 0 = 2 . A time dependence in e i t is assumed in the present analysis and suppressed throughout.

Asymptotic solution to the integral via the modified method of steepest descents
The integral in Eq.
(1) is separated in two parts, E 1z =K (I 1 + I 2 ), with K= μ 0 /2 , and I 1 and I 2 respectively corresponding to the first and second term inside the brackets. We first consider the integral I 1 and impose 1 and 2 to have positive imaginary parts in order for the field to be bounded at infinity. In the corresponding first quadrant of the complex k-plane, we discard the branch cut running from the point k=k 2 since its contribution -characterized by an attenuation factor exp Im k 2 { } x ( ) that drops to nearly zero within a quarter-wavelength distance x -is negligible for the optical frequencies of interest compared to that of the branch cut running from k=k 1 . There is also a pole singularity in the denominator of I 1 located at ±k p = ±k 0 p where p = 1 2 1 + 2 ( ) defines the SPP's "effective permittivity". The proper pole ( +k p ) lies very close to the branch point k 1 and we will see later on that this aspect entails particular considerations. We make the simplification 1 =1 such that k 1 =k 0 and perform the successive transformations k = k 0 sin , x = R 2 sin and z + h ( ) = R 2 cos , where R 2 is the distance from the mirror image of the dipole in region 2 to the observation point. These procedures enable one to express I 1 as the integral over an angular spectrum of plane waves: The complex variable = + i defines the angle between the direction of propagation and the x-axis. The saddle-point is found by setting the derivative of the argument in exponential equal to zero, thus yielding = where is real. The location of the pole in the -plane is given by k p = k 0 sin p . As a numerical example, at the excitation wavelength =852nm and . As shown on Fig. 2, the original integration contour C is then deformed into the steepest-descent contour (SDC) whose path, cos ( ) cosh = 1 , is shifted to pass through the saddle-point at = 2 where the highest accuracy is assigned to observation points along the surface. Since the pole p has been crossed by the path in the process, and is positioned below the SDC and very close to the saddle-point (see Fig. 2), the pole is captured and its contribution must be accounted for with a residue [12,14]. To carry out the steepest-descent integration in Eq.
(2), we first perform the change of variable , such that the integral I 1 is rewritten: where f ( ) = g ( ) h ( ) and: Because a pole singularity lies very close to the saddle-point = 2, the classical steepestdescent method cannot be applied directly. To circumvent this problem, the procedure that we will employ here is a modified steepest-descent technique described in [12]. In this approach we separate the pole and analytical contributions, respectively f P ( ) and f A ( ) , by subtracting out the pole term from the main integrand: where we keep only the first term and neglect the remaining higher even-order correction We separate the integral I 1 in two parts, I 1 = I P + I A , where I P and I A respectively denote the pole and analytical contributions from f P ( ) and f A ( ) . We first consider the analytical part: Since the function f A 0 ( ) has no dependence in , it can be removed from the integrand.
Upon evaluating the remaining kernel, Next we compute the pole contribution: where P ( ) = d e 2 2 p ( ) + and = k 0 R 2 2 . Special care must be taken to evaluate P ( ) since it can be defined by different functions depending on whether the pole is located in the lower-half (Im{ p }<0) or upper-half (Im{ p }>0) -plane. We can find a solution valid over the entire -plane by first evaluating P ( ) for Im{ p }>0 and then assuming Im{ p }<0 in the resulting integral expression. In which case we obtain t with t = i p , such that: where: U R 2 , The complete air-side normal E-field is written E 1z = K I P + I A + I I + I D ( ) . Upon adding the terms I A and I I together and performing some algebraic manipulations, one obtains the following general closed-form expression: Equations (15) having essentially a free-space cylindrical nature. The boundary wave lies in the geometrical shadow so as to compensate for the discontinuity in the geometrical-optics field across the planar interface [15]. We refrain from using the denomination "creeping wave", which was previously chosen in [6], to identify this radiative wave since it generally refers to a surface mode propagating in the geometrical shadow around convex surfaces. The contributions of the SPP and BW respectively correspond to the first and second terms in Eq. (19). In that

Direct dipole field
Reflected field SPP field Boundary wave field expression, the "mismatch parameter" k = k p k 0 ( ) defines the difference between the evanescent and free-space wavevectors. We normalize the total surface wave field, which we denote E sw , by taking out the factor μ 0 Ae i from Eqs. (17)-(18) and moving the x root from the denominator to the numerator in order to eliminate any sign ambiguity when x < 0: The + superscript in Eq. (19) indicates that the SW originates on the right of the source dipole and propagates in the +x direction. An identical SW is excited in the -x direction whose field is the complex conjugate: E sw x ] , which is accurate for k 0 x>2. This outcome is not incidental because the steepest-descent method -which is asymptotically exact -was used to obtain the solution. Thus in principle, the BW could be written in "exact" form with the

Fresnel diffraction effects and the asymptotic propagation regime
The inspection of Eq. (19) reveals that the SPP excited through diffraction at the interface is not pure but rather modulated by a slowly oscillating envelope U(x) owing to Fresnel diffraction effects [Eq. (12) evaluated at = 2]: The first term on the right-hand side of Eq. (20) defines the classical SPP mode whereas the second erfc term describes the fringe pattern generated by the process of coupling incident homogeneous light into the evanescent mode via scattering. As a note, Eq. (20) can equivalently be expressed using the conventional complex Fresnel integral, ) . As shown on Fig. 3, the modulation function Re{U(x)} evolves from U(0)=0.5 at the origin before reaching the first peak at some distance x 1 and then oscillates onward with a large distinctive period L = 2 k to lower amplitudes before asymptotically growing towards ever higher amplitudes due to the positive imaginary component inside the argument of erfc. The location of the m-th peak from the origin is predicted by x m = 2m 5 4 ( ) k with m=1,2,3…; which yields x 1 =20.8μm, x 5 =242.2μm and x 9 =463.6μm for the example of Fig. 3. The complex modulated SPP field is superposed on the monotonously decaying BW. As shown on the same figure, the sum of the two co-propagating surface waves produces a much weaker modulation, also with period L, of the total amplitude. Clearly, this beating behavior -more discernible at higher optical frequencies -is related to the wavevector mismatch between the evanescent and homogeneous modes, for when k 0 in the PEC limit (perfect electric conductor) all fringing effects disappear. The most salient feature is the first peak at x 1 = 3 4 k ( ) which may represent the only detectable Fresnel modulation of the total field in practice due to the typically low visibility of the following fringes. To our best knowledge, this phenomenon of field pattern created by the interference of the two co-propagating surface waves -one of which is the evanescent SPP and the other a radiative wave -both originating from the same source, has not yet been reported in the literature.

Characterization of the field at the origin and in the near-zone propagation regime
The first point away from the origin where the real part of the oscillating BW field takes a zero value is approximately located at x= /8, yielding: E bw x /8 the real part is well-behaved; while for x< /8 it rapidly diverges towards negative values to become singular at the limit x 0 + . By contrast, the SPP mode is well defined at x=0: These previous remarks suggest that the SW field at the origin can be reasonably approximated by x 1 The BW's amplitude scales as 1/(4 | k| x) 1/2 in accordance with the 1/x 1/2 damping predicted in [6]; while the SPP modulus follows Re{U(x)} exp(-Im{k p }x). The "critical distance" where the moduli of the SPP and BW coincide is approximately located at x c =4 /(67| k|). For the conditions =852nm and 2 = 33.22+1.17i, the predicted value is x c =1.64μm. For propagation lengths below the critical distance (0<x<x c ) the BW provides the dominant contribution whereas the long-range SPP is the main vector thereafter (Fig. 4). The 1/e intensity decay lengths of the SPP and BW are respectively L spp =1/(2 Im{k p }) and L bw =e/(4 | k|). For the same previous conditions, the BW's characteristic decay length L bw =1.9μm is much shorter compared to the SPP's L spp =122μm. As a consequence, the effective value (k sw =2 / sw ) of the composite SW's wavevector (or effective index: n sw = k sw /k 0 ) is significantly affected by the rapid decay of the BW inside the near-zone regime (0<x<2L bw ) as shown in Fig. 4. The behavior of the near field depicted by the analytical solution in the source origin's vicinity, brings a clear and cohesive physical explanation to the reported accounts of transient phenomena in the first few wavelengths of propagation [3,4,7,10].

Dipolar model of the near-field interactions between nano-objects
As indicated earlier, when light is incident on a nano-object (nano-slit or nano-groove) with dimensions smaller than the wavelength of light, the scattered far-field radiation from the nano-object can be considered as originating from an infinitesimal horizontal electric dipole (HED) located on the structure. This assumption has been rigorously validated in Green's tensor numerical calculations by Lévêque et al. [10]. Hence in principle, many diffraction problems involving nanostructured surfaces can be modeled by replacing the individual nanoobjects with point-sized dipoles. Indeed we demonstrate in this Section that recently reported experiments using groove-slit and double-slit configurations are accurately modeled by placing radiating linear HEDs at the corresponding sites of the surface nano-defects.

Groove-slit transmission
We first consider the groove-slit experiment performed by Gay et al. [3] and theoretically investigated by Lalanne and Hugonin [6]. The main interactions in this setup are described by the simple model illustrated in the diagram below (Fig. 5) where the groove-slit center-tocenter separation distance is controlled by the variable d. An incident x-polarized plane wave of amplitude E i is scattered on both 100nm-wide groove and slit as HEDs. The fraction of incident light that is coupled by the groove into a z-polarized +x-directed SW is determined by the multiplication factor E sw + (d) where is the amplitudecoupling coefficient and E sw + (d) describes the field propagation of the excited SW. The diffraction-launched SW then impinges on the adjacent nano-object (i.e. slit) where it is partially reflected and re-radiated. The re-radiated x-polarized component E x sw that is coupled into the slit is defined by the factor E sw + (0), where backconversion reciprocity (SW-toradiation and vice-versa) is assumed for the coupling coefficient and the factor E sw + (0) accounts for the generation of a new SW along the slit's left-wall. Indeed, as demonstrated in [9] a new SW originates from the slit's top-left corner and whose initial value, ( ) e i 2 , introduces a /2 phase-shift in magnitude. This intrinsic /2 phase-shift between the SW generated at the groove and the directly incident light is consistent with earlier experimental [3,4] and theoretical studies [10]. The x-polarized light normally incident on the slit similarly generates two SWs: one in the forward and another in the backward x-direction. As described in [9], the SW launched from the groove does not directly interfere with the normally incident light since they have orthogonal polarizations; instead it is the re-radiated x-polarized component E x sw that interferes with the transmitted component E 0 in the slit. The complete process of interference at the slit, as described by this simple model, is expressed by E 0 = t 0 E i respectively denote the x-polarized transmitted components relating to the left incident SW and the normal incident plane wave. The total transmitted intensity I is a function of groove-slit distance d, and is normalized with that without adjacent groove ( I 0 = E 0 2 ) such as performed in [3,16]. Figure 6 shows the strong correlation between the experimental data [3], the fully-vectorial numerical calculations [16]  Comparison between experimental (circles) [3], numerical (dashed line) [16], and the analytical (solid line) results with =852nm and 2 = -33.22+1.17i.

Double-slit near field
The field just above the surface between two subwavelength slits milled in a gold film, separated by 2d=10.44μm and illuminated with =974.3nm normal incident TM-polarized light, was recently measured through a scanning near-field optical microscope (SNOM) by Aigouy and co-workers [7]. Neglecting the tangential component of the E-field at the interface, the patterns recorded from fluorescence emission are expected to scale with E z 4 .
Modeling once more each nano-object (i.e. the slits) as a linear HED, the resulting near-field intensity pattern arises from the interference of two counter-propagating SWs, E z x ( ) = E sw where the origin of the x-axis is located at the half slit-to-slit separation distance. In Fig. 7 the near field between the slits, a E z x ( ) 4 + b , is plotted in the vertical axis, where b=0.33 is the background illumination offset taken from the original data and a=0.031 is the best-fit gain factor. The decay, phase and pattern periodicity predicted by the analytical model, with 1 =1 and 2 = 44.05+3.24i, closely match the experimental data.
Clearly, the dipolar model is again fully consistent with the near-field structure depicted in real-world experiments. It is worth to emphasize that the two studies [3,7] reported in this Section both assess the SW field in the near-zone propagation regime (0<x<2L bw ). The corresponding results yielded for either investigation by the analytical dipolar model indicate that the asymptotic solution is fairly accurate inside this transient near-zone as expected from the accuracy range k 0 x>2, or stated alternatively, x> / (see Subsection 3.1).

Conclusion
In summary, we provide a rigorous closed-form description of the surface wave generated via the diffraction of a horizontal Hertzian dipole over the air-metal interface by solving the corresponding Sommerfeld integral with a modified method of steepest descents. The asymptotic solution -accurate for distances x> / from the origin -demonstrates that the total surface near field is composed of two distinct components as previously evidenced in [6]: a long-range surface plasmon polariton (SPP) evanescent wave and a short-range boundary wave (BW) with free-space cylindrical properties. The dynamics of both constituent waves are fully revealed through simple analytic expressions, and appropriate parameters are defined to quantitatively and continuously describe their properties across the entire propagation range: near and asymptotic. Moreover, our calculations predict that the wavevector mismatch of the two co-propagating surface modes creates a weak periodic beating of the total field amplitude, which is noticeable at relatively large distances from the source and high optical frequencies. The closed-form solution is further validated via comparison between a dipolar model and recent experimental data, which demonstrates excellent quantitative agreement. In the process we show that the generation of surface waves by diffraction -and their ensuing interactions between one-dimensional subwavelength-sized nano-defects along the metallodielectric interface -are conveniently and accurately modeled with linear Hertzian dipoles. Hence, the theoretical formalism presented in this contribution can be used for the characterization and engineering of the near-field interactions in plasmonic and nano-optical devices while alleviating the reliance on time-consuming computer simulations.