We studied the evolution of the Extraordinary Electromagnetic Transmission (EET) through subwavelength hole arrays versus hole size. Here, we show that for large holes EET vanishes and is replaced by another unusual transmission. A specific hole size is found where all the characteristics of the EET vanish and where most usual models fail to describe the transmission except full 3D simulations. The transition between these two domains is characterized by the discontinuity of parameters describing the transmission, in particular the resonance frequency. This transition exhibits a first order phase transition like behavior.
© 2008 Optical Society of America
In recent years, the discovery of an unexpected Extraordinary Electromagnetic Transmission (EET) through arrays of subwavelength holes  has generated numerous experimental and theoretical works. The EET shows an enhancement of several order of magnitude  with respect to usual hole theory [3, 4], and can even exceed the surface ratio occupied by the holes. Since this demonstration, many experiments have been made both to characterize and model this EET, in the optical [5, 6, 7, 8, 9, 10, 11, 12], infrared [13, 14] and terahertz ranges [15, 16, 17, 18, 19, 20]. In the terahertz domain, EET is related to Surface Plasmon Polaritons . These experiments studied the evolution of the EET versus hole diameter and shape [16, 9, 22, 23], lattice geometry [1, 24, 25], film thickness , etc.
Several models have been designed to describe this EET. Most of them are based either on surface plasmons or dynamical diffraction effects. Furthermore, these models are often completed by full three-dimensional (3D) simulations, either Finite Element Modeling (FEM) or Finite Difference Time Domain (FDTD) method.
Here we experimentally studied in the terahertz domain the continuous evolution of the EET through arrays of subwavelength holes versus hole size. We show that for large holes (but still subwavelength dimensions) EET vanishes and is replaced by a transmission of different kind. Most interestingly, the transition between these two domains is characterized by the discontinuity of many parameters describing the transmission, in particular the resonance frequency, and also exhibits a first order phase transition-like behavior. These results have been confronted with several models and simulations, and full 3D FEM simulation was found to be the only one calculation describing the discontinuity. After presenting the experimental results and data analysis, we shall then discuss the origin of the discontinuity and the necessity to use full 3D calculation.
2. Experiments and results
The long-wavelength terahertz domain offers the advantage of a precise control of the geometry. It allows an accurate design of hole shape and dimensions of arrays of subwavelength holes. The realization of ultra-thin (λ/50) nickel plates removes possible wave-guiding Fabry-Pérot resonances. Furthermore, no plasma resonance is present in the terahertz domain, contrary to the visible one. One must notice that a unique set of experimental conditions can be achieved in the terahertz range. Here, the arrays may at the same time have a thickness h very small compared to the wavelength (h≈8µm in our experiments) and large compared to the skin depth δ of the electromagnetic field (δ≪1µm). Thus the arrays may be considered to be almost purely two dimensional (D≫h≫δ where D is the size of the holes). On the contrary in the visible, the skin depth is relatively larger, since it scales as √λ. The system is 3D in the visible range (D≈h≫δ), and our results may not be reproducible in the visible range since such conditions are hardly obtained.
The structures used in the present experiment are made of free-standing electroformed ultrathin nickel plates, with precision of design better than 1 µm (figure 1A). Series of twodimensional arrays of square and round subwavelength holes (54 arrays of 16 by 16 holes) have been designed and analyzed. Terahertz spectra are recorded using standard terahertz time-domain spectroscopy . Broadband linearly polarized subpicosecond single cycle pulses of terahertz radiation are generated and coherently detected by illuminating photoconductive antennas with two synchronized femtosecond laser pulses (figure 1B). The sample is positioned on a 10mm circular hole, in the linearly polarized, frequency independent, 4.8mm-waist (1/e in amplitude) Gaussian THz beam. Numerical Fourier transform of the time-domain signals with and without the sample gives access to the transmission spectrum of the subwavelength structures. The dynamics of the EET is recorded over 240 ps, yielding a raw 3GHz frequency precision after numerical Fourier transform, with 104 signal to noise ratio in 300ms acquisition time. The resonance frequency is obtained by fitting the resonance with a parabolic function taking into account 20 points over the resonance. Therefore the precision is better than the single point precision, and below 1GHz.
Characteristic spectra are presented in figure 2. For small holes, typical features of 2D lattice EET can be observed (figure 2, colored lines) with sharp maxima, asymmetric resonance and resonance frequencies scaling as √2 . A resonance frequency shift versus hole size is also observed. For very tiny holes, first and second resonances agree well with theoretical resonance frequencies given by Bloch wave model at 0.5 and 0.707 THz, respectively , for a lattice period of L=600 µm. However for larger holes the spectra exhibit unusual features (figure 2, black line). The spectra exhibit large, symmetric resonances with resonance frequencies remaining constant with the hole size, and scaling as a factor 2. These are not characteristics of the EET anymore, even though the transmission energy still exceeds the sum of transmission of all the holes. All resonances not scaling as integers have disappeared. It then appears that the spectral features originate from a complex interplay between the geometrical resonance of the array and the shape and size of the holes. The array periodicity provides approximatively the resonance positions (L=600µm gives ν=0.5 THz as first resonance). As previously found by , it can be observed that the increase of the size of the holes shifts the resonance frequencies toward lower frequencies. The size of the individual holes acts as a high pass filtering of the incident wave and this filtering interacts with the lattice structure. Then, the shape of the holes finely further tunes the resonance frequency. We observe that two types of transmissions are possible through these arrays depending on the hole size. For small sizes, the plates behave according to EET characteristics, whereas for large holes, the light differently diffracts on the array.
3. Quantitative study
We quantitatively studied the generation and propagation of EET by analyzing the terahertz transmission spectra of the series of free-standing subwavelength 2D lattices with increasing quantity of metal, while the periodicity remained constant. In order to have unitless parameters so that arrays of different shape could be compared, a filling parameter p is defined as the ratio between the surface of the metal Sm and the surface of the holes Sa in the lattice,
Then for square holes of width a, p=L 2/a 2-1 and for round holes of diameter d, with L the lattice period. The evolution of frequency ν of the first resonance is given by the normalized frequency shift
where ν 0 is the limit resonance frequency for tiny holes, given by Bloch theory as . Furthermore, we define U as the total electromagnetic energy transmitted through the plate
where At and Ao are the measured amplitudes of the electric fields with and without the sample, respectively. We also define a quality factor by . It describes how strongly a small change in the quantity of metal affects the transmission characteristic of the plates. The last parameter ΔW is the full width at half maximum of the first resonance. We focused on the observation and evolution of these parameters, which provide fundamental insights on the nature of the generation and propagation of the EET.
The evolution of Δν, Q and ΔW versus p are given by figure 3 for square and round holes with a lattice period L=600µm. Each point (filling parameter) is provided by a corresponding free-standing plate. For each hole shape, two domains are clearly distinguishable, limited by a discontinuity of all the parameters, and at the same filling parameter pc. However, pc is different for square and round holes: pc=1.03±0.01 for square holes and pc=1.37±0.01 for round holes. Transmission spectra just before and after the transition filling parameter for round holes are shown in figure 4.More precisely, Δν is constant below pc and corresponds to the first order of diffraction. Across pc, Δν jumps to lower values, and then continuously decreases toward 0. The jump is small, respectively 5 and 19GHz for square and round holes, but easily observable and is also supported by the slopes of the curves. Evolution of ΔW is similar, except that it not constant before discontinuity, in particular for round holes. The behavior of Q is more complex. First, the evolution in each domain is non monotonic. The discontinuity is clearly visible for square holes, less important for round holes. All the curves with the same hole shape show discontinuity for the same value of pc. We emphasize the point that pc is different for square and round holes. It is due to different limit conditions between adjacent interacting holes. As a consequence, the distribution of distance (then the quantity of available metal) is different between square and round holes, much sharper for round holes. And the ability for the electromagnetic wave to couple through surface waves is modified.
Experiments have shown an unusual transmission transition through the arrays. We tested several models in order to determine which domains could be described and under which validity conditions (the list of model used is clearly not exhaustive). Since none of the models gave satisfactory results near the discontinuity domain, we also carried out full 3D ab initio finite element simulations for directly solving Maxwell’s equations.
1. A model widely used in the microwave range has been developed by C.C. Chen in 1973 . This model takes the finite thickness (2l 1=8 µm in our experiments) into account. It is based on the superposition theorem, and decomposes the problem into symmetric and antisymmetric plane wave excitations. Both electric and magnetic fields, on both sides of the holes, are expanded into a set of Floquet modes Φpqr, where p and q are the spatial modes, and r denotes TE or TM mode . Each mode has a propagation constant γpq, and a modal wave admittance ξfpqr. Inside the holes, the waves are expressed as waveguide modes Ψmn. The amplitude of the electric field is given by A 00r. By matching the boundary conditions on both sides of the holes,
where E t is the transverse electric field, and Fmnr and Ymnr are the modal coefficients of the waveguide modes.
This model neglects the near field interactions and makes the hypothesis that the electromagnetic wave can be broken into waveguide modes in a subwavelength hole. The calculation are the solid green lines on figure 5A. This model gives a general good shape evolution for the three parameters. It fails rapidly for p<pc, does not catch the discontinuities and Δν diverges rapidly for p≈pc.
2. The Fano model [29, 30, 20] was originally derived to describe autoionization phenomena. Its extension to EET distinguishes two interfering contributions: a continuum non resonant one associated with the direct scattering of the incident field by the subwavelength holes |ψE〉, and a discrete resonant one related to EET excitation |ϕ〉. The transmission is then described by the interaction of a continuum with an isolated resonant level. This model provides an interesting description of the transmission through subwavelength arrays. Although it does not directly deal withMaxwell’s equations, it allows a new approach in modeling the interaction of light with the subwavelength hole arrays.
The total transmission is given by 
where α is a numerical constant, q 2 is proportional to the ratio of the transition probabilities to the isolated level, and to the portion of the continuum that is not interacting. Eφ is the energy of the isolated level and Hilb(VE) is the Hilbert transform of VE=〈ψE|V|ϕ〉, where V is the interaction potential between the isolated level and a level of the continuum with energy E. Contrary to the approximation of the Fano model previously used in EET , we re-introduced the original dependence of VE with respect to E . We assumed a parabolic potential for V which provides a Gaussian wave function for the fundamental level of |ϕ〉, and leads to . Here, A and B are parameters of the extended model . The calculations are the solid red lines on figure 5A. This model gives good fits for high p values, but begins to differ from experimental curves near pc. This model does not catch the discontinuities.
3. In the Fourier Modal Method (FMM) [32, 33, 34], the incident electromagnetic field is decomposed into a modal basis. All modes, even the evanescent ones, are solved through an eigenvalue problem in the Fourier space, which provides the projection of the qth mode wave vector γq along z axis (eigenvalues) and the field distribution Φq(x,y) (eigenvectors). The field components around the hole area are written as
where r=x or y, uq is the amplitude of the upward decaying modal fields, and dq is the amplitude of the downward decaying modal fields. This model also partially neglects near field interactions and is based on modal decomposition. The precision of the FMM and its ability to describe experimental results depend in principle on the limit of the development. In describing our experimental data, we show simulation with a 30th order development. We also checked the convergence of the algorithm for higher order calculation. We found less than 1% variation in field transmission. Surprisingly, this model is also unable to describe the discontinuity between the two kinds of transmission. It is possible that part of the near field interactions are not entirely described, as well as the complete decomposition of the field on the surface of the metal. The calculations are the dashed red lines on figure 5B. This model gives better result for square holes than for round holes.
4. A model based on surface plasmon scattering [35, 36] was also considered. An incoming electric field E in of wave vector k in is decomposed near the holes into a two dimensional surface wave that is scattered away from the hole as a spherical waves with and εm the dielectric constant of the metal. The total surface plasmon field is evaluated as the sum of all fields emitted from each holes,
with e SP is a unitary polarization vector, δ n is the nth source from the center of the array. This model is based on the assumption that surface plasmons are excited by the plane wave, or that the arrays generate surface plasmon like interactions . Near field interactions are neglected. The calculations are the dashed green lines on figure 5B. This model is in good agreement with experimental results for p>pc. Results are still good near pc. But it does not catch the discontinuities and fails to describe experimental results for p<pc.
5. Finally, as none of the considered models were able to fit the discontinuities we performed full 3D numerical simulations. We carried out a direct local resolution of Maxwell’s equations through a full 3D ab initio Finite Element Method (FEM) analysis [20, 38] of the electric field propagating through the array. This method allows the calculation of the transmitted THz electromagnetic field and takes into account the nearfield effects on the array. To reduce the size of the simulation box, we used a unitary cell with one hole in its center, with adequate boundary conditions. The precision of the simulations are controlled by progressively reducing the adaptive mesh size, in particular in the subwavelength holes. Typical mesh dimensions are λ/700 in the holes and λ/50 outside, yielding to a total precision on the transmitted electric field of 0.1%. The relative permittivity of nickel is ε=-9.7×103+1.1×105 i, and relative permeability is 100 (in Gaussian units) [39, 40]. The calculations are the solid blue lines on figure 5B. The FEM analysis is the only one to fit the three discontinuities on square and round holes. The simulations give the general evolution for Δν and ΔW with small shifts, which are the consequence of multiple reflections of the electromagnetic waves inside the simulation box. The evolution of Q is also well described by the simulation.
We shall now discuss here the origin of the observed discontinuity. For small apertures, the transmission properties are given by a complex interplay between waves diffracted by the holes, and waves at the surface of the metal. When the aperture size increases, less metal remains to guide the coupling between holes, up to a point where the surface waves can no more couple the holes and where a transition between two coupling regime appears. An interesting parallel can be drawn with waveguides, even though the plate thickness is negligible compared to the wavelength in our experiment. The frequency cut-off of a hollow waveguide provides the limit wavelength above which the wave can no more propagate. It is given by the cross section of the waveguide. For circular aperture, the wavelength cut-off of the TE10 mode is λc=2.61d and is equal to 575µm at pc=1.37 (d=220µm). For symmetry reasons, the corresponding mode for square apertures is TE11. This is justified by the invariance of the shape of the transmission spectra with respect to the incident polarization angle we experimentally observed, as well as in . Here, the cut-off is λc=√2a=595µm with pc=1.03 (a=421µm). Both cutoff wavelengths for square and round holes are very close to the lattice period of 600µm. Therefore, the transition seems to arise when the quantity of metal goes below the quantity given by the corresponding cut-off wavelength. At this point, the number of degrees of freedom of the surface waves diminishes. The propagation at the surface of the metal being prohibited, the system evolves from 2D to 1D characteristics, as hinted by the shapes of the spectra before and after pc (see fig. 4): 2D as the resonance frequencies scale as √2 for p>pc and 1D as they scale as integers for p<pc.
We know discuss the reasons why all the models but FEM fit for the region of enhanced transmission but fail to reproduce the discontinuity. Existing models of EET agree well with our experimental results for small holes, but none shows any discontinuity. One approach is the question of the developments of the Maxwell’s equations solution in such theories under our experimental conditions, and then of the involved scalar or vectorial approximations. Most generally these models are derived from Green’s theorem which gives the field ψ(x) behind the aperture and from Kirchhoff’s approximation . Assuming small holes compared to the wavelength , one obtains
where ψ 0 is the field just before the hole, and R=x-x′. Even though this approximation gives mostly very good results in Dirichlet (or Neumann) boundary conditions, it is known  that it does not hold in mixed boundary conditions, where both ψ and ∇′ψ(x′) are fixed. Here, the plate thickness is negligible compared to the wavelength, resulting in a discontinuity of ψ in the hole and then in mixed boundary conditions. To obtain the exact boundary conditions, one then must apply Maxwell’s equations. Otherwise, subsequent calculation of the diffracted fields may contain errors. Similar demonstration can be done with vectorial theory. However, one must emphasize that Kirchhoff’s approximation gives good results for small holes, showing that the boundary conditions are essentially given by Dirichlet conditions. Near the transition, the wave propagation at the surface of the metal is strongly attenuated, and do not anymore only contribute to the boundary conditions. The fields inside the holes also play a role and the boundary conditions are then mixed. Therefore only full 3D resolution of Maxwell’s equations can provide a satisfactory result.
Two conclusions can be drawn from these results. Most models are able to fit part of the dynamics of the EET. For p<pc, none of the models used here is able to describe the evolution of the three parameters nor do they reproduce the discontinuity of the parameters. Finally only finite element programming is able to describe the discontinuities in the evolution of Δν, Q and ΔW.
Bethe theory [3, 4, 42] describes the transmission through one subwavelength hole. Considering no interaction between the holes of the array, the transmitted spectra through the arrays may be fitted for p<pc using the sum of one dimensional array of individual Bethe transmissions
with N the number of holes, r j the position of the center of hole j, r′j the coordinates inside hole j, and H o and E o the magnetic and the electric field, respectively. If we add an effective spatial period corresponding to the first resonant frequency experimentally detected, we obtain spectra fitting reasonably well the experimental ones (figure 6).
Finally, a surprising parallel can be drawn between the characteristics of the transmission through the subwavelength lattices, and a first order phase transition. Indeed, the frequency of the first resonance (or its width), as well as the quality factor exhibit their discontinuity at the same value of filling parameter. Assimilating the frequency shift Δν to the order parameter of the system consisting in the electromagnetic field coupled to the subwavelength lattice, and the quality factor Q to the heat capacity, one obtains that the system displays a discontinuity of both its order parameter and heat capacity at the same value pc. This is literally the definition of a first order phase transition [43, 44]. This parallel could prove to be helpful in conceiving a more general description of the complex interaction between electromagnetic waves and arrays of subwavelength holes, thanks to the important theoretical background of phase transition theories .
In this paper we showed that the characteristics of the transmission of an electromagnetic field through arrays of subwavelength holes strongly depend on the quantity of metal of the lattices. For small holes, all the characteristics of the EET can be found and several models may partially model the propagation of the electromagnetic waves. For larger holes the EET vanishes, and a new transmission takes place. This new transmission is characterized by symmetric resonances, with frequencies scaling as integers, and an electromagnetic transmission that exceeds the sum of each hole transmission. Finally, the transition between the two transmissions exhibits similar characteristics to first order phase transitions.
We thank Daniel R. Grischkowsky and R. Alan Cheville for the generous donation of the THz antenna used for this work.
References and links
1. T. W. Ebbesen, H. J. Lezec, H. F. Ghaemi, T. Thio, and P. A. Wolff, “Extraordinary optical transmission through subwavelength hole arrays,” Nature 391, 667–668 (1998). [CrossRef]
3. H. A. Bethe, “Theory of diffraction by small apertures,” Phys. Rev. 66, 163–182 (1944). [CrossRef]
4. C. J. Bouwkamp, “Diffraction Theory,” Rep. Prog. Phys. 17, 35–100 (1954). [CrossRef]
6. T. J. Kim, T. Thio, T. W. Ebbesen, D. E. Grupp, and H. J. Lezec, “Control of optical transmission through metals perforated with subwavelength hole arrays,” Opt. Lett. 24, 256–258 (1999). [CrossRef]
7. L. Martin-Moreno, F. J. Garcia-Vidal, H. J. Lezec, K. M. Pellerin, T. Thio, J. B. Pendry, and T. W. Ebbesen, “Theory of extraordinary optical transmission through subwavelength hole arrays,” Phys. Rev. Lett. 86, 1114–1117 (2001). [CrossRef] [PubMed]
8. S. C. Hohng, Y. C. Yoon, D. S. Kim, V. Malyarchuk, R. Müller, Ch. Lienau, J. W. Park, K. H. Yoo, J. Kim, H. Y. Ryu, and Q. H. Park, “Light emission from the shadows: surface plasmon nano-optics at near and far fields,” Appl. Phys. Lett. 81, 3239–3241 (2002). [CrossRef]
9. R. Gordon, A. G. Brolo, A. McKinnon, A. Rajora, B. Leathem, and K. L. Kavanagh, “Strong polarization in the optical transmission through elliptical nanohole arrays,” Phys. Rev. Lett. 92, 037401 (2004). [CrossRef] [PubMed]
10. A. Degiron, H. J. Lezec, W. L. Barnes, and T.W. Ebbesen, “Effects of hole depth on enhanced light transmission through subwavelength hole arrays,” Appl. Phys. Lett. 81, 4327–4330 (2002). [CrossRef]
12. E. Devaux, T. W. Ebbesen, J.-C. Weeber, and A. Dereux, “Launching and decoupling surface plasmons via micro-gratings,” Appl. Phys. Lett. 83, 4936–4939 (2003). [CrossRef]
14. Y.-H. Ye and J.-Y. Zhang, “Middle-infrared transmission enhancement through periodically perforated metal films,” Appl. Phys. Lett. 84, 2977–2980 (2004). [CrossRef]
15. J. Gomez Rivas, C. Schotsch, P. Haring Bolivar, and H. Kurz, “Enhanced transmission of THz radiation through subwavelength apertures,” Phys. Rev. B 68, 201306(R) (2003).
17. F. Miyamaru and M. Hangyo, “Finite size effect of transmission property for metal hole arrays in subterahertz region,” Appl. Phys. Lett. 84, 2742–2745 (2004). [CrossRef]
18. J. Gomez Rivas, C. Janke, P. H. Bolivar, and H. Kurz, “Transmission of THz radiation through InSb gratings of subwavelength apertures,” Opt. Express 13, 847–859 (2005), http://www.opticsinfobase.org/abstract.cfm?URI=oe-13-3-847. [CrossRef]
19. H. Cao and A. Nahata, “Resonantly enhanced transmission of terahertz radiation through a periodic array of subwavelength apertures,” Opt. Express 12, 1004–1010 (2004), http://www.opticsinfobase.org/abstract.cfm?URI=oe-12-6-1004 [CrossRef] [PubMed]
20. J.-B. Masson and G. Gallot, “Coupling between surface plasmons in subwavelength hole arrays,” Phys. Rev. B 73, 121401(R) (2006).
22. T. Thio, H. F. Ghaemi, H. J. Lezec, P. A. Wolff, and T. W. Ebbesen, “Surface-plasmon-enhanced transmission through hole arrays in Cr films,” J. Opt. Soc. Am. B 16, 1743–1748 (1999). [CrossRef]
23. A. K. Azad, Y. Zhao, and W. Zhang, “Transmission properties of terahertz pulses through an ultrathin subwavelength silicon hole array,” Appl. Phys. Lett. 86, 141102 (2005). [CrossRef]
24. F. Przybilla, C. Genet, and T.W. Ebbesen, “Enhanced transmission through Penrose subwavelength hole arrays,” Appl. Phys. Lett. 89, 121115 (2006). [CrossRef]
26. D. Grischkowsky, S. R. Keiding, M. van Exter, and Ch. Fattinger, “Far-infrared time-domain spectroscopy with terahertz beams of dielectrics and semiconductors,” J. Opt. Soc. Am. B 7, 2006–2015 (1990).
27. C. C. Chen, “Transmission of microwave through perforated flat plates of finite thickness,” IEEE Trans. Microwave Theory Technol. 21, 1–6 (1973). [CrossRef]
28. J.A. Besley, N.N. Akhamediev, and P.D. Miller, “Periodic optical wavequides: exact Floquet theory and spectral properties,” Studies in Applied Mathematics 101, 343–355 (1998). [CrossRef]
29. U. Fano, “Effects of configuration interaction on intensities and phase shifts,” Phys. Rev. 124, 1866–1875 (1961). [CrossRef]
30. C. Genet, M.P. van Exter, and J.P. Woerdman, “Fano-Type interpretation of red shifts and red tails in hole array transmission spectra,” Opt. commun. 225, 331–336 (2003). [CrossRef]
31. J.-B. Masson, A. Podzorov, and G. Gallot, “Generalized parabolic Fano model of extraordinary electromagnetic transmission in subwavelength hole arrays,” submitted.
32. K. J. Klein Koerkamp, S. Enoch, F. B. Segerink, N. F. van Hulst, and L. Kuipers, “Strong influence of hole shape on extraordinary transmission through periodic arrays of subwavelength apertures,” Phys. Rev. Lett. 92, 183901 (2004). [CrossRef] [PubMed]
33. L. F. Li, “New formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. Soc. Am. A 14, 2758–2767 (1997). [CrossRef]
34. C.M. Soukoulis, “Photonic Crystals and Light Localization in the 21st Century,” NATO Science Series Vol. 563 (Kluwer, Dordrecht, 2001).
35. J.-Y. Laluet, E. Devaux, C. Genet, T. W. Ebbesen, J.-C. Weeber, and A. Dereux, “Optimization of surface plasmons launching from subwavelength hole arrays: modelling and experiments,” Opt. Express 15, 3488–3495 (2007), http://www.opticsinfobase.org/abstract.cfm?URI=oe-15-6-3488. [CrossRef] [PubMed]
36. E. Devaux, T.W. Ebbesen, J.-C Weeber, and A. Dereux, “Launching and decoupling surface plasmons via micro-gratings,” Appl. Phys. Lett. 83, 4936–4939 (2003). [CrossRef]
38. Comsol Multiphysics, Comsol Inc., Burlington, MA.
39. M. A. Ordal, L. L. Long, R. J. Bell, S. E. Bell, R. R. Bell, J. R.W. Alexander, and C. A Ward, “Optical-properties of the metals Al, Co, Cu, Au, Fe, Pb, Ni, Pd, Pt, Ag, Ti, and W in the infrared and far infrared,” Appl. Opt. 22, 1099–1119 (1983). [CrossRef] [PubMed]
40. Magnetic Properties of Metals, Landolt-Bornstein, Group III: condensed matter, Springer-Verlag, Berlin (1986).
41. J. D. Jackson, Classical Electrodynamics, 3rd edition, John Wiley & Sons (1999).
42. M. Born and E. Wolf, Principle of Optics, 7th edition, Cambridge University Press (2001).
43. F. Schawbl, Statistical Mechanics, Springer (2000).
44. H.E. Stanley, Introduction to phase transitions and critical phenomena, Oxford science publications (1971).
45. J.-B. Masson and G. Gallot, “Experimental evidence of percolation phase transition in surface plasmons generation,” ArXiv:cond-mat/0611280v1 (2006).