Complex Bloch-modes calculation of plasmonic crystal slabs by means of finite elements method

Giuseppe Parisi; Pierfrancesco Zilio; Filippo Romanato

doi:10.1364/OE.20.016690

1. Introduction

Periodic metallic nanostructures have received great attention in recent years due to their surprising optical properties. The interaction of electromagnetic fields with the free electron gas of metals gives rise to particular electromagnetic field waves which are known as Surface Plasmon Polaritons (SPPs) [1]. Due to their properties of extreme light confinement and sensitivity to the dielectric environment, SPPs have found lots of interesting application in the fields of sensing [2], thin film photovoltaics [3], near field imaging [4], subwavelength waveguides [5] etc.

In particular, periodically nanostructured metal-dielectric interfaces have been subject of many researches [1,4] due to their unexpected optical properties such as Extraordinary Optical Transmission [6–10] and negative refraction [11]. Surface Plasmons propagation across such structures exhibits analogous properties of light propagation in dielectric photonic crystals [4]. A complicated system of band gaps has indeed been calculated and measured using different metallic gratings [4]. Shallow surface features can be assigned an effective refractive index describing their interaction with SPPs, thus enabling a complete analogy with 2D photonic crystals. The description of plasmonic Bloch waves propagating in a plasmonic crystal slab, however, is a much more challenging task with respect to the description of photonic Bloch waves in fully periodic dielectric systems.

First of all, in photonic crystals analyses, typically transparent non dispersive materials are considered [12]. In these cases the Helmholtz eigenvalue equation reduces to a generalized linear eigenvalue equation, which is readily numerically treated with standard linear solvers. The band structure of the photonic crystal, ω(k), is then obtained as a function of the crystal momentum k. However, when material dispersion plays a crucial role, as in the case of plasmonic crystals, the resulting eigenvalue equation is nonlinear [13]. Besides time consuming methods involving non-linear iterative solvers [14], other methods have been proposed to solve such kind of problems, which are based on the solution of Helmholtz’s equation considering as eigenvalue the wave vector instead of the frequency (obtaining therefore the bands as k(ω)) [15–19]. Generally these methods end up with a much more tractable quadratic eigenvalue equation in k, which can be solved by proper linearization procedures [17]. Moreover the imaginary part of k as well as the real one is naturally retrieved. In a recent work Fietz et al. [18] showed a modal analysis method based on reformulation of the eigenvalue equation in weak form, which enables to obtain a quadratic eigenvalue equation in k. The weak formulation finds a natural solution in the frame of the Finite Elements Method, which inherently handles weak forms of partial differential equations. This modal method, originally developed by Hiett et al. [16], demonstrated to be very powerful since it allows dealing with completely general materials, dispersive, lossy and possibly anisotropic.

The analysis of Photonic Crystal Slabs (PCS), however, requires proper handling not only of truly bound modes but also of leaky modes. This term addresses eigenmodes whose crystal momentum lies inside the light cone. In these cases the modes can couple to waves propagating in the semi-infinite half spaces surrounding the slab, resulting in radiative losses [12,13]. Simulation of open space boundaries is a tricky task both for methods based on discretization of real space (like Finite Elements and Finite Differences) as well as for methods which discretize the wave vectors space (like Plane Wave Expansion method). A common approach introduces a fictitious periodicity in the direction normal to the slab (Super-cell approach [13]). The period is chosen sufficiently large in order to decouple the bound modes of the slabs. However, the artificial periodicity produces many spurious modes and moreover it perturbs the physical leaky eigenmodes profile. Plane Wave Expansion Method has been extended by Shi et al. [15] including Perfectly Matched Layers (PMLs) which correctly absorb leakage radiation. These layers are characterized by slowly space varying relative permittivity and permeability, properly designed in order to minimize reflection of plane waves impinging on them at arbitrary incidence. The method can deal with photonic crystal slabs with arbitrarily dispersive materials but only in case of absence of losses, excluding therefore the important class of plasmonic crystal slabs. Other methods based on FDTD can handle the lossy dispersive problem but are computationally very expensive [15]. To our knowledge, a complete and efficient numerical treatment of metallic planar crystals still lacks.

In this paper we calculate the Bloch band structure of dispersive lossy photonic crystal slabs adopting the weak formulation of the Helmholtz’s equation and discretizing the system by means of Finite Elements Method. In addition to the analysis presented in [18], we include PMLs within the unit cell domain in order to properly deal with leaky modes radiation. We show how the technique is effective in analyzing plasmonic crystal slabs considering first a test example, a simple mono-periodic sinusoidal metal-dielectric interface. This structure turns out to be a useful test and allows to properly set PMLs parameters. Then we apply the method to an example of particular physical interest for its properties of Extraordinary Optical Transmission (EOT), i.e. a bi-periodic array of square holes in a metal film. Although it has been widely investigated in recent years [7], few works considered the relationship between the optical Bloch eigenmodes of this structure and its EOT properties [9]. This is mainly because retrieving the complete complex band structure in presence of strong leakage radiation still remains a challenging task, and it is the typical case of plasmonic crystals with lattice constant comparable with the impinging light wavelength and with deep modulation of the metal-dielectric interface. The modal analysis method we propose comes in useful in clarifying the relationship between optical response of the structure and periodicity-induced resonant modes.

