Theory on the scattering of light and surface plasmon polaritons by arrays of holes and dimples in a metal film

The scattering of light and surface plasmon polaritons (SPPs) by finite arrays of either holes or dimples in a metal film is treated theoretically. A modal expansion formalism, capable of handling real metals with up to thousands of indentations, is presented. Computations based on this method demonstrate that a single hole scatters a significant fraction of incoming light into SPPs. It is also observed that holes and dimples scatter SPPs into light with similar efficiencies, provided the depth of the dimple is larger than its radius. Finally, it is shown that in arrays the normalized-to-area emittances in the out-of-plane and SPP channels present different dependences with the number of holes.


Introduction
Surface plasmon polaritons (SPPs) are electromagnetic modes bound to a metaldielectric interface [1]. Their unique optical properties have, among other applications, the potential for designing highly integrated photonic circuits with length scales much smaller than those currently achieved. Such circuits should convert light into SPPs, propagate them to logic elements where they are to be processed and ultimately converted back into light, as discussed in [2].
The theoretical analysis of these processes is hampered by the difficulty of computing the optical properties of nanostructured metals. Standard techniques capable of providing virtually exact results to Maxwell's equations can nowadays only treat either very small systems or systems with a high degree of symmetry. In this way, optical transmission through a single hole has been studied with the Multiple Multipole Method [3], Finite Difference Time Domain (FDTD) simulations [4], the Green Dyadic approach [5,6,7], and special methods devised for treating circular holes [8] or 1D systems (slits) [9]. Optical transmission through periodic infinite arrays of holes have been studied with FDTD [10,11]. Nevertheless, approximate methods are required when computing either the optical transmission through finite arrays of holes, or the scattering coefficients of SPPs impinging on finite structures. Regarding the scattering properties of SPPs, several approaches have been developed in order to study collections of 1D scatterers (i.e. with translational symmetry in one direction) [12,13,14,15]. In the 2D case (holes, protrusions, dimples...) the available theoretical formalisms are based on the coupling of dipoles, an approximation that is only valid when the dimensions of the scatterers are much smaller than the wavelength (specially when, as usual, the polarizability of the scatterer is represented by its quasi-static value) [16,17,18].
Our first main goal is to present a method for treating the electromagnetic properties of up to thousands of indentations in a real metal. We apply this formalism for computing (i) how much energy goes into the different SPP channels when an array of holes is back illuminated (see figure 1(a)), and (ii) the scattering coefficients of SPPs by finite arrays of 2D indentations (see figure 1(b)).
The indentations can be either holes or dimples (i.e. holes closed at one end), with arbitrary size and shape, placed at arbitrary positions. Our approach is based on a modal expansion of the fields in different spatial regions. It is a non-trivial extension from the formalism previously developed for dealing with either 1D indentations in a real metal [19] or 2D indentations in a perfect electrical conductor (PEC), where EM fields can not penetrate [20]. The model we shall present has already been applied to computing the plasmon coupling efficiency of a single defect (both with circular and rectangular shape) [21]. Furthermore, it was used for calculating the optical transmittance through different systems: a single rectangular hole [22], a periodic array of rectangular holes [23], and a finite array of circular holes [24]. Excellent agreement was obtained in all these cases with both experimental results [21,24] and FDTD simulations [22,23,24] (when available). However, no detailed account of its derivation has been presented before.
Our second main goal in this paper is to analyze the results on launching and decoupling of SPP by metal gratings reported in [25]. These experiments motivate our choice of geometrical parameters.
The paper is organized as follows. Section 2 is devoted to presenting the theoretical framework used in the manuscript. In order to facilitate the reading, in Section 2 we give an overview of the method and define the quantities studied throughout the paper. The derivation of the formalism as well as several useful expressions are given in the appendixes. In section 3 we present the optical response of a single circular hole (illuminated by either a plane wave or a SPP), while in section 4 we look at the scattering coefficients for arrays of circular holes.

