1. Introduction

The study of the optical properties of Photonic Crystals (PCs) has relied on the generation of Bloch-mode photonic band structures. Information such as the existence of forbidden band gaps, phase and group velocity, group velocity dispersion, photonic densities of states, propagation anisotropy, etc. may be conveniently obtained by inspection or manipulation of the computed propagation bands [1, 2]. A commonly adopted approach for the calculation of PC band structures using frequency-domain methods consists in obtaining the allowed propagation frequencies ω for specific Bloch k-vectors as eigenvalues of an electromagnetic wave equation defined in a unit-cell of the periodic medium. If the PC incorporates only dispersionless materials (i.e., the dielectric constant ε(ω) constant for all frequencies), the problem is reduced to a linear generalized eigenvalue system of the form A_kx_k=ω ² B_kx_k [1–4].

An important class of metamaterials composed of a periodic arrangement of metallic inclusions in a dielectric host, or Subwavelength Plasmonic Crystals (SPCs) [5–7], are particularly useful in the UV and optical frequency ranges, where the permittivity of the constitutive metal inclusions displays a plasma-like dispersion. In [5], such structures were shown to display a negative refractive index propagation band in a subwavelength regime (~λ/10), and simulations in fact predicted the possibility of subwavelength-resolution focusing. In this case, the full variation of the metallic permittivities must be taken into account, which results in a nonliear equation of the form A_k(ω)x_k=ω ² B_k(ω)x_k, or, more generally, T_k(ω)x_k=0. Solutions to this type of problem are generally obtained through iterative algorithms -normally extensions of algorithms for linear problems [8], or developments from Newton’s Method [9]- in which successive approximations of ω are computed within a chosen convergence domain. Successful and efficient convergence usually requires a good initial guess; in addition matrix T(ω) and, depending on the method, its derivatives, must be recalculated at several frequencies within the search domain, causing the solution process to be time-consuming and computationally intensive. Hence, for this reason, band structure diagrams are not abundant in the metamaterials and plasmonic crystals literature.

Here, we revisit the Finite Element Method (FEM) formulation for the calculation of 2D photonic crystal bands developed by Hiett et al. [4]. This method is shown to yield a quadratic eigenvalue equation in the Bloch k-vector magnitude; frequency in this case is a parameter, so that material dispersion is readily taken into account. While the quadratic eigenvalue equation is nonlinear, it is a more tractable problem than the general nonlinear case above, and can yield solutions more efficiently. In addition, this formulation inherently yields bands of purely imaginary- and complex-wavevector Bloch modes, which may be particularly hard to obtain with nonlinear search routines. Imaginary k modes are generally available at the photonic bandgap frequencies and play an important role in representing the evanescent field inside a finite or semi-infinite PC slab under external excitation. This has been suggested in [10], in which an equivalent formulation to that presented here, based on the plane-wave expansion method [3], was used to calculate imaginary k bandgap modes in PCs composed of nondispersive, purely dielectric materials. On the other hand, modes with complex, rather than purely imaginary, k-vectors may also exist within, and even outside, bandgap regions that have considerable influence in the transmission and reflection transfer functions of finite plasmonic crystals.

In Sec. 2, we derive the FEM formulation for the plasmonic crystal eigenvalue problem following similar steps as [4]. This formulation is implemented in a simple way with the COMSOL Finite Elements package [11] and used to produce the photonic bands presented in Sec. 3, for a square-lattice SPC displaying a negative refraction band [5,6]. In these references, the negative refraction band was calculated only along one of the high-symmetry crystal directions; in the present article, the photonic bands are obtained over the entire irreducible first Brillouin zone, such that the isotropy of the negative refraction band near k=0 is revealed. A comparison between scattering data from a finite SPC and its band diagram is also performed, in which the role of the multiple existing Bloch modes (including those with imaginary and complex wavevectors) is determined. In particular, we report on the existence of resonances in transmission through a finite crystal slab that stem from the coexistence of two SPC bands in the same frequency range, being similar in origin and character to Fano resonances [12]. In Sec. 4 we present conclusions.

2. Finite-element method