The work is organized as follows. In section 2 and 3 we will recall the weak form FEM formulation of the eigenvalue problem, pointing out also how PMLs can be implemented in this solution scheme. In section 4 we will apply the technique to the determination of the complex Bloch-bands of a simple two-dimensional metallic structure, i.e. a sinusoidal metal dielectric interface. Finally in section 5 the method will be used to calculate bound and leaky modes of a three-dimensional bi-periodic plasmonic crystal slab. In particular it will be pointed out how the method succeeds in reproducing the leaky modes field patterns. In both sections, 4 and 5, the modal results are compared to scattering FEM simulations, which, in turn, serve as test tool for the accuracy and reliability of the method.

2. The FEM simulation

The aim of the present section is to briefly revisit the FEM weak formulation [17,18] of the following Helmholtz’s eigenvalue equation

\nabla \times (\hat{p} \nabla \times V) - ω^{2} \hat{q} V = 0,

where

V = E, \hat{p} = 1 / \hat{μ} (r, ω), \hat{q} = \hat{ε} (r, ω)

for the electric field formulation and

V = H, \hat{p} = 1 / \hat{ε} (r, ω), \hat{q} = \hat{μ} (r, ω)

for the magnetic field formulation. Here

\hat{ε} = \hat{ε} (r, ω)

and

\hat{μ} = \hat{μ} (r, ω)

are respectively the frequency-dependent permittivity and permeability tensors of the plasmonic crystal. They are supposed to be periodic in x, y and z direction. The problem reduces to find Bloch-waves solutions of the form [12,13]

V (r; k) = e^{- i k \cdot r} u (r) = e^{- i (k_{x} x + k_{y} y + k_{z} z)} u (r),

where k is the Bloch-vector. The Finite Elements method relies on the so-called weak formulation of Eq. (1), which consists in annulling the following residues

R_{k} (v,u) = \int_{Ω} v \cdot [\nabla \times (\hat{p} \nabla \times V (r; k)) - ω^{2} \hat{q} V (r; k)] d^{3} r,

where v is a weight function and the integration is over the unitary cell volume [20]. After inserting expression (2) in Eq. (3) and some straightforward algebra, Eq. (3) yields to

\begin{array}{l} \int_{Ω} d^{3} r (- v \cdot [k \times \hat{p} (k \times u)] - i v \cdot [k \times \hat{p} (\nabla \times u)] - i (\nabla \times v) \cdot \hat{p} (k \times u)) + \\ + \int_{Ω} d^{3} r ((\nabla \times v) \cdot \hat{p} (\nabla \times u) - \frac{ω^{2}}{c^{2}} v \cdot \hat{q} u) + \int_{\partial Ω} d A v \cdot [\hat{n} \times \hat{p} (- i k \times u+ \nabla \times u)] = 0, \end{array}

In order to turn Eq. (4) into an eigenvalue problem, the three degrees of freedom that comprise the Bloch wave-vector k must be reduced to one; i.e by fixing a particular k-vector direction, $\hat{λ}$ . In this case Eq. (4) turns out to be a quadratic eigenvalue equation in the k-vector projection along the chosen direction, i.e. $λ = k \cdot \hat{λ}$ . Following the usual FEM discretization procedure [20], Eq. (4) is then turned into a matrix equation [17,18,20] and opportunely linearized.

The presented weak formulation of the eigenvalue problem is inherently handled by several Finite Elements software packages. All the models and examples in this paper have been performed and calculated with COMSOL Multiphysics, which allows the user to specify a custom field equation to be solved. In particular the software automatically turns the input Eq. (4) into an algebraic linear equation system. For more details about FEM we remand to Ref [20,21].

3. The perfect matched layer (PML) implementation

A photonic crystal slab (PCS) is a photonic structure with two-dimensional periodicity and finite thickness along the third spatial dimension [12], as depicted in Fig. 1 (green structure). As mentioned before, the modal analysis of a PCS requires proper truncation of the computational domain in the direction normal to the slab. Such a truncation may be performed by introducing the so-called Perfectly Matched Layers at the upper and lower boundaries of the unit cell (the violet domains in Fig. 1). Their function is to absorb radiation coming from the slab, at any frequency and angle of incidence.

Fig. 1 Scheme of photonic crystal slab (green) with PMLs (violet) truncating the cladding domains.

Download Full Size | PDF

As proposed by Sacks et al. [22,23], PMLs can be treated as uniaxial anisotropic absorbers whose permittivity and permeability tensors are specified according to the following relations:

\hat{ε} = ε \hat{Λ}, \hat{μ} = μ \hat{Λ},

where, in case of absorption in the z-direction,

\hat{Λ} = (\begin{matrix} c & 0 & 0 \\ 0 & c & 0 \\ 0 & 0 & c^{- 1} \end{matrix}),

