Complex k band diagrams of 3D metamaterial/photonic crystals.

Chris Fietz; Yaroslav Urzhumov; Gennady Shvets

doi:10.1364/OE.19.019027

1. Introduction

The numerical simulation of electromagnetic fields inside both metamaterial and photonic crystals is an important tool for analyzing these periodic structures. In particular, the eigenmodes of crystals, defined as freely propagating waves not coupled to external currents, are often of the most interest. The conventional method [2,3] of numerically solving for crystal eigenmodes is to define the geometry of the unit cell of the crystal of interest and the differential equation that the fields must obey in this geometry and then impose Bloch periodic boundary conditions. The partial differential equation problem is then discretized using one of the many standard methods (finite element, finite integral, finite difference, etc.) thereby turning it into an algebraic eigenvalue problem with a finite number of degrees of freedom and the frequency ω as the eigenvalue. This finite sized eigenvalue problem is then solved numerically. An important detail of this method is that the Bloch wavenumber k is chosen beforehand, and the frequency is then computed as a function of the wavenumber, yielding the dispersion curves ω = ω(k). This is the most commonly used method for calculating dispersion curves of the electromagnetic waves propagating in photonic crystals or in closely related metamaterial crystals.

There are however, many instances where it is more convenient to specify the frequency ω and solve for the wavenumber as a function of frequency: k = k(ω). At least five such instances can be identified. First, metamaterials often contain dispersive materials such as metals, where the dielectric function strongly depends on the frequency of the wave. In this case, the eigenfrequency problem becomes a nonlinear eigenvalue problem and must be solved iteratively [4]. In contrast, when solving for the wavenumber as a function of frequency, the resulting eigenvalue problem only needs to be solved once. Second, to describe the microscopic electromagnetic response of a metamaterial inclusion many researchers use tabulated data for the permittivity, such as the data found in Ref. [5]. These tabulated material parameters are functions of a real valued ω and should not be used for ω(k) eigenvalue simulations where ω is often complex valued. Third, it is often useful to solve for the wavenumber as the eigenvalue because of the information contained in the complex wavenumber including the decay lengths of the electromagnetic modes (either due to finite dissipation or because of the evanescent nature of the mode) and the figure of merit [6] of negative index modes. Fourth, in the majority of experiments electromagnetic fields inside metamaterial/photonic crystals are excited by external sources producing time-harmonic fields with real valued frequencies. A complex wavenumber eigenvalue simulation provides the correct field distribution in the photonic crystal relevant to such an experiment. Finally, this approach provides a very natural way of calculating the so-called isofrequency surfaces corresponding to ω(k) = const, where ω is real. Isofrequency diagrams are fundamentally important for predicting wave refraction at the PC interfaces [7] and for calculating the density of states [8].

There are several previously published methods on calculating k(ω) dispersion curves including variations of the plane-wave expansion method [9, 10] and diagonalizing the crystal transfer matrix [11, 12]. A method of solving for complex wavenumber dispersion curves using the FEM has been proposed [13, 14] but only for 2D crystals. The benefits of a complex wavenumber 2D FEM are becoming better appreciated and use of this method is becoming more common [15–19]. It is therefore timely to generalize this method of Complex Wavenumber Eigenvalue Simulations (CWES) from two to three dimension, which is the object of this paper. This 3D complex wavenumber eigenvalue simulation was recently used as part of a metamaterial homogenization procedure [20].

The basic theory behind solving for complex wavenumber eigenvalues using the FEM discretization is explained in Sec. 2. The underlying field equations for the electric and magnetic fields and the necessary boundary conditions are discussed. In Sec. 3 we apply this method to a 3D photonic crystal consisting of non-dispersive dielectric spheres. The photonic band structure for electromagnetic waves propagating both parallel and obliquely to the principal axes is calculated, and different types of modes (transverse and longitudinal) are identified. This simple model [21] can be solved with the commercial FEM software package COMSOL Multiphysics. In Sec. 4 we calculate the dispersion curves for a negative index fishnet metamaterial [22] and demonstrate the existence of four distinct negative index modes. We also calculate the two-dimensional isofrequency contours for the least damped transverse mode and demonstrate from the shape of the isofrequency contours that this transverse mode is hyperbolic.

2. The finite element eigenvalue problem

2.1. The field equation

In this section we present the FEM formulation for solving for the magnetic field. Electromagnetic wave propagation is described by the Maxwell equations which can be rearranged into a wave equation for either the electric field E or the magnetic field H. The wave equation for the magnetic field is

\nabla \times (\frac{1}{ɛ} \nabla \times H) - μ \frac{ω^{2}}{c^{2}} H = 0 .

Here ɛ(x) and μ(x) are the microscopic permittivity and permeability of the metamaterial/photonic crystal of interest. Due to the periodic nature of the crystal both are are assumed to be scalar functions periodic in the crystal lattice. According to Bloch’s theorem [2, 3] the magnetic field can be represented as the product of a periodic function and an exponential factor

H (x) = u (x) exp [i (ω t - k \cdot x)],

where ω is the frequency of the wave and k is the wavevector of the Bloch-Floquet wave. u(x) is a vector function which is periodic in the crystal lattice. By inserting Eq. (2) into Eq. (1) we obtain an equivalent field equation for u