We consider a periodic 2D medium with electromagnetic waves propagating in the xy plane. Solutions toMaxwell’s equations are classified as either TE (H=H_z ẑ, where H is the magnetic field) or TM (E=E_z ẑ, where E is the electric field), with the respective z-components obeying the wave equation

For TE waves, ϕ=H_z, p=1/ε, q=1 and ε is the electric permittivity; for TM waves, ϕ=H_z, p=1, q=ε. In Eq.(1), ω/c is the vacuum angular wavenumber.

The medium is periodic, and hence ε(r+T)=ε(r), where T is a lattice translation vector. From Bloch’s theorem we consider solutions ϕ(r)=u(r)exp(-i k·r), with u(r+T)=u(r) and k in the first Brillouin zone. Along with Eq.(1), this leads to:

Following the Galerkin procedure to reduce Eq.(1) to a weighted-residuals expression [13], we multiply that equation on both sides by a weight function w(r), and note that

and

to arrive at

This expression is integrated over a closed domain Γ with boundary δ _Γ that covers one unit cell of the periodic medium (the integration domain may actually extend over an integer number of unit cells). Making use of the Divergence theorem, the following weak-form expression is obtained:

In Eq.(6), n̂ is the outward normal unit vector to the boundary δ _Γ. Since u is periodic, the line integral vanishes. The resulting integrodifferential equation may be transformed into matrix format by following the usual FEM discretization procedure [4, 13]: the domain Γ is divided into several triangular subdomains (elements) in which locally supported expansion functions are defined; u is expanded in terms of such functions within each element; w is taken to be each one of the local expansion functions inside each element; and the material parameters p and q are allowed to be constant inside each element. Then the following quadratic matrix eigenvalue equation in k results:

Here, u is a vector containing the (complex) coefficients of the expansion of u, and matrices A through D may be individually related to each integral in Eq.(6) by inspection. The results in Sec. 3 are all calculated using second-order Lagrange elements. Explicit expressions for the matrices in this case can be found in [13]. The eigenvalue equation may be solved at a fixed frequency ω (and thus fixed ε(ω)) and k-vector direction. The eigenvalue itself is the k-vector amplitude, k. Due to the periodicity of the system, Re{k} is expected to be periodic in k-space. Note that, for real ε, all matrices are Hermitian, so that the eigenvalues k are either purely real or are complex conjugate pairs (k,k̄) [14]. The most common way of solving the Quadratic Eigenvalue Problem is by linearization, which results in a (linear) system twice the original size [10,14]. Other possibilities are to use nonlinear iterative methods such as Nonlinear Inverse Iteration or Rayleigh Coefficient Iteration [8].

3. Rectangular array plasmonic crystal

Here, the TE polarization (magnetic field perpendicular to the crystal plane) band structure of the plasmonic crystal studied in [5, 6], consisting of a 2D rectangular array of infinite metallic cylinders in air, is analyzed. The metal permittivity was considered to follow the Drude model where ε=1-ω ² _p/ω(ω-iω_c), with ω_p the plasma frequency and ω_c the damping frequency related to the mean electron collision rate. The square unit cell has a lattice parameter a=c/ω_p and the cylinder radii are R=0.45a. All of our results below are given in terms of the normalized frequency Ω=ω/ω_p.

For comparison, the transmission (t) and reflection (r) coefficients for plane-waves incident on a finite (i.e., with N of unit cells) crystal slab were obtained using the same technique as described in [15]. The incident plane wavew is in the Γ-X direction of the square-lattice Brillouin zone. Furthermone, an effective index n _eff is calculated using expressions from the scattering-parameter technique of Smith et al. [15, 16], i.e.:

with

and

These expressions model a 1D Fabry-Pérot cavity for waves propagating with a single wavevector at each frequency. This in principle precludes their direct application in situations when multiple waves associated with different k-vectors exist simultaneously - for instance in systems displaying Fano-type resonances as shown below.

Figure 1 shows the calculated band diagram in the Γ-X direction of the square lattice irreducible Brillouin zone, with ω_c=0, in terms of both k_r and k_i (k=k_r+ik_i). Modes for which k_i=0 are marked with crosses, while complex-k modes are marked with circles and triangles. Duplicates of such eigensolutions are also produced by the formulation, with k′_r=k_r±2m·π/a (k_r in the first B.Z., m=1,2, …) and identical k_i. Such solutions only arise because they satisfy the established periodic boundary conditions at the edges of the unit-cell; they do not participate in the Bloch-mode expansion of the total field in the crystal. Notice as well that since the calculated bands are for a crystal with zero damping losses (w_c=0), the existence of complex eigenvalues is completely unrelated to power dissipation.