being c = α - iβ, with α equal to the real part of the adjacent medium relative permittivity and β > 0. This condition, without any further restriction on β, can assure perfect absorption of any plane wave incident upon the boundary of an infinitely thick PML, independently of frequency or incidence angle [22]. As was noticed [23], however, the discrete approximation of fields and material parameters results in a spurious impedance loading at the interface between PMLs and physical domains, and significant reflections are found in presence of constant β. In other words, after discretization, PMLs are still absorbing materials and waves that propagate across them are still attenuated, but the boundary between PML and regular medium is no longer reflectionless. This mismatch problem in the discretized space can be tempered by using spatially varying material parameters [23]:

β = \bar{σ} \frac{{| (z \mp | z_{0} |) |}^{n}}{L^{n}} .

The plus (minus) sign refers to PMLs located in the lower (upper) half space. In Eq. (7)

\bar{σ} = σ / ω

, and

σ

, L and n are constant parameters. L is the PML thickness while

\bar{σ}

is the wave attenuation rate within the PML. Condition (7) assures a smooth increasing of the damping rate of waves incoming into the absorbers and significantly improves PML absorption. Optimal PML performance, however, requires a careful optimization of parameters σ, L and n. The problem will be addressed in next section. We notice that, in general, the β parameter depends on ω and therefore the PML is a dispersive material. This does not represent a problem since, as mentioned before, in the present modal method frequency is fixed.

We notice that PMLs are not the only option in order to properly truncate open boundaries. Absorbing Boundary Conditions (ABC) [20] can also be used to absorb plane wave leaky radiation. A detail comparison between effectiveness of PMLs and ABCs in the present method is out of the scope of the present work. However, we notice that if the leakage radiation is not in the form of a single plane wave, ABCs are expected to be not suitable.

4. 2D sinusoidal grating and PML test

We first start with a 2D example (Fig. 2(a) ) which serves as a guideline to set the PMLs parameters, optimizing them in order to minimize reflections. The unit cell consists of a sinusoidal metal-dielectric interface in the x-z plane, infinitely extended in y-direction (out of plane). In the z-direction the unitary cell is limited above by a PML layer and below by the bulk metal. The period (d) in the x-direction is set to 600nm while the peak-valley amplitude of the sinusoid (a) is increased from 0 to 150nm. We adopt silver as metal, its permittivity being taken from Palik [24]. Periodic boundary conditions are set in both x and z directions. In z direction this condition is not mandatory because of the presence of the PML. This choice is made in order to reduce to zero the last term in Eq. (4). Since the PMLs are expected to absorb all radiation impinging on them, the particular boundary condition set at the end of the PML layer actually does not influence the resulting field distribution.

Fig. 2 (a) Scheme of unit cell for the sinusoidal plasmonic grating; (b) Variation of a sample eigenvalue (λ_{L = d} = 0.0905356 at ω = 1.7∙10¹⁵Hz) as a function of PML thickness L and σ, ∆λ is defined as ∆λ_m = λ_{L = (m + 1)d} - λ_{L = md} (m>0), scale is in decibels.

Download Full Size | PDF

In the case of mono-periodic gratings, as in the considered example, one is typically interested in studying surface plasmon modes whose k-vector is parallel to the grating vector, $G = (2 π / d) \hat{x} .$ In this case the weak form of Helmholtz’s equation, Eq. (4), reduces to the following simplified form

\int_{Ω} d^{2} r [λ^{2} \frac{u v}{ε_{33}} + λ \frac{i}{ε_{33}} (\frac{\partial u}{\partial x} v - u \frac{\partial v}{\partial x}) + \frac{\partial u}{\partial z} \frac{\partial v}{\partial z} \frac{1}{ε_{11}} + \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} \frac{1}{ε_{33}} - {(\frac{ω}{c})}^{2} u v μ_{22}] = 0.

where λ = k_x, ε_ii and μ_ii are the diagonal components of the relative permittivity and permeability tensors of the medium involved (both physical domains and PMLs).

The need of terminating the computational domain with proper boundary condition, in any case, results in the generation of undesired spurious modes, even in presence of PMLs. Most of these nonphysical modes look like guided modes within the PML domains or between slab and absorbing boundary. A common way to discriminate physical leaky or bound modes from PML modes is to employ an average field intensity based filter [15]. The physical modes, even leaky modes, are expected to be those ones mostly confined in proximity of the plasmonic slab and rapidly decaying within the PML domain. We calculated the following quantity for all eigenmodes

κ = \frac{{〈 {| H |}^{2} 〉}_{(1)}}{{〈 {| H |}^{2} 〉}_{(2)}}

where <|H|²>₍₁₎, <|H|²>₍₂₎ are the averages of the squared H field norm in a region far from the grating in which all evanescent fields are vanished and within a thin air layer (tens of nm) close to the metal surface respectively, i.e. regions marked with (1) and (2) in Fig. 2(a) respectively. We found out that κ values ranging from 0 to 20% are typically obtained in case of physical bound or leaky eigenmodes, whereas higher values are found for nonphysical PML–related modes. We used therefore the quantity κ to filter out non-physical modes.