Modal expansion formalism
We consider a set of indentations (either holes or dimples) of arbitrary shape, and arbitrarily placed in a planar metal film (infinite in the x-y plane and having finite thickness h). The system can be divided in the three regions shown in figure 1. Region I and region III are dielectric semi-spaces characterized by the real dielectric constant ǫ 1 and ǫ 3 , respectively. Region II represents the corrugated metal film with a wavelength-dependent dielectric function ǫ M . Holes or dimples could also be filled with a dielectric ǫ 2 . We assume that the system is illuminated by EM wavefields coming from region I.
We expand the EM fields on the eigenmodes of each region, and match them at the boundaries. The finite dielectric constant of the metal is taken into account by using Surface Impedance Boundary Conditions (SIBCs) [26]. Roughly speaking, the results obtained with the SIBCs can be understood as a second order Taylor expansion in z s = 1/ √ ǫ M . The first order term (z s = 0) is the result obtained within the PEC approximation, and was reported in [20]. The use of SIBCs does not simply represent a quantitative improvement over results obtained with the PEC approximation: it also allows for the computation of scattering of SPPs (which do not exist in flat PEC interfaces but already appear when SIBCs are employed). Notice however that SIBC can not describe tunneling of EM fields directly across the metal. Therefore our method is only applicable to metal thickness larger than 2-3 skin depths, when direct tunneling is negligible. A cautionary remark: when dimples are considered, this restriction in thickness should be applied from the bottom of the dimple rather than from the metal surface.
After some algebraic manipulations we end up with the following coupled system of equations for E α and E ′ α , which are essentially the modal amplitudes of the electric field at the input and output interfaces of the indentations, respectively (α is an index labeling each waveguide mode in each indentation considered in the calculation).
Details of the derivation, as well as the expressions for the different quantities involved can be found in Appendix A. Let us just mention that the ingredients required for the calculation are the spatial dependence of waveguide modes in each indentation, and their propagation constants. These modes are analytically known Optical properties of circular hole in a gold film illuminated by a normal-incident plane wave, for different hole radii r. The metal thickness is fixed to h=150. The system under study is schematically represented in the inset to panel (b), illustrating the two relevant scattering channels: radiation and launching of SPP. Panel (a) shows the spectral dependence of the normalized- for some hole geometries; otherwise they can be numerically computed solving a 2D problem. In this section we just give the physical meaning of the different quantities in (1). I α represents the direct illumination over waveguide mode α. It is proportional to the overlapping of the incident electric field with the mode in the indentation. The rest of the terms take into account the (self-consistent) wandering of fields between the indentations. Σ α is related to the bouncing back and forth of a given waveguide mode inside an indentation, due to the discontinuity that the waveguide mode faces at the end of the indentation. The main difference between holes and dimples is the presence of G V α , which only appears for holes, and is related to the EM field on one side of the hole due to the presence of an EM field on the other side. The coupling between different waveguide modes is given by the "propagator" G αβ . This takes into account that the EM field emitted by each point within object β can be "collected" by the object α. The propagator G ′ γν differs from G αβ in the constituent parameters only, i.e. G ′ γν is a function of ǫ 3 while G αβ depends on ǫ 1 . Notice that this formalism is reminiscent of the Green Dyadic approach [27]. In that method the EM fields must be computed in the volume inside the indentations. In contrast, our method provides the EM fields everywhere in terms of the fields at the openings of the indentations. This very compact representation allows for the treatment of the optical properties of systems involving a large number of indentations, at a relatively low computational cost [24].
The price we pay is that the modal expansion formalism is approximate and relies on the SIBCs, while the Green Dyadic method is virtually exact. The implementation of modal expansion calculations usually faces another source of inaccuracy. The required overlaps between plane waves and waveguide modes are only known for some waveguide geometries defined in a PEC. Nevertheless with some minor corrections, waveguide modes in a PEC can still yield good agreement with experimental results (and exact calculations, when available). For instance, using the propagation constant of the waveguide in a real metal greatly improves the prediction of the spectral position of resonances. Another improvement is to enlarge the hole side, so as to simulate the real penetration of the field at the lateral wall defining the waveguide. Doing so, our method has provided good agreement with the FDTD simulation for the optical transmission through a single rectangular hole [22]. In that case, the hole was enlarged as proposed in [28]. In this paper we will only consider circular indentations, for which the propagation constant of waveguide modes in a real metal can be easily computed [29]. The expressions for the waveguide modes are borrowed from those of a circular waveguide in the PEC [29,30]. The radius of the hole is phenomenologically enlarged by one skin depth. Such enlargement provides the best agreement with FDTD simulations for an infinite periodic array of holes [24].
As said before, once the coefficients E α and E ′ α are known, it is possible to obtain the fields everywhere. In particular, the energy power scattered into different channels can be computed by integration of the relevant time-averaged components of the Poynting vector. We define W z as the total radiated power crossing a fictitious plane placed at constant z in region I. The definition of W II z and W ′ z is the same as W z , except the appropriate fictitious plane is placed in region II and III, respectively. Notice that we are maintaining the notation used in the system of equations, where primed quantities refer to region III. The values of these quantities in terms of the set E α and E ′ α can be found in Appendix C. In a real metal, EM energy power can also leave the system as a SPP wave. We define W spp and W ′ spp as the energy power in the SPP channel that crosses an ideal cylinder with axis parallel to the z direction, placed in either region I or III, respectively. Details and expressions for these quantities can be found in Appendix D.
We shall present results in terms of normalized emittances J, defined as the total power W divided by both the incident power impinging onto the area covered by indentations in region I, W inc , and the total number of holes in the system N . In other words, for each scattering channel we define a corresponding emittance as J = W/(W inc N ), independently of whether the incident field is a SPP or a radiation field. However, when illumination is via a SPP, W inc is defined as the total power crossing the infinite imaginary strip perpendicular to the metal surface, whose base coincides with the maximum geometrical cross section of the collection of indentations. For a lossy system the total scattered power in the SPP channel W spp (D.6) is a function of the observation point, R, at the metal surface. As plasmon's decay with in-plane distance is well known, it is possible to define the total scattered power in the plasmon channel as if the plasmon were excited at R = 0. We denote this quantity by Ω spp (D.7), which satisfies W spp = Ω spp exp(−2|Im[k spp ]|R) (D.6). Nevertheless, the comparison of efficiencies into the SPP channel between systems of different sizes is not trivial in lossy systems. We need to treat separately the case of back-illuminated and SPP illuminated arrays.
When comparing efficiencies into the SPP channel for back-illuminated arrays of different sizes at fixed R, larger systems may have emitters of SPPs closer to that reference distance R, leading to exponentially larger values of Ω spp than those of smaller arrays. In order to avoid this geometric artifact in the definition of Ω spp , we analyze the SPP field at a fixed distance R e from the edge of the array. If we measure distances from the center of the array: R = R e + L sys , where we define L sys as the distance from the center of the array to its edge. Then the normalized SPP emittance The "cross-section" σ spp is the quantity we focus on throughout the paper, as it provides the normalized intensity into the SPP channel on a circle passing through the edge of the array.
When the array is illuminated by a SPP, however, using the same convention may . Spectral dependence of both σ ′ spp (a) and J ′ z (b) for a "horizontal" linear array of circular holes with radius r=125 nm, h=150 nm, and period P=760 nm, for different array sizes. The system under study, schematically represented in the inset, is back-illuminated by a plane wave.
lead to spurious size dependences in the "cross-section" σ spp . It is not convenient to place the origin on the center of the array since the SPP wave may not even reach the center for large arrays and strong SPP scattering. Hence, we place the origin at the center of the first column of indentations encountered by the incident SPP and analyze the response of different arrays at the same distance R from this origin, irrespectively of the array size. Therefore, in this case, we have J spp = σ spp exp(−2|Im[k spp ]|R), with Let us stress that the previous conventions are chosen in order to obtain a crosssection into the SPP channel which allows for the comparison between arrays of different sizes; however, these conventions do not affect neither the total flux W spp nor the flux normalized to both area and number of holes J spp .
Let us close this section by reminding that current conservation imposes relations that are useful for checking the correctness of both the derived expressions and the computer codes. For instance, if absorption is neglected, the identity J II z = J ′ z + J ′ spp should be fulfilled, i.e the fraction of the incident energy traversing the hole, J II z , must be either radiated into freely propagating waves, J ′ z , or scattered into SPP at the output side of the metal film, J ′ spp . On the input side we have a similar relation, which now depends on whether the incident field is a plane wave or a SPP. Although all results presented in this paper include absorption, we have checked that for a lossless metal our code provides current conservation up to a relative error of 0.1% of the incident flux.