Same as in the case of completely dielectric photonic crystals, real-k modes extend over the entire crystal, with a phase evolution given by the k(Ω) dispersion and with no amplitude variation apart from that of the periodic Bloch envelope. In Fig. 1, purely real eigenvalues are found in the ranges Ω≈0.2-0.36 and Ω≈0.54-0.6.

It has been suggested that imaginary k modes play an important role in representing the decaying field inside a finite or semi-infinite PC upon wave incidence at bandgap frequencies [10]. The inference that the imaginary part of the Bloch wavevector relates to the field decay rate in this situation is supported by our analysis. In Fig. 1, non-real eigenvalues may be either purely imaginary (k_r=0), or complex. Modes with purely imaginary k (found in the intervals Ω=0.43-0.52 and Ω=0.6-0.61) have no phase variation across a unit-cell, and are non-propagating. Complex k modes, on the other hand, present a phase variation concomitant with an exponential amplitude decay (or growth), which translate into interference effects in scattering from a finite crystal. In Fig. 1, complex k modes with |k_r| <π/a are found between Ω=0.36 and 0.43, and between Ω=0.52 and 0.54. Bands with k_i marked with triangles in Fig. 1(c) have real parts at the Brillouin zone boundaries, k_r=±π/a, such that the phase varies in full cycles across a unit cell. These solutions will henceforth be referred to as zone-boundary modes.

Figure 2(a) shows the transmission and reflection coefficients for a crystal with N=4 and ω_c/ω_p=0.001. In Figs. 2(b) and (c), the band structure with k_r and k_i on the horizontal axis are presented once again. The superimposed thick curves are the effective k-vectors k _eff=Ω·n _eff, with effective index n _eff obtained from from Eq. (8). A clear correspondence exists between the high-reflectivity band in the frequency range Ω=0.36-0.54 and the exclusive existence of non-real-k modes in Figs. 2(b) and 2(c). The correspondence between the the purely real band at Ω<≈ 0.36 and k _eff is also apparent - in this case, disregarding high-decay-rate k_r=±π modes, only single propagating modes are available at each frequency, which is an assumption of the n _eff expression, as pointed out above. In the range Ω≈0.43-0.56, k _eff coincides almost exactly with the lowest-decay-rate purely imaginary Bloch band. We infer that the waves travelling through the slab suffer an exponential attenuation given by the lower-amplitude imaginary Bloch k. In addition, waves related to the larger k_i modes either are not as strongly excited, or simply decay too quickly to be meaningful in the determination of k _eff.

 

Fig. 1. (a)Square-lattice plasmonic crystal band structure as function of k_r. (b) Band structure as function of k_i. Crosses indicate bands with k_i=0; circles, k_i≠0; and triangles, k_r=±π/a, k_i≠0

Download Full Size | PDF

 

Fig. 2. (a)Amplitude-squared transmission and reflection coefficients, (b) Band structure as function of k_r and (c) Band structure as function of k_i for a.

Download Full Size | PDF

