Photonic band structures of periodic arrays of pores in a metallic host: tight-binding beyond the quasistatic approximation

We have calculated the photonic band structures of metallic inverse opals and of periodic linear chains of spherical pores in a metallic host, below a plasma frequency $\omega_{\text{p}}$. In both cases, we use a tight-binding approximation, assuming a Drude dielectric function for the metallic component, but without making the quasistatic approximation. The tight-binding modes are linear combinations of the single-cavity transverse magnetic (TM) modes. For the inverse-opal structures, the lowest modes are analogous to those constructed from the three degenerate atomic p-states in fcc crystals. For the linear chains, in the limit of small spheres compared to a wavelength, the results are the"inverse"of the dispersion relation for metal spheres in an insulating host, as calculated by Brongersma \textit{et al.} [Phys.\ Rev.\ B \textbf{62}, R16356 (2000)]. Because the electromagnetic fields of these modes decay exponentially in the metal, there are no radiative losses, in contrast to the case of arrays of metallic spheres in air. We suggest that this tight-binding approach to photonic band structures of such metallic inverse materials may be a useful approach for studying photonic crystals containing metallic components, even beyond the quasistatic approximation.


I. INTRODUCTION
The photonic band structures of composite materials have been studied extensively. Such band structures are defined by the relation between frequency ω and Bloch vector k in media in which the dielectric constant is a periodic function of position. A major reason for such interest is the possibility of producing photonic band gaps, i.e., frequency regions, extending through all k-space, where electromagnetic waves cannot propagate through the medium. Such media have many potentially valuable applications, including possible use as filters and in films with rejection-wavelength tuning. [1] In systems with a complete photonic band gap, the spontaneous emission of atoms with level splitting within the gap can be strongly suppressed. [2] Since light cannot travel through the photonic band gap materials (Bragg diffracted backwards), one of their applications can be a complete control over wasteful spontaneous emission in unwanted directions when a device, such as a laser, is embedded inside a 3D photonic crystal. [3] 2D photonic crystals can be used as optical microcavities, microresonators, [4] waveguides, [5] lasers, [6] or fibers [7] while 1D photonic crystals can be used as Bragg gratings or optical switches. [8] The photonic band structure of a range of materials has been studied using a plane wave expansion method. Typically, the method converges easily when the dielectric function is everywhere real, but more slowly, or not at all, when the dielectric function has a negative real part, as occurs when one component is metallic. For example, McGurn et al. [9] used this method to calculate the photonic band structure of a square lattice of metal cylinders in two dimensions (2D) and of an fcc lattice of metal spheres embedded in vacuum in 3D. They found that that method converged well when the filling fraction f (i.e., volume fraction of metal spheres or cylinders) sat-isfied f ≤ 0.1%.
Kuzmiak et al. [10] used the same method to calculate the photonic band structures for 2D metal cylinders in a square or triangular lattice in vacuum. For low f and ω > ω p , the calculated photonic band structures are just slightly perturbed versions of the dispersion curves for electromagnetic waves in vacuum. However, for ω < ω p and H-polarized waves (magnetic field H parallel to the cylinders), they obtained many nearly flat bands for ω < ω p ; these bands were found to converge very slowly with increasing numbers of plane waves. They later extended this work to systems with dissipation. [11] To describe dispersive and absorptive materials, they used a complex, position-dependent form of dielectric function. They also introduced a standard linearization technique to solve the resulting nonlinear eigenvalue problem.
Zabel et al. [12] extended the plane wave method to treat periodic composites with anisotropic dielectric functions. In particular, they studied the photonic band structures of a periodic array of anisotropic dielectric spheres embedded in air. They found that the anisotropy split degenerate bands, and narrowed or even closed the band gaps. Much further work on anisotropic photonic materials has been carried out since this paper (see, e.g., Ref. [2]).
A different type of periodic metal-insulator composite is a periodic arrangement of metallic spheres in an insulating host. Brongersma et al. [13] studied the dispersion relation for coupled plasmon modes in such a linear chain of equally spaced metal nanoparticles, using a near-field electromagnetic (EM) interaction between the particles in the dipole limit. They also studied the transport of EM energy around the corners and through tee junctions of the nanoparticle chain-array.
Park and Stroud [14] also studied the surface-plasmon dispersion relations for a chain of metallic nanoparticles in an isotropic medium. They used a generalized tight-binding calculations, including all multipoles. This approach is more exact than the previous point-dipole calculation, [13] in a quasistatic limit, but still leaves out non-quasistatic effects associated with radiative damping (i.e., effects associated with the non-vanishing of ∇ × E, where E is the electric field. They calculated the lowest bands as well as many higher bands and compared their results with those in Ref. [13].
Weber and Ford [15] have shown that all calculations within the quasistatic approximation omit important interactions between transverse plasmon waves and free photon modes, even if the interparticle separation is small compared to the wavelength of light. Thus, most quasistatic calculations need to have certain corrections included at particular values of the wave vector.
Recently, Gaillot et al. [16] have studied the photonic band structures of another type of structure, a so-called inverse opal structure. This structure is an fcc lattice of void spheres in a host of another material. Such a structure can be prepared, e.g., starting from an opal structure made of spheres of a convenient substance, infiltrating it with another material, then dissolving away the spheres. In the work of Ref. [16], the photonic band structure of Si inverse opal was calculated as a function of the infiltrated volume fraction f of air voids using three-dimensional finite difference time domain (3D FDTD) method. It was found that for certain values of f , a complete band gap opens up between the eighth and ninth bands.
In the present work, first we study the photonic band structure of an inverse opal structure, such as that investigated in Ref. [16], but instead of dielectric materials such as Si, we consider metals as the infiltrated materials. Thus, the material we study is also the inverse of the fcc array of metal spheres studied by McGurn et al. [9] Such metallic inverse opal structures have recently become of great interest, because it has been found that Pb inverse opals exhibit superconductivity. [17] These workers have studied the response of these materials to an applied magnetic field, and have found a highly non-monotonic fractional flux penetration into the Pb spheres as a function of the applied field.
As a second example, we study the photonic band structure of a linear chain of nanopores in a metallic medium. This is an inverse structure of a linear chain of metallic nanospheres, of which the dispersion relation is given in Ref. [13]. As anticipated, we get a kind of "inverse image" of the dispersion relation found by Ref. [13] in our system.
For both types of structures, our primary method for studying the photonic band structures below the plasma frequency ω p is a tight-binding approximation which is valid even in the non-quasistatic regime. Because the analogs of the tight-binding atomic states decay exponentially in the metallic host medium, the resulting tight-binding waves do not lose energy radiatively, as do the corresponding waves along one-dimensional chains of metallic nanoparticles in air. Furthermore, because the modes are expanded in "atomic" states rather than plane waves, there is no convergence problem as there can be in the plane wave case.
The remainder of this paper is organized as follows. In Section II, we first present the formalism for calculating the transverse magnetic (TM) and transverse electric (TE) modes of a single spherical cavity in a metallic host. We then describe the method for calculating the photonic band structures of metallic inverse opals and of linear chains of nanopores in a metallic host, using a simple tight-binding approach for ω < ω p . In Section III, we give the numerical results for the TM and TE modes of a single cavity and those of the tight-binding method for the metal inverse opals and the linear chain of nanopores. Section IV presents a summary and discussion.

II. FORMALISM
In this section, we present a summary of the equations determining the band structure of a photonic crystal containing a metallic component with Drude dielectric function ǫ(ω) = 1 − ω 2 p /ω 2 and an insulating component of dielectric constant unity. The insulating component is assumed to be present in the form of identical spherical cavities of radius R. We first write down the equations for the TM and TE modes of a spherical cavity in a Drude metal. Then, we present a tight-binding method for ω < ω p .

A. Spherical Cavity
As a preliminary to calculating the photonic band structure, we first discuss the modes of a single spherical cavity in a Drude metal host. We begin with the TM modes of the cavity, then the TE modes.

TM Modes
It is convenient to describe the modes of the embedded cavity in terms of the B field. To that end, we combine the two homogeneous Maxwell equations to obtain a single equation for B: Here, we have expressed the position-and frequencydependent dielectric function ǫ(x, ω) as 1/ǫ(x, ω) = θ(x)/(1 − ω 2 p /ω 2 ) + 1 − θ(x), where the step function θ(x) = 1 inside the metallic region and θ(x) = 0 elsewhere. Multiplying this equation by ω 2 − ω 2 p and simplifying, we obtain Thus, inside the spherical void, we have while inside the metal, For a spherical void within a metallic host, it is convenient to solve these equations in spherical coordinates. It is readily found that, for the TM modes, the nonvanishing components of the solutions for B and E for Eq. (5) and Eq. (2) are [18] where k = ω/c, j ℓ and n ℓ are spherical Bessel functions, and the subscripts φ, r, and θ denote components of the corresponding fields in spherical coordinates. Likewise, the solutions of Eqs. (6) and (2) within the metal will be where k ′ = ω 2 − ω 2 p /c. The requirements that the normal displacement D and tangential E should be continuous at R gives the two conditions and ∂u ℓ (r) where R is the radius of the spherical cavity. Since the fields at the center of the void sphere must be finite, we also have where we have normalized the solution so that the coefficient A ℓ = 1. From Eq. (9) we have whereas, from Eq. (10), we get The coefficients C ℓ and D ℓ can then be determined from these boundary conditions. The results for ω < ω p , can be obtained by making the substitution k ′ → ik ′ , with k ′ real. In this case, the radial component of Eq. (6) takes the form which is the modified Bessel equation. The solutions of Eq. (14) are the modified spherical Bessel functions i ℓ and k ℓ (note that this k ℓ is different from the wave vectors k and k ′ ). In this case, Eq. (8) becomes where k ′ = ω 2 p − ω 2 /c. In addition, Eqs. (10), (12), and (13) are transformed, respectively, into and It is of interest to consider the specific case of a spherical cavity in an infinite medium. In this case, C ℓ = 0 because i ℓ (x) diverges at large x. As a result, Eqs. (15), (17), and (18) and From Eq. (20) we have D ℓ = j ℓ (kR)/k ℓ (k ′ R), and hence Eq. (21) becomes We can readily obtain the asymptotic forms of the solutions when kR ≪ 1 and k ′ R ≪ 1. In this case j ℓ (kR) and In this limit, Eq. (22), after some algebra, reduces to simply Since The largest value, ω = 2/3ω p , occurs at ℓ = 1 and the limiting value for large ℓ is ω = ω p / √ 2.

TE Modes
For the TE mode, inside the spherical void, we have whereas inside the metal, we have We now use these equations to calculate E φ , B r , and B θ . From Eq. (25) we get where k = ω/c. The solutions of Eq. (26) are where k ′ = ω 2 − ω 2 p /c.
Since normal B and tangential H should be continuous on the boundaries, we obtain the conditions as in the TM case, and Since the fields must be finite at the center of the void sphere, we can choose we also take the coefficient while, from Eq. (30), we get The corresponding equation for ω < ω p , can again be obtained by the transformation k ′ → ik ′ . The TE modes for ω < ω p using the modified spherical Bessel functions i ℓ (x) and k ℓ (x) satisfy Eqs. (14), (15), (30), and (17). The only changes are in Eq. (18), which becomes In an infinite medium, these conditions become, from Eq. , (35) If we consider the asymptotic forms of the solutions when kR ≪ 1 and k ′ R ≪ 1 as we did for the TM modes, Eq. (35) simplifies to which gives ℓ = −1/2. Since ℓ must be a positive integer, we see that there are no eigenvalues for TE modes in the limit kR ≪ 1 and k ′ R ≪ 1.

B. Tight-Binding Approach to Modes for ω < ωp
We now turn from describing the single-cavity modes to a discussion of the band structure for a periodic array of such cavities. In conventional periodic solids, the tightbinding method is very useful in treating narrow bands.
In what follows, we try to suggest an analogous tightbinding approach for the lowest set of TM modes in a periodic lattice of spherical cavities in a metallic host, in the frequency range ω < ω p . We apply the resulting method, first, to an fcc lattice of pores, and then to a linear chain of spherical pores in a metallic host.
Even though these are TM modes, it is convenient to describe them now in terms of their electric fields. We denote the electric field of the λth mode by E λ (x). This field satisfies is the "Hamiltonian" of this system. Since O is a Hermitian operator, the eigenstates corresponding to unequal eigenvalues ω 2 λ /c 2 and ω 2 µ /c 2 are orthogonal and may be chosen to be orthonormal. (The orthogonality may also be proved directly by integration by parts.) The orthonormality relation is Since E λ (x) is real for ω < ω p , the complex conjugation is, in fact, unnecessary. In Sec. II A 1, our paper already gives the equations determining the electric and magnetic fields of isolated TM modes for a spherical cavity. The lowest set corresponds to ℓ = 1, and there should be three of these. For a spherical cavity, all three are degenerate, i.e., all three have the same eigenfrequencies. Even though the three modes have equal frequencies, one can always choose an orthonormal set, with electric fields E 1 , E 2 , and E 3 satisfying the orthonormality relation in Eq. (38).
In order to obtain the tight-binding band structure built from these three modes, we need to calculate matrix elements of the form corresponding to two single-cavity modes associated with different cavities centered at the origin and at R. Here, O is the "Hamiltonian" of the system as defined implicitly in Eq. (37). Next, we introduce normalized Bloch states associated with the three ℓ = 1 single-cavity modes. In order to do this, we first make the usual tight-binding assumption that the "atomic" states corresponding to different cavities are orthogonal: This orthogonality of states on different cavities is reasonable since the fields fall off exponentially with separation. The orthonormal Bloch states then take the form where k is a Bloch vector, and the R's are the Bravais lattice vectors. In writing Eq. (41), we have assumed that there are N identical spherical cavities, and that the Bloch states satisfy the usual periodic boundary conditions of Born-von Karman type. We also introduce the elements of the "Hamiltonian" matrix We can then obtain the frequencies ω(k) by diagonalizing a 3 × 3 matrix as follows: where ω at is the eigenvalue of a single-cavity mode. The solutions to these equations give the three p-bands for a periodic lattice of cavities in a metallic host. This procedure is analogous to that used in the well-known procedure for obtaining tight-binding bands from three degenerate p-bands in the electronic structure of conventional solids (see, for example, Ref. [19]). We briefly comment on the connection between this approach and that used by earlier workers. [13,14] In this work, the authors treat wave propagation along a chain of metallic nanoparticles. They use the tight-binding approximation, as we do, but in the quasistatic approximation in which one assumes that ∇ × E = 0. This approximation is reasonable when both the particle radii and the interparticle separations are small compared to a wavelength, but is not accurate in other circumstances. Furthermore, even in the small-particle and small-separation regime, this approximation still fails to account for the radiation which occurs at certain wave numbers and frequencies. The present approach would generalize this tight-binding method to (a) three dimensions as well as one; (b) pore modes instead of small particle modes; and most importantly (c) larger pores and larger interparticle separations, via extension beyond the quasistatic approximation.
Next, we discuss the numerical evaluation of the required matrix elements, Eq. (39). The relevant electric fields are given in this paper, but in spherical coordinates. It is not difficult to convert these into Cartesian coordinates. The operator O is just a little trickier. We first since E β is an eigenstate of O R with an eigenvalue ω 2 at /c 2 . But since we are assuming that the overlap integral between "atomic" electric field states centered on different sites vanishes, the term involving O R does not contribute to the matrix element M α,β , which is therefore just given by We can also write where is a step function which is unity inside the cavity centered at R ′ and is zero otherwise. A reasonable approximation to Eq. (46) might be to include just R ′ = 0. In this case, we finally will get where the integral runs just over the cavity centered at the origin. As a further approximation, we can just replace E β (x−R) by the value of this function at the origin, i.e., E β (−R). Then this field can be taken outside the integral and we just have where once again the integral runs over the cavity centered at the origin. Next, we attempt to calculate the relevant quantities needed to solve for this matrix element. In order to use the tight-binding approach we will need to normalize the individual eigenstates E α . Therefore, we will begin by obtaining this normalization. For ℓ = 1, the u ℓ (r)'s are r times spherical Bessel functions. We write this field as where we have used the relation P 1 (cos θ) = cos θ, and introduced the normalization constant C 1 , which will be determined below. Similarly, where we use P 1 1 (cos θ) = − sin θ. For r > R, we have and We will need the integrals of the Cartesian components of the field over the volume of the sphere centered at the origin. Let us assume we are considering the z mode, i.e., the one for which θ refers to the angle from the z axis. Then the symmetry of the problem shows that only the z component of the electric field will have a nonzero integral. Also, we have that E z,in (r, θ) = E r,in cos θ − E θ,in sin θ. (54) Thus, after a little algebra, we find that the integral of this field over the volume of the cavity is Next, we work out the coefficient D 1 . It is determined by the boundary conditions at r = R. These conditions are that D r and E θ should be continuous at r = R. These two conditions determine not only the value of D 1 but also the allowed frequency. After a bit of algebra, we find that The allowed value of ω is given by Eq. (22). Finally, we need the normalization constant C 1 . We choose this so that the integral of the square of the electric field for a single-cavity mode should be normalized to unity. This condition may be written If we write and then we can express the normalization condition as Therefore, we can now write out an explicit expression for the matrix element M α,β (R) given in Eq. (49). For the αth mode, the integral of E α over the volume of a cavity is a vector in the αth direction. To evaluate Eq. (49), we need the component of the αth mode in the βth direction at a position R. Let us first consider the zth mode (α = z). We can use Eqs. (52) and (53) to rewrite this field in Cartesian coordinates with the additional equations cos θ = z r , We just substitute these expressions back into Eqs. (52) and (53) to get the Cartesian components of the field for a mode parallel to the z axis. For the mode parallel to the x axis, we just permute the coordinates cyclically: z → x, x → y, and y → z. Similarly, for the y modes, we make the permutation (x, y, z) → (z, x, y).
Using these results, we should be able to compute all the elements in the tight-binding matrix and hence obtain the band structure for the photonic p-bands in the tightbinding approximation, in either one or three dimensions.

III. NUMERICAL RESULTS
For the inverse opals we arbitrarily assume a lattice constant d = 500 √ 2 nm, and a void sphere radius R = 150 nm as in Fig. 1(a). This choice is the same as that of Ref. [17], where the Pb inverse opal has this lattice constant. Since the volume of the primitive unit cell is v c = d 3 /4, this corresponds to a void volume fraction f = 0.160. For the linear chain of nanopores (see below) this d is the separation between two nanopores and R is the radius of a nanopore as in Fig. 1(b).
The metallic dielectric functions we assume for the inverse opals and linear chain of nanopores are of the usual Drude form, where ω p is the plasma frequency of the conduction electrons. ǫ(ω) < 0 when ω < ω p , while ǫ(ω) > 0 when ω > ω p . Our calculations are thus carried out assuming that the Drude relaxation time τ → ∞. For a metal in its normal state, ω 2 p = 4πne 2 /m, where n is the conduction electron density and m is the electron mass. Note that with this choice of dielectric function, the entire band structure can be expressed in scaled form. That is, the scaled frequency ωd/c is a function only of the scaled wave vector kd, and the band structures are parameterized by the two constants ω p d/c and f for the case of inverse opals.
Since we are considering void spheres in inverse opals and linear chains of nanopores, it is of interest to consider electromagnetic wave modes in a single cavity, which could be considered a single "atom" of the void lattice. We show only results for ω < ω p , since these are the results most relevant to possible narrow-band photonic states in the inverse opal structure. Our results for ω < ω p for an isolated spherical cavity in an infinite medium, and when kR ≪ 1 and k ′ R ≪ 1 are given in Table I. These two inequalities are reasonable for our inverse opal system parameters d = 500 √ 2nm, R = 150nm, and ω p d/c = 1, because The (modified) spherical Bessel functions in Eq. (22) are extremely close to the ω axis for ℓ > 5, so that it is difficult to get eigenfrequencies for ℓ > 5 in the isolated spherical cavity. However the eigenfrequencies continue to exist even for ℓ > 5 when kR ≪ 1 and k ′ R ≪ 1.  The solutions to Eq. (35) do not exist for ω < ω p with ω p d/c = 1. This fact is consistent with that the eigenvalues for ω < ω p do not exist for TE modes when kR ≪ 1 and k ′ R ≪ 1.
Assuming ω p d/c = 1 and using ω at d/(2πc) = 0.1296 for ℓ = 1 in an infinite medium, we get the tight-binding results in Fig. 2. This figure shows three separate bands in the X-U -L region and X-W -K region as expected for the p-bands. The bandwidth is relatively small as M α,β (R)d 2 ∼ 0.001, which proves the general relation between the bandwidth and the overlap integral. [19] All three bands are degenerate at k = 0 (the Γ point). In addition, there is a double degeneracy when k is directed along either a cube axis (Γ-X) or a cube body diagonal (Γ-L), the higher (concave upward) bands being degenerate in both cases. The lower two bands have a band gap at the U point, and these bands cross at the W point.
Next, we turn to the band structure of a periodic linear chain of spherical nanopores in a Drude metal host. For this linear chain, the Bravais lattice vectors are R = d(0, 0, ±n), where ±n is the nth nearest-neighbor, d is the separation between two nanopores and we assume that the chain is directed along the z axis. We can calculate the tight-binding band structure including as many sets of neighbors ±n as we wish. To compare our results with those in Ref. [13], we first use their parameters, R = 25 nm and d = 75 nm, together with their overlap parameter ω 1 = 1.4 × 10 15 rad/s. These combine to give ω p d/c = 0.35. The "atomic" frequency is found by solving Eq. (22) and gives ω at d/(2πc) = 0.0454 for ℓ = 1 in an infinite medium. Our resulting tightbinding dispersion relations are shown in Fig. 3 with only nearest-neighbors included. Note that our frequencies are given in unit of 2πc/d while the results of Ref. [13] are not scaled. Our results are exactly the inverse images of theirs -that is, we would get their curves (to within a constant of proportionality) if we reflect our curves through the horizontal line of the atomic level, and the transverse (T) branches are twofold degenerate as are theirs, while the longitudinal (L) branch is nondegenerate. Our eigenfrequency for a single-cavity ω at corresponds to their resonance frequency ω 0 . As we increase the number of nearest-neighbors (nn's) included, the separation between the L and T branches increases at the zone center but decreases at the zone boundary, as shown in Fig. 4; the same trend is seen in Fig. 1 in Ref. [13]. The sum also converges quickly, so there is little difference between the dispersion relation including through the next-nearest-neighbors and that including through the 5th nearest-neighbors.
We have carried out similar calculations using other values of the parameter ω p d/c, namely 1.0, 2.0, and 5.0. Such calculations are possible here because our calculations are non-quasistatic, so that the overlap integral between neighboring spheres falls off exponentially with separation. The results are given in Figs. 5, 7, and 9, respectively. The corresponding results including more overlap integrals are shown in Figs. 6, 8, and 10, respectively. It is also striking that, as ω p d/c increases in going from Fig. 3 to Figs. 5, 7, and 9, the ratio r LT of the width of the L band to that of the T band steadily decreases. In Fig. 3, r LT > 1, in Fig. 9, r LT < 1, while in Fig. 7 (for which ω p d/c = 2.0), r LT ∼ 1.
One could also say that, except for an overall scale factor, Fig. 9 looks like an inverted image of Fig. 3 about the horizontal line of ω at . The dispersion relations for the intermediate value ω p d/c = 2.0 has nearly perfect symmetry about the horizontal line of ω at (the T branches are nearly reflections of the L branch about the horizontal line of ω at ), as in Fig. 7. For the nn case, the T and L bands cross at ±π/(2d), as can be seen in Figs. 3, 5, 7, and 9. When further neighbors are included, they cross at smaller values than |π/(2d)|, as can be seen in Figs. 6, 8, and 10, but the crossing points get closer to ±π/(2d) as ω p increases. Also, the effects of including further neighbors become smaller as ω p increases; they are smallest at ω p d/c = 5.0, as can be seen in Fig. 10. Next we consider values of R/d other than 1/3, but still keeping the same value of ω 1 = 1.4 × 10 15 rad/s (i.e., ω p d/c = 0.35). For a smaller R/d = 0.25, the variation of the band energies with k becomes smaller, as seen in Fig. 11, than it is in Fig. 3, but the crossing points between the L and T branches still occur at ±π/(2d). This behavior becomes clearer when the results for several values of R/d are plotted together as in Fig. 12. As R/d increases, the variation of the band energies with k, and the separation between the L and T branches at both the zone center and zone boundary, increase, but the L and T branches still cross at ±π/(2d). If we include more neighbors up to fifth nearest-neighbors, but consider only up to R/d = 0.4, we get the dispersion relations shown in Fig. 13. These show the same trends as in Fig. 12, except that the band crossing points occur at values of |k| slightly less than |π/(2d)|. Furthermore, the separation between the L and T bands increases slightly at k = 0, but decreases slightly at k = ±π/d.

IV. DISCUSSION
In this work we have calculated the photonic band structures of metal inverse opals and of a linear chain of spherical voids in a metallic host for frequencies below ω p , when ℓ = 1 using a tight-binding approximation. In both cases, we include only the ℓ = 1 "atomic" states of the voids. As a possible point of comparison, we have also computed the same band structures using the asymptotic forms of the spherical and modified spherical Bessel functions for small void radius. In this asymptotic region, there are only TM modes. The results for the linear chain of voids can be considered as the "inversions" of those in Ref. [13], in the sense discussed earlier. In other words, if we reflect our L and T branches with respect to the atomic energy level, we would get their L and T modes. Although we did not discuss this approach, we did attempt to use the plane wave expansion method to calculate the band structure for the inverse opals, similarly to Refs. [9] and [10]. Just as found in those papers, the photonic bands for modes below ω p depend on the number of plane waves included in the expansion and on the type of field, B or E, used in the expansion. Furthermore, this plane wave expansion method gives a large number of flat bands below ω p , which are difficult to interpret physically. Because of this problem, and because of the apparent non-convergence of this approach with the number of plane waves, we do not present these results here. By contrast, the band structures above ω p , when calculated using the plane wave expansion, varied smoothly with k and converged well with the number of plane waves included.
In calculating the tight-binding band structure for the linear chain (and the inverse opal structure), one should, in principle, include all the neighbors. But in practice, for the linear chain, it is sufficient to include only up to the fifth nearest-neighbors. This calculation is easily carried out, since the matrix in Eq. (39) is already diagonal in x, y, z and the sum converges quickly. In fact, even the inclusion of neighbors beyond the first two sets changes the band structure very little. The smallness of the further neighbor effect is particularly apparent when ω p d/c = 5.0, as in Fig. 10.
Although we studied the 3D and 1D lattices, we have not investigated a 2D lattice of spherical pores in a metallic host. For the 2D case, the matrix element M α,β (R) can again be readily calculated using our tight-binding approximation. It is expected to have some nonzero offdiagonal elements in addition to the diagonal elements. Thus band structures somewhat similar to the 3D band structures shown in Fig. 2 are also expected in the 2D case.
In the quasistatic case, for metal grains in air, when R/d is greater than about 0.4, it becomes important to include more than just ℓ = 1 as in Ref. [14]. Inclusion of such higher ℓ's might be rather difficult in the present dynamical case, though it would be possible in the quasistatic limit for 1D chains of spherical nanopores.
In the present work, we have considered only the case of one cavity per primitive cell and ℓ = 1. It would be of a great interest to consider multiple cavities per unit cell. Of course, in this case, the dimension of the matrix in Eq. (39) will increase and the Bloch states in Eq. (41) will acquire an additional index. It should be straightforward to extend the present work to such a case, which would make an interesting subject for future work.
The surface effect of the metal nanopores on the eigenstates becomes prominent as the radius of a pore (R) or the ratio of radius to center-to-center separation (R/d) increases. Shockley surface states form on metal surfaces depending on the type of metal, such as the Fermi wavevector k F or the Fermi energy ǫ F . These surface states will interact with the eigenstates, resulting in the change of eigenvalues. However we think these effects are negligible, especially in the quasistatic limit since in the asymptotic limit, these oscillatory interactions take the form [20] E asym pair (a) ∝ where a is the interatomic separation, Θ is the effective interaction phase shift, k F the Fermi wavevector of the isotropic surface state, and ǫ F the Fermi energy. The proportionality constant gives the consequences of scattering into bulk states. [20] The very slow a −2 decay of these interactions allows them to play a role at large separations, but the overall magnitude is small due to the sinusoidal term.
In summary, we have described a tight-binding method for calculating the photonic band structure of a periodic composite of spherical pores in a metallic host, and have applied it to both 1D and 3D systems. The method is fully dynamical, and is not limited to very small pores. The method does not have the convergence problems found when the magnetic or electric field is expanded in plane waves. Furthermore, there are no radiation losses to consider, unlike the complementary case of small metal particles in an insulating host, because the fields associated with these modes outside the pores are exponentially decaying. Thus, this method may be useful for a variety of periodic metal-insulator composites. It would be of interest to compare these calculations to experiments on such materials.