\frac{k^{2}}{ɛ} u - \frac{k}{ɛ} (k \cdot u) - i k \times (\frac{1}{ɛ} \nabla \times u) - i \nabla \times (\frac{1}{ɛ} k \times u) + \nabla \times (\frac{1}{ɛ} \nabla \times u) - μ \frac{ω^{2}}{c^{2}} u = 0,

which can be interpreted as an eigenvalue problem and solved for the Bloch wavenumber k as the eigenvalue. The spatial profile of the eigenmode u(x) is also recovered providing the magnetic field profile according to Eq. (2) and the electric field profile according to

E (x) = \frac{1}{i ɛ ω / c} \nabla \times H = \frac{1}{i ɛ ω / c} (- i k \times u + \nabla \times u) exp [i (ω t - k \cdot x)]

2.2. The finite element model

There are several commercial FEM software programs (COMSOL Multiphysics by COMSOL, HFSS by Ansys, Vector Fields Opera by Cobham Technical Services, etc.) that are available for modeling metamaterial/photonic crystals. These commercial software packages provide a convenient graphical user interface for defining a crystal’s geometry, meshing the computational domain, and visualizing the electromagnetic fields. This allows for models to be quickly developed and tested. Of the many commercial FEM codes available, the authors of this paper are only aware of one (COMSOL Multiphysics) that allows the user to specify the field equation to be solved. The simulation examples and results presented here were obtained using COMSOL. The FEM model used for the photonic crystal described in Sec. 3 is available for download [21]. In what follows, only the most essential features of the FEM approach are reviewed; more detailed treatments can be found elsewhere [23–25].

The FEM is based on setting the integral of a so-called weak expression over the domain of interest to zero. Doing so ensures the field equation is satisfied and also creates boundary conditions. The weak expression corresponding to Eq. (3) is

\begin{array}{l} F_{H} (v, u) = \frac{k^{2}}{ɛ} v \cdot u - \frac{1}{ɛ} (k \cdot v) (k \cdot u) - i \frac{1}{ɛ} v \cdot [k \times (\nabla \times u)] - i (\nabla \times v) \cdot \frac{1}{ɛ} (k \times u) \\ + (\nabla \times v) \cdot \frac{1}{ɛ} (\nabla \times u) - μ \frac{ω^{2}}{c^{2}} v \cdot u, \end{array}

where v(x) is a test function. When the integral of the weak expression over the unit cell Ω of the crystal is set to zero, integrating by parts gives us two separate integrals (Eq. (6)). The first integral enforces the field equation. The second integral is over the boundary of the domain and represents a natural boundary condition [23, 24],

\begin{array}{l} 0 & = & \int_{Ω} d^{3} x F_{H} (v, u) \\ = & \int_{Ω} d^{3} x v \cdot [- \frac{1}{ɛ} k \times (k \times u) - i \frac{1}{ɛ} k \times (\nabla \times u) - i \nabla \times (\frac{1}{ɛ} k \times u) \\ + \nabla \times (\frac{1}{ɛ} \nabla \times u) - μ \frac{ω^{2}}{c^{2}} u] \\ + \oint_{\partial Ω} d A v \cdot [\hat{n} \times \frac{1}{ɛ} (- i k \times u + \nabla \times u)], \end{array}

where n̂ is the vector normal to the boundary. On an external boundary, the natural boundary condition enforced by the integral in Eq. (6) over the boundary ∂Ω forces the expression n̂ × (−ik × u + ∇ × u)/ɛ to be equal to zero. Recalling Eq. (4) we note that this simply enforces the boundary condition n̂ × E = 0. This is known as the perfect electric conductor or PEC boundary condition. This natural boundary condition is the default if no other boundary condition is enforced. On an internal boundary within the unit cell the surface integrals over each side of the boundary must be equal to each other. The effect is that the tangential components of the electric field must be continuous across the internal boundary or n̂ × E ⁺ = n̂ × E ⁻ where E ⁺ and E ⁻ are the electric fields on opposite sides of the internal boundary.

The periodicity of u is enforced by imposing periodic boundary conditions on the exterior boundaries of the unit cell. In COMSOL, these periodic boundary conditions override the natural boundary condition. However, if a PEC boundary condition is desired inside the unit cell (e.g., on the surface of a metal inclusion) this can be accomplished by removing the subdomain representing the metal inclusion. Now only the exterior side of the metal boundary remains and the natural boundary condition forces the tangential electric fields to zero on this boundary.

If a perfect magnetic conductor or PMC boundary condition (n̂ × H) is desired while solving for the magnetic field, this can be enforced with constraints [23] on the tangential magnetic field on the boundary.

In order to turn Eq. (6) into an eigenvalue problem, the three degrees of freedom that comprise the Bloch wavevector k must be reduced to one by restricting two degrees of freedom. This is accomplished by setting k = k ₀ + λ k̂ _n where λ will be the eigenvalue solved for, k ₀ is an offset vector and k̂ _n is a unit vector (k̂ _n · k̂ _n = 1). Note that since there are many possible ways to restrict these two degrees of freedom, the resulting complex k(ω) vs. ω dispersion curves are in general not unique. The FEM turns the weak form and accompanying boundary conditions into an algebraic problem, in this case a quadratic eigenvalue problem [26]:

A \vec{u} + λ B \vec{u} + λ^{2} C \vec{u} = 0,

where A, B and C are N × N matrices and u⃗ is an N × 1 solution vector. N is the number of degrees of freedom of the discretized system. Terms in the weak form (Eq. (5)) that are zero, first and second order in λ contribute to the A, B and C matrices respectively. This quadratic eigenvalue problem can be recast in the form of a generalized eigenvalue equation

(\begin{matrix} A & B \\ 0 & 1 \end{matrix}) (\begin{matrix} \vec{u} \\ λ \vec{u} \end{matrix}) = λ (\begin{matrix} 0 & - C \\ 1 & 0 \end{matrix}) (\begin{matrix} \vec{u} \\ λ \vec{u} \end{matrix}) .

When using COMSOL to solve the FEM problem, this linearization is performed automatically during the solution stage.

2.3. Electric field formulation

The previous discussion focused on solving for the magnetic field H or rather the periodic function u equal to the magnetic field with the exponential Bloch factor removed. This is especially convenient when an inclusion requires a PEC boundary condition since that is the natural boundary condition when solving for H. However, solving for the electric field is very similar to solving for the magnetic field. The wave equation for the electric field for a free wave is

\nabla \times (\frac{1}{μ} \nabla \times E) - ɛ \frac{ω^{2}}{c^{2}} E = 0 .

As before, we replace the electric field with a periodic vector field times an exponential factor

E (x) = u (x) exp [i (ω t - k \cdot x)],

producing the new field equation

\frac{k^{2}}{μ} u - \frac{k}{μ} (k \cdot u) - i k \times (\frac{1}{μ} \nabla \times u) - i \nabla \times (\frac{1}{μ} k \times u) + \nabla \times (\frac{1}{μ} \nabla \times u) - ɛ \frac{ω^{2}}{c^{2}} u = 0.

The corresponding weak form for this field equation is

\begin{array}{l} F_{E} (v, u) & = & \frac{k^{2}}{μ} v \cdot u - \frac{1}{μ} (k \cdot v) (k \cdot u) - i \frac{1}{μ} v \cdot [k \times (\nabla \times u)] - i (\nabla \times v) \cdot \frac{1}{μ} (k \times u) \\ + (\nabla \times v) \cdot \frac{1}{μ} (\nabla \times u) - ɛ \frac{ω^{2}}{c^{2}} v \cdot u, \end{array}

which is equivalent to Eq. (5) if ɛ and μ are interchanged. Integrating this weak form over the crystal unit cell by parts and setting its value to zero again gives produces two integrals, one enforcing the field equation and a surface integral enforcing the boundary condition n̂ × H = 0. Thus the PMC boundary condition is the natural boundary condition when solving for the electric field.

3. Example: dielectric photonic crystal

For a demonstration of the CWES method of calculating complex k dispersion curves we use a simple photonic crystal as an example. The unit cell, pictured in Fig. 1, is a cube with a dielectric sphere at the center surrounded by vacuum. The sphere has a radius of 0.3a, where a is the lattice constant of the cubic array, and a permittivity of ɛ = 5 – i0.01. Arrays of spheres as well as arrays of point dipoles are both systems that have been studied extensively by metamaterial researchers [27–30].

Fig. 1 Complex k dispersion curves and field profiles for eigenmodes of the photonic crystal pictured in the inset assuming k ₀ = 0 and k̂ _n = x̂. (a) Real part of k_x(ω) for a transversely polarized mode and a diagram of the crystal unit cell. (b) Imaginary part of k_x(ω) for a transversely polarized mode and a field profile for the ẑ polarized transverse mode. There are two transverse modes, ŷ and ẑ electrically polarized, which are degenerate. (c) Real part of k_x(ω) for two longitudinally polarized modes and a field profile for the magnetic longitudinal mode. (d) Imaginary part of k_x(ω) for two longitudinally polarized modes and a field profile for the electric longitudinal mode. The longitudinal mode with the passband near ω = 4.5c/a is magnetically polarized in the x̂ direction and the longitudinal mode with the passband near ω = 5c/a is electrically polarized in the x̂ direction. The longitudinal modes correspond to Mie’s dipole resonances. For all dispersion curves the dotted lines are the result of a conventional ω(k) eigenvalue simulation. For all field profiles the frequency is ω = 2c/a with arrows representing D_y and D_z and color representing D_x.

Download Full Size | PDF

As mentioned in Sec. 2.2, it is necessary to restrict two of the three degrees of freedom of the Bloch wavevector k. There are many possible ways to do this. As the first example, we calculate the dispersion curves corresponding to the propagation along a principal axis. To simulate this we set k ₀ = 0 and k̂ _n = x̂. The results of this eigenvalue simulation for the frequency range 1c/a ≤ ω ≤ 5.5c/a are plotted in Fig. 1 as ω vs. k_x = x̂ · k = λ.

