OSA's Digital Library

Optics Express

Optics Express

  • Editor: Michael Duncan
  • Vol. 14, Iss. 7 — Apr. 3, 2006
  • pp: 2552–2572
« Show journal navigation

Extraordinary transmission through 1, 2 and 3 holes in a perfect conductor, modelled by a mode expansion technique

J.M. Brok and H.P. Urbach  »View Author Affiliations


Optics Express, Vol. 14, Issue 7, pp. 2552-2572 (2006)
http://dx.doi.org/10.1364/OE.14.002552


View Full Text Article

Acrobat PDF (360 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

We discuss a mode expansion technique to rigorously model the diffraction from three-dimensional pits and holes in a perfectly conducting layer with finite thickness. On the basis of our simulations we predict extraordinary transmission through a single hole, caused by the Fabry-Perot effect inside the hole. Furthermore, we study the fundamental building block for extraordinary transmission through hole arrays: two and three holes. Coupled electromagnetic surface waves, the perfect conductor equivalent of a surface plasmon, are found to play a key role in the mutual interaction between two or three holes.

© 2006 Optical Society of America

1. Introduction

The scattering from a single sub-wavelength hole in a conducting metal layer was already modelled in the 1940s and 1950s [1

01. H.A. Bethe, “Theory of diffraction by small holes,” Phys. Rev. 66, 163 (1944). [CrossRef]

, 2

02. J. Meixner and W. Andrejewski, “Strenge Theorie der Beugung ebener elektromagnetischer Wellen and der vol-lkommen leitenden Kreisscheibe und an der kreisformigen Offnung im vollkommen leitenden ebenen Schirm,” Annalen der Physik 7, 157–168 (1950). [CrossRef]

, 3

03. C. Flammer, “The vector wave function solution of the diffraction of electromagnetic waves by circular disks and apertures. I. Oblate spheroidal vector wave functions,” J. Appl. Phys. 24, 1218–1223 (1953). [CrossRef]

, 4

04. C. Flammer, “The vector wave function solution of the diffraction of electromagnetic waves by circular disks and apertures. II. The diffraction problems,” J. Appl. Phys. 24, 1224–1231 (1953). [CrossRef]

, 5

05. C.J. Bouwkamp, “Diffraction theory,” Reports on progress in physics 17, 35–100 (1954). [CrossRef]

]. Apart from the assumption that the metal is perfectly conducting, the metal layer was also assumed to be infinitely thin; these assumptions made it possible to calculate the diffraction analytically. However, with respect to enhanced transmission through sub-wavelength hole arrays, as reported in Ref. [6

06. T.W. Ebbesen, H.J. Lezec, H.F. Ghaemi, T. Thio, and P.A. Wolff, “Extraordinary optical transmission through sub-wavelength hole arrays,” Nature 391, 667–669 (1998). [CrossRef]

], these assumptions are not valid. Furthermore, it is believed that the excitation of surface plasmons by neighboring holes plays an important role in the enhanced transmission.

However, the specific case of two holes and their mutual interaction has not received much attention. This is remarkable, because it is a fundamental system and building block of a hole array. Yet, especially when the two holes are a wavelength or more apart, as the corresponding computational domain is then also large, this system is not easy to calculate.

For two-dimensional structures, such as infinite slits or circular symmetric setups, or for periodic systems in which the computational domain is only as large as one unit cell, current computer power suffices to obtain acceptable accuracy. For a finite number of three-dimensional holes that are far apart, however, a clever computational scheme is needed to prevent computation times of days or even weeks.

Fig. 1. Problem under consideration. Multiple rectangular holes in a perfectly conducting layer with finite thickness.

We find that extraordinary transmission occurs for a single hole when the lowest order waveguide mode is just above cut-off. This is due to the Fabry-Perot resonance of this lowest order mode inside the waveguide. Furthermore, we calculate the transmission through two and three holes. In order to understand the mutual interaction between two and three holes, we normalize by the transmisson through an identical but single hole. In this way we are able to isolate the effect that the presence of the second (and third) hole has on the transmission through the first. By changing the polarization of the incident light and by replacing one of the holes by a pit, we are able to show that coupled electromagnetic surface waves cause the enhanced and decreased transmission.

2. Problem definition and system parameters

Let (x,y,z) be a rectangular Cartesian coordinate system. Perpendicular to the z-axis we have a perfectly conducting layer with finite thickness D. In this layer a finite number of rectangular holes and pits are present. See Fig. 1. A hole is a rectangular cylinder that is open on both sides and is as long as the thickness of the layer; a pit has an open end at one side, either at z = D/2 or at z = -D/2, and a depth dp < D. The sub- or superscript p denotes the number of the pit or hole. The lengths in the x- and y-direction are Lxp and Lyp , respectively. The cross-section of a hole or a pit p is given by Ω p = {(x,y)|x0p<x<x0p+Lxp ,y0p<y<y0p+Lyp }. The halfspaces z > D/2 and z < -D/2 are filled with homogeneous dielectrics with index of refraction nu and nl , respectively. Every hole and pit is filled with a homegeneous dielectric with index of refraction np . The corresponding relative permeabilities are ϵu = nu2, ϵl = nl2 and ϵp =np2. The magnetic permeability is μ 0 everywhere.

A monochromatic incident field can originate from above and/or from below. The wavelength of the field in free space is given by λ. The local wavelengths are λu = λ/nu , λl = λ/nl and λp = λ/np . The corresponding wave vectors are ku = 2π/λu , kl = 2π/λl and kp = 2π/λp . The harmonic time dependence of the electromagnetic field is given by the factor exp(-iωt), with ω < 0, which will be omitted throughout.

3. Mode expansions

In section 3.2, we describe the electromagnetic field above and below the conducting layer. This field is expanded in propagating and evanescent plane waves, that are characterized by their polarization and their direction of propagation.

3.1. Inside the holes and pits

Solving Maxwell’s equations inside the pits and holes means finding solutions of the scalar Helmholtz equation for every Cartesian component of the electromagnetic field. The boundary conditions imply that, at a perfect conductor, the tangential electric field as well as the normal magnetic field vanish. We then find solutions that are called waveguide modes. These are propagating or evanescent in the z-direction:

[Eα(r)Hα(r)]=[EαxyHαxy]e±zz,
(1)

with propagation constant γz given by [21

21. Here and henceforth the square root of a complex number z is defined such that for real z > 0 we have √z > 0 and √-z = +iz, with the branch cut along the negative real axis.

]:

γZ=kp2γx2γy2,
(2)

where γx and γy determine the spatial behaviour in x and y-direction:

γx=mxπLxp,γy=myπLyp,
(3)

with mx and my integers. The bold subscript α= (α 1,α 2,α 3,α 4) is a multi-index that describes four discrete variables: α 1 (or p) denotes the pit number, α 2 indicates the polarization (TE or TM), α 3 is determined by mx and my and α 4 specifies whether the mode is travelling upwards or downwards.

Because the matching conditions at the interfaces z = ±D/2 only involve the x- and y- component of the fields, it is convenient to introduce the following notation:

exyzîZ×[îZ×Exyz]=(ExEy0),
(4a)
hxyzîZ×Hxyz=(HyHx0),
(4b)

with î z the unit vector in the z-direction. In this way, the lower case e and h are the rotated transverse components of the electric and magnetic field. Furthermore, we split the transverse components of the modes into a real part that depends on x and y and a complex part that depends on z:

[eαxyzhαxyz]=υᾱxy[ηα(Z)ζα(Z)],
(5)

where the subscript α ̄ = (α 1,α 2,α 3) and thus the transverse vectorfield vα ̄ do not depend on the direction of propagation of the mode.

We normalize the real parts of the modes by [22

22. J.D. Jackson, Classical Electrodynamics, 3rd ed. (Wiley, New York, 1999).

]:

υᾱυᾱΩp∫∫Ωp[υᾱ,xxyυᾱ,xxy*+υᾱ,yxyυᾱ,yxy*]dxdy=1,
(6)

where the superscript * denotes complex conjugation [23

23. Although να ̄(x,y) is real, we include the conjugation for consistency of the notation.

]. Furthermore, the modes are orthogonal such that for different modes α ̄ and α ̄′:

υᾱυα′̄Ωp=0,ifᾱᾱ
(7)

Note that the time averaged Poynting vector in the z-direction of a mode is given by:

Sα,z=12Re(Eα,xHα,y*Eα,yHα,x*)=12Re(ηαζα*)(υᾱ,xυᾱ,x*+υᾱ,yυᾱ,y*),
(8)

and, hence, the scalar product 〈να ̄|να ̄Ωp a of a mode α ̄ with itself is proportional to the flow of energy of this mode through a plane of constant z.

For a full listing of the waveguide modes, see Appendix A. The mode functions (Eα,H α) are complete in the following sense: any time harmonic electromagnetic field with frequency ωsatisfying the source-free Maxwell equations inside the holes and pits can be expressed as a linear combination of these mode functions. Hence, for z between D/2 and -D/2, we have:

[Epit(r)Hpit(r)]=αaα[Eα(r)Hα(r)],
(9)

for some expansion coefficients aα that will be determined by matching the field inside the holes and pits to the field above and below the conducting layer.

3.2. Above and below the layer

The total electric and magnetic field above and below the conducting layer consist of the (known) incident field, its corresponding reflected field (that results from the incident field when the conducting layer does not contain any pits or holes) and the scattered field (that results from the presence of the pits and holes):

[E(r)H(r)]=[Ei(r)Hi(r)]+[Er(r)Hr(r)]+[Es(r)Hs(r)].
(10)

The reflected field can easily be calculated from the incident field and we therefore consider the sum of the incident and reflected field to be known. The scattered field can be written as:

[Es(r)Hs(r)]=β1bβ[Eβ(r)Hβ(r)]dβ2,
(11)

where the coefficients bβ are still to be determined and (Eβ,H β are plane waves with wave vector k u above the layer and k l below the layer. The transverse components (kx ,ky ) of the wave vector are real and the z-component is given by:

kzu=+ku2kx2ky2,
(12a)
kz=k2kx2ky2.
(12b)

The sign before the square root follows from the assumed time dependence exp(-iωt) and from the fact that the scattered field propagates away from the conducting layer. The subscript β = (β 1,β 2) is a short notation for the polarization (β 1) and the x- and y-component of the wave vector (β 2 = (kx ,ky )). The polarization can either be S or P. S-polarized means that the z-component of the electric field is zero (and thus corresponds to TE polarization inside the holes and pits), while for P-polarization the z-component of the magnetic field is zero (TM polarization). Note that the integral ∫dβ 2 is a shorthand notation for ∫∫dkx dky .

We use Eq. (4) to obtain the transverse components of the plane waves and, as before, we split these into a part that depends on x and y and a part that depends on z:

[eβxyzhβxyz]=υβxy[ηβ(z)ζβ(z)].
(13)

These are given in Appendix B.

Analogous to the normalization of the waveguide mode functions, we normalize νβ such that:

υβυβ2∫∫2[υβ,x(x,y)υβ,x(x,y)*+υβ,y(x,y)υβ,y(x,y)*]dxdy=δβ1β1δ(β2β2).
(14)

Note that the integration is over an entire plane. The first δ is the Kronecker delta and the second is the two-dimensional Dirac delta function:

δ(β2)=14π2ei(kxx+kyy)dxdy,β2=kxky.
(15)

Because of the use of rotated transverse electric and magnetic fields, we have the following convenient relations between the transverse components of the electric and the magnetic field of the plane waves:

hβxyz=kzωμ0eβxyz,β1=S,
(16a)
hβxyz=ωεε0kzeβxyzβ1=P.
(16b)

Furthermore, every time-harmonic solution of Maxwell’s equations with frequency ω in the half spaces z > D/2 and z < -D/2, that propagates away from the conducting layer can be expanded in terms of the plane waves (Eβ,Hβ). In particular, we have, for the scattered transverse electric field:

esxyz=β1bβeβxyz2,
(17)

where, as before, the integration over β 2 is a short-hand notation for integrating over kx and ky .

By using Eq. (14) and (16) we define an operator A that works on any two-dimensional vectorfield f : ℝ2 → ℂ2:

A(f)kzωμ0fυβ2S2υβ2s2+ωεε0kzfυβ2P2υβ2Pdβ2,
(18)

where the superscript S or P naturally means that β 1 = S or β 1 = P. This operator is basically the integral version of the operator kωμ0× × that can be applied to the electric field of a plane wave to calculate the corresponding magnetic field. Please note that the factor ωεε0kz in the second integral is singular for kx2 + ky2 = ω 2 ϵϵ 0 μ 0. This is, of course, the fingerprint of the coupled electromagnetic surface wave. Although the integrand is integrable, in the numerical implementation prudence is necessary.

In any plane z is constant and in particular for z = ±D/2, the scattered transverse magnetic field can now be expressed in terms of the electric field:

hs(x,y,±D2)̄=A[es(x,y,±D2)].
(19)

This equation holds for all (x,y) with -∞ < x,y < ∞.

4. Matching at the interfaces

At the interfaces z = ±D/2, we have the following relations for the tangential electric and the tangential magnetic field:

epit=ei+er+es,xy,z=±D2,
(23a)
hpit=hi+hr+hs,xyα1Ωα1,z=±D2,
(20b)

Here, Ωα1 is, as before, the cross-section of the pit or hole that is denoted by index α 1. In Eq. (20a), because the layer is perfectly conducting, the sum of the incident and reflected tangential electric field vanishes at z = ±D/2, hence:

epit=es,xy,z=±D2.
(21)

Using this together with Eq. (19), we have for Eq. (20b):

hpit=hi+hr+A(epit),
(22)

which is valid for all (x,yD/2) within the holes and pits. The waveguide modes that constitute e pit and h pit vanish outside the pits and holes, as indicated by the rectangle function ∏ in Eq. (31) and (32) in Appendix B.

In order to obtain a system of equations that is suitable for numerical implementation, we project this equation on the function ν α¯ by using the scalar product defined in Eq. (34):

α4aαςα(±D2)αaαηα(±D2)A(vᾱ)vᾱΩp=hi+hrvᾱΩp,
(23)

where the summation over α 4 is a summation overthe two directions of propagation (α 4 is not contained in α ̄). This equation is valid for all α ̄′ hence for all α 1 (counting the number of holes and pits), for all α 2 (TE and TM polarization) and for all α 3 (the mode numbers, mx and my ). Consequently, solving the system of Eq. (23) for all α ̄′ and for z = ±D/2 gives the waveguide mode expansion coefficients a α. Note that the term on the right acts as a source term. The factor 〈A (ν α ̄)|ν α ̄′〉Ωp is called the interaction integral. Physically speaking, it describes the interaction of a waveguide mode α ̄, via the scattered plane waves through operator A, with another mode α ̄′. In Appendix C we will discuss some of its properties and a method to compute it numerically.

To obtain an expression for the scattered field, we use Eq. (14) to project Eq. (21) for the tangential electric field on the plane wave mode function e β:

bβ=αaαeαeβ.
(24)

In this way, the scattered field can be expressed in the amplitudes of the modes of the pits and holes:

[Es(r)Hs(r)]=β1αaαeαeβ[Eβ(r)Eβ(r)]2.
(25)

It can be shown that this integrand is integrable everywhere, except, possibly, at the edges of the holes and pits. The occurrence of infinite fields near infinitely sharp, conducting wedges is well-known [24

24. J. van Bladel, Singular Electromagnetic Fields and Sources, 1st ed. (Clarendon Press, Oxford, 1991).

]. For a protruding, right angle, perfectly conducting wedge, the field components perpendicular to the sharp edge may become infinite like r -1/3, where r is the distance to the edge. The field components parallel to the edge remain finite. Furthermore, the charge density always remains finite. At an intruding wedge, like the inner part of a waveguide, all field components remain finite.

With Eq. (23) and (24) we have formulated our three-dimensional vectorial scattering problem as a linear system for the amplitudes of the waveguide modes only. Since these modes are parametrized by two parameters (γx and γy ), we have thus reduced the three-dimensional scattering problem to a two-dimensional numerical problem.

5. Numerical implementation

Of course, when implementing our diffraction problem into a computer code we will have to truncate the infinite series of waveguide modes. For large γx and γy , and depending on the size Lxp and Lyp , the mode of concern will be evanescent in the z-direction. For large imaginary γz , the mode will only penetrate into the hole or pit a very small distance. It is therefore reasonable to expect that only the modes with a small imaginary γz will contribute to the total result. Roberts [17

17. A. Roberts, “Electromagnetic theory of diffraction by a circular aperture in a thick, perfectly conducting screen,” J. Opt. Soc. Am. A 4, 1970–1983 (1987). [CrossRef]

] used 168 modes, while García-Vidal and coworkers [15

15. F.J. Garcia-Vidal, E. Moreno, J.A. Porto, and L. Martin-Moreno, “Transmission of light through a single rectangular hole,” Phys. Rev. Lett. 95, 103,901 (2005). [CrossRef] [PubMed]

] used only one. To compare the results of two calculations, one with a number N and the second with a smaller number Ñ, we define the following measure:

FNN˜=∫∫∫Vp(ε0εpEN˜EN2+μ0HN˜HN2)dxdydz∫∫∫Vp(ε0εpEN˜2+μ0HN˜2dxdydz,
(26)
Fig. 2. Relative error FNN ̃, with N = 2600 as a function of the number of waveguide modes (Ñ = 120,440,960,1680). Setup is a single hole with a perpendicular incident, linearly polarized plane wave. As a reference, 1/N is also plotted.
Fig. 3. Line scan of the absolute value of the x-component of the electric field in a plane with constant z, through the center of the hole, for various numbers of waveguide modes. Setup is a single hole,Lx =Ly =D = λ/5, with a perpendicular incident, linearly polarized plane wave.

where the integration is over the volume Vp of the hole (or pit). Here, (E Ñ,H Ñ) is the elec tromagnetic field inside the hole for which the series are truncated after Ñ waveguide modes and (E N , H N ) is the electromagnetic field inside the hole obtained by truncation after N waveguide modes. Hence, this measure corresponds to the error in the energy.

Fig. 2 shows this error as a function of the number of unknowns for a few typical setups. It is clear that only several hundreds of unknowns per hole or pit are enough to model our three-dimensional problem accurately, provided that the thickness of the layer is not too small.

Fig. 3 shows a line scan of the electric field inside a single hole, for various numbers of waveguide modes. Fig. 3(a) shows the field at the entrance of the hole, whereas Fig. 3(b) shows the field in the middle of the hole. Large numbers of waveguide modes are needed to show the singular behaviour at the rim of the hole, as is clear from Fig. 3(a). Inside the hole, where the fields are smooth and regular, the difference with respect to the electric field between 120 waveguide modes and as much as 2600 waveguide modes is not much more than 1 percent. Hence, if the singular behaviour of the fields dominates the solution, as is the case when the conducting layer is thin as compared to the wavelength, our mode expansion technique is probably not the most suitable method.

A small system of equations with only several hundreds of unknowns is solved on a regular desk top computer in only a few seconds. Consequently, most computing time is spent on calculating the interaction integral, which takes a few hours. As discussed in Appendix C, these integrals contain an exponential factor that oscillates violently when mode α ̄ and α ̄′ live in pits or holes that are far apart. Moreover, the integral contains the factor 1/kz that is singular on the circle given by kx2 + ky2 = k 2. As stated before, this is the fingerprint of the coupled electromagnetic surface wave. The integrand is still integrable, but a careful implementation is required.

However, because the interaction integral only involves the plane z = +D/2, it does not depend on the following important parameters: the thickness D of the conducting layer; the index of refraction np inside the pit or hole; whether the scatterer is a pit or a hole and, in case of a pit, its depth dp . Consequently, once the interaction integrals are calculated for a certain setup, we can vary these parameters with negligible computational effort. This is a great advantage of our method. The possibility to construct a library of calculated interaction integrals is also beneficial.

6. Extraordinary transmission

In this section we discuss our first results. Calculations were done for single as well as multiple holes and pits. For all calculations, we took into account a number of 440 waveguide modes, such that the error in the energy is less than 1 percent. All holes and pits are square (Lx = Ly = L) and the index of refraction above and below the layer as well as inside the pits and holes is taken to be unity.

In the following we (among other things) calulate the energy flux through a hole for various setups. This energy flux through a plane with z = z 0 is calculated directly from the coefficients of the waveguide modes in the following way:

∫∫ΩpSzdxdy=ᾱα4α˜412Re[aᾱα4aᾱα˜4*ηαα4(z0)ςαα˜4(z0)*],
(27)

where α 4 denotes the direction of propagation of the waveguide mode. Note that two waveguide modes that have an opposite direction of propagation but that are otherwise identical together produce a non-zero energy flux.

6.1. Extraordinary transmission through a single hole

Fig. 4 shows the energy flux through a single hole as a function of the layer thickness. The energy flux is normalized by the energy that is incident on the area of the hole (the scalar optics normalization). Results are shown for various sizes of the hole. For sizes where all waveguide modes are below or at cut-off, the amount of energy coming through the hole decreases exponentially with the layer thickness, as expected. When the lowest order modes are just above cut-off a strong modulation of the energy flux with layer thickness is seen. The period of this modulation is half of the effective wavelength (2π/γz ) of the propagating mode, indicating that the interference of this mode with its own reflection is responsible for the increased and decreased transmission. It follows from Fig. 4 that if the lowest order mode is just above cut-off, extraordinary transmission of a factor of about 1.5 seems possible. If the size of the hole is increased further, more and more modes are propagating and the normalized energy flux decreases below unity. Going from L = λ to L = 2λ, the energy flux increases to just below unity. For large holes, one expects an energy flux of unity, of course. The energy flux that is shown, is calculated directly from the coefficients found for the waveguide modes and hence, it is not necessarily the energy that will travel along the z-axis and, possibly, arrive at a far field detector. However, because of the perfect conductor assumption, none of this energy is absorbed.

Fig. 4. Energy flux through single hole as a function of layer thickness, normalized by the energy that is incident on the area of the hole. Incident field is a perpendicular, linearly polarized plane wave. The hole is square: Lx = Ly = L with different values of L for every curve, as listed in the legend. The number at the right of each curve is the number of waveguide modes above cut-off.
Fig. 5. Polar plot of the near field scattering from a single square pit with a depth d = λ/4. For different sizes, the Poynting vector in the radial direction is shown, at a half circle with radius λ, with its center coinciding with the center of the pit, at z = D/2. Black line is for the (y,z)-plane, gray for the (x,z)-plane. Incident field is a perpendicular plane wave, with its electric field linearly polarized along the x-direction. The radial scale is arbitrary, but equal for all three figures.

Fig. 5 gives information on the direction along which most energy is scattered. It shows polar plots of the energy, scattered from a single pit with depth λ/4. Shown is the Poynting vector in the radial direction along a half circle with radius of one wavelength, hence the scattering in the near field. The scattering in the plane in which the incident electric field is polarized can be non-zero along the interface, because the corresponding electric field then points in the z-direction. The scattering in the perpendicular plane can not be along the z-direction, for the tangential field at a perfect conductor must be zero. For a pit with size L = λ/5 the scattering is like that of a dipole, whereas for larger pits the scattering is mainly along the optical axis (the z-axis).

Fig. 6. Energy flux through a hole, normalized by the energy flux through an identical single hole, as a function of the distance between two holes (or one hole and one pit). Pits and holes are all square, with Lx = Ly = L = λ/4. Incident field is a linearly polarized plane wave, with polarization as stated at the top of the figures. Incidence is always perpendicular, except for the dashed and dotted lines in Fig. 6(b). Here, kxi ≠ 0 means that the plane of incidence is the (x,z)-plane. Fig. 6(d) shows cross-sections of the four used geometries. The bold arrows denote where the energy flux is calculated. The thickness of the layer is λ/2 and the depth of the pits is λ/4.

6.2. Extraordinary transmission through multiple holes

In Fig. 6, we show the effect that the presence of a second and third hole or pit have on the energy flux, as calculated with equation (27), through the first hole. As a function of the distance between the centers of the two or three scatterers, the energy flux through one hole is calculated. This energy flux is normalized by the energy flux through an identical, but single hole, also calculated with our method. With this special normalization, we are able to isolate the effect that the presence of the second and third scatterer have on the transmission through the first. The incident field is again a linearly polarized plane wave. To distinguish the two basic directions of polarization we define a reference plane through the line that connects the centers of the scatterers and the z-axis. The electric field of the incident plane wave is either directed perpendicular to this plane or else parallel to this plane. Fig. 6(a) shows the normalized energy flux through one of two holes and through one of three holes for perpendicular polarization. A modulation of the energy flux as a function of distance between the centers of the holes only occurs for holes that are less than two wavelengths apart. We believe that enhanced or decreased transmission in this case is caused by the coupling of evanescent fields scattered from one hole to the other and/or by polarization rotation at the corners of the hole. The solid lines in Fig. 6(b) show the same calculations, but now for parallel polarization. The modulation of the energy flux is now also present for large distances between the holes. Its period is equal to the wavelength of the incident field. Its amplitude for three holes is twice the amplitude for two holes. Furthermore, the amplitude is proportional to the inverse of the distance between the centers of the holes, as expected for a cylindrical wave. If the propagation direction of the incident plane wave is slightly tilted (dashed and dotted line in Fig. 6(b)) then a phase shift occurs that is equal to the delay that the incident field experiences in reaching the farthest hole as compared to the nearest hole. Fig. 6(c) shows, for the incident field polarized parallel, the normalized energy flux through a hole in the presence of a pit. This pit has its open end at the upper side (dashed, black line) or at the lower side (solid, gray line). The field is incident from above. The modulation period for the latter is half that of the first. Furthermore, the amplitude for both cases is much smaller than for the case of two holes (solid, black line).

Fig. 7. Comparison between the scalar optics normalization (solid line) and the single hole normalization (dashed line). The vertical axis shows the energy flux through one of two holes, the horizontal axis the distance between the centers of the two holes. The holes are square, with L = 0.49λ (black) and L = 0.6λ (gray). The thickness of the conducting layer is λ/2.
Fig. 8. Cartesian components of the Poynting vector in a plane of constant z, at a distance of λ/20 below the metal layer. Gray scale is in arbitrary units, the same for all figures. The layer contains two holes with Lx = Ly = λ/4 and the thickness of the layer is λ/2. A parallel polarized plane wave is incident from above, perpendicular incidence. Top figure corresponds to a maximum in energy throughput, lower figure to a minimum. See arrows in Fig. 6(b).

Fig. 7 shows a comparison of the usual scalar optics normalization, where the transmission is normalized by the energy that is incident on the area of the hole, with the normalization used in Fig. 6. The setup is two square holes, at varying distance of each other. The incident field is a parallel incident plane wave. Shown is, again, the transmission through one of these two holes. The transmission, when normalized by the transmission through an identical but single hole, always varies around unity. Hence, the presence of the second hole can lead to an increase as well as a decrease in transmission. However, when the scalar optics normalization is used, we already know from Fig. 4 that the layer thickness has an influence on the transmission. This influence is of another nature. For holes that are so small that all modes are evanescent (shown in black), the transmission decreases exponentially with layer thickness. When the size of the holes is such that at least one mode is propagating (shown in gray) the Fabry-Perot effect can cause enhanced transmission for a single hole.

We think that, for the cases in which the incident electric field is polarized perpendicular to the reference plane, the enhanced and reduced transmission that occurs when the single hole normalization is used, is mainly caused by waves with kz = 0 that are scattered along the metal surface. These scattered waves cause a periodicity of a wavelength when two scatterers at the same side of the metal layer contribute and a periodicity of half the wavelength when there is only one source, as is the case when one hole is accompanied by a pit with its open end at the non-illuminated side. In this case, the surface wave that is excited at the exit of the hole travels to the pit. It there excites another surface wave that travels back to the hole and interferes. This results in a phase shift that corresponds to twice the distance between the hole and the pit. See also Ref. [9

09. H.F. Schouten, N. Kuzmin, G. Dubois, T.D. Visser, G. Gbur, P.F.A. Alkemade, H. Blok, G.W. Hooft, D. Lenstra, and E.R. Eliel, “Plasmon-assisted two-slit transmission: Young’s experiment revisited,” Phys. Rev. Lett. 94, 053,901 (2005). [CrossRef] [PubMed]

]. The doubling of the amplitude for three holes as compared to two holes also fits nicely in the picture of the interfering surface waves.

For a perfect conductor, a wave propagating along the surface with its wave vector equal to the wave vector of the incident light, is the analogue of a surface plasmon. It is often stated that perfect metals do not support surface plasmons. However, we believe that the expression surface plasmon (or surface plasmon polariton) is confusing. The ending -on suggests a sort of localization. For a surface plasmon, as described by, for example, Raether [25

25. H. Raether, Surface Plasmons on Smooth and Rough Surfaces and on Gratings, 1st ed. (Springer-Verlag, Berlin, 1988).

], this means that it is bound to the interface between the metal and the dielectric. The derivation by Raether of the surface plasmon wave vector (which is only valid for a flat interface) is also valid for an interface between a perfect conductor and a dielectric. Then, the plasmon wave vector is just the wave vector of the light. The penetration depth inside the metal is zero and the charges inside the metal oscillate only in the plane of the interface between the metal and the dielectric. The plasmon - or, better, coupled electromagnetic surface wave - then has a constant field strength in the half space above the metal and is not bound to the surface. The existence and the physical nature of the phenomenon, however, are the same for both a conductor with finite conductivity and an idealized perfect conductor. The only difference is the finite decay length and absorption caused by finite conductivity.

Fig. 8 shows the three cartesian components of the Poynting vector in a plane that is λ/20 below the metal layer. The metal layer contains two square holes. The small arrows in Fig. 6(b) indicate the data points for which we calculated the near field as shown in the upper and lower figure. Hence, the incident electric field is polarized in the x-direction. This means that the wave along the surface is mainly propagating energy in the x-direction. For the two setups, the difference is indeed largest for the x-component of the Poynting vector. Note that, near the holes, especially in the upper figure, the z-component of the Poynting vector points towards the metal layer.

7. Conclusion

We have described a mode expansion method to rigorously calculate the diffraction by a perfectly conducting layer with finite thickness, that contains rectangular pits and/or holes. The electromagnetic field above and below the conducting layer is written as an integral over plane waves. These plane waves can be S- or P-polarized and they can be propagating or evanescent. The field inside the pits and holes is expanded into waveguide modes, that can have TE or TM polarization and can also be propagating or evanescent. A system of equations is derived by matching the tangential field components at the interfaces. As unknowns, this system only contains the expansion coefficients of the waveguide modes and is therefore very small. For each pit or hole with a size that is of the order of the wavelength of the incident light, about 400 waveguide modes are sufficient for an accuracy in the energy of less than one percent. Once the system of equations of a particular diffraction geometry is composed, important parameters such as the thickness of the conducting layer and the index of refraction inside the pit or hole can be varied with negligible computational effort.

We have shown results on extraordinary transmission through single as well as multiple holes. For a single hole, with transverse sizes such that the lowest order waveguide mode is just above cut-off, a strong modulation of the transmitted energy as a function of layer thickness was obtained. This is due to the interference of this lowest order mode with its own reflection (the Fabry-Perot effect). For two or three holes, a strong modulation in the transmitted energy was found as a function of the distance between neighbouring holes, but only for one polarization of the incident field. Since we normalize the transmitted energy by the energy transmitted by a single, identical hole, the displayed increase and decrease is only due to the presence of the second (and third) hole. The fact that this modulation is hardly present when the incident field is polarized in the perpendicular direction clearly points to a plasmon mechanism. Actually, the term coupled electromagnetic surface wave would be more appropriate, as the boundedness of this kind of wave depends on the conductivity of the metal.

Except for absorption, all relevant physics is present in the model. This means that, apart from the microwave and terahertz frequency regimes, the results found can also be of use in the optical domain, as long as the penetration depth inside the metal is not too large as compared to the thickness of the metal layer and the distances between the holes or pits.

A. The waveguide modes

We here give a complete listing of the waveguide modes. The rotated transverse components of the modes are split into a transverse part and a part that depends on z:

[eαxyzhαxyz]=υᾱxy[ηα(z)ζα(z)].
(28)

With:

Γα=γx2+γy2,
(29)

and

Λα={2(LxpLyp)12,2(LxpLyp)12,ifmx0andmy0,ifmx=0ormy=0,
(30)

we first define the following auxiliary functions:

υᾱx̄pȳp={ΛαΓαx̄pȳp[γycos(γyx̄p)sin(γyȳp)γxsin(γxx̄p)cos(γyȳp)],α2=TE,ΛαΓα(x̄p,ȳp)[γxcos(γxx̄ p)sin(γyȳ p)γysin(γxx̄ p)cos(γyȳ p)],α2=TM,
(31)
ϑᾱx̄ pȳp={iΛαΓαx̄pȳpcos(γxx̄p)cos(γyȳp),α2=TE,iΛαΓαx̄pȳpsin(γxx̄p)sin(γpȳp),α2=TM,
(32)

where we have introduced local coordinates for every pit or hole: px - x0p , ȳp = y - y0p. Furthermore, the function ∏( p ,ȳ p ) is a rectangle function that indicates that the mode functions are identical to zero outside the cross-sectional area of the p-th hole:

[H(x̄p)H(x̄pLxp)][H(ȳp)H(ȳpLyp)],
(33)

where H(x) is the Heaviside step function. By using (31) it is easy to see that two different modes are orthogonal in the sense:

υᾱυᾱΩp∫∫Ωp[υᾱ,x(x,y)υᾱ,xxy*+υᾱ,y(x,y)υᾱy(x,y)*]dxdy=0,ᾱᾱ.
(34)

Hence, {ν α¯ } is an orthonormal system with respect to the scalar product (34).

The following auxiliary functions are needed for the z-dependent parts η α and ζα in Eq. (5):

fα(z)={exp[z(zz1p)],exp[z(zz2p)],ikpγz1cos(γzz),isin(γzz),γzkpε,α4=,γzkpε,α4=+,γzkp<ε,α4=,γzkp<ε,α4=+,
(35)
gα(z)={exp[z(zz1p)],exp[z(zz2p)],kpγz1sin(γzz),cos(γzz),γzkpε,α4=,γzkpε,αa=+,γzkp<ε,α4=,γzkp<ε,α4=+,
(36)

The constants z1p and z2p are the z-coordinates of the upper and lower end of the pit or hole, respectively. Hence, for holes we have z1p = D/2 and z2p = -D/2, for pits we either have z1p = D/2 and zp 2 = D/2 - dp or we have zp 1 = -D/2 + dp and zp2 = -D/2. We have introduced the constant factors like exp(iγz z1p) to make sure that, for imaginairy γz , the moduli of the exponents are always equal to or smaller than unity. In Eq. (35) and (36) we use the sine and cosine functions instead of the exponential function if the propagation constant γz is very small as compared to kp , with smallness parameter ϵ [26

26. Although the ± for α4 now does not have anything to do with the propagation direction, for consistency, we stick to this notation.

]. If we would not do this, for γz = 0 for some mode, we would miss the mode function that is linear in z and our set of modes would not be complete [27

27. The reader might wonder why we do not use the sine and cosine form always instead of using the exponential form only for |γz /kp | < ϵ. However, for large and purely imaginary γz , the functions cos(γzz) and sin(γzz) increase exponentially with |z|, which is not very convenient for numerical implementation.

]. In our implementation, we take ϵ = 10-5. For the z-dependent parts we now have:

ηα(z)={ωμ0gα(z),γzμ0εpε0fα(z),α2=TE,α2=TM,
(37)
ζα(z)={γzfα(z),kpgα(z),α2=TE,α2=TM.
(38)

Note that the modes are propagating in the z-direction if Γα < kp , while for Γα > kp the modes are evanescent. For a square pit or hole (Lx = Ly ) with Lx < λp /2 all modes are evanescent. Finally, the (rotated) transverse components of the modes are then given by:

[eαxyzhαxyz]=υᾱxy[ηα(z)ζα(z)],
(39)

and the longitudinal components:

Eαxyz={0,μ0εpε0ϑᾱxy(z),α2=TE,α2=TM,
(40a)
Eα,zxyz={ϑᾱxy(z)0,α2=TE,α2=TM.
(40b)

We note here that the normalization of the waveguide modes involves only ν α ̄ , which is the part of the transverse field that does not depend on z. This means that the above z-dependent part is only defined up to a constant. We have chosen this constant such that for both TE and TM polarization the waveguide modes have the same order of magnitude.

B. The plane waves above and below the layer

The transverse components of the plane waves E β,H β are divided into a part that depends on x and y and a part that depends on z:

[eβxyzhβxyz]=υβxy[ηβ(z)ζβ(z)].
(41)

We will give a listing of these functions here. With:

Γβ=kx2+ky2,
(42)

and

Λβ=12π,
(43)

we define the following functions that are independent of z:

υβxy={ΛβΓβei(kxx+kyy)(kykx),β1=S,kx2+ky2>0,Λβ(01),β1=S,kx2+ky2=0,ΛβΓβei(kxx+kyy)(kxky)β1=P,kx2+ky2>0,Λβ(10),β1=P,kx2+ky2>0,
(44)
ϑβxy=ΛβΓβei(kxx+kyy),
(45)

For the z-dependent part we have the following auxiliary function:

fβ(z)={exp[ikzu(zD2)],z>D2,exp[ikzl(z+D2)],z<D2,
(46)

and the actual z-dependent parts:

ηβ(z)={ωμ0fβ(z),β1=S,kzμ0εε0fβ(z),β1=P,
(47)
ζβ(z)={kzfβ(z),β1=S,kfβ(z),β1=P.
(48)

Note that k, kz and ϵ in the above equations are either ku , kzu and ϵu or kl , kzl and ϵl , depending on z. The z-components of the plane waves are given by:

Eβ,z(x,y,z)={0,β1=S,μ0εε0ϑβ(x,y)fβ(z),β1=P,
(49a)
Hβ,zxyz={ϑβxyfβ(z)β1=S,0,β1=P.
(49b)

C. The interaction integral

In this appendix we elaborate on the interaction integral that describes the interaction of waveguide mode α ̄ , via the scattered field through operator A , with another waveguide mode α ̄′. We will first write out the integral and then we will discuss the numerical implementation.

C.1. The interaction integral

The interaction integral is given by:

A(υᾱ)υᾱΩp=kzωμ0υαυβ2SΩpυαυβ2SΩp*dβ2+ωεε0kzυαυβ2PΩpυαυβ2PΩp*dβ2.
(50)

Please recall that β 2 = (kx ,ky ) and that ∫dβ 2 = ∫∫dkx dky . We consider one of the scalar products, with β= (β 1,β 2) and β 1 = S,P:

υᾱυβΩp=x0px0p+Lxpy0py0p+Lypυᾱxyυβxy*dxdy,
ei(kxx0p+kyy0p)0Lpx0Lpyυᾱx̄̄pȳpυβx̄pȳp*d̄xpd̄x̄p,
(51)

where we have changed to local coordinates. This final double integral is in fact a Fourier integral that can be calculated analytically:

Fᾱβ1(kx,ky)0Lxp0Lypυᾱ(x̄p,ȳp)υβ̄(x̄p,ȳp)*dx̄pdx̄p,
=ΛαΛβΓαΓβ{γykycmxp(kx)smyp(ky)+γxkxsmxp(kx)cmyp(ky),α2=TE,β1=S,γxkycmxp(kx)smyp(ky)γykxsmxp(kx)cmyp(ky)=0,α2=TE,β1=S,γykxcmxp(kx)smyp(ky)γxkysmxp(kx)cmyp(ky),α2=TE,β1=P,γxkxcmxp(kx)smyp(ky)+γykysmxp(kx)cmyp(ky),α2=TE,β1=P,
(52)

where the functions cmjp and smjp with subscript j = x,y are given by:

cmjp(kj)0Ljpcos(γjz)eikjzdz={ikjγj2kj2[1(1)mjeikjLj],kj±γj,12Lj,kj=±γj,kj0,Lj,kj=γj=0,
(53a)
smjp(kj)0Ljpsin(γjz)eikjzdz={γjγj2kj2[1(1)mjeikjLj],kj±γj,12iLj,kj±γj,kj0,0,kj=γj=0.
(53b)

Note that F α ̄ β1(kx ,ky ) is zero for all (kx ,ky ) when the waveguide mode is TM polarized and the plane wave is S-polarized. Hence, these polarizations do not interact. For the interaction integral we now have:

Fig. 9. The division of the integration area into 12 domains. Not on scale.
A(υᾱ)υᾱΩp̄=kzωμ0ei(kxΔx+kyΔy)FᾱS(kx,ky)Fᾱ'S(kx,ky)*dkxdky+ωεε0kzei(kxΔx+kyΔy)FᾱP(kx,ky)FᾱP(kx,ky)*dkxdky,
(54)

with Δ x = xp0 - x0p and Δ y = yp 0 - y0p, the distance between pit p′ and p. From the above equation it is clear that this double integral is difficult for two reasons. First, the factor exp[i(kx Δ x + ky Δ y )] oscillates violently when pit P′ and p are far apart. Second, the factor kz1 is singular for kx2+ky2 = k 2. The integrand is still integrable, but care has to be taken.

C.2. Numerical computation of the interaction integrals

Since we have to calculate a lot of interaction integrals, we have to find a method that is both fast and accurate. To select the best numerical integration routine we concentrate on the two most difficult parts of the interaction integral: the circle in the (kx ,ky )-plane where the square root kz = (k 2 - kx2 - ky2)1/2 is equal to zero and the exponential factor exp(ikx Δ x + iky Δ y ) that oscillates violently when the two holes or pits under consideration are far apart. The square root term would be best tackled with polar coordinates, whereas the exponential term would be easier to integrate with cartesian coordinates. To overcome this problem, we split the integration area in 12 domains. See Fig. 9. From symmetry properties of the integrand, it follows that it suffices to integrate over half the (kx ,ky )-plane [28

28. The interaction integral also has some other properties that can save a lot of computation time. These properties are outside the scope of this paper.

]. Furthermore, in the domains 1, 6 and 11 that are situated within the circle kx2+ky2 = k 2, only the real parts need to be calculated. For the other domains, we only need the imaginary parts:

A(υᾱ)υᾱΩp=NNinmfRe(IN)+NNoutmfiIm(IN),
(55)

with mf a multiplication factor and ,N a set of domain numbers:

forα1=α1:mf=4,Nin={1,11},Nout={2,3,4,5,12},
(56a)
forα1α1:mf=2,Nin={1,6,11},Nout={2,3,4,5,7,8,9,10,12},
(56b)

and with the double integral over one of the domains. Here, I α ̄α ̄, is shorthand notation for the integrand of the interaction integral.

IN=∫∫domainNIαᾱ̄(kx,ky)dkxdky.
(57)

The two half rings, domains 11 and 12, are called the inner d-ring [(k - δ)2k 2 + k 2k 2] and the outer d-ring [(k 2 <kx2+ky2 ≤ (k+δ)2]. The value of δ is chosen such, that the number of oscillations of the exponential factor inside the rings is small. However, the d-rings should be wide enough to contain the steepest part of the square root factor (typically larger than or equal to k/10). Within the two d-rings, we choose a polar coordinate system. Furthermore, we apply a substitution to get rid of the square root singularity. We then use a standard, adaptive quadrature routine from the NAG foundation toolbox for Matlab, D01FCF [29

29. This routine is a multi-dimensional adaptive quadrature over a hyper-rectangle. For a description, see the internet link: http://www.nag.com/nagware/mt/doc/d01fcf.html.

].

Regarding the domains 1 to 10, we assume that the square root factor is sufficiently flat. In these areas, we split the integrand in a slowly varying part and the (possibly) quickly oscillating exponential factor. Domains 1, 2, 6 and 7 are the areas that are bounded by the inner and outer d-ring. We approximate this boundary linearly on a cartesian grid. The slowly varying part is also approximated linearly, on a cartesian grid. This linear approximation times the exponential factor can now be integrated exactly. The domains 3, 4, 5 and 8, 9, 10 are rectangular areas, bounded by K 1 and k max. Here, we approximate the slowly varying part parabolically and we integrate exactly this parabolic approximation times the exponent.

Acknowledgements

This research was supported by the Dutch Technology Foundation STW.

References and links

01.

H.A. Bethe, “Theory of diffraction by small holes,” Phys. Rev. 66, 163 (1944). [CrossRef]

02.

J. Meixner and W. Andrejewski, “Strenge Theorie der Beugung ebener elektromagnetischer Wellen and der vol-lkommen leitenden Kreisscheibe und an der kreisformigen Offnung im vollkommen leitenden ebenen Schirm,” Annalen der Physik 7, 157–168 (1950). [CrossRef]

03.

C. Flammer, “The vector wave function solution of the diffraction of electromagnetic waves by circular disks and apertures. I. Oblate spheroidal vector wave functions,” J. Appl. Phys. 24, 1218–1223 (1953). [CrossRef]

04.

C. Flammer, “The vector wave function solution of the diffraction of electromagnetic waves by circular disks and apertures. II. The diffraction problems,” J. Appl. Phys. 24, 1224–1231 (1953). [CrossRef]

05.

C.J. Bouwkamp, “Diffraction theory,” Reports on progress in physics 17, 35–100 (1954). [CrossRef]

06.

T.W. Ebbesen, H.J. Lezec, H.F. Ghaemi, T. Thio, and P.A. Wolff, “Extraordinary optical transmission through sub-wavelength hole arrays,” Nature 391, 667–669 (1998). [CrossRef]

07.

H.F. Schouten, T.D. Visser, D. Lenstra, and H. Blok, “Light transmission through a subwavelength slit: Waveg-uiding and optical vortices,” Phys. Rev. E 67, 036,608 (2003). [CrossRef]

08.

J. Bravo-Abad, L. Martin-Moreno, and F.J. Garcia-Vidal, “Transmission properties of a single metallic slit: From the subwavelength regime to the geometrical-optics limit,” Phys. Rev. E 69, 026,601 (2004). [CrossRef]

09.

H.F. Schouten, N. Kuzmin, G. Dubois, T.D. Visser, G. Gbur, P.F.A. Alkemade, H. Blok, G.W. Hooft, D. Lenstra, and E.R. Eliel, “Plasmon-assisted two-slit transmission: Young’s experiment revisited,” Phys. Rev. Lett. 94, 053,901 (2005). [CrossRef] [PubMed]

10.

M.G. Moharam and T.K. Gaylord, “Rigorous coupled-wave analysis of metallic surface-relief gratings,” J. Opt. Soc. Am. A 3, 1780–1787 (1986). [CrossRef]

11.

J.A. Porto, F.J. Garcia-Vidal, and J.B. Pendry, “Transmission resonances on metallic gratings with very narrow slits,” Phys. Rev. Lett. 83, 2845–2848 (1999). [CrossRef]

12.

Q. Cao and P. Lalanne, “Negative role of surface plasmons in the transmission of metallic gratings with very narrow slits,” Phys. Rev. Lett. 88, 057,403 (2002). [CrossRef] [PubMed]

13.

L. Martin-Moreno and F.J. Garcia-Vidal, “Optical transmission through circular hole arrays in optically thick metal films,” Opt. Express 12, 3619–3628 (2004). [CrossRef]

14.

S.H. Chang, S.K. Gray, and G.C. Schatz, “Surface plasmon generation and light transmission by isolated nanoholes and arrays of nanoholes in thin metal films,” Opt. Express 13, 3150–3165 (2005). [CrossRef] [PubMed]

15.

F.J. Garcia-Vidal, E. Moreno, J.A. Porto, and L. Martin-Moreno, “Transmission of light through a single rectangular hole,” Phys. Rev. Lett. 95, 103,901 (2005). [CrossRef] [PubMed]

16.

F.J. Garcia de Abajo, “Light transmission through a single cylindrical hole in a metallic film,” Opt. Express 10, 1475–1484 (2002). [PubMed]

17.

A. Roberts, “Electromagnetic theory of diffraction by a circular aperture in a thick, perfectly conducting screen,” J. Opt. Soc. Am. A 4, 1970–1983 (1987). [CrossRef]

18.

J.M. Brok and H.P Urbach, “A mode expansion technique for rigorously calculating the scattering from 3D subwavelength structures in optical recording,” J. Mod. Opt. 51, 2059–2077 (2004). [CrossRef]

19.

J.B. Pendry, L. Martin-Moreno, and F.J. Garcia-Vidal, “Mimicking surface plasmons with structured surfaces,” Science 305, 847 (2004). [CrossRef] [PubMed]

20.

A.P Hibbins, B.R. Evans, and J.R. Sambles, “Experimental verification of designer surface plasmons,” Science 308, 670 (2005). [CrossRef] [PubMed]

21.

Here and henceforth the square root of a complex number z is defined such that for real z > 0 we have √z > 0 and √-z = +iz, with the branch cut along the negative real axis.

22.

J.D. Jackson, Classical Electrodynamics, 3rd ed. (Wiley, New York, 1999).

23.

Although να ̄(x,y) is real, we include the conjugation for consistency of the notation.

24.

J. van Bladel, Singular Electromagnetic Fields and Sources, 1st ed. (Clarendon Press, Oxford, 1991).

25.

H. Raether, Surface Plasmons on Smooth and Rough Surfaces and on Gratings, 1st ed. (Springer-Verlag, Berlin, 1988).

26.

Although the ± for α4 now does not have anything to do with the propagation direction, for consistency, we stick to this notation.

27.

The reader might wonder why we do not use the sine and cosine form always instead of using the exponential form only for |γz /kp | < ϵ. However, for large and purely imaginary γz , the functions cos(γzz) and sin(γzz) increase exponentially with |z|, which is not very convenient for numerical implementation.

28.

The interaction integral also has some other properties that can save a lot of computation time. These properties are outside the scope of this paper.

29.

This routine is a multi-dimensional adaptive quadrature over a hyper-rectangle. For a description, see the internet link: http://www.nag.com/nagware/mt/doc/d01fcf.html.

OCIS Codes
(050.1220) Diffraction and gratings : Apertures
(050.1960) Diffraction and gratings : Diffraction theory

ToC Category:
Diffraction and Gratings

History
Original Manuscript: January 27, 2006
Revised Manuscript: March 16, 2006
Manuscript Accepted: March 16, 2006
Published: April 3, 2006

Citation
J. M. Brok and H. P. Urbach, "Extraordinary transmission through 1, 2 and 3 holes in a perfect conductor, modelled by a mode expansion technique," Opt. Express 14, 2552-2572 (2006)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-14-7-2552


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. H.A. Bethe, "Theory of diffraction by small holes," Phys. Rev. 66, 163 (1944). [CrossRef]
  2. J. Meixner and W. Andrejewski, "Strenge Theorie der Beugung ebener elektromagnetischer Wellen and der vollkommen leitenden Kreisscheibe und an der kreisformigen Offnung im vollkommen leitenden ebenen Schirm," Annalen der Physik 7, 157-168 (1950). [CrossRef]
  3. C. Flammer, "The vector wave function solution of the diffraction of electromagnetic waves by circular disks and apertures. I. Oblate spheroidal vector wave functions," J. Appl. Phys. 24, 1218-1223 (1953). [CrossRef]
  4. C. Flammer, "The vector wave function solution of the diffraction of electromagnetic waves by circular disks and apertures. II. The diffraction problems," J. Appl. Phys. 24, 1224-1231 (1953). [CrossRef]
  5. C.J. Bouwkamp, "Diffraction theory," Reports on progress in physics 17, 35-100 (1954). [CrossRef]
  6. T.W. Ebbesen, H.J. Lezec, H.F. Ghaemi, T. Thio, and P.A. Wolff, "Extraordinary optical transmission through sub-wavelength hole arrays," Nature 391, 667-669 (1998). [CrossRef]
  7. H.F. Schouten, T.D. Visser, D. Lenstra, and H. Blok, "Light transmission through a subwavelength slit: Waveguiding and optical vortices," Phys. Rev. E 67, 036,608 (2003). [CrossRef]
  8. J. Bravo-Abad, L. Martin-Moreno, and F.J. Garcia-Vidal, "Transmission properties of a single metallic slit: From the subwavelength regime to the geometrical-optics limit," Phys. Rev. E 69, 026,601 (2004). [CrossRef]
  9. H.F. Schouten, N. Kuzmin, G. Dubois, T.D. Visser, G. Gbur, P.F.A. Alkemade, H. Blok, G.W. Hooft, D. Lenstra, and E.R. Eliel, "Plasmon-assisted two-slit transmission: Young’s experiment revisited," Phys. Rev. Lett. 94, 053,901 (2005). [CrossRef] [PubMed]
  10. M.G. Moharam and T.K. Gaylord, "Rigorous coupled-wave analysis of metallic surface-relief gratings," J. Opt. Soc. Am. A 3, 1780-1787 (1986). [CrossRef]
  11. J.A. Porto, F.J. Garcia-Vidal, and J.B. Pendry, "Transmission resonances on metallic gratings with very narrow slits," Phys. Rev. Lett. 83, 2845-2848 (1999). [CrossRef]
  12. Q. Cao and P. Lalanne, "Negative role of surface plasmons in the transmission of metallic gratings with very narrow slits," Phys. Rev. Lett. 88, 057,403 (2002). [CrossRef] [PubMed]
  13. L. Martin-Moreno and F.J. Garcia-Vidal, "Optical transmission through circular hole arrays in optically thick metal films," Opt. Express 12, 3619-3628 (2004). [CrossRef]
  14. S.H. Chang, S.K. Gray, and G.C. Schatz, "Surface plasmon generation and light transmission by isolated nanoholes and arrays of nanoholes in thin metal films," Opt. Express 13, 3150-3165 (2005). [CrossRef] [PubMed]
  15. F.J. Garcia-Vidal, E. Moreno, J.A. Porto, and L. Martin-Moreno, "Transmission of light through a single rectangular hole," Phys. Rev. Lett. 95, 103,901 (2005). [CrossRef] [PubMed]
  16. F.J.G.I. de Abajo, "Light transmission through a single cylindrical hole in a metallic film," Opt. Express 10, 1475-1484 (2002). [PubMed]
  17. A. Roberts, "Electromagnetic theory of diffraction by a circular aperture in a thick, perfectly conducting screen," J. Opt. Soc. Am. A 4, 1970-1983 (1987). [CrossRef]
  18. J.M. Brok and H.P. Urbach, "A mode expansion technique for rigorously calculating the scattering from 3D subwavelength structures in optical recording," J. Mod. Opt. 51, 2059-2077 (2004). [CrossRef]
  19. J.B. Pendry, L. Martin-Moreno, and F.J. Garcia-Vidal, "Mimicking surface plasmons with structured surfaces," Science 305, 847 (2004). [CrossRef] [PubMed]
  20. A.P. Hibbins, B.R. Evans, and J.R. Sambles, "Experimental verification of designer surface plasmons," Science 308, 670 (2005). [CrossRef] [PubMed]
  21. Here and henceforth the square root of a complex number z is defined such that for real z > 0 we have √z > 0 and √z = +i√z, with the branch cut along the negative real axis.
  22. J.D. Jackson, Classical Electrodynamics, 3rd ed. (Wiley, New York, 1999).
  23. Although υ¯α(x,y) is real, we include the conjugation for consistency of the notation.
  24. J. van Bladel, Singular Electromagnetic Fields and Sources, 1st ed. (Clarendon Press, Oxford, 1991).
  25. H. Raether, Surface Plasmons on Smooth and Rough Surfaces and on Gratings, 1st ed. (Springer-Verlag, Berlin, 1988).
  26. <other>. Although the ± for α4 now does not have anything to do with the propagation direction, for consistency, we stick to this notation. </other>
  27. The reader might wonder why we do not use the sine and cosine form always instead of using the exponential form only for | γz/kp|< ε. However, for large and purely imaginary γz, the functions cos (γzz) and sin (γzz) increase exponentially with |z|, which is not very convenient for numerical implementation.
  28. The interaction integral also has some other properties that can save a lot of computation time. These properties are outside the scope of this paper.
  29. This routine is a multi-dimensional adaptive quadrature over a hyper-rectangle. For a description, see the internet link: http://www.nag.com/nagware/mt/doc/d01fcf.html.

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited