Complex k band diagrams of 3D metamaterial/photonic crystals

A finite element method (FEM) for solving the complex valued k({\omega}) vs. {\omega} dispersion curve of a 3D metamaterial/photonic crystal system is presented. This 3D method is a generalization of a previously reported 2D eigenvalue method. This method is particularly convenient for analyzing periodic systems containing dispersive (e.g., plasmonic) materials, for computing isofrequency surfaces in the k-space, and for calculating the decay length of the evanescent waves. Two specific examples are considered: a photonic crystal comprised of dielectric spheres and a plasmonic fishnet structure. Hybridization and avoided crossings between Mie resonances and propagating modes are numerically demonstrated. Negative index propagation of four electromagnetic modes distinguished by their symmetry is predicted for the plasmonic fishnets. By calculating the isofrequency contours, we also demonstrate that the fishnet structure is a hyperbolic medium.


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 [3,4] 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 four such instances can be identified. First, metamaterials often contain dispersive materials such as metals, where the dielectric function strongly depends on the frequency on the wave. In this case, the eigenfrequency problem becomes a nonlinear eigenvalue problem and must be solved iteratively [5]. In contrast, when solving for the wavenumber as a function of frequency, the resulting eigenvalue problem only needs to be solved once. Second, 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. Third, 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. Fourth, 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 [1,2] but only for 2D crystals. The benefits of the complex wavenumber 2D FEM are becoming better appreciated and use of this method is becoming more common [13,14,15,16,17]. 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 [18].
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. In Sec. 4 we calculate the dispersion curves for a negative index fishnet metamaterial [19] 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.

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 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 [3,4] the magnetic field can be represented as the product of a periodic function and an exponential factor 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 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

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 COM-SOL. In what follows, only the most essential features of the FEM approach are reviewed; more detailed treatments can be found elsewhere [20,21,22]. 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 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 [20,21], wheren 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 expressionn × (−ik × u + ∇ × u) /ε to be equal to zero. Recalling Eq. (4) we note that this simply enforces the boundary conditionn × 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 orn × 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 [20] 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 0 + λk n where λ will be the eigenvalue solved for, k 0 is an offset vector andk n is a unit vector (k n ·k n = 1) that defines the direction of the complex wavenumber eigenvalue λ . The FEM turns the weak form and accompanying boundary conditions into an algebraic problem, in this case a quadratic eigenvalue problem [23]: 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 When using COMSOL to solve the FEM problem, this linearization is performed automatically during the solution stage.

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 As before, we replace the electric field with a periodic vector field times an exponential factor producing the new field equation The corresponding weak form for this field equation is 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 conditionn × H = 0. Thus the PMC boundary condition is the natural boundary condition when solving for the electric field.

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. 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 = 0 andk 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. 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 thex direction making it a magnetic bulk plasmon. The longitudinal mode with the passband near ω = 5c/a is electrically polarized in thex 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 [24].
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. [1,23] 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 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 0 = ω/c sin(π/6)ŷ andk n =x and solve the resulting eigenvalue problem. The resulting photonic band structure is plotted in Fig. 2.
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 [23,1].
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.

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 [19] demonstrated that the so-called fishnet metamaterial supports a negative index eigenmode for near infrared wavelengths of about λ 0 ≈ 1.7µm. The dimensions and composition of the unit cell taken from Ref. [19] are shown in Fig. 3. The fishnet metamaterial is made of alternating layers of Ag and MgF 2 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 ε Ag = 1 − ω 2 p /(ω(ω − iγ)) with the same parameters cited in Ref. [19]: ω p = 9eV and γ = 0.054eV. It should be noted that in Ref. [19] 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 2 , no values were provided in the original paper so we used an oscillator model of the dielectric permittivity from Ref. [25].
The negative index mode reported in Ref. [19] propagates in theẑ direction. Therefore, we have applied CWES to the specific case of k 0 = 0 andk n =ẑ. The resulting dispersion curves Fig. 3. The fishnet metamaterial from Ref. [19]. 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 2 with a thickness of 50nm. The widths of the the crisscrossing fishnet strips are b 1 = 265nm and b 2 = 565nm. are plotted in Fig. 4 for the 100T Hz < ω/2π < 300T Hz 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. [19] is labeled E x to indicate that it is transversely polarized with the electric field pointing in thex 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 200T Hz. 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 thê y direction and the mode labeled E z is longitudinally polarized with the electric field pointing in theẑ direction.
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. [19] 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 150T Hz < ω/2π < 200T Hz frequency range. Somewhat unexpectedly, the mode experimentally observed in Ref. [19] (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 [19], which could be attributed to fabrication imperfections.
The existence of multiple negative index waves in the fishnet structure has important experimental implications. For some frequencies (e.g., at 150T Hz) 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 [19] only excited the E x mode by using the light normally incident onto the  Fig. 3. In addition to the transverse mode electrically polarized in thê x direction labeled E x which was identified in Ref. [19] 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 2 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 2 strips. Arrows represent in-plane electric field and color represents the H y field. For all field profiles the frequency is 175T Hz and each region is labeled Ag or MgF 2 according the the material of the region. Unlabeled regions are vacuum.
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 [19] would give rise to two distinct beams. The existence of the additional strongly  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 220T Hz indicates that it is no longer a negative index mode at these higher frequencies.
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 [26,27]. Strong spatial dispersion is not unexpected because the fishnet's transverse period is about λ /2. Recent theoretical [28] and experimental [29] 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. [19]. This is done by introducing the plasmonic parameter T p [30], 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 175T Hz 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 [29,28].
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 [31]. Traditionally, isofrequency diagrams are drawn using a conventional ω(k) eigenvalue simulation [3]. 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 [28] 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 0 = k yŷ andk 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 150T Hz to 180T Hz. The resulting eigenmodes can be be excited by a plane wave incident on an interface between vacuum and the fishnet crystal parallel with thex-ŷ 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.
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 170T Hz 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 = ω/c sin(π/6) and a frequency of 170T Hz 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 [32]. 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 thex-ŷ vacuum/fishnet interface is expected according to Fig. 6 despite the fact that a negative index eigenmode is excited. This does not contradict the experiment by Valentine et al. in Ref. [19], 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.

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.