Before extracting the dispersion curves we performed a PML test in order to appropriately set the constants σ, L and n in the PML model, Eq. (7). In Fig. 2(b) we consider a sample eigenvalue λ and report its variation as a function of PML thickness L, i.e. ∆λ_m = λ_{L = (m + 1)d} - λ_{L = md} (m>0), for several σ values, keeping fixed n, ω and mesh element density. As can be seen, PML performances improve with increasing PML thickness and decreasing σ (∆λ converges to a small constant value related to the numerical precision). This can be explained as follows. The variation in $λ$ is related to non-zero reflections at the interface between PML and air. As previously mentioned, PML is perfectly reflectionless only when solving the exact wave equation. Reflections keep low as long as the discretization is a good approximation of the exponentially decaying wave within the PML, provided the PML is thick enough to completely absorb the incoming wave. In order to increase the accuracy of the approximation, two methods are possible [25,26]: (a) for a fixed PML thickness, increasing mesh element density in the PML or, (b), at a fixed mesh element density, decreasing σ parameter, in order to turn on wave absorption more gradually, and increasing PML thickness. Clearly option (b) works since, as long as the PML material is slowly varying, the wave decay is more diluted in space and is better resolved by the given mesh density. Of course both options (a) and (b) require increasing the number of elements within the PML domain. An acceptable compromise between accuracy and computational cost was found taking the PML thickness L equal to 4d and the σ value close to ω, keeping the mesh density the same as in the air domain.

With regard to the n parameter, it was shown for example in [23,27], that a simple quadratic or cubic turn-on of the PML absorption usually produces negligible reflections for a PML only few wavelengths thick.

Once set the tensor parameters in Eq. (7) we performed the modal analysis of the structure in Fig. 2(a). In Fig. 3(a) we report a reflectance map obtained by a scattering FEM simulation of a grating with amplitude 75nm. This kind of maps is commonly used in literature to deduce the real part of SPP modes dispersions by looking at the reflectance dips [1]. Black lines are the real part of the SPP modes calculated with the modal analysis. As can be seen, the modal analysis correctly predicts the dips observed in the reflectance maps, which correspond to coupling of impinging light to SPP modes.

Fig. 3 (a) Reflectance map compared with the calculated dispersion curves (black lines), green dashed line is the light line; (b) Comparison between bands obtained with (black) and without (red) PML domains (in both cases the tolerance parameter κ was set to 0.1, and periodicity in z direction is 5μm).

Download Full Size | PDF

For a better understanding of the numerical improvement due to the PML insertion into the unitary cell, we report in Fig. 3(b) a comparison between the modal curves calculated with and without the PML addiction. In both cases, a fictitious periodicity of 5μm in z-direction is set. In absence of PML, the modal curves present a discontinuous behavior. The fictitious periodicity in z-direction unavoidably restricts the allowed spectrum of possible frequencies for a given momentum Re(k_x) to a discrete set. This discretization occurs even in presence of PMLs, but is much more refined, tending to a continuum of states in the limit of infinite PML mesh resolution (which mimics the perfect open boundary condition).

In Fig. 4 we report the real and imaginary parts of the mode dispersions for increasing grating amplitudes. As well known [1], a flat metal-dielectric interface sustains a SPP wave whose dispersion is given by $k_{S P P, F} = k_{0} {(ε_{m} ε_{d} / (ε_{m} + ε_{d}))}^{1 / 2}$ , where ε_m,d are respectively the metal and dielectric relative permittivities, while k₀ is the vacuum wave vector. When a shallow periodic perturbation is introduced the mode dispersion remains close to the flat case. The only remarkable difference is observed at the crossing point of the flat SPP dispersions at the edge of first Brillouin zone (k_x = π/d), where an energy gap appears with increasing grating amplitude [28–30]. As can be seen the modal analysis correctly reproduces the effect. Actually it can be noted that the method converges to a solution also at frequencies within the band gap and we should more correctly term this frequency range as a gaplike region [31]. In particular the band, at frequencies close to the lower gap edge, rises steeply through a continuum of states toward the upper gap edge. The gap edges are actually not well defined, especially for higher amplitude gratings. This behavior was already noticed in similar modal analyses performed with real frequency and complex Bloch-wave-vector [32,33]. The frequency gap however is more evident looking at the imaginary part of the mode dispersions (Fig. 4(b)). At frequencies corresponding to the gap the imaginary part assumes great values, reaching a maximum at the center of the gap. The phenomenon becomes more evident with increasing grating depth, as expected. Within the gap the modes have a great dissipation which prevents them to propagate in the direction orthogonal to the grooves.

Fig. 4 Real part (a) and imaginary part (b) of the modes varying the sinusoidal grating amplitude a: a = 0 (flat case, black), 30nm (cyan), 50nm (magenta), 75nm (green), 100nm (red), 150nm (blue). The dashed black curve in (a) is the light line.

Download Full Size | PDF

For frequencies above the upper gap edge (ω > 1.7∙10¹⁵Hz) the imaginary part of the Bloch-wave-vector has values much greater than its corresponding values at frequencies below the gap. This is because the folded plasmonic band enters the light cone and the radiative coupling to propagating waves in the upper air half space becomes possible.

In order to have an independent check of the reliability of imaginary part of mode eigenvalues, we considered the following FEM model (see Fig. 5(a) ). An SPP is launched on a flat metal-air interface at the left boundary of the model. The metal interface then continues with a sinusoidal profile. The SPP wave, reaching the sinusoidal grooves, is partially reflected and partially couples to the grating Bloch-mode. Once the SPP Bloch mode is correctly excited, it propagates along the grating with its own complex propagation constant. In particular, from the field profile in Fig. 5(a) we can extract the attenuation constant, α, by means of exponential fit, obtaining therefore the imaginary part of the propagation constant as Im(λ) = 1/(2α) [1].

Fig. 5 (a) z-component of the magnetic field obtained by exciting the proper mode by mean of excitation boundary conditions. Frequency is ω = 3.5∙10¹⁵Hz, d = 600nm and a = 50nm. (b) Comparison between the imaginary part of the modes calculated with illumination (blue curve) and by modal analysis (magenta curve) respectively at fixed amplitude a = 50nm.

Download Full Size | PDF

We performed the calculation for different frequencies at fixed period d = 600nm and amplitude a = 50nm. Results are reported in Fig. 5(b). We see that there is a good agreement between the direct modal calculation and the indirect method based on SPP excitation for a wide spectral range of our interest (up to ω ≈4∙10¹⁵Hz). For frequencies above ω ≈4∙10¹⁵Hz a growing mismatch is observed since the considered PML setup is not optimal any more for too short wavelengths. The deviation may be reduced by tuning again the PML parameters.

5. 2-D periodic array

In this section we apply the modal analysis to the plasmonic crystal slab depicted in Fig. 1, consisting of an array of squared nano-holes with sizes a_x = a_y = a milled in a thin silver slab. Period and silver thickness are fixed to d = 940nm and h = 200nm respectively. Two PMLs are introduced in the unit cell in order to absorb the leakage radiation propagating toward open space. They are set at distance z₀ = ± 1470nm from the slab and their parameters are set according to previous section.

Holey thin metal films exhibit very interesting plasmonic properties of extraordinary optical transmission which have been extensively studied in a vast literature [6–11]. However, what is usually omitted in literature is the detailed visualization of both the real and imaginary dispersions of the SPP Bloch modes of the considered structures. This is due the strong leakage radiation damping that affects the modes for large holes sizes, which is hardly handled by standard numerical techniques.

In order to verify the effectiveness of the method in presence of both low and high radiative losses, we carried out two simulations, the first for a small hole size, a = 250nm, the second for a wider hole, a = 500nm. We focus on modes along x-direction, $k = λ {\hat{k}}_{x}$ . Figure 6 reports the transmittance maps obtained by FEM scattering simulations for both hole sizes and for TM and TE polarized impinging plane waves. The superimposed black lines are the real parts of the SPP Bloch modes of the structure calculated with the modal analysis. As can be seen, Bloch-modes dispersions directly correlate to transmittance features of the structure. This strict connection has been noticed in several works [6–9]. The calculated bands approximately follow the dispersions given by the well known grating coupling relation

ω {(k_{x}, k_{y})}_{m, n} = \frac{1}{n_{e f f}} \sqrt{{(Re (k_{x}) + \frac{2 π m}{d})}^{2} + {(Re (k_{y}) + \frac{2 π n}{d})}^{2}},

where n_eff is the real part of effective index of the SPP mode propagating on a flat silver-air interface [1] and m, n are integers. We consider dispersions along the x-direction by setting Re(k_y) = 0. The white and red lines in Fig. 6(a,c,e) denote the (m,n) = (

\pm

1,0) and the (0,

\pm

1) flat SPP dispersion curves respectively, according to Eq. (10). Figure 6(e), 6(f) report details of Fig. 6(a). As is seen the dispersions found with the modal analysis more precisely account for the transmittance features observed.

Fig. 6 Maps of TM and TE transmittance through arrays of square nano-holes (period d = 940nm) in a 200nm thick Silver film in air compared with the calculated Bloch-modes dispersion curves (black lines, labeled as (|m|,|n|) ^± ); (a) and (b) refer to holes with size a = 250 nm, (c) and (d) to a = 500nm; (e) and (f) are zooms of the regions marked with 1 and 2 respectively in (a). Black dashed line marks the light line, white and red solid lines mark the flat SPP dispersions, ( ± 1,0) and (0, ± 1) respectively. Figure 6(f) contains also a comparison between the x-component of the magnetic field profile in the z-y plane obtained with modal analysis and illumination respectively.

Download Full Size | PDF

The modes in general split into two categories: antisymmetric (lower frequency, labeled as (|m|,|n|)^-) and symmetric modes (higher frequency, labeled as (|m|,|n|)⁺), depending on the distribution of the dominant magnetic field component with respect to the z = 0. This is a consequence of the mirror symmetry of the system with respect to the z = 0 plane [12]. The comparison between modal analysis and direct illumination gives important information about the modes which can be excited with the two types of illumination. We see that the (1,0) ^± modes can be excited only in presence of TM illumination and both correspond to transmittance peaks (see Fig. 6(e)). On the other hand, different (0,1) ^± modes can be excited both in TM and TE illumination. In particular, in Fig. 6(f), only one transmittance peak is clearly visible corresponding to the antisymmetric TM (0,1)^- mode (see the comparison between H_x in the y-z plane obtained with modal analysis and illumination at frequency of 2.022∙10¹⁵Hz).