For clarity, we have only plotted the three least evanescent (i.e. possessing the smallest values of Im(k_z)) eigenmodes. These three k(ω) dispersion curves are plotted as solid lines. For comparison, we have also ploted a conventional ω(k) dispersion curve as dotted line. The three eigenmodes in Fig. 1 are described as either transverse or longitudinal according to their polarization. The symmetry of the dispersion curves is such that for every solution k(ω) there is also the solution –k(ω) indicating that this is a reciprocal crystal. The dispersion curves for the transverse modes plotted in Fig. 1 in fact represent two polarization-degenerate modes because of the symmetry of the crystal. The longitudinal mode with the passband near ω = 4.5c/a is magnetically polarized in the x̂ direction making it a magnetic bulk plasmon. The longitudinal mode with the passband near ω = 5c/a is electrically polarized in the x̂ direction and is an electric bulk plasmon. The field profiles of both longitudinal modes indicate that the passbands correspond to Mie’s resonances of the dielectric sphere [31].

The transverse mode dispersion curve has a band in the approximate frequency range 4.6c/a < ω < 4.8c/a with a large value of Im(k_x), indicating it is an evanescent band, but a Re(k_x) that is equal to neither 0 nor π/a as is typical of ω(k) dispersion curves. As described in Refs. [13, 26] for a quadratic eigenvalue problem with hermitian matrices (corresponding to a lossless crystal) the eigenvalues must always be real or come in complex-conjugate pairs. The dielectric photonic crystal under consideration has very low loss, so this lossless condition approximately holds true for the dispersion curves in Fig. 1. The transverse band in the 4.6c/a < ω < 4.8c/a frequency band is one half of a complex conjugate pair, the other half is a transverse doubly degenerate mode not shown here. At the frequency of ω ≈ 4.8c/a the two modes that make up this complex conjugate pair both enter a passband and split, the plotted mode going to the Γ point and the unplotted mode going to the band edge (this unplotted mode corresponds to the dotted lines in Fig. 1 from the ω(k) simulation). Note that in this passband there are two pairs of propagating doubly polarization degenerate modes or four propagating modes in total.

The transverse eigenmodes plotted in Fig. 1 can be excited by a plane wave normally incident onto the vacuum/photonic crystal interface if the interface is parallel to the y-z plane. The longitudinally polarized modes could not be excited by a normally incident wave without the aid of a coupling device at the interface. If the incident beam of light is not normal to the interface, if for example the incident beam has a wavenumber laying in the x-y plane but at a 30° angle from normal then a different set of eigenmodes will be excited at the interface. To simulate these excited eigenmodes we set k ₀ = ω/csin(π/6)ŷ and k̂ _n = x̂ and solve the resulting eigenvalue problem. The resulting photonic band structure is plotted in Fig. 2.

Fig. 2 Complex wavenumber dispersion curves and field profiles for eigenmodes of the photonic crystal pictured in Fig. 1(a) assuming k ₀ = ω/csin(π/6)ŷ and k̂ _n = x̂. Modes excited by p or s polarized incident light are plotted with solid or dashed lines respectively. (a) Real part of k_x(ω) for two transverse hybrid modes and a expanded view of the avoided crossing in Re(k_x) space. (b) Imaginary part of k_x(ω) for two transverse hybrid modes, an expanded view of the avoided crossing in Im(k_x) space (plotting the same modes as the expanded view in Re(k_x) space), and a field profile for the E_z polarized transverse hybrid mode. (c) Real part of k_x(ω) for two longitudinal hybrid modes and a field profile for the magnetic longitudinal hybrid mode. (d) Imaginary part of k_x(ω) for two longitudinal hybrid modes and a field profile for the electric longitudinal hybrid mode. The magnetic longitudinal hybrid mode is excited by s polarized incident light and the electric longitudinal hybrid mode is excited by p polarized incident light. For all field profiles the frequency is ω = 2c/a with arrows representing D_y and D_z and color representing D_x.

Download Full Size | PDF

The eigenmodes in Fig. 2 are roughly split into predominantly-transverse (transverse hybrid) and predominantly-longitudinal (longitudinal hybrid) modes. The hybridization between the transverse and longitudinal modes is caused by the finite symmetry-breaking k_y. Both the transverse and longitudinal hybrid modes in Fig. 2 can be characterized by the polarization of the incident light that couples to them. For a p polarized incident beam (electric field in the x–y plane) the p polarized eigenmode is excited (plotted in Fig. 2 with solid lines) and for an s polarized incident beam (electric field in the ẑ direction) the s polarized eigenmode is excited (plotted in Fig. 2 with dashed lines).

At the frequency ω ≈ 4c/a the transverse hybrid modes and the longitudinal hybrids modes appear to cross in a propagating band. An expanded view of this region in Fig. 2 plotting both transverse and longitudinal hybrid modes shows that the apparent crossing actually occurs in a band gap. Viewed in complex k_x space it is clear that this is actually an avoided crossing and in the band gap the transverse and longitudinal hybrid modes form complex conjugate pairs [13, 26].

We see that even for a simple photonic crystal the CWES method of calculating the dispersion curve produces rich and complex results. In particular, it is not possible to solve for evanescent eigenmodes using the conventional ω(k) method for calculating dispersion curves.

4. Example: negative index plasmonic fishnet