The k _eff resonances observed in the interval Ω ≈ 0.53-0.6, are associated with Fabry-Pérot resonances of waves related to the purely real band within this range. It is apparent, however, that the complete transmission transfer function also involves the coexisting purely imaginary band. Evidence for this is that the imaginary part of k _eff follows the k_i band (see Fig. 2(c)), rather than tending to zero. Furthermore, the transmission transfer functions for the N=4 structure in Fig. 3(b) can be roughly approximated by a function of the form t(Ω)=α(Ω)+i·β (Ω)·t _FP(Ω), with t _FP(Ω) the Fabry-Pérot cavity transmission coefficient. Examples of this may be seen in Fig. 3(a), in which the transmission coefficient amplitude |t| for 4- and 3-unit cell structures are compared to those obtained from the expression above, under the following conditions: the Fabry-Pérot cavity single-pass phase varied as k_r ·N · a, with k_r from the purely real band (N is the number of unit cells); the cavity reflection coefficient was taken as r=0.92; β=0.5; the term α(Ω) was made to vary as exp(-k_i·N·a), where k_i is from the purely imaginary band; the total transmission was normalized to its maximum value. General features such as sharp peaks (related to resonances of t _FP) and alternating maxima and minima in-between sharp resonances (related to minima of t _FP with phases of respectively odd and even multiples of π/2), are similar in both numerical and analytical curves, and for both cavity lengths. Thus, the terms α(Ω) and β (Ω) may be related to the amplitudes of two waves at the output side of the slab: one which is transmitted without a phase change, corresponding to the purely imaginary k-vector Bloch mode, and a second that propagates dispersively given by the real-k band, leading to the Fabry-Pérot resonances. In this case, both purely real and imaginary bands have relevant contributions to the field within the finite slab and thus a strong influence on k _eff. Conversely, n _eff obtained from Eq.(8) does not allow one to discern the propagation characteristics of the different SPC modes. Notice as well that from Fig. 3(b), transmission resonances are quite sharp with null damping losses (ω_c=0), with peaks approaching unity. Introduction of even small damping (ω_c=0.001ω_p) causes strong degradation of the resonance quality factor. The asymmetric shape of the transmission transfer function around resonant peaks is similar to that of Fano absorption lineshapes, which in quantum systems arise from the interference of a discrete state with a continuum [12]. The transmission resonances presented here have a similar origin to these, namely the interference between waves forming the Fabry-Pérot resonant peaks and the continuum formed by the evanescent waves.

 

Fig. 3. (a) Comparison between numerically calculated transmission coefficient amplitudes for 3- and 4-unit-cell crystal slabs (thin lines) and the analytical expression |t(Ω)|=A·| exp(-k_i(Ω)·N·a)+i·β· t _FP(k_i(Ω)·N·a)| (thick lines), with k_i(Ω) and k_r(Ω) taken from the purely real and imaginary Bloch bands respectively. (b)Amplitude and phase of the transmission transfer function for the 4-unit-cell crystal with ω_c/ω_p=0.001 (thick line) and ω_c=0 (thin line).

Download Full Size | PDF

Finally, it is important to point out that the upper real-k bands discussed above are negative refraction bands, with group velocity v_g of opposite sign to the phase velocity ω/k; however, as shown in [5,6], the smallest number of rods that can be meaningfully analyzed in the search for optical magnetism is four, since each rod must be surrounded by other rods to create a magnetically active resonance.

In the range Ω≈0.36-0.43, the transmitted power is extremely low, yet the transfer function presents a series of anti-resonances corresponding to the observed peaks in Im{k _eff} with concomitant step-transitions in Re{k _eff}. At the same time, disregarding the resonances, both real and imaginary parts of the effective wavevector closely follow the corresponding photonic bands. Hence, this behavior has an origin similar to the resonances discussed above, given the existence of two evanescent Bloch bands, as seen in Figs. 1(a) and 1(b).

We next analyze the effects or small damping losses in the calculated band structure. In Fig. 4(a) and 4(b), the real and imaginary k bands for a lossless and a lossy structure with ω_c/ω_p=0.001 are compared. A noticeable difference between the two cases is the disappearance of the the small bandgap centered at Ω ≈ 0.36. Whereas in the lossless case the purely real k_r-bands reach the Brillouin zone boundary at gap frequencies, in the lossy structure the corresponding bands do not reach k_r=π/a, turning around at some k _max <π/a where the group velocity defined as v_g=dω/dk_r becomes infinite. Notice that in the lossless case the group velocity as defined above is infinite throughout the photonic bandgap, since k_r=π/a everywhere in this range. It is apparent that the band turn-around observed in the lossy crystal case is caused by a less-efficient Bragg wave-interference in the presence of damping losses. This effect was also observed in [17], for photonic crystals containing polaritonic inclusions. It must be pointed out that within the bandgap frequency region, the group velocity as defined above does not correspond to the velocity with which an incident wavepacket would traverse a finite length of the SPC, so that v_g>c does not imply superluminal propagation [18].

 

Fig. 4. (a) Band structure for (a)k_r and (b)k_i. Dots correspond to ω_c=0, crosses to ω_c/ω_p=0.001.