In case of large holes, we note that the (1,0)^- mode strongly deviates from the flat SPP dispersion curves, whereas the (1,0)⁺ is almost unperturbed. It is also evident how the (1,0)^- one is correlated to a much higher transmittance peak than the (1,0)⁺ (Fig. 6(c)).

In Fig. 7(a) and 7(b) we report respectively the fields H_y and |E_x| in the x-z plane, for the two TM modes observed at frequency ω = 1.6∙10¹⁵Hz with a = 250nm (I-II) and a = 500nm (III-IV). As can be expected, for small hole sizes, the modes are well confined close to the metal slab. The mode confinement decreases in case of a = 500nm and a stronger leakage radiation is observed in the cladding regions. The presence of the PMLs, however, makes it possible to properly absorb this radiation and allows reconstructing the correct mode field profiles. This is clearly seen in Fig. 7(b). In fact, the almost uniform |E_x| in the region far from the slab indicates that no waves reflected from the PMLs are present. Comparing Fig. 7(b) I and III with 7(b) II and IV, we observe that the electric field intensity is mainly concentrated within the hole in the anti-symmetric modes, while the field has a node at the plane z = 0 in the symmetric case.

Fig. 7 (a) H_y and (b) |E_x| fields of the TM modes found at frequency of 1.6∙10¹⁵Hz for a = 250nm (I,II) and a = 500nm (III,IV). The calculated eigenvalues in the four cases are respectively k_x/(π/d) = 0.3828 + 0.0019i, 0.3878 + 0.0011i, 0.2792 + 0.1251i, 0.3829 + 0.0025i. Modes are classified as antisymmetric (I,III) and symmetric (II,IV) according to the symmetry of H_y field with respect to the z = 0 plane. PML size is fixed at 2d for convenience.

Download Full Size | PDF

As just pointed out in [7–11], peaks in optical transmission through a holey metal film can be described in terms of symmetric and antisymmetric coupling of plasmonic modes between the two horizontal metal-dielectric interfaces of the slab. In case of symmetric coupling, SPPs at the two horizontal metal-dielectric interfaces are weakly coupled via evanescent fields inside the hole. The mechanism depends mainly on the periodicity of the structure and is the main EOT channel for arrays of small holes. In the antisymmetric coupling, instead, single-interface SPPs are strongly coupled via a Fabry-Pérot resonance inside the hole. This resonance depends more on single hole characteristics (size and depth) rather than on the periodicity of the structure, and turns out to be the dominant EOT mechanism for arrays of large holes [10]. We refer to [7–11] for a detailed description of these EOT mechanisms.

This simplified picture for EOT is confirmed by our modal analysis, provided we replace the horizontal-SPP-coupling description with a Bloch-modes description. As a matter of fact, we find a good matching between SPP Bloch-modes of the structure and EOT peaks for the 250nm-holes array (Fig. 6(a), 6(b)). This is in agreement with the fact that, for small holes, the nature of the EOT is closely related to the periodicity of the structure. In case of large holes, on the other hand, it can be noticed that the TM bands do not exactly follow the wide transmittance maxima observed. Moreover, in presence of TE illumination, no Bloch modes are correlated to the transmittance maximum found around ω = 1.62∙10¹⁵Hz. This suggest that the main EOT peaks observed are not strictly correlated to the periodicity of the structure but are mostly related to single-hole Fabry-Pérot resonances [9,10] and therefore they are reasonably not completely caught by an in-plane modal analysis. By contrast, we note that the symmetric modes, both in TE and TM illumination, are only weakly perturbed by the increased holes size and exactly match the corresponding transmittance peaks visible in the maps. This persistent matching confirms that these transmittance features are strictly related to the excitation of a periodicity-induced Bloch modes of the structure.

Finally we report in Fig. 8 the real and imaginary parts of the complete band structure within the first Brillouin zone along the Г-X crystal direction for both cases at a = 250 nm and a = 500nm. The colored bands refer to the TM symmetric (red) and anti-symmetric (blue) modes. The most striking effect of increasing the hole size is found looking at the imaginary parts of the modes. In particular, as can be expected, the antisymmetric mode, being related to the vertical Fabry-Pérot resonance of the structure, is the most sensitive to the hole size and shows a huge increase of imaginary part. At frequencies close to 1∙10¹⁵ Hz the small frequency gap found for a = 250 nm is much wider at a = 500 nm. The hole size variation corresponds to a variation of the metal filling fraction along the mode propagation direction and acts on the gap-size as the metal amplitude variation did in the sinusoidal grating example discussed in section 4. On the other hand, the imaginary part of the symmetric mode is almost the same for a = 250 nm and 500 nm.

Fig. 8 Real and imaginary parts of the modes for a = 250nm (a,b) and a = 500nm (c,d). Blue and red dots represent the TM (1,0)^- and TM (1,0)⁺ modes respectively, black dots represent the TE and TM (0,1) ^± modes. Insets in (a) report the dominant magnetic field component profile (colorscale) and the (H) field (arrows) respectively in a x-y plane laying 10 nm above the slab at frequency of 2.05x10¹⁵Hz for the (|m|,|n|)^- modes.