The second example highlighting the versatility of the CWES method is a plasmonic negative index metamaterial (NIM) shown in Fig. 3. Because this metamaterial contains dispersive (plasmonic) components, it is even more convenient to use this method because the dielectric permittivities of metals are tabulated for real-valued frequencies. Recent experiments [22] demonstrated that the so-called fishnet metamaterial supports a negative index eigenmode for near infrared wavelengths of about λ ₀ ≈ 1.7μ m. The dimensions and composition of the unit cell taken from Ref. [22] are shown in Fig. 3. The fishnet metamaterial is made of alternating layers of Ag and MgF₂ with thicknesses of 30nm and 50nm respectively. This layered structure was milled with a focused ion beam into crisscrossing strips with widths of 265nm and 565nm. The crystal lattice constants are a_x = a_y = 860nm and a_z = 80nm. For the dielectric response of the Ag we used the Drude model $ɛ_{A g} = 1 - ω_{p}^{2} / (ω (ω - i γ))$ with the same parameters cited in Ref. [22]: ω_p = 9eV and γ = 0.054eV. It should be noted that in Ref. [22] the scattering frequency γ was increased by a factor of three from that of standard bulk Ag to account for additional scattering losses. For the dielectric response of MgF₂, no values were provided in the original paper so we used an oscillator model of the dielectric permittivity from Ref. [32].

Fig. 3 The fishnet metamaterial from Ref. [22]. The lattice constants of the unit cell are a_x = a_y = 860nm and a_z = 80nm. The fishnet is made up of alternating layers of Ag with a thickness of 30nm and MgF₂ with a thickness of 50nm. The widths of the the crisscrossing fishnet strips are b ₁ = 265nm and b ₂ = 565nm.

Download Full Size | PDF

The negative index mode reported in Ref. [22] propagates in the ẑ direction. Therefore, we have applied CWES to the specific case of k ₀ = 0 and k̂ _n = ẑ. The resulting dispersion curves are plotted in Fig. 4 for the 100THz < ω/2π < 300THz frequency range. For clarity we have only plotted the four eigenmodes that have the smallest values of Im(k_z) and therefore the longest propagation lengths. The negative index mode found in Ref. [22] is labeled E_x to indicate that it is transversely polarized with the electric field pointing in the x̂ direction. We can see that it is a negative index mode because (according to the convention where fields are Bloch periodic with the exponential factor exp[i(ωt – k · x)]) a negative index mode should have a wavenumber whose real and imaginary parts have the same sign. In Fig. 4, that describes the E_x mode for frequencies below 200THz. In addition, we observe three other negative index modes, though they generally have a shorter propagation lengths than the E_x mode. Of these three modes, the two labeled E_y are transversely polarized with the electric field pointing in the ŷ direction and the mode labeled E_z is longitudinally polarized with the electric field pointing in the ẑ direction.

Fig. 4 TOP: Real (a) and imaginary (b) parts of k_z(ω) for several eigenmodes of the fishnet metamaterial shown in Fig. 3. In addition to the transverse mode electrically polarized in the x̂ direction labeled E_x which was identified in Ref. [22] we see two other transverse modes electrically polarized in the ŷ direction labeled E_y and a longitudinal mode electrically polarized in the ẑ direction labeled E_z. MIDDLE: Field profiles for the transverse E_x mode (c) and the longitudinal E_z mode (d) on a cross-section laying on the x-y plane in the middle of the MgF₂ layer. Arrows represent in-plane electric field and color represents the E_z field. BOTTOM: (e) Field profile of the transverse E_x mode on a cross-section laying on the x-z plane halfway between two thin Ag-MgF ₂ strips. Arrows represent in-plane electric field and color represents the H_y field. For all field profiles the frequency is 175THz and each region is labeled Ag or MgF ₂ according the the material of the region. Unlabeled regions are vacuum.

Download Full Size | PDF

The figure of merit (FOM) is a number commonly used to quantify the quality of a negative index material. It is often defined [6] as the ratio between the real and imaginary parts of the index of refraction. For eigenmodes of the fishnet crystal with real ω and a complex k pointing in the ẑ direction this is equivalent to the ratio between the real and imaginary parts of k_z or FOM≡ Re(k_z)/Im(k_z). Without knowledge of the complex wavenumber of an eigenmode of a metamaterial, the figure of merit must be calculated indirectly. For example, in Ref. [22] the FOM is estimated numerically from simulated transmission through a fishnet sample. With knowledge of the complex wavenumbers of the crystal eigenmodes, it is possible to calculate the FOM directly. The FOM for each of the four eigenmodes identified in Fig. 4 is plotted in Fig. 5(a), along with the propagation lengths plotted in Fig. 5(b). A positive FOM corresponds to the real and imaginary parts of k_z having the same sign, therefore indicating a negative index eigenmode. We see again that all four modes are in fact negative index modes in a fairly broad 150THz < ω/2π < 200THz frequency range. Somewhat unexpectedly, the mode experimentally observed in Ref. [22] (labelled E_x) has the lowest FOM despite having the longest propagation length. The theoretically predicted FOM≈ 7.5 of the E_x mode is a factor of 2 higher than the experimentally measured FOM [22], which could be attributed to fabrication imperfections.

Fig. 5 The FOM (a) and propagation length (b) for the four fishnet eigenmodes identified in Fig. 4. Note that though the mode labeled E_x has the largest propagation length it also has the smallest FOM. The negative value of FOM for the E_x mode for frequencies above 220THz indicates that it is no longer a negative index mode at these higher frequencies.

