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

Kwangmoo Kim; D. Stroud

doi:10.1364/OE.21.019834

1. 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 three-dimensional (3D) photonic crystal [3]. Two-dimensional (2D) photonic crystals can be used as optical microcavities, microresonators [4], waveguides [5], lasers [6], or fibers [7] while one-dimensional (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 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) satisfied 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., [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. Many other workers have carried extensive work, both theoretical and experimental, on plasmonic waves in 1D and occasionally 2D arrays of metallic nanoparticles in an insulating host [14–27]. Park and Stroud [28] 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, but still in the quasistatic approximation where it is assumed that ∇ × E = 0, where E is the electric field.

Weber and Ford [29] 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.

Later, Gaillot et al. [30] 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 [30], the photonic band structure of Si inverse opal was calculated as a function of the infiltrated volume fraction f of air voids using 3D finite difference time domain (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 [30], 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 [31]. 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 [13].

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 1D chains of metallic nanoparticles in air. The absence of radiative decay has been previously noted in a multiple-scattering calculation of the band structure of a periodic array of pores in a host metal [32]. 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.

Wave propagation through void networks is of interest, in part, because of the special formalism needed to treat it. First, it is straightforward to go beyond the quasistatic regime (in which the electric field is assumed to be curl-free), as we discuss below. Thus, calculations can readily be carried out even for voids which are not much smaller than the wavelength of light. Furthermore, in contrast to waves propagating along chains or other periodic arrays of metal particles, there are no radiative losses, because the waves in the host region (i.e., in this case, the metallic region) are exponentially decaying.

The remainder of this paper is organized as follows. In Section 2, 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 3, 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 4 presents a summary and discussion.

2. 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 - ω_{p}^{2} / ω^{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.

2.1. 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.

2.1.1. 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 $\nabla \times E = \frac{i ω}{c} B$ and $\nabla \times B = - \frac{i ω}{c} ε E$ to obtain a single equation for B. We express the position- and frequency-dependent dielectric function ε(x, ω) as $1 / ε (x, ω) = θ (x) / (1 - ω_{p}^{2} / ω^{2}) + 1 - θ (x)$ , where the step function θ(x) = 1 inside the metallic region and θ(x) = 0 elsewhere. Then, after a little algebra, we obtain

[ω^{2} - ω_{p}^{2} (1 - θ (x))] \nabla \times (\nabla \times B) = \frac{ω^{2}}{c^{2}} (ω^{2} - ω_{p}^{2}) B .

This expression gives rise to different equations inside and outside the void.

For a spherical void within a metallic host, these equations are conveniently solved in spherical coordinates. For the TM modes, the solutions for B and E are given in the standard literature [33]. The coefficients of the solutions inside and outside the cavity can then be determined from the boundary conditions at r = R. The allowed frequencies for ω < ω_p are then found to satisfy

{k^{'}}^{2} [j_{ℓ} (k R) + k R j_{ℓ}^{'} (k R)] = - k^{2} \frac{j_{ℓ} (k R)}{k_{ℓ} (k^{'} R)} [k_{ℓ} (k^{'} R) + k^{'} R k_{ℓ}^{'} (k^{'} R)],

where k = ω/c and

k^{'} = {(ω_{p}^{2} - ω^{2})}^{1 / 2} / c

are the wave vectors inside and outside the void, j_ℓ is the spherical Bessel function, and k_ℓ is the modified spherical Bessel function (note that this k_ℓ is different from the wave vectors k and k′). In the limit kR ≫ 1 and k′R ≪ 1, we can readily obtain the asymptotic forms of the solutions to Eq. (2). The result is k′²(ℓ + 1) = k²ℓ, or, for a Drude metal,

ω^{2} = \frac{ℓ + 1}{2 ℓ + 1} ω_{p}^{2} .

The largest value,

ω = \sqrt{2 / 3} ω_{p}

, occurs at ℓ = 1 and the limiting value for large ℓ is

ω = ω_{p} / \sqrt{2}

.

2.1.2. TE Modes

A similar procedure for the electric field of the TE modes leads to the self-consistency condition

j_{ℓ} (k R) + k R j_{ℓ}^{'} (k R) = \frac{j_{ℓ} (k R)}{k_{ℓ} (k^{'} R)} [k_{ℓ} (k^{'} R) + k^{'} R k_{ℓ}^{'} (k^{'} R)] .

It is readily found that, in the asymptotic regime when kR ≪ 1 and k′R ≪ 1, there are no allowed eigenvalues for the TE modes.

2.2. 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 tight-binding method is very useful in treating narrow bands. In what follows, we try to suggest an analogous tight-binding 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

\nabla \times (\nabla \times E_{λ} (x)) + \frac{ω_{p}^{2} θ (x)}{c^{2}} E_{λ} (x) \equiv 𝒪 E_{λ} (x) = \frac{ω_{λ}^{2}}{c^{2}} E_{λ} (x),

where

𝒪 = \nabla \times (\nabla \times) + (ω_{p}^{2} / c^{2}) θ (x)

is the “Hamiltonian” of this system. Since 𝒪 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

\int E_{λ}^{*} (x) \cdot E_{μ} (x) d x = δ_{λ, μ} .

Since E_λ(x) is real for ω < ω_p, the complex conjugation is, in fact, unnecessary.

In Sec. 2.1.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₁, E₂, and E₃ satisfying the orthonormality relation in Eq. (6).

In order to obtain the tight-binding band structure built from these three modes, we need to calculate matrix elements of the form

M_{α, β} (R) = \int E_{α}^{*} (x) \cdot 𝒪 E_{β} (x - R) d x,

corresponding to two single-cavity modes associated with different cavities centered at the origin and at R. Here, 𝒪 is the “Hamiltonian” of the system as defined implicitly in Eq. (5).

Next, we introduce normalized Bloch states associated with the three ℓ = 1 single-cavity modes. In order to do this, we first make the standard tight-binding assumption that the “atomic” states corresponding to different cavities are orthogonal:

\int E_{λ}^{*} (x - R) \cdot E_{μ} (x - R^{'}) d x = δ_{λ, μ} δ_{R, R^{'}} .

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

E_{k, λ} (x) = N^{- 1 / 2} \sum_{R} e^{i k \cdot R} E_{λ} (x - R),

where k is a Bloch vector, and the R’s are the Bravais lattice vectors. In writing Eq. (9), 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

M_{λ, μ} (k) = \sum_{R} e^{i k \cdot R} M_{λ, μ} (R) .

We can then obtain the frequencies ω(k) by diagonalizing a 3 × 3 matrix as follows:

det | M_{λ, μ} (k) - (\frac{ω^{2} (k)}{c^{2}} - \frac{ω_{at}^{2}}{c^{2}}) δ_{λ, μ} | = 0,

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, e.g., [34]).

We briefly comment on the connection between this approach and that used by earlier workers [13,28]. 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) 3D as well as 1D; (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. (7). 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 𝒪 is just a little trickier. We first note that 𝒪 = 𝒪_R + 𝒪′, where 𝒪_R is the single-cavity operator: 𝒪_R = ∇ × (∇×) if x is inside the Rth cavity and $𝒪_{R} = \nabla \times (\nabla \times) + ω_{p}^{2} / c^{2}$ otherwise. Now we also have

𝒪_{R} E_{β} (x - R^{'}) = \frac{ω_{at}^{2}}{c^{2}} E_{β} (x - R^{'}),

since E_β is an eigenstate of 𝒪_R with an eigenvalue

ω_{at}^{2} / c^{2}

.

But since we are assuming that the overlap integral between “atomic” electric field states centered on different sites vanishes, the term involving 𝒪_R does not contribute to the matrix element M_α,β, which is therefore just given by

M_{α, β} (R) = \int E_{α} (x) \cdot 𝒪^{'} E_{β} (x - R) d x .

We can also write

𝒪^{'} = \frac{ω_{p}^{2}}{c^{2}} \sum_{R^{'}} θ_{R^{'}} (x),

where

θ_{R^{'}} (x) = θ (x - R^{'}),

is a step function which is unity inside the cavity centered at R′ and is zero otherwise.

A reasonable approximation to Eq. (14) might be to include just R′ = 0. In this case, we finally will get

M_{α, β} (R) ~ \frac{ω_{p}^{2}}{c^{2}} \int E_{α} (x) \cdot E_{β} (x - R) d x,

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

M_{α, β} (R) ~ \frac{ω_{p}^{2}}{c^{2}} E_{β} (- R) \cdot \int E_{α} (x) d x,

where once again the integral runs over the cavity centered at the origin.

The calculation of this matrix element from this integral expression is straightforward. The Cartesian components of the normalized eigenfunctions E_α(x) are readily calculated from the solutions discussed in Sec. 2.1.1 and the integral in Eq. (17). In fact, it turns out that all the integrals entering the matrix element can be obtained analytically. Given this matrix element, the computation of the full tight-binding band structure of the three p-bands is also straightforward.

3. Numerical Results

Before showing our numerical results, we first point out that, for a metal with a Drude dielectric function, the photonic band structures can be entirely expressed in terms of suitable dimensionless parameters. Specifically, for the inverse opal structure, the two relevant dimensionless parameters are R/d, the ratio of the void radius to the fcc lattice constant, and ω_pd/c. Of course, other combinations of these parameters would serve equally well. Given these parameters, the scaled frequencies ω/ω_p are functions of the scaled wave vector kd. In what follows, we consider only those scaled units.

For the inverse opals we arbitrarily assume that the ratio of void sphere radius R to the nearest neighbor distance $d / \sqrt{2}$ is 3/10, or $R / d = 3 / (10 \sqrt{2})$ , corresponding to a sphere filling fraction f = 0.160. In our calculations, we also arbitrarily use the value ω_pd/c = 1. For typical metallic values of ω_p, this would correspond to d of order 20 nm. For the linear chain of nanopores (see below), we use d to denote the separation between the centers of two adjacent nanopores and R to denote the radius of a nanopore as in Fig. 1(b); we take the ratio R/d = 1/3 at first, then change it later.

Fig. 1 Schematic diagram for (a) an inverse opal structure with a lattice constant d and a void sphere radius R; (b) a linear chain of nanopores with a pore separation d and a nanopore radius R.

Download Full Size | PDF

Our band structures for the inverse opals are expressed in terms of the standard notation for k values at symmetry points in the Brillouin zone. These are Γ = (0, 0, 0), X = (2π/d)(0, 0, 1), U = (2π/d)(1/4, 1/4, 1), L = (2π/d)(1/2, 1/2, 1/2), W = (2π/d)(1/2, 0, 1), and K = (2π/d)(3/4, 0, 3/4).

The metallic dielectric functions we assume for the inverse opals and linear chain of nanopores are of the usual Drude form,

ε (ω) = 1 - \frac{ω_{p}^{2}}{ω^{2}},

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,

ω_{p}^{2} = 4 π n e^{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 of ω/ω_p.

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 1. These two inequalities are reasonable for the choice of “inverse opal” system parameters $R / d = 3 / (10 \sqrt{2})$ and ω_pd/c = 1, because

\begin{matrix} k R & = & \frac{ω}{c} R < \frac{ω_{p}}{c} R = \frac{ω_{p} d}{c} \frac{R}{d} = \frac{3}{10 \sqrt{2}} = 0.2121, \\ k^{'} R & = & \frac{\sqrt{ω_{p}^{2} - ω^{2}}}{c} R = \sqrt{{(\frac{ω_{p} R}{c})}^{2} - {(\frac{ω R}{2})}^{2}} = \sqrt{{(\frac{ω_{p} d}{c} \frac{R}{d})}^{2} - {(k R)}^{2}} \\ = & \sqrt{{(\frac{3}{10 \sqrt{2}})}^{2} - {(k R)}^{2}} = \sqrt{0.045 - {(k R)}^{2}} < \sqrt{0.045} = 0.2121 . \end{matrix}

The (modified) spherical Bessel functions in Eq. (2) 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.

Table 1. TM mode frequencies ω′ = ωd/(2πc), where ω < ω_p and ω_pd/c = 1, calculated for an isolated spherical cavity (“Infinite medium”) and those when both kR ≪ 1 and k′R ≪ 1. The (modified) spherical Bessel functions are extremely close to the ω′ axis for ℓ > 5, so that it is difficult to get eigenfrequencies for ℓ > 5 in the isolated spherical cavity. However this does not happen when kR ≪ 1 and k′R ≪ 1.

View Table

The solutions to Eq. (4) do not exist for ω < ω_p with ω_pd/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.

For our fcc calculations, we calculate the band structure including only the 12 nearest-neighbors of the cavity at the origin. Thus R = (d/2)(±1, ±1, 0), (d/2)(±1, ∓1, 0), (d/2)(±1, 0, ±1), (d/2)(±1, 0, ∓1), (d/2)(0, ±1, ±1), and (d/2)(0, ±1, ∓1). Assuming ω_pd/c = 1.0 and using ω_atd/(2πc) = 0.1296 (ω_at = 0.8143ω_p) 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, which behave as expected for the p-bands. The bandwidth is relatively small as M_α,β(R)d² ∼ 0.001, which proves the general relation between the bandwidth and the overlap integral [34]. 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.

Fig. 2 Tight-binding inverse opal band structure for ω < ω_p with $R / d = 3 / (10 \sqrt{2})$ and ω_pd/c = 1.0, using ω_atd/(2πc) = 0.1296 (ω_at = 0.8143ω_p) for ℓ = 1 in an infinite medium. The horizontal dotted line represents the “atomic” level.

Download Full Size | PDF

In Fig. 3, we plot the frequency of the triply degenerate Γ point as a function of the number of nearest neighbor shells included, up to seven shells. It can be seen that the frequency changes significantly over this range of shells, though we find, of course, that the triple degeneracy is retained (since it is required by symmetry). The band structure will certainly converge extremely well for a sufficient number of neighbor shells, because the hopping integral will eventually fall off exponentially with separation, with an inverse decay length ${(ω_{p}^{2} - ω^{2})}^{1 / 2} / c$ . For the present case, ω_pd/c = 1.0, so this exponential decay does not fully set in until a fairly large number of shells is included. For a larger value of ω_pd/c, the eigenfrequencies will converge much more quickly with number of shells. This problem is worse in 3D than in 1D (see below), because the magnitude of the hopping integral, as a function of separation r, varies as 1/r³ for small ω_pr/c, while the number of terms in each shell increases as r². Convergence is assured, however, with a sufficient number of shells, because of the exponential decay which sets in at large r. As is seen below, the convergence is much faster in 1D.

Fig. 3 Dependence of the triply degenerate frequency at Γ on the number of nearest neighbor shells included in the tight-binding calculation, for the inverse-opal calculation shown in Fig. 2. Up to seven shells are included.

Download Full Size | PDF

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. We use the ratio R/d = 1/3 and the dimensionless parameter ω_pd/c = 0.35. These are arbitrarily chosen to be the same as used in [13]. The “atomic” frequency is found by solving Eq. (2) and gives ω_atd/(2πc) = 0.0454 (ω_at = 0.8150ω_p) for ℓ = 1 in an infinite medium. Our resulting tight-binding dispersion relations are shown by open triangles in Fig. 4 with only nearest-neighbors included. The transverse (T) branches are twofold degenerate, while the longitudinal (L) branch is non-degenerate. 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 sum also converges quickly, so there is only a slight difference between the dispersion relation including through the next-nearest-neighbors and that including through the 5th nearest-neighbors. However, if one includes more than nearest-neighbor overlap, the L and T branches no longer cross exactly at k = ±π/(2d). [32] considered the long-wavelength tight-binding limit of their multiple-scattering dispersion relation for a periodic linear chain of dielectric cavities in a metallic host, and also found that the L and T branches of the ℓ = 1 dispersion relations crossed at k = ±π/(2d) if only nearest-neighbor hopping is included.

Fig. 4 Tight-binding results of a periodic chain of nanopores in a Drude metal host, for ω < ω_p. We take R/d = 1/3 and ω_pd/c = 0.35, using ω_atd/(2πc) = 0.0454 (ω_at = 0.8150ω_p) for ℓ = 1 in an infinite medium. The horizontal dotted line represents the “atomic” level. Three different numbers of neighbors are included: nearest-neighbors (nn’s), next-nearest-neighbors (nnn’s), and fifth-nearest-neighbors (5nn’s). In this and the following two plots, “L” and “T” denote the longitudinal and transverse branches, respectively.

Download Full Size | PDF

To check the convergence of the number of nn’s included, we have calculated the ω at k = 0 since the difference between different numbers of nn’s is the most evident there. For ω_pd/c = 1.0 and ω_atd/(2πc) = 0.1291 (ω_at = 0.8112ω_p), we find that both the L and T frequencies change by less than 0.3% in going to the 5th nn shell, and are unchanged to within 0.05% thereafter, up to 10 nn shells. This convergence is quicker than in 3D and can be readily seen in the ω(k) plots.

We have carried out similar calculations using other values of the parameter ω_pd/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, and the corresponding results including more overlap integrals, for a typical example, ω_pd/c = 5.0, are shown in Fig. 5, since the results for ω_pd/c = 1.0 and 2.0 are similar to Fig. 4 except for the increase of ω along the y-axis and ω_at. It is also striking that, as ω_pd/c increases in going from Fig. 4 to 5, the ratio r_LT of the width of the L band to that of the T band decreases. In Fig. 4, r_LT > 1, while in Fig. 5, r_LT < 1.

Fig. 5 Same as Fig. 4, except ω_pd/c = 5.0 and ω_atd/(2πc) = 0.5691 (ω_at = 0.7152ω_p).

Download Full Size | PDF

One could also say that, except for an overall scale factor, Fig. 5 looks like an inverted image of Fig. 4 about the horizontal line of ω_at. For the nn case, the T and L bands cross at ±π/(2d), while they cross at smaller values than |π/(2d)| when further neighbors are included, as can be seen in Figs. 4 and 5, but the crossing points get closer to ±π/(2d) as ω_pd/c increases. Also, the effects of including further neighbors become smaller as ω_pd/c increases; they are smallest at ω_pd/c = 5.0, as can be seen in Fig. 5.

Next we consider values of R/d other than 1/3, but still keeping the same value of ω_pd/c = 0.35 and including up to the fifth nearest-neighbors. For a smaller R/d = 0.25, the variation of the band energies with k becomes smaller, as seen by open triangles in Fig. 6, than it is in Fig. 4, but the crossing points between the L and T branches still occur at values of |k| slightly less than |π/(2d)|. This behavior becomes clearer when the results for more values of R/d are plotted together as in Fig. 6. 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 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 compared to the results with only nn’s included. We show only R/d up to 0.4 in this Figure because, in the quasistatic limit, there is evidence that for larger values of R/d the dispersion relations are significantly modified by higher values of ℓ[28].

Fig. 6 Plotting together three different results for ω < ω_p, all with ω_pd/c = 0.35, but with different (R/d)’s: R/d = 0.25 and ω_atd/(2πc) = 0.04546 (ω_at = 0.8161ω_p); R/d = 0.33 and ω_atd/(2πc) = 0.04544 (ω_at = 0.8157ω_p); R/d = 0.40 and ω_atd/(2πc) = 0.04543 (ω_at = 0.8156ω_p), with inclusion of up to the fifth nearest-neighbors. For the larger values of R/d, it may be necessary to include more than just ℓ = 1.

Download Full Size | PDF

4. 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 somewhat resemble those of [13] for a chain of metallic spheres in an insulating host, except that the L branch lies above the doubly degenerate T branch.

It is of interest to compare our work to that of other workers. In particular, Lidorikis et al. [35] have previously developed a method based on a linear combination of “atomic” orbitals to treat propagation of photonic waves in a periodic composite. Their method gives an excellent account of the photonic band structure of a periodic array of parallel dielectric cylinders in a host of a smaller dielectric constant, and agrees very well with that computed by expanding the fields in plane waves. However, our work deals with a very different system from [35]. We treat a metal dielectric composite, in which the dielectric is a network of voids and the metal has a frequency-dependent dielectric function which is negative at the frequencies of interest. By contrast, [35] deals with a composite in which both components have frequency-independent, positive dielectric constants. As a result, our tight-binding bands are linear combinations of the individual plasmons associated with the voids, whereas those of [35] come from Mie resonances of the 2D dielectric cylinders. We believe that these can never be very narrow if the cylinders have a positive dielectric constant, and therefore, their bands are not plasmonic in character as ours are.

A somewhat different method, based on multiple scattering theory, has been used to treat 1D chains of nanocavities in a metallic host in [32], and a 3D extension of this multiple-scattering method has been used to treat an fcc crystal of silicon spheres in a metallic host in [36]. Although the approach in these two papers is quite different from the tight-binding method used in our calculations, it is still of interest to compare results from the two methods where possible. For example, we can compare our tight-binding results for the inverse opal structure [Fig. 2] along the Γ − X direction in the Brillouin zone with those presented in [36] in their Fig. 1(a). A quantitative comparison is not possible, because our calculations and those of [36] are carried out for a quite different fcc lattice constant, pore volume fraction, and dielectric constant of the material in the pore space (ε = 11.9 in [36], ε = 1 in our calculations). Nevertheless, the doubly degenerate bands in our calculation appear to have the same general shape as the corresponding bands shown in their Fig. 1(a), that is, an increase with increasing k_z starting from the Γ point followed by a flattening of these bands as the point X is approached. The non-degenerate band along Γ − X in our calculations falls monotonically starting from the Γ point (lower band along the line Γ − X in our Fig. 2), whereas the corresponding band shown in their Fig. 1(a) seems to be nearly flat. We tentatively attribute this difference primarily to the difference in the parameters of the two calculations.

In tight-binding calculations, it is important to ascertain how sensitive the results are to the number n of nearest neighbor shells included in the calculations. In order to answer this question, we show in Fig. 3 the calculated energies at Γ as functions of n, for the inverse opal structure. As mentioned earlier, the frequency will certainly converge extremely well for a sufficient number of neighbor shells, because the hopping integral will eventually fall off exponentially with separation. This convergence is faster in 1D than in 3D because the magnitude of the hopping integral varies as 1/r³ for small ω_pr/c, while the number of terms in each shell increases as r².

In the case of 1D bands, [32] shows that their multiple-scattering formalism, in the limit ka ≪ 1 (where a is the cavity radius), leads to a dispersion relation of the tight-binding form ω = ω₀ + ω₁ cos(kd), where d is distance between the sphere centers. This result is obtained provided that one neglects interactions other than between nearest neighbor cavities, and also that one disregards interactions between the lowest (ℓ = 1) plasmon band and all the higher bands. In our present tight-binding calculation, we do not need to assume ka ≪ 1, and we are also able to go beyond nearest neighbor hopping, though we do include only the ℓ = 1 bands. Thus, our approach is somewhat different from that of [32], though it gives similar results in certain parameter regimes. While we do not include modes with ℓ > 1, our approach could readily be extended to do so. Finally, our bands do not cross exactly at |k|d = π/2 when we include more than one shell of neighbors in the tight-binding calculation; so our dispersion relations are not exactly of the form ω₀ + ω₁ cos(kd) found in [32]. Thus, in short, our 1D results bear some similarities to those of [32], but are obtained in a different way and are calculated in somewhat different regimes.

Next, we briefly discuss the non-orthogonality of our basis functions for the TM “atomic” modes corresponding to different spherical voids. The electric fields of these modes are given in [33]. This non-orthogonality also arises in the usual tight-binding method as applied to electronic states in solids (see [34]), when one uses atomic states as basis functions. In that case, the common procedure is to neglect the overlap integral between states centered on different atoms. In our case, we have neglected the overlap integral

S_{λ μ} (R) = \int E_{λ}^{*} (x) \cdot E_{μ} (x - R) d x

for modes centered on voids separated by a nonzero lattice vector R. This neglect is reasonable, for our problem, because the basis functions decay exponentially into the metal, and thus the overlap should be small.

If one wishes to correct for the non-orthogonality, there is a well-established procedure for doing so. It has been discussed in the electronic case, e.g., in [37] and [38], based on a transformation originally developed by Löwdin [39]. For the present problem, one would obtain a Hamiltonian H with an orthonormal basis from the the original Hamiltonian ℋ with a nonorthogonal basis by means of the transformation

H = {(I + S)}^{- 1 / 2} ℋ {(I + S)}^{- 1 / 2},

where S is a matrix whose elements are given by Eq. (20), and I is the identity matrix. To carry out this procedure, one would simply need to calculate the overlap matrix elements S_λμ (R) and execute the transformation in Eq. (21). We believe that the resulting corrections to the band structure would be small, because, as in the analogous electronic problem for narrow tight-binding bands, the off-diagonal elements of the overlap matrix are expected to be small for exponentially decaying basis functions, compared to those of the original Hamiltonian ℋ.

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 [28]. Inclusion of such higher ℓ’s might be rather difficult in the present dynamical case, though it would be straightforward in the quasistatic limit for 1D chains of spherical nanopores.

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.

Acknowledgments

This work was supported by the National Science Foundation through the Materials Research Science and Engineering Center at The Ohio State University ( DMR-0820414), and by Department of Energy Grant No. DE-FG02-07ER46424. All of the calculations using plane wave expansions were carried out on the P4 Cluster at the Ohio Supercomputer Center, with the help of a grant of time. We also thank Korea Institute for Advanced Study for providing computing resources (KIAS Center for Advanced Computation Abacus) for this work.

References and links

1. J. Q. Xia, Y. R. Ying, and S. H. Foulger, “Electric-Field-Induced Rejection-Wavelength Tuning of Photonic-Bandgap Composites,” Adv. Mater. 17, 2463–2467 (2005) [CrossRef] .

2. Kurt Busch and Sajeev John, “Liquid-Crystal Photonic-Band-Gap Materials: The Tunable Electromagnetic Vacuum,” Phys. Rev. Lett. 83, 967–970 (1999) [CrossRef] .

3. Eli Yablonovitch, “Inhibited Spontaneous Emission in Solid-State Physics and Electronics,” Phys. Rev. Lett. 58, 2059–2062 (1987) [CrossRef] [PubMed] .

4. A. Scherer, O. Painter, B. D’Urso, R. Lee, and A. Yariv, “InGaAsP photonic band gap crystal membrane microresonators,” J. Vac. Sci. Technol. B 16, 3906 (1998) [CrossRef] .

5. Attila Mekis, J. C. Chen, I. Kurland, Shanhui Fan, Pierre R. Villeneuve, and J. D. Joannopoulos, “High Transmission through Sharp Bends in Photonic Crystal Waveguides, ” Phys. Rev. Lett. 77, 3787–3790 (1996) [CrossRef] [PubMed] .

6. O. Painter, R. K. Lee, A. Scherer, A. Yariv, J. D. O’Brien, P. D. Dapkus, and I. Kim, “Two-Dimensional Photonic Band-Gap Defect Mode Laser,” Science 284, 1819–1821 (1999) [CrossRef] [PubMed] .

7. F. Benabid, J. C. Knight, G. Antonopoulos, and P. St. J. Russell, “Stimulated Raman Scattering in Hydrogen-Filled Hollow-Core Photonic Crystal Fiber,” Science 298, 399–402 (2002) [CrossRef] [PubMed] .

8. Y. Cao, J. O. Schenk, and M. A. Fiddy, “Third order nonlinear effect near a degenerate band edge,” Opt. Photonics Lett. 1, 1–7 (2008) [CrossRef] .

9. Arthur R. McGurn and Alexei A. Maradudin, “Photonic band structures of two- and three-dimensional periodic metal or semiconductor arrays,” Phys. Rev. B 48, 17576–17579 (1993) [CrossRef] .

10. V. Kuzmiak, A. A. Maradudin, and F. Pincemin, “Photonic band structures of two-dimensional systems containing metallic components,” Phys. Rev. B 50, 16835–16844 (1994) [CrossRef] .

11. V. Kuzmiak and A. A. Maradudin, “Photonic band structures of one- and two-dimensional periodic systems with metallic components in the presence of dissipation,” Phys. Rev. B 55, 7427–7444 (1997) [CrossRef] .

12. I. H. H. Zabel and D. Stroud, “Photonic band structures of optically anisotropic periodic arrays,” Phys. Rev. B 48, 5004–5012 (1993) [CrossRef] .

13. Mark L. Brongersma, John W. Hartman, and Harry A. Atwater, “Electromagnetic energy transfer and switching in nanoparticle chain arrays below the diffraction limit,” Phys. Rev. B 62, R16356–R16359 (2000) [CrossRef] .

14. Stefan A. Maier, Pieter G. Kik, Harry A. Atwater, Sheffer Meltzer, Elad Harel, Bruce E. Koel, and Ari A.G. Requicha, “Local detection of electromagnetic energy transport below the diffraction limit in metal nanoparticle plasmon waveguides,” Nat. Mater. 2, 229–232 (2003) [CrossRef] [PubMed] .

15. S. A. Maier, M. L. Brongersma, P. G. Kik, S. Meltzer, A. A. G. Requicha, and H. A. Atwater, “Plasmonics–A Route to Nanoscale Optical Devices,” Adv. Mater. 13, 1501–1505 (2001) [CrossRef] .

16. Prashant K. Jain, Wenyu Huang, and Mostafa A. El-Sayed, “On the Universal Scaling Behavior of the Distance Decay of Plasmon Coupling in Metal Nanoparticle Pairs: A Plasmon Ruler Equation,” Nano Lett. 7, 2080–2088 (2007) [CrossRef] .

17. LinLin Zhao, K. Lance Kelly, and George C. Schatz, “The Extinction Spectra of Silver Nanoparticle Arrays: Influence of Array Structure on Plasmon Resonance Wavelength and Width,” J. Phys. Chem. B 107, 7343–7350 (2003) [CrossRef] .

18. Q.-H. Wei, K.-H. Su, S. Durant, and X. Zhang, “Plasmon Resonance of Finite One-Dimensional Au Nanoparticle Chains,” Nano Lett. 4, 1067–1071 (2004) [CrossRef] .

19. L. A. Sweatlock, S. A. Maier, H. A. Atwater, J. J. Penninkhof, and A. Polman, “Highly confined electromagnetic fields in arrays of strongly coupled Ag nanoparticles,” Phys. Rev. B 71, 235408 (2005) [CrossRef] .

20. Shengli Zou and George C. Schatz, “Narrow plasmonic/photonic extinction and scattering line shapes for one and two dimensional silver nanoparticle arrays,” J. Chem. Phys. 121, 12606–12612 (2004) [CrossRef] [PubMed] .

21. A. F. Koenderink and Albert Polman, “Complex response and polariton-like dispersion splitting in periodic metal nanoparticle chains,” Phys. Rev. B 74, 033402 (2006) [CrossRef] .

22. Andrea Alù and Nader Engheta, “Theory of linear chains of metamaterial/plasmonic particles as subdiffraction optical nanotransmission lines,” Phys. Rev. B 74, 205436 (2006) [CrossRef] .

23. Alexandre Bouhelier, Renaud Bachelot, Jin Seo Im, Gary P. Wiederrecht, Gilles Lerondel, Sergei Kostcheev, and Pascal Royer, “Electromagnetic Interactions in Plasmonic Nanoparticle Arrays,” J. Phys. Chem. B 109, 3195–3198 (2005) [CrossRef] .

24. Vadim A. Markel and Andrey K. Sarychev, “Propagation of surface plasmons in ordered and disordered chains of metal nanospheres,” Phys. Rev. B 75, 085426 (2007) [CrossRef] .

25. Yu-Rong Zhen, Kin Hung Fung, and C. T. Chan, “Collective plasmonic modes in two-dimensional periodic arrays of metal nanoparticles,” Phys. Rev. B 78, 035419 (2008) [CrossRef] .

26. W. Jacak, J. Krasnyj, J. Jacak, A. Chepok, L. Jacak, W. Donderowicz, D. Z. Hu, and D. M. Schaadt, “Undamped collective surface plasmon oscillations along metallic nanosphere chains,” J. Appl. Phys. 108, 084304 (2010) [CrossRef] .

27. Matthew D. Arnold, Martin G. Blaber, Michael J. Ford, and Nadine Harris, “Universal scaling of local plasmons in chains of metal spheres,” Opt. Express 18, 7528–7542 (2010) [CrossRef] [PubMed] .

28. Sung Yong Park and David Stroud, “Surface-plasmon dispersion relations in chains of metallic nanoparticles: An exact quasistatic calculation,” Phys. Rev. B 69, 125418 (2004) [CrossRef] .

29. W. H. Weber and G. W. Ford, “Propagation of optical excitations by dipolar interactions in metal nanoparticle chains,” Phys. Rev. B 70, 125429 (2004) [CrossRef] .

30. D. Gaillot, T. Yamashita, and C. J. Summers, “Photonic band gaps in highly conformal inverse-opal based photonic crystals,” Phys. Rev. B 72, 205109 (2005) [CrossRef] .

31. Ali E. Aliev, Sergey B. Lee, Anvar A. Zakhidov, and Ray H. Baughman, “Superconductivity in Pb inverse opal,” Physica C 453, 15–23 (2007) [CrossRef] .

32. G. Gantzounis and N. Stefanou, “Cavity-plasmon waveguides: Multiple scattering calculations of dispersion in weakly coupled dielectric nanocavities in a metallic host material,” Phys. Rev. B 74, 085102 (2006) [CrossRef] .

33. See, e.g.,J. D. Jackson, “Earth and Ionosphere as a Resonant Cavity: Schumann Resonances,” in Classical Electrodynamics, 3rd ed. (Wiley, New York, 1999), pp. 374–376.

34. See, e.g.,N. W. Ashcroft and N. D. Mermin, “Problem 2. Tight-Binding p-Bands in Cubic Crystals” & “General remarks on the tight-binding method,” in Solid State Physics (Saunders College Publishing, Orlando, 1976), pp. 189–190& pp. 184–185.

35. E. Lidorikis, M. M. Sigalas, E. N. Economou, and C. M. Soukoulis, “Tight-Binding Parametrization for Photonic Band Gap Materials,” Phys. Rev. Lett. 81, 1405–1408 (1998) [CrossRef] .

36. N. Stefanou, A. Modinos, and V. Yannopapas, “Optical transparency of mesoporous metals,” Solid State Commun. 118, 69–73 (2001) [CrossRef] .

37. B. A. McKinnon and T. C. Choy, “Significance of nonorthogonality in tight-binding models,” Phys. Rev. B 52, 14531–14538 (1995) [CrossRef] .

38. Madhu Menon and K. R. Subbaswamy, “Transferable nonorthogonal tight-binding scheme for silicon,” Phys. Rev. B 50, 11577–11582 (1994) [CrossRef] .

39. Per-Olov Löwdin, “On the NonOrthogonality Problem Connected with the Use of Atomic Wave Functions in the Theory of Molecules and Crystals,” J. Chem. Phys. 18, 365–375 (1950) [CrossRef] .

	Infinite medium	kR ≪ 1, k′R ≪ 1

ℓ = 1	0.1296	0.1299
ℓ = 2	0.1232	0.1233
ℓ = 3	0.1203	0.1203
ℓ = 4	0.1186	0.1186
ℓ = 5	0.1178	0.1175

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

Abstract

1. Introduction