Single defect
In this section we shall study the optical properties of a single circular hole in a real metal, when illuminated by either a plane wave or a SPP. We shall analyze both the out-of-plane radiation and how much energy goes into the plasmon channel.
Let us first consider back-illumination of the hole by a plane wave, impinging at normal incidence. The optical transmittance through a single hole perforated in a real metal has been largely studied both experimentally [21,24,31,32,33] and theoretically [3,4,8,22,34]. In general, the transmittance is characterized by Fabry-Pérot peaks, one of them being very broad and appearing close to the cutoff wavelength of the fundamental waveguide mode [22,34]. The previously cited theoretical works have looked at the spatial dependence of scattered fields, finding that a hole launches SPPs. However, to the best of our knowledge, no computation has yet been performed on the efficiency of the SPP launching by a back-illuminated hole (Ref. [21] considered the coupling into the SPP channel at the interface of incidence). We have performed such calculation, analysing the scattered EM arising from the SPP pole in the Green dyadic (see Appendix D). This analysis is possible even in the presence of absorption, although in lossy metals the actual EM fields at the metal surface decay with the distance from the hole due to the dissipation of SPP energy into heat. Figure 2 renders the (normalized-to-area) fraction of incident energy that is scattered by a circular hole into either the SPP (σ ′ spp , panel (a)) or radiation (J ′ z , panel (b)) channels. Both quantities are computed in the region of transmission. The metal considered is gold, whose dielectric constant is fitted to Palik's data [35]. The film thickness is h=150 nm, and different hole radii are studied. For the sake of simplicity we assume a freestanding metal film since the peak investigated in the experiment concerns the metal-air interface. As figure 2 shows, both σ ′ spp and J ′ z present non-monotonous spectral dependencies, due to a broad resonance appearing close to the cutoff wavelength of the fundamental waveguide mode, λ c . For Gold, we obtain λ c (r=125 nm) = 589 nm, λ c (r=175 nm) = 732 nm, and λ c (r=200 nm) = 811 nm. Notice that spectra are plotted from λ = 600 nm, because far from λ < 600 nm results are not reliable because the SIBC approximation breaks down(i.e ǫ m (λ) ≪ 1 is not longer fulfilled in that regime). However, in the considered spectral window, as much as 30% of the energy impinging into the holes can be converted into SPPs. Figure 2 (c) shows that the ratio σ ′ spp /J ′ z is a smooth radius-dependent function of wavelength. For the parameters considered, this ratio is of the order of 0.3-0.8, being larger for smaller holes. In the long-wavelength limit (λ ≫ r) this ratio can be worked out analytically; we find that σ ′ spp /J ′ z ∼ |z s | ∼ λ −1 for a metal represented by a Drude dielectric constant. This decrease of the coupling to SPPs with wavelength originates from the weakening of the SPP confinement to the surface, which translates into a weaker coupling to indentations.
In our formalism the ratio σ ′ spp /J ′ z is independent of metal thickness. However, both σ ′ spp and J ′ z depend on h, as shown in figure 3 for a single hole with r = 125. σ ′ spp and J ′ z decay with both h and λ in the spectral window considered, where all fields inside the hole are evanescent.
Let us consider now the decoupling of SPPs into radiation, after they have been scattered by either a hole or a dimple (see scheme in figure 1(b)). To the best of our knowledge, this problem is virtually unexplored. We are only aware of a study on SPP scattering by a shallow dimple with a Gaussian profile [36]. Figure 4 shows the spectral dependence of the energy power radiated into the region of incidence, J z , when a SPP impinges into either a circular hole (a) or a circular dimple (b). In both cases the radius is r =125 nm. Several values of h are considered (h is the metal thicknesses in the case of a hole, and the depth of the indentation in the case of a dimple). The inset of figure 4(b) renders the ratio J z /σ spp (notice that this is the inverse of the quantity shown in figure 2(c)). This ratio depends on the radius of the indentation, but does not depend neither on h nor  Figure 6. (Color). Spectral dependence of both σ ′ spp (a) and J ′ z (b) for a "vertical" linear array of circular holes with radius r=125 nm, h=150 nm, and period P=760 nm, for different array sizes. The system under study, schematically represented in the inset, is back-illuminated by a plane wave.
on whether the indentation is a hole or a dimple. For the considered case this ratio is always larger than unity, implying that a single hole scatters better into radiation than back into SPPs. Figure 4 also shows that a hole and a dimple couples SPPs into radiation with similar intensity, provided the indentation is deep enough. Dimples may provide higher values of J z than holes, due to the lack of radiation into region III. However, shallow holes radiate more than dimples. This is trivial when h → 0 since the dimple then disappears. We find that, in the subwavelength regime, h ≈ r is a rough rule of thumb for the depth at which the scattering of SPPs into radiation is the same for dimples and holes of the same size. For smaller h, conversion of SPPs into radiation is substantially smaller for a dimple than for a hole. Notice though that this "transition" depth depends on the evanescent decay of the fields inside the indentation, and therefore on wavelength. For instance, a hole illuminated by a SPP at λ =800 nm radiates 35% more than a dimple, when both have h =150 nm.