Download Full Size | PDF

As mentioned above, in the lossless case, the matrices involved are Hermitian, such that eigenvalues are either purely real or come in complex conjugate pairs. The latter case is observed in Ω=0.36-0.43, where two degenerate modes at each k_r exist. The degeneracy is evident from the opening of bifurcations in k_i at Ω≈ 0.36 and splitting of k_i bands in Ω=0.36-0.43 due to small damping losses. The singled-out band in Fig. 4(a) is a perturbed version of a zone-boundary band that exists at Ω⪷0.365 in the lossless case. Above Ω⪷0.365, the latter ceases to be a zone-boundary band, and becomes a degenerate complex-k band. This transition is evident in the departure of the perturbed band from k_r=π/a and its merging with a second complex-k band around Ω=0.36.

Finally, to show that the method may be conveniently applied to the generation of bands in arbitrary directions in the k-vector space, Figs. 5 and 6 display, respectively, the low-frequency and two high-frequency purely real-k bands. These diagrams were obtained simply by changing the direction of the unit k-vector in Eq. (6).

 

Fig. 5. Low-frequency, purely real bands of the square-lattice plasmonic crystal over the first Brillouin zone.

Download Full Size | PDF

Figure 5 shows that the lowest purely real band is highly isotropic at lower frequencies, however becomes anisotropic as the photonic bandgap is approached. Figure 6 shows slices of the upper real bands, which in reality form a closed surface. The negative group velocity band Ω=0.57-0.6 is highly isotropic near its top, thus making the square array of circular rods a perfectly isotropic negative-index metamaterial for the range of frequencies Ω=0.58-0.6.

 

Fig. 6. High-frequency, purely real bands of the square-lattice plasmonic crystal over the first Brillouin zone.

Download Full Size | PDF

On the other hand, the longitudinal-wave, or Bulk Plasmon, band near Ω=0.0.61 predicted in [5] displays very strong anisotropy near the Γ-point. This anisotropy is apparent in Fig. 7, where the Bulk Plasmon band was plotted for the Γ-M (0°) and Γ-X (45°) directions.

 

Fig. 7. Bulk Plasmon band along Γ-M (0°) and Γ-X (45°) directions.

Download Full Size | PDF

It is important to notice that, in [5] and [6], the origin of the negative refraction band above was related to the hybridization of quasi-static multipole resonances of the plasmonic rods in close proximity to each other. Over the negative-index band, the lattice constant is ~λ/10, where λ is the vacuum wavelength. Thus we believe that negative refraction occurs in the longwavelength rather than in the Bragg diffraction regime.

4. Conclusions

The Finite-Element formulation presented in Section 2 is convenient for the calculation o photonic band structures of periodic metamaterials composed of dispersive materials. The formulation yields a quadratic eigenvalue equation system that may be easily implemented with commercial, generic Finite-Element packages such as COMSOLMultiphysics, and that may be more efficiently solved than the general nonlinear eigenvalue problem that is normally solved in the metamaterials community. In addition, calculation of bands in arbitrary directions in k-space simply involves determination of the components of the unit k-vector, which is convenient for the analysis of medium anisotropy.

An important advantage associated with solving the quadratic eigenvalue problem is that it yields Bloch eigenmodes with not only real but also complex wavevectors, which, as suggested in [10], may be associated with many effects observed in the scattering of waves from finite crystals, especially at, but not limited to, bandgap frequency ranges. A clear indication of this in fact is the existence of the Fano-type resonances in transmission as shown in Section 3.

It is clear from the analysis in Section 3 that the effective index n _eff determined from Eq.(8) via scattering from finite crystals will not completely coincide with effective indices of particular Bloch modes when multiple such modes exist that influence the scattering process at similar levels. At the same time, just the Bloch band structure by itself does not give an immediately clear picture of wave scattering and therefore should be used with some attention in predicting the scattering response of finite structures. At any rate, both methods may be used in conjunction to form a more complete understanding of the propagation characteristics of a metamaterial.

Finally, it must be pointed out that alternative formulations to the present finite element-based one presented in this article may also lead to quadratic eigenvalue equations, see for instance the plane-wave expansion formulation used in [10].