Download Full Size | PDF

The existence of multiple negative index waves in the fishnet structure has important experimental implications. For some frequencies (e.g., at 150THz) the propagation length of the longitudinal mode (E_z) is as long as that of the primary transverse mode (E_x). While the original experiments [22] only excited the E_x mode by using the light normally incident onto the x–y vacuum/fishnet interface, one can envision exciting both E_x and E_z modes using obliquely incident light. For example, if p-polarized light is obliquely incident with the incident wavevector laying in the x-z plane (electric field laying in the x-z plane), then both E_x and E_z modes will be launched into the fishnet with different phase velocities. Their refraction by the fishnet prism [22] would give rise to two distinct beams. The existence of the additional strongly dispersive longitudinal modes (bulk plasmons) also suggests that the fishnet is a metamaterial with significant spatial dispersion that can have a strong effect on the surface waves at the vacuum/fishnet interface [33, 34]. Strong spatial dispersion is not unexpected because the fishnet’s transverse period is about λ/2.

Recent theoretical [35] and experimental [36] studies demonstrated that negative index propagation in fishnet structures is not limited to the optical frequency range and can, in fact, be observed with microwaves. The main physical distinction (excluding their appropriately scaled sizes) between microwave and optical structures is that the metals in the former can be accurately described as PECs. On the other hand, in the optical range metals are described as being ”plasmonic”, which means that optical field penetration into the metal is significant. Our simulation can quantify how important the plasmonic properties of the metal are for the fishnet structure from Ref. [22]. This is done by introducing the plasmonic parameter T_p [37], which is the number that characterizes the plasmonic nature of a metamaterial and it is defined as the ratio of the kinetic energy of the plasmonic electrons to the magnetic energy in the unit cell of a crystal. Strong plasmonic effects and the importance of electrostatic resonances imply T_p ≫ 1. We find that the plasmonic parameter for the mode labelled E_x in Fig. 4 at a frequency of 175THz is T_p = 0.24. Being less than one indicates this mode is not predominately plasmonic in nature because the electromagnetic fields do not significantly penetrate into the Ag. This is consistent with a recent demonstration of a negative index mode in a PEC fishnet structure [35, 36].

Finally, we describe another application of the CWES method: calculation of isofrequency ω(k) = const contours in metamaterial crystals. From the early days of metamaterials reseach, isofrequency diagrams have been used as a simple visual tool for studying refraction at the vacuum/metamaterial interfaces, especially in the context of negative index propagation and negative refraction [38]. Traditionally, isofrequency diagrams are drawn using a conventional ω(k) eigenvalue simulation [2]. This provides no information about propagation lengths which are important in lossy metamaterials such as the fishnet. The conventional approach is also highly laborious for plasmonic fishnets because of the strong frequency dependence of the dielectric permittivities of metals. Other semi-analytic techniques used for analyzing wave propagation in highly symmetric PEC-based fishnets [35] do not apply to plasmonic fishnets of interest.

Figure 6 shows an isofrequency diagram calculated using the CWES approach. A single isofrequency contour was obtain by fixing the real frequency ω, setting k ₀ = k_y ŷ and k̂ _n = ẑ, and then scanning k_y from –π/a_y to π/a_y. The eigenvalue in this simulation is a complex valued k_z. This procedure was repeated for several values of ω from the 150THz to 180THz. The resulting eigenmodes can be be excited by a plane wave incident on an interface between vacuum and the fishnet crystal parallel with the x̂-ŷ plane if the wavevector of the incident wave is confined to the ŷ-ẑ plane. Because the ŷ component of the incident wavevector is real valued, the ŷ component of the wave excited in the fishnet crystal must also be real-valued. k_y = sinθω/c where θ is the incidence angle with respect to the normal z. The eigenmode excited in the fishnet crystal decays in the ẑ direction as indicated by the imaginary part of k_z plotted in Fig. 6.

Fig. 6 (a) Isofrequency contours for vacuum plotted with respect to k_z and k_y (b) Isofrequency contours for the E_x polarized eigenmode of the fishnet crystal plotted with respect to Re(k_z) and k_y. (c) Isofrequency contours for the E_x polarized eigenmode of the fishnet crystal plotted with respect to Im(k_z) and k_y. The black arrows indicate the direction of the phase velocity and the red arrows indicate the direction of the group velocity at a frequency of 170THz and k_y = ω/csin(π/6).

Download Full Size | PDF

The isofrequency contours in Fig. 6 are hyperbolic in appearance. We can study refraction at the interface between vacuum and the fishnet crystal by comparing the isofrequency contours of the fishnet to those of vacuum, which are also shown in Fig. 6. The vacuum isofrequency contours are circular. As can be seen in Fig. 6 for frequency of 170THz and with an incident angle of 30° the component of the incident wavevector tangential to the interface (k_y) must be matched to the eigenmodes of the fishnet. For k_y = ω/csin(π/6) and a frequency of 170THz the isofrequency diagram shows two eigenmodes with the correct value of k_y and equal and opposite values of k_z. We find the correct mode by calculating the group velocity v_g = ∂ω/∂ Re(k) which is by definition normal to the isofrequency contours. In the absence of anomalous dispersion the group velocity indicates the direction of energy flow in the crystal [39]. Choosing the correct solution requires us to choose the solution that has energy flowing in the positive ẑ direction (i.e. away from the interface). This selects the solution with a negative Re(k_z). This is a negative index mode in the sense that in the ẑ direction the phase velocity and group velocity have opposite signs.