Download Full Size | PDF

Black dots in Fig. 8(a) mark the TE and TM (0,1) ^± modes. Their imaginary parts present strong increases in correspondence of the edge of the first Brillouin zone around the frequency of 2.2∙10¹⁵ Hz (Fig. 8(b) inset). This behavior indicates the presence of gaps similar to the ones observed for the TM modes. In the insets of Fig. 8(a) we report the dominant magnetic field component profiles of the (|m|,|n|)^- modes for a = 250 nm in a x-y plane laying 10 nm above the slab at frequency of 2.05∙10¹⁵ Hz. We see that the TM (1,0)^- mode is a transversal mode with the magnetic field polarized in the y-direction, while both the TE and TM (0,1)^- modes are longitudinal modes with the magnetic field mostly polarized in the x-direction of propagation.

Looking at the real dispersions of the modes it can be noticed a band bending at frequencies around ω = 2∙10¹⁵ Hz near k_x = 0. Correspondingly, a divergence in the imaginary parts is observed. Similar deviations are typical of real-frequency (and complex propagation constant) eigenvalue methods. They are not physical and were already reported elsewhere [18,31,33].

6. Conclusions

In conclusion we presented a full vectorial Finite Elements based numerical method for the modal analysis of photonic crystal slabs in presence of dispersive lossy materials. In particular, the important class of plasmonic crystal slabs of arbitrarily complex geometries can be handled. The method relies on the reformulation of Helmholtz eigenvalue equation in weak form, including Perfectly Matched Layers in the unit cell in order to simulate open boundaries. Results were firstly obtained in the simple case of 2D sinusoidal metal-dielectric interface and were compared to scattering FEM simulation. Good agreement was found, provided that the PML parameters are properly set. Our results prove that PML implementation allows to effectively study leaky modes, characteristic features of photonic crystal slabs, thus enabling the reconstruction of the correct radiative eigenmode profile.

The method was then applied to the more complex case of periodic arrays of holes in a silver metal film, enabling to investigate the role that plasmonic Bloch modes play in the Extraordinary Optical Transmission phenomenon presented by the structure. By comparing Bloch-modes dispersions to transmittance maps, indeed, it was possible to discriminate between optical features mainly related to periodicity from those mainly dependent upon single hole characteristics (with a weak array interaction). A study as a function of the hole size, moreover, revealed that in the case of small holes the EOT phenomenon is mainly correlated to Bloch modes, while for large hole arrays it is mainly influenced by single hole resonances.

Acknowledgments

This work has been supported by SPLENDID project (Surface Plasmonics for Enhanced Nano-Detectors and Innovative Devices, CARIPARO foundation) and PLATFORMS project of Padova University.

References and links

1. H. Rather, Surface Plasmons (Springer-Verlag, 1988).

2. J. N. Anker, W. P. Hall, O. Lyandres, N. C. Shah, J. Zhao, and R. P. Van Duyne, “Biosensing with plasmonic nanosensors,” Nat. Mater. 7(6), 442–453 (2008). [CrossRef] [PubMed]

3. H. A. Atwater and A. Polman, “Plasmonics for improved photovoltaic devices,” Nat. Mater. 9(3), 205–213 (2010). [CrossRef] [PubMed]

4. A. V. Zayats, I. I. Smolyaninov, and A. A. Maradudin, “Nano-optics of surface plasmon polaritons,” Phys. Rep. 408(3-4), 131–314 (2005). [CrossRef]

5. M. L. Brongersma and P. G. Kik, Surface Plasmon Nanophotonics, (Springer, 2007).

6. F. J. Garcia de Abajo, “Colloquium: Light scattering by particle and hole arrays,” Rev. Mod. Phys. 79(4), 1267–1290 (2007). [CrossRef]

7. F. J. Garcia-Vidal, L. Martin-Moreno, T. W. Ebbesen, and L. Kuipers, “Light passing through subwavelength apertures,” Rev. Mod. Phys. 82(1), 729–787 (2010). [CrossRef]

8. H. Liu and P. Lalanne, “Microscopic theory of the extraordinary optical transmission,” Nature 452(7188), 728–731 (2008). [CrossRef] [PubMed]

9. P. Lalanne, J. C. Rodier, and J. P. Hugonin, “Surface plasmons of metallic surfaces perforated by nano-hole arrays,” J. Opt. A, Pure Appl. Opt. 7(8), 422–426 (2005). [CrossRef]

10. J. Bravo-Abad, L. Martín-Moreno, F. J. García-Vidal, E. Hendry, and J. Gómez Rivas, “Transmission of light through periodic arrays of square holes: From a metallic wire mesh to an array of tiny holes,” Phys. Rev. B 76(24), 241102 (2007). [CrossRef]

11. A. Mary, S. G. Rodrigo, L. Martin-Moreno, and F. J. Garcia-Vidal, “Plasmonic metamaterials based on holey metallic films,” J. Phys. Condens. Matter 20(30), 304215 (2008).

12. J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and D. Robert, Meade, Photonic Crystals - Molding the Flow of Light, 2nd ed. (Princeton University Press, 2007).

