1. Introduction

The scattering from a single sub-wavelength hole in a conducting metal layer was already modelled in the 1940s and 1950s [1, 2, 3, 4, 5]. 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], 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.

To model this collective effect and to determine its importance, a lot of geometries have been studied and a lot of methods have been used. To name just a few: Infinitely long, single slits [7, 8] and double slits [9] in metals with finite conductivity were modelled with the Green’s function approach. Coupled-wave analysis has been used to study gratings of slits (see, for example, [10, 11, 12]). Gratings that are periodic in two dimensions were studied with a modal method [13] and with the finite-difference time-domain method [14]. The modal method [15] and the boundary element method [16] were used for a single rectangular hole in a perfect conductor.

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.

Download Full Size | PDF

We present a rigorous modal method, similar to the methods in Refs. [17, 15, 18], that turns our three-dimensional diffraction problem into a two-dimensional numerical problem. Using this method, we calculate the transmission through one, two and three holes and we determine the influence that a second and third hole have on the transmission through the first. Parameters such as size of the hole, distance between two holes and thickness of the conducting layer are varied.

For the modal method that we use, the perfect conductor assumption is essential. Hence, our method is of quantitative value when the imaginary part of the index of refraction is so large, that the skin depth of the metal is small compared to the typical length scales of the geometry. This is for example the case in the terahertz and microwave frequency regimes, where also a lot of research is done in sub-wavelength hole arrays [19, 20]. However, we argue that, apart from absorption of energy in the metal, all relevant physics is present in our model. That also includes the surface plasmon equivalent of a perfect conductor: an electromagnetic surface wave. Thus, the understanding we gain from our results will also help us to understand extraordinary transmission in the optical regime.

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.

This paper is organized as follows. In the next section, we will define the problem under consideration. Then, we will describe the field above and below the layer, as well as inside the holes by mode expansions. Matching these expressions gives the system of equations that needs to be solved. After a few words on the numerical implementation, we will show results on single as well as multiple holes and pits.

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 d_p < D. The sub- or superscript p denotes the number of the pit or hole. The lengths in the x- and y-direction are $L_x^p$ and $L_y^p$ , respectively. The cross-section of a hole or a pit p is given by Ω _p = {(x,y)|$x_0^p$<x<$x_0^p$+$L_x^p$ ,$y_0^p$<y<$y_0^p$+$L_y^p$ }. The halfspaces z > D/2 and z < -D/2 are filled with homogeneous dielectrics with index of refraction n_u and n_l , respectively. Every hole and pit is filled with a homegeneous dielectric with index of refraction n_p . The corresponding relative permeabilities are ϵ_u = $n_u^2$, ϵl = $n_l^2$ and ϵ_p =$n_p^2$. The magnetic permeability is μ ₀ 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 = λ/n_u , λ_l = λ/n_l and λ_p = λ/n_p . The corresponding wave vectors are k_u = 2π/λ_u , k_l = 2π/λ_l and k_p = 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 each pit or hole, the electromagnetic field is expanded in a set of propagating or evanescent waveguide modes. In the next section, we describe these modes, that are characterized by the geometry of the pit or hole they live in, by their polarization, by their spatial frequency and by their direction of propagation.

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:

with propagation constant γ_z given by [21]:

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

with m_x and m_y integers. The bold subscript α= (α ₁,α ₂,α ₃,α ₄) is a multi-index that describes four discrete variables: α ₁ (or p) denotes the pit number, α ₂ indicates the polarization (TE or TM), α ₃ is determined by m_x and m_y and α ₄ 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:

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:

where the subscript α ̄ = (α ₁,α ₂,α ₃) 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]:

where the superscript ^* denotes complex conjugation [23]. Furthermore, the modes are orthogonal such that for different modes α ̄ and α ̄′:

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

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:

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):

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:

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 (k_x ,k_y ) of the wave vector are real and the z-component is given by:

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 β = (β ₁,β ₂) is a short notation for the polarization (β ₁) and the x- and y-component of the wave vector (β ₂ = (k_x ,k_y )). 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β ₂ is a shorthand notation for ∫∫dk_x dk_y .

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:

These are given in Appendix B.

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

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:

Analogous to the scalar product for the waveguide mode functions, the scalar product of two identical plane waves 〈ν_β|ν_β〉_ℝ² is related to the flow of energy of the plane wave through a plane of constant z.

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:

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:

where, as before, the integration over β ₂ is a short-hand notation for integrating over k_x and k_y .

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

where the superscript S or P naturally means that β ₁ = S or β ₁ = P. This operator is basically the integral version of the operator $\frac{k}{{\mathrm{\omega \mu }}_0}\times$ × that can be applied to the electric field of a plane wave to calculate the corresponding magnetic field. Please note that the factor $\frac{\mathrm{\omega \epsilon }{\epsilon }_0}{k_z}$ in the second integral is singular for $k_x^2$ + $k_y^2$ = ω ² ϵϵ ₀ μ ₀. 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:

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:

Here, Ω_α1 is, as before, the cross-section of the pit or hole that is denoted by index α ₁. 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:

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

which is valid for all (x,y,±D/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 ν _{$\overline{\alpha }$ ′} by using the scalar product defined in Eq. (34):

where the summation over α ₄ is a summation overthe two directions of propagation (α ₄ is not contained in α ̄). This equation is valid for all α ̄′ hence for all α ₁ (counting the number of holes and pits), for all α ₂ (TE and TM polarization) and for all α ₃ (the mode numbers, m_x and m_y ). 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 _β′:

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

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]. 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 $L_x^p$ and $L_y^p$ , 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] used 168 modes, while García-Vidal and coworkers [15] 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:

 

Fig. 2. Relative error F_NN ̃, 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.

Download Full Size | PDF

 

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,L_x =L_y =D = λ/5, with a perpendicular incident, linearly polarized plane wave.

Download Full Size | PDF

where the integration is over the volume V_p 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/k_z that is singular on the circle given by $k_x^2$ + $k_y^2$ = k ². 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 n_p inside the pit or hole; whether the scatterer is a pit or a hole and, in case of a pit, its depth d_p . 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 (L_x = L_y = 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 ₀ is calculated directly from the coefficients of the waveguide modes in the following way:

where α ₄ 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: L_x = L_y = 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.

Download Full Size | PDF

 

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.

Download Full Size | PDF

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 L_x = L_y = 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, $k_x^i$ ≠ 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.

Download Full Size | PDF

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.

Download Full Size | PDF

 

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 L_x = L_y = λ/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).

Download Full Size | PDF

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 k_z = 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]. 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], 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:

With:

and

we first define the following auxiliary functions:

where we have introduced local coordinates for every pit or hole: x̄_p ≡ x - $x_0^p$ , ȳ_p = y - $y_0^p$. Furthermore, the function ∏(x̄ _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:

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:

Hence, {ν _{$\overline{\alpha }$} } 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):

The constants $z_1^p$ and $z_2^p$ are the z-coordinates of the upper and lower end of the pit or hole, respectively. Hence, for holes we have $z_1^p$ = D/2 and $z_2^p$ = -D/2, for pits we either have $z_1^p$ = D/2 and z_{p 2} = D/2 - d_p or we have z_{p 1} = -D/2 + d_p and $z_p^2$ = -D/2. We have introduced the constant factors like exp(iγz $z_1^p$) 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 k_p , with smallness parameter ϵ [26]. 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]. In our implementation, we take ϵ = 10^-5. For the z-dependent parts we now have:

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

and the longitudinal components:

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:

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

and

we define the following functions that are independent of z:

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

and the actual z-dependent parts:

Note that k, k_z and ϵ in the above equations are either k_u , $k_z^u$ and ϵ_u or k_l , $k_z^l$ and ϵ_l , depending on z. The z-components of the plane waves are given by:

Note that we did not take special precautions for the case that k_z = 0. For the waveguide modes in the previous section, we made sure that, when it happens that γ_z = 0, the set of mode functions is still complete. The plane waves in the upper and lower half spaces, however, form a continuous set, parametrized by -∞ < k_x ,k_y < ∞. The plane waves with k_z = 0 are only a set of measure zero in the space of all plane waves and are therefore irrelevant for the completeness.

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:

Please recall that β ₂ = (k_x ,k_y ) and that ∫dβ ₂ = ∫∫dk_x dk_y . We consider one of the scalar products, with β= (β ₁,β ₂) and β ₁ = S,P:

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

where the functions $c_{mj}^p$ and $s_{mj}^p$ with subscript j = x,y are given by:

Note that F _{α ̄} ^β₁(k_x ,k_y ) is zero for all (k_x ,k_y ) 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.

Download Full Size | PDF

with Δ _x = x^p ′₀ - $x_0^p$ and Δ _y = y^p′ ₀ - $y_0^p$, 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(k_x Δ _x + k_y Δ _y )] oscillates violently when pit P′ and p are far apart. Second, the factor $k_z^{-1}$ is singular for $k_x^2$+$k_y^2$ = k ². 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 (k_x ,k_y )-plane where the square root k_z = (k ² - $k_x^2$ - $k_y^2$)^1/2 is equal to zero and the exponential factor exp(ik_x Δ _x + ik_y Δ _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 (k_x ,k_y )-plane [28]. Furthermore, in the domains 1, 6 and 11 that are situated within the circle $k_x^2$+$k_y^2$ = k ², only the real parts need to be calculated. For the other domains, we only need the imaginary parts:

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

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

The two half rings, domains 11 and 12, are called the inner d-ring [(k - δ)² ≤k ² + k ² ≤ k ²] and the outer d-ring [(k ² <$k_x^2$+$k_y^2$ ≤ (k+δ)²]. 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].

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 ₁ and k _max. Here, we approximate the slowly varying part parabolically and we integrate exactly this parabolic approximation times the exponent.