However, because the shape of the isofrequency contours is hyperbolic, the phase and group velocities in the ŷ direction have the same sign. Therefore, positive refraction at the x-–ŷ vacuum/fishnet interface is expected according to Fig. 6 despite the fact that a negative index eigen-mode is excited. This does not contradict the experiment by Valentine et al. in Ref. [22], where the fishnet structure was cut in the shape of a prism. In that experiment, the wave propagation through the fishnet was entirely in the ẑ direction enabling the group and phase velocities to be antiparallel. Negative refraction indeed occurs at a vacuum-fishnet interface tilted with respect to the principal axis of the crystal.

Finally, though we calculated the isofrequency contours in Fig. 6 by varying k_y and solving for k_z, we could have alternately varied the angle of k by setting k = λcos(θ)ẑ +λsin(θ)ŷ, varying the angle θ and solving for λ as the eigenvalue. This would have produced similar but slightly different isofrequency contours. This is an example of how the complex k dispersion curves in this paper are not unique and depend on how one restricts the three degrees of freedom k_x, k_y and k_z to one degree of freedom λ.

5. Conclusion

In conclusion, we have presented a three-dimensional realization of the Complex Wavenumber Eigenvalue Simulation (CWES) approach to calculating photonic band structure of metamaterial/photonic crystals. A detailed implementation of CWES using FEM discretization is described. The CWES approach was applied to two periodic photonic structures: (a) photonic crystal comprised of dielectric spheres, and (b) plasmonic fishnet metamaterial supporting negative refractive index waves. For case (a) we have used the results of CWES to identify both transverse and longitudinal modes and investigated their coupling that gives rise to avoided crossings and bandgap formation. For case (b) we have identified, for the first time, four negative-index modes of the fishnet structure and computed hyperbolic isofrequency surfaces for the least-damped transverse mode. The sphere photonic crystal model is available for [21].

Acknowledgments

We would like to acknowledge the developers of MUMPS [40], ARPACK [41,42], PETSC [43] and SLEPC [44]. Support of the staff of the Texas Advanced Computing Center is gratefully acknowledged. This work was supported by the Air Force Research Laboratory and the Office of Naval Research.

References and links

2. J. D. Joannopoulos, R. D. Meade, and J. N. Winn, Photonic Crystals: Molding the Flow of Light, 2nd ed. (Princeton University Press, 2008).

3. K. Sakoda, Optical Properties of Photonic Crystals, 2nd ed. (Springer, 2004).

4. G. Shvets and Y. A. Urzhumov, “Engineering the electromagnetic properties of periodic nanostructures using electrostatic resonances,” Phys. Rev. Lett. 93, 243902 (2004). [CrossRef]

5. E. D. Palik, Handbook of Optical Constants of Solids (Acedemic Press, 1985), vol. 1.

6. S. Zhang, W. Fan, K. J. Malloy, S. R. J. Brueck, N. C. Panoiu, and R. M. Osgood, “Near-infrared double negative metamaterials,” Opt. Express 613, 4922 (2005). [CrossRef]

7. M. Notomi, “Theory of light propagation in strongly modulated photonic crystals: refractionlike behavior in the vicinity of the photonic band gap,” Phys. Rev. B 62, 10696 (2000). [CrossRef]

8. R. C. McPhedran, L. C. Botten, J. McOrist, A. A. Asatryan, C. M. de Sterke, and N. A. Nicorovici, “Density of states functions for photonic crystals,” Phys. Rev. E 69, 016609 (2004). [CrossRef]

9. T. Suzuki and P. L. Yu, “Tunneling in photonic band strucures,” J. Opt. Soc. Am. B 12, 804 (1995). [CrossRef]

10. E. Istrate, A. A. Green, and E. H. Sargent, “Behavior of light at photonic crystal interfaces,” Phys. Rev. B 71, 195122 (2005).

11. A. J. Ward, J. B. Pendry, and W. J. Stewart, “Photonic dispersion surfaces,” J. Phys.: Condens. Matter 7, 2217–2224 (1995). [CrossRef]

12. Z. Y. Li and L. L. Lin, “Photonic band structure solved by a plane-wave–based transfer-matrix method,” Phys. Rev. E 67, 046607 (2003). [CrossRef]

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

14. C. Engström, C. Hafner, and K. Schmidt, “Computations of lossy bloch waves in two-dimensional photonic crystals,” J. Comput. Theor. Nanosci. 6, 1–9 (2009). [CrossRef]

15. X. Zhang, M. Davanco, Y. Urzhumov, and G. Shvets, “A subwavelength near-infrared negative index material,” Appl. Phys. Lett. 94, 131107 (2009). [CrossRef]

16. X. Zhang, M. Davanco, K. Miller, T. W. J. C. Wu, C. Fietz, D. Korobkin, X. Li, G. Shvets, and S. R. Forrest, “Interferometric characterization of a subwavelength near-infrared negative index metamaterial,” Opt. Express 18, 17788–17795 (2010). [CrossRef] [PubMed]