13. K. Sakoda, Optical Properties of Photonic Crystal, ser. Optical Sciences. (Springer, New York, 2001).

14. A. Ruhe, “Algorithms for the nonlinear eigenvalue problem,” SIAM J. Numer. Anal. 10(4), 674–689 (1973). [CrossRef]

15. S. Shi, C. Chen, and D. Prather, “Revised plane wave method for dispersive material and its application to band structure calculation of photonic crystal slabs,” Appl. Phys. Lett. 86(4), 043104 (2005). [CrossRef]

16. B. P. Hiett, J. M. Generowicz, S. J. Cox, M. Molinari, D. H. Beckett, and K. S. Thomas, “Application of finite element methods to photonic crystal modelling,” IEE Proc. Sci. Meas. Technol. 149(5), 293–296 (2002). [CrossRef]

17. M. Davanco, Y. Urzhumov, and G. Shvets, “The complex Bloch bands of a 2D plasmonic crystal displaying isotropic negative refraction,” Opt. Express 15(15), 9681–9691 (2007). [CrossRef] [PubMed]

18. C. Fietz, Y. Urzhumov, and G. Shvets, “Complex k band diagrams of 3D metamaterial/photonic crystals,” Opt. Express 19(20), 19027–19041 (2011). [CrossRef] [PubMed]

19. G. Shvets and Y. A. Urzhumov, “Electric and magnetic properties of sub-wavelength plasmonic crystals,” J. Opt. A, Pure Appl. Opt. 7(2), S23–S31 (2005). [CrossRef]

20. J. Jin, The Finite Element Method in Electromagnetics, 2nd ed. (John Wiley & Sons, Inc., 2002).

21. F. Tisseur and K. Meerbergen, “The quadratic eigenvalue problem,” SIAM Rev. 43(2), 235–286 (2001). [CrossRef]

22. Z. S. Sacks, D. M. Kingsland, R. Lee, and J. F. Lee, “A perfectly matched anisotropic absorber for use as an absorbing boundary condition,” IEEE Trans. Antenn. Propag. 43(12), 1460–1463 (1995). [CrossRef]

23. S. D. Gedney, “An Anisotropic Perfectly Matched Layer-Absorbing Medium for the truncation of FDTD Lattices,” IEEE Trans. Antenn. Propag. 44(12), 1630–1639 (1996). [CrossRef]

24. D. Edward, Palik, Handbook of optical constants of solids, (Academic Press, 1985).

25. G. Steven, Johnson, “Notes on Perfectly Matched Layers (PMLs),” Lecture notes Massachusetts Institute of Technology Massachusetts (2008), 1–16 (2007).

26. S. G. Johnson, P. Bienstman, M. A. Skorobogatiy, M. Ibanescu, E. Lidorikis, and J. D. Joannopoulos, “Adiabatic theorem and continuous coupled-mode theory for efficient taper transitions in photonic crystals,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 66(6), 066608 (2002). [CrossRef] [PubMed]

27. M. Imran and G. Andrew, Kirk, “Accurate determination of quality factor of axisymmetric resonators by use of perfectly mathed layer,” arXiv:1101.1984v2 [physics.optics], 1–7(2011).

28. W. L. Barnes, T. W. Preist, S. C. Kitson, and J. R. Sambles, “Physical origin of photonic energy gaps in the propagation of surface plasmons on gratings,” Phys. Rev. B Condens. Matter 54(9), 6227–6244 (1996). [CrossRef] [PubMed]

29. P. Vincent and M. Neviere, “Corrugated dielectric waveguides: A numerical Study of the Second-Order Stop bands,” Appl. Phys. (Berl.) 20(4), 345–351 (1979). [CrossRef]

30. P. Lalanne, J. P. Hugonin, and P. Chavel, “Optical properties of deep lamellar grating: a coupled Bloch-mode insight,” J. Lightwave Technol. 24(6), 2442–2449 (2006). [CrossRef]

31. K. C. Huang, E. Lidorikis, X. Jiang, J. D. Joannopoulos, K. A. Nelson, P. Bienstman, and S. Fan, “Nature of lossy Bloch states in polaritonic photonic crystals,” Phys. Rev. B 69(19), 195111 (2004). [CrossRef]

32. A. Kolomenskii, S. Peng, J. Hembd, A. Kolomenski, J. Noel, J. Strohaber, W. Teizer, and H. Schuessler, “Interaction and spectral gaps of surface plasmon modes in gold nano-structures,” Opt. Express 19(7), 6587–6598 (2011). [CrossRef] [PubMed]

33. Y. Ding and R. Magnusson, “Band gaps and leaky-wave effects in resonant photonic-crystal waveguides,” Opt. Express 15(2), 680–694 (2007). [CrossRef] [PubMed]

Complex Bloch-modes calculation of plasmonic crystal slabs by means of finite elements method

Abstract

1. Introduction

2. The FEM simulation

3. The perfect matched layer (PML) implementation

4. 2D sinusoidal grating and PML test

5. 2-D periodic array

6. Conclusions

Acknowledgments

References and links

Cited By

Figures (8)

Equations (10)

Optics Express