Array of circular holes
We shall now investigate scattering by arrays of circular holes of both a normal incident plane wave and a SPP.
Let us first focus on the transmittance and the launching of SPP via backillumination of the array. In the perfect conductor approximation, the normalizedto-area optical transmission through linear arrays of holes shows resonances of the same order as those in square arrays [20]. This occurs if the array is lined up parallel to the incident electric field. In contrast, the normalized-to-area optical transmission through an array of holes perpendicularly oriented to the incident field is very similar to that of a single hole. Here we consider these configurations and analyze whether the same phenomenology holds for the launching of SPPs.
Take a set of holes with radius r=125 nm, separated by a distance P = 760 nm in a metal layer of thickness h=150 nm (again, the geometrical parameters of [25]). parallel to the incident electric field ("horizontal" array), a perpendicularly aligned one ("vertical" array), and a square array, respectively. In all three figures, panel (a) shows the spectral dependence of the normalized-to-area SPP emittance, σ ′ spp , while panel (b) renders the normalized-to-area transmittance J ′ z , for a different number of holes.
Clearly, holes in the "vertical" configuration ( figure 6) are virtually independent. Conversely scattering of light by "horizontal" arrays (figure 5) shows resonances due to interaction between holes. These resonances also appear in the square arrays (figure 7) ‡. Notably, the wavelengths of SPP emittance peaks appear slightly redshifted with respect to the periodicity. Actually, the position of the SPP emittance peak is also slightly redshifted with respect to the wavelength of SPPs of the flat metal interface with in-plane wavevector k spp = 2π/P (this wavelength is λ = 777 nm for gold and period P = 760 nm). This redshift reflects the fact that the EM field spends some time inside the indentations, slightly modifying the resonant conditions obtained by considering geometrical arguments only.
At resonance, the normalized-to-area transmittance J ′ z is larger in the square array than in the "horizontal" one. Additionally, the normalized-to-area emittance σ ′ spp has already reached saturation for arrays with 20x20 holes, while it has not saturated yet for "horizontal" arrays with N = 30 holes. This different behavior with number of holes can be understood by analyzing the angular dependence of the power scattered into SPP, σ ′ spp (λ, θ) (θ = 0 coincides with the direction of the incident electric field). Panel (a) of figure 8 renders a contour plot σ ′ spp (λ, θ) for an isolated circular hole in Au (r=125 nm, h=150 nm). As the figure illustrates, a hole launches SPPs mainly in the directions parallel and antiparallel to the incident electric field. This is because a SPP is partly a longitudinal wave, with its in-plane wavevector parallel to the in-plane electric field. The SPP launched by a "horizontal" array of holes at angles that are not close to either θ = 0 or θ = 180 degree are not further scattered ‡ Notice that, in order to better resolve the details close to the peaks, the spectral window shown here is smaller than the one used when presenting results for the single hole . Angular and wavelength dependence of the emittance of SPPs, σ ′ spp (λ, θ), by a back-illuminated set of holes. The incident field is a plane wave with electric field pointing along one of the cartesian axis of the array, direction from which θ is measured (in degree). Panel (a): single circular hole, panel (b): 10 × 10 square array of circular holes with P =760 nm. In both cases r =150 nm and h =150 nm. Panel (c) renders the curves of maxima (heavy line) and minima (thin lines) intensity for a diffraction grating of point emitters [37]. and thus contribute to σ ′ spp . In a square array, however, those "off-axis" SPPs which are launched by holes close to the center of the array can be scattered by other holes. This scattering can be either into SPPs or into radiation. Such process increases J ′ z and decreases σ ′ spp , and explains the different size dependencies of "horizontal" and square arrays.
The angular dependence of the SPP emittance for a square array of 10x10 holes, is shown in figure 8(b). The hole radius and metal thickness are those of the single hole case of figure 8(a), while the period is P = 760 nm. The angular dependence of the SPP emittance by the array can be understood as the one due to the single hole, modulated by the interference pattern characteristic of a diffraction grating [37]. We recall that, for light impinging normally on a finite "vertical" diffraction grating composed of N point emitters, intensity maxima appear at angles θ given by the relation λ/d = 2 sin 2 (θ/2)n/N where both n and n/N are integers. Angles of intensity minima satisfy the same equation with n integer and n/N non-integer. Similarly, finite "horizontal" gratings satisfy λ/d = 2| sin θ|n/2N , with the same conditions as the previous case governing the presence of maxima and minima in terms of n/N . The maxima of both horizontal and vertical gratings are represented by heavy lines in figure 8(c), while the minima are represented by thin lines. The emittance of SPP by a hole array ( figure 8(b)) combines both the pattern of maxima and minima observed in figure 8(c) and the single hole emittance ( figure 8(a)).
The influence of hole geometry on the SPP and radiation efficiencies in arrays is illustrated in figure 9, which renders the emittance spectra for a "horizontal" linear array of 30 circular holes with period P=760 nm as function of both its radius r and metal thickness h. The hole array with r=125 nm and h=150 nm, already studied in figure 5, is compared with both an array of smaller h and the same r and an array of larger r and the same h. Reducing the film thickness to h=100 we increase the total , and the ratio σ ′ spp /J ′ z for a "horizontal" linear array of 30 circular holes with period P=760 nm and varying radius r and metal thickness h. The system under study is backilluminated by a plane wave. transmission through the system increasing the amount of energy scattered in both channels, σ ′ spp in figure 9(a) and J ′ z in figure 9(b), keeping the the same radio as for h=150, see figure 9(c). On the other hand, when the circle radius is increased to r=150 nm the transmittance of the system is favored against σ ′ spp , leading to a lower relative ratio. Notice, however, that these properties may present different dependences with number of holes for different hole sizes. This point deserve further investigation which exceeds the scope of the present paper.
Let us now focus on the scattering of an incident SPP by a "horizontal" linear array of indentations. As seen, the properties of a square array of indentations are very similar to those of "horizontal" arrays. It has been experimentally found [25] that square arrays of holes (with radius r =125 nm drilled in a gold film 150nm thick) decouple SPPs into radiation, while dimples with the same radius and 50 nm depth do not radiate any measurable EM signal. In both cases the period of the array was P =760 nm. Figure 10 presents the spectral dependence of the normalized out-of-plane radiation, J z , for "horizontal" arrays (holes: panel (a), dimples: panel (b)). The same geometrical parameters as in the experiments were used. Additionally, in figure 10(b) we present results for dimples with the same depth as that of the holes considered in the experiment. For the geometrical parameters and spectral window considered, this figure shows that, J z develops a resonance, which appears slightly redshifted from the period. As figure 10 shows, arrays formed with either holes or dimples with h=150 nm scatter SPP into radiation with similar efficiencies. However, in line with the experimental results, the computed intensity of this process for arrays of dimples with h=50 nm is around five times smaller, at the wavelength considered in the experiment (λ = 800 nm). In any case, notice that for these geometrical parameters, for which scattering by a single indentation is weak, the normalized out-of-plane radiation is of the same order for arrays and for single indentations (results for the single hole of figure 4 are included in figure 10 in order to facilitate the comparison).
As we have seen, "horizontal" arrays of holes produce resonant behavior in both the launching of SPPs and the radiated energy. This occurs independently of whether the array is illuminated with a plane wave or by a SPP. However, all these efficiencies present very different dependencies with the size of the array. In order to illustrate this point we consider "horizontal" arrays of holes, with r=125 nm and P =760 nm, drilled in a gold film of thickness h =150 nm. Figure 11(a) renders the evolution of the maximum value of both σ ′ spp and J ′ z with the number of holes (at the resonance around λ ≈ 790nm), when the array is back-illuminated by a plane wave. This figure completes the analysis of the peaks in figure 5 considering arrays of larger sizes. For the geometrical parameters analyzed, σ ′ spp reaches a maximum for arrays of about 35 holes. In contrast, the saturation of J ′ z has not yet been reached even for arrays of 150 holes. In Figure 11(b) we plot the size dependence of both the maximum of J z and the maximum of σ spp when the "horizontal" array is illuminated by a SPP. We recall that these quantities are normalized to the total number of holes. In this case, for large arrays, both J z and σ spp decrease with the size of the array. This is because part of the incident SPP current is radiated as it progresses along the array, therefore weakening the illumination of the holes at the further end of the array. The values of J z are even smaller in "horizontal" arrays than in single holes. Consequently, "vertical" arrays (which, as we have seen behave as a collection of independent holes) are a better choice for out-coupling of SPPs than "horizontal" arrays.
As result of the hole-hole interaction, we find that in a hole array at resonance Jz > σ spp , in contrast with the result for the single hole. This behavior can be understood in terms of simple physical arguments. In-plane emittance for a single hole occurs essentially along the forward-backward direction, see figure 8(a). Along this direction, at resonance, the EM fields of the individual holes are in phase, interfering constructively. Therefore, the total EM field along these directions of favored SPP emission increase with number of holes. On the other hand, the out-of-plane single hole emission (not shown here) is practically isotropic. Therefore, the interference of the individual holes is either constructive or destructive, depending on the optical path to a given observation point. As a result of having both constructive and destructive interference, the total contribution to this scattering channel is smaller than for the in-plane radiation. It is important to emphasize that the emittances reported in figure 11 are measured at resonance. This explains why the figure considers N 5: for smaller N the peak is practically not yet developed, see figure 5 and figure 10.