17. C. Fietz and G. Shvets, “Current-driven metamaterial homogenization,” Phys. B 405, 2930–2934 (2010). [CrossRef]

18. M. Conforti and M. Guasoni, “Dispersive properties of linear chains of lossy metal nanoparticles,” J. Opt. Soc. Am. B 27, 1576–1582 (2010). [CrossRef]

19. C. Fietz and G. Shvets, “Homogenization theory for simple metamaterials modeled as one-dimensional arrays of thin polarizable sheets,” Phys. Rev. B 82, 205128 (2010). [CrossRef]

20. A. Pors, I. Tsukerman, and S. I. Bozhevolnyi, “Effective constitutive parameters of plasmonic metamaterials: a rigorous approach,” arXiv:1104.2972v1(2011).

21. Both .m and .mph files for the sphere photonic crystal model can be acquired by contacting Chris Fietz at fietz.chris@gmail.com.

22. J. Valentine, S. Zhang, T. Zentgraf, E. Ulin-Avila, D. A. Genov, and G. Bartal, “Three-dimensional optical metamaterial with a negative refractive index,” Science 455, 376–380 (2008).

23. W. B. J. Zimmerman, Process Modelling and Simulation with Finite Element Methods (World Scientific, 2004).

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

25. P. P. Silvester and R. L. Ferrari, Finite Elements for Electrical Engineers3rd ed. (Cambridge University Press, 1996).

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

27. C. L. Holloway, E. F. Kuester, J. Baker-Jarvis, and P. Kabos, “A double negative (DNG) composite medium composed of magnetodielectric spherical particles embedded in a matrix,” IEEE Trans. Antennas Propag. 51, 2596–2603 (2003). [CrossRef]

28. R. A. Shore and A. D. Yaghjian, “Traveling waves on two- and three-dimensional periodic arrays of lossless scatterers,” Radio Sci. 42, RS6S21 (2007). [CrossRef]

29. C. R. Simovski, “Bloch material parameters of magneto-dielectric metamaterials and the concept of Bloch lattices,” Metamaterials 1, 62–80 (2007). [CrossRef]

30. A. Alú, “First-principle homogenization theory for periodic metamaterial arrays,” arXiv:1012.1351v2(2011).

31. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (John Wiley & Sons, 1998). [CrossRef]

32. J. M. Siqueiros, R. Machorro, and L. E. Regalado, “Determination of the optical constants of MgF2 and ZnS from spectrophotometric measurements and the classical oscillator method,” Appl. Opt. 27, 2549–2553 (1988). [CrossRef] [PubMed]

33. P. Halevi, Spatial Dispersion in Solids and Plasmas (Elsevier Science Publishers, 1992).

34. M. A. Shapiro, G. Shvets, J. R. Sirigiri, and R. J. Temkin, “Spatial dispersion in metamaterials with negative dielectric permittivity and its effect on surface waves,” Opt. Lett. 31, 2051 (2006). [CrossRef] [PubMed]

35. L. Jelinek, R. Marques, and J. Machac, “Fishnet metamaterials - rules for refraction and limits of homogenization,” Opt. Express 18, 17940–17949 (2010). [CrossRef] [PubMed]

36. M. Navarro-Ca, M. Beruete, M. Sorolla, and I. Campillo, “Negative refraction in a prism made of stacked subwavelength hole arrays,” Opt. Express 16, 560–566 (2008). [CrossRef]

37. Y. A. Urzhumov and G. Shvets, “Optical magnetism and negative refraction in plasmonic metamaterials,” Solid State Commun. 146, 208–220 (2008). [CrossRef]

38. C. Luo, S. G. Johnson, J. D. Joannopoulos, and J. B. Pendry, “All-angle negative refraction without negative effective index,” Phys. Rev. B 65, 201104 (2002).

39. L. Brillouin, Wave Propagation and Group Velocity (Academic Press, 1960).

40. MUltifrontal Massively Parallel Solver (MUMPS 4.9.2) User’s guide (Lyon, 2009).

41. R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK Users’ Guide, Solution of Large-Scale Eigenvalue Problems by Implicitly Restarted Arnoldi Methods (SIAM, Philadelphia). [PubMed] [PubMed]

42. K. J. Maschhoff and D. C. Soerensen, “PARPACK: An efficient portable large scale eigenvalue package for distributed memory parallel architectures,” Lect. Notes Comp. Sci. 1184, 478–486 (1996). [CrossRef]

43. S. Balay, J. Brown, K. Buschelman, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, B. F. Smith, and H. Zhang, “PETSc users manual,” Tech. Rep. ANL-95/11 - Revision 3.1, Argonne National Laboratory (2010).

44. J. E. Roman, E. Romero, and A. Tomas, “SLEPc users manual,” Tech. Rep. DSIC-II/24/02 - Revision 3.1, D. Sistemas Informáticos y Computación, Universidad Politécnica de Valencia (2010).

Complex k band diagrams of 3D metamaterial/photonic crystals.

Abstract

1. Introduction

2. The finite element eigenvalue problem

2.1. The field equation

2.2. The finite element model

2.3. Electric field formulation

3. Example: dielectric photonic crystal

4. Example: negative index plasmonic fishnet

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (6)

Equations (12)

Optics Express