Summary
We have presented a theoretical formalism for dealing with optical properties of sets of indentations in an optically thick metal film. The indentations may have arbitrary shape and can be arbitrarily placed at any of the two metal interfaces. The formalism is based on a modal expansion of the fields in the different regions of space. Fields are matched by Surface Impedance Boundary Conditions. These boundary conditions take into account the dielectric constant of a real metal and can therefore be used to approximately describe surface plasmon polaritons. Hence, our formalism can be used for analyzing optical problems like transmission through hole arrays, scattering of SPP by sets of indentations and launching of SPP by back-illuminated arrays, among others. Our systematic study of the optical properties of both circular isolated indentations and arrays of indentantions in gold has found that: • Close to the cutoff condition, a back-illuminated hole scatters a significant fraction of light into SPPs in the region of transmission (≈ 30%). • A hole and a dimple of the same size scatter SPPs into radiation with comparable efficency if, roughly speaking, the depth of the indentation is not smaller that its radius. Otherwise, holes couple SPP into radiation more efficiently. • Linear arrays oriented along the direction of the impinging electric field exhibit resonances similar to those appearing in square arrays of holes. This holds in the coupling into both radiation and SPPs. Furthermore, the appearance of these resonances is independent of whether the illumination is via a plane wave or a SPP.
• Linear arrays oriented perpendicularly to the direction of the impinging electric field essentially behave as collection of independent holes. • Illuminating by a SPP or a plane wave results in a substantially different behavior of the emittance as function of the array size. • The SPP channel saturates for smaller number of holes than the radiation channel, when a hole array is illuminated by a plane wave.
We solve Maxwell's equations (written in CGS units) for the electric, E, and magnetic, H, fields in a system comprising a planar metal film corrugated with a set of indentations. The film is in the x-y plane and has a finite thickness h. The indentations can be either holes or dimples (i.e. holes that do not pierce through the film), with arbitrary shape and placed at arbitrary positions at both input and output interfaces. We assume the system lies in a rectangular supercell, with lattice parameters L x and L y along the x and y axes, respectively. This supercell may be real (if we are considering a bona fide periodic system) or artificial, if the number of defects is finite. In the latter case, the limit L x , L y → ∞ must be taken. The space within the supercell is divided into the three regions shown in figure 1. Region I and region III are dielectric semi-spaces characterized by the real dielectric constant ǫ 1 and ǫ 3 , respectively. Region II (which extends from z = 0 to z = h) represents the corrugated metal film with a wavelength-dependent dielectric function ǫ M . Holes or dimples could be also filled with a dielectric ǫ 2 . We assume the EM energy incident on the metal layer is coming from region I (see figure 1).
In order to take advantage of known solutions of Maxwell's equations the EM fields of the constituent media are represented in a convenient basis. We use Dirac's notation, such that r |E = E(r ) = (E x (r ), E y (r )) t , where t stands for transposition and r = (x, y).
In region I, see figure 1, the fields are expanded into an infinite set of plane waves with parallel wavevector k = (k x , k y ) and polarization σ. We consider that the incident radiation (labeled by the superscript 0) is either a plane wave or a SPP. Arbitrary illumination can be easily considered by decomposing first the incident wavepacket into plane waves, and subsequently integrating the optical response of the array to each plane wave. The fields can be written in terms of the reflection amplitudes ρ k σ as where u z is unitary vector directed along z-direction and the summation runs over wavevectors of the form k = k 0 + K R , K R being a vector of the supercell reciprocal lattice. The real space representation of the plane waves is given by where we assume a k σ-independent normalization k σ|k σ = 1. The electric and magnetic fields in (A.1) are related through the admittance Y k p = k ω ǫ 1 /k z and Y k s = k z /k ω where k ω = 2π/λ (λ is the wavelength of the incident radiation), and k 2 + k 2 z = ǫ 1 k 2 ω . In the transmission region (region III in figure 1, where z ≥ h) the EM fields are also expanded in plane waves, and expressed in terms of the transmission amplitudes t k σ Note that all quantities in region III are primed in order to distinguish them from those in region I. Waveguide modes are the most natural choice for expanding the EM field inside the indentations. The main ingredients needed in the present formalism are both the propagation constants and the EM fields of waveguide modes in all indentations. In general, these quantities are attained by solving an EM problem with reduced (cylindrical) symmetry. For holes with a given symmetry (circular, rectangular...) these modes are known analytically, if the surrounding metal is considered a perfect electrical conductor. In such favorable situation, the penetration of the EM field into the real metal can be taken into account by phenomenologically enlarging the hole.
Writing the fields as function of the expansion coefficients C α and D α , we have where α labels a certain waveguide mode inside a certain aperture. In what follows we shall refer to a given element of the set {α} as an object. The propagation constant along the z direction of eigenmode α is given by k αz , while Y α = k αz /k ω and Y α = ǫ 2 k ω /k αz are the admittances for TE and TM modes, respectively. In order to take into account the dielectric properties of a real metal, we impose surface impedance boundary conditions (SIBCs) [26]. SIBCs establish the following relation between the tangential components of the EM fields at the horizontal interfaces: E = z s H × u n , where u n is the vector normal to the surface and pointing to the interior of the metal and z s = 1/ √ ǫ M . The previous continuity equation must be satisfied at each interface. Projecting over the different plane waves, at both z = 0 and z = h, we get where f ± α = 1 ± z s Y α and δ ij is the Kronecker's delta. We have also defined the quantities which (approximately) represent the modal amplitudes of the electric field at the input and output interfaces of the indentations, respectively. The overlapping integral, k σ|α , of the plane wave and the waveguide mode are analytically known for simple geometries of the indentations. For circular holes these expressions can be found in [30,38].
Let us now impose the continuity of the tangential component of the magnetic field. As this boundary condition only holds at the openings, we must project over waveguide modes, which form a complete set in that region. Plugging ρ k σ and t k σ (A.5) into the projection over waveguide modes, we end up with the following set of for a dimple. (A.8) and G ′ γν differs from G αβ only in the constituent parameters, i.e., G ′ γν is function of ǫ 3 while G αβ depends on ǫ 1 .
The interpretation of these quantities is as follows: I α takes into account the direct illumination over object α, Σ α is related to the bouncing back and forth of the waveguide fields inside the indentations, G V α reflects the coupling between EM fields at the two sides of a given indentation, while G αβ and G ′ αβ are propagators that couple "objects" at the same side of the metal film.
If the system is periodic, with periodicity defined by L x and L y , G αβ is calculated through the discrete sum over Bragg diffraction modes defined above. For non-periodic systems, L x and L y define a fictitious super-cell and the limit L x , L y → ∞ must be taken. This can be done analytically, simply replacing all sums in diffraction modes by integrals in k , i.e. k σ → L x L y (2π) −2 σ d 2 k . Notice that the diverging factor L x L y cancels out with the normalization factor of the plane waves.
Once the self-consistent E α , E ′ α are found, (A.5) gives the expansion coefficients ρ k σ and t k σ , which in turn can be used to obtain the EM fields everywhere (Appendix B is devoted to derive the expressions for the EM fields in region III).
So far we have focused on in-plane components of the EM-fields. The z-components of the field can be readily computed using Maxwell equations.
For z s = 0 (perfect conductor case) the expressions presented in this Appendix simplify to those derived in [20]. Notice however, that in [20] the definitions of I α , Σ α , G v , G αβ , and G ′ αβ contain an additional multiplicative factor ı. As explained in [39], this extra factor is useful for drawing a parallelism with the "tight-binding" method employed in solid state physics, but it will be omitted here.
We also assume that the rest of the integrand, F (k ), is practically constant in the region of the upper complex semi-plane close to this SPP pole, and that we can replace the exact result of the integral by its residue, i.e.
Special care must be taken when integrating the Bessel function, as this function diverges for large arguments both in the lower and upper complex semi-planes. This problem can be solved by using the identity 2J n (x) = H (1) n + H (2) n , where H (1) n and H (2) n are the n-th order Hankel functions of first and second kind, respectively. The part of the argument containing H (1) n must be integrated in the upper complex semiplane, while H (2) n has to be integrated in the lower semi-plane. Only the integration over the upper complex semi-plane encloses the SPP pole. The integral in the lower semi-plane vanishes identically. As a result, and using the asymptotic expression H where γ = iǫ 2 k 2 ω z s e −iπ/4 [2π] −1/2 and Y spp = −z −1 s is the admittance evaluated at k spp . Comparison with numerical results shows that, the previous expression is an excellent approximation in the optical regime for good metals (like Au and Ag) already for R 2λ, despite being a long-distance asymptotic result.
Up to this point we have assumed a single indentation is placed at the origin of the system of reference . We now consider the case of an arbitrary number of indentations, located at arbitrary positions R α = |R α | u α . In order to compute the long-distance asymptotic expression, at a given observation point R, of the SPP field launched by the indentations, we only take into account terms of order O(R α /R). Then |R − R α − r ′ | ≈ R − R α · u R − r ′ · u R , where u R is a unitary vector directed from the indentation to the point R and r ′ is the local coordinate inside the hole. The SPP EM fields can be computed after integrating the propagator over the different areas occupied by the indentations, as indicated in (B.1).
After performing the integration, we obtain for the non-vanishing components of the SPP field: where V = −k spp /k spp z . The function f (u R ) gives the angular dependence of the total launched SPP field, and accounts for the interference of the SPPs launched by all the different waveguide modes in all the different indentations where the angular dependence of the SPP field scattered by a given normalized waveguide mode (with in-plane electric field given by E α (r ′ )) is provided by where Ω spp = ǫ 2 |z s |k 2 andf 2 = (2π) −1 2π 0 |f (u R ) | 2 dθ. This is the result integrated over all angles. Clearly, the quantity |f (u R )| 2 determines the angular distribution of the SPP scattered power.
For a lossless material z s = iIm[z s ], Im[k spp ] = 0, and Re[k spp z ] = 0. Therefore W spp is independent of R; otherwise W spp decays as function of R. For lossy systems Ω spp still provides the total energy scattered into the plasmon channel; the dependence with distance to the scatterers of W spp follows the decay of the energy stored in the SPP, due to absorption by the metal.
Finally, we adopt the following normalization for the scattering problem related to the out-coupling of a SPP into propagating radiation. We assume a SPP propagating at the metal-dielectric interface is incident on the system composed of either hole or dimples, at an angle θ 0 taken from the x-axis, see figure 1(b). The incident EM fields are given by E inc = e ik spp ·r (cos θ 0 , sin θ 0 , V ) t /Y spp , H inc = e ik spp ·r (− sin θ 0 , cos θ 0 , 0) t , (D.8) from which we derive the time-averaged Poynting vector along an arbitrary radial direction r, i.e. S inc r = e −2|Im[k spp ·r]| Re[k spp · u r ]/(ǫk ω ). We define the incident power, W spp inc , integrating S inc r in the rectangular area at x = 0, which extends along the side of the defect in the y-direction, and along z > 0 in the z-direction. For a SPP impinging onto a single circular indentation (with radius r) we obtain