Novel approach to photonic bands with frequency-dependent dielectric constants

Kazuaki Sakoda; Jun Kawamata

doi:10.1364/OE.3.000012

1. Introduction

Photonic crystals, which have spatially periodic dielectric structures, can control the nature of radiation field.[1–3] They can create photonic band gaps, or the frequency ranges in which no eigenmodes of radiation field exist. Then the fundamental optical properties of atoms and molecules such as the rate of spontaneous emission may change drastically when they are embedded in the photonic crystals. In addition, localized eigenmodes of radiation field may appear in the photonic band gaps when appropriate defects are introduced to the regular dielectric structures of the photonic crystals.[4, 5] Because the localized modes have high quality factors and show extremely narrow transmission spectra, they can be used as good resonators and optical filters. Moreover, as was pointed out by Yablonovitch,[6] the single-mode light-emitting diodes that utilize the emission through the localized modes may have such properties as good temporal and spatial coherence, low noise, high efficiency, and high modulation rate. Then the micro-fabrication of the photonic crystals with semiconducting materials have been tried energetically.

When we analyze and design the photonic crystals made from semiconductors that work in visible and near infrared frequency ranges, we have to take into account the large frequency dependence of their dielectric constants. However, most of the theoretical analyses of the photonic band structures have been made for the frequency-independent case. An exception is the work reported by Kuzmiak et al.[7] They calculated the photonic band structure, or the dispersion relation, of a two-dimensional (2D) array of metallic cylinders with a frequency-dependent dielectric constant of Drude type by means of the plane-wave expansion method. They showed that the wave equation derived from the Maxwell’s equations leads to a well-defined eigenvalue problem for this particular case when the electric field of the eigenmode is parallel to the cylinder axis (E polarization) and the plane-wave expansion method gives a good convergence. But we have to devise another method when the wave equation does not lead to a well-defined eigenvalue problem, and then the good convergence of the plane-wave expansion is not necessarily guaranteed.

Then, in this paper, we would like to present a new and general method to calculate the dispersion relations of the photonic crystals with frequency-dependent dielectric constants. Our method is based on the numerical simulation of the radiation process of a virtual oscillating dipole moment embedded in the photonic crystal. Our algorithm is suitable to vector processing and the CPU time can be reduced by a considerable amount by using a vector processor. This method was first formulated and applied to the problem of mid-gap eigenmodes localized at point [8] and line defects.[9] In Sec. 2, we will extend this method in order to treat frequency-dependent dielectric constants and to calculate the dispersion relations. Then it will be applied to the particular problem of Kuzmiak et al, [7] i.e., the photonic band structure of the 2D square array of metallic cylinders for E polarization. The numerical results will be presented and compared with those of Kuzmiak et al in Sec. 3.

2. Theory

First, let us summarize our method briefly. Originally we considered the radiation process of the virtual oscillating dipole moment P _d(r, t) that was embedded in the photonic crystal with a structural defect:

P_{d} (r, t) = μ δ (r - r_{0}) \exp (- iωt),

where μ, and r ₀ are the magnitude and the position of the dipole moment, respectively, ω is the angular frequency of the oscillation, and δ is Dirac’s delta function. Then we showed by means of a Green’s-function method[10] that the induced electric field, E(r, t), and the electromagnetic energy radiated in a unit time, U, are given as follows.[8]

E (r, t) ≃ - \frac{2 π ω_{d} {μ \cdot E_{d}^{*} (r_{0})} E_{d} (r) \exp (- iωt)}{V (ω - ω_{d} + iγ)},

U ≃ \frac{π ω_{d}^{2} γ {∣ μ \cdot E_{d} (r_{0}) ∣}^{2}}{V {{(ω - ω_{d})}^{2} + γ^{2}}},

where ω_d and E _d(r) denote the eigen-angular frequency and the eigenfunction of the relevant defect mode, γ is its decay rate, V is the volume of the photonic crystal. E _d(r) is normalized as follows.

\int_{V} ϵ (r) {∣ E_{d} (r) ∣}^{2} d r = V .

Here ϵ(r) is the position-dependent dielectric constant of the photonic crystal. When we derived Eqs. (2) and (3), we assumed that ω was close to ω_d that was isolated in the photonic band gap and neglected the contribution from all other eigenmodes. If we can evaluate the frequency dependence of U, we can obtain the eigen-angular frequency as the resonance frequency [see Eq. (3)], and then the induced electric field is proportional to the eigenfunction [see Eq. (2)]. Then we discretized the wave equation derived from the Maxwell’s equations to obtain a difference equation and solved it numerically. We could obtain ω_d as the resonance frequency as we had expected, and the calculated values were very close to observed ones. [8, 9]

Now, we extend this method in order to treat the frequency dependence of the dielectric constants and to calculate the dispersion relations. We begin with the following two equations.

\nabla \times E (r, t) = - \frac{1}{c} \frac{\partial}{\partial t} H (r, t),

\nabla \times H (r, t) = \frac{1}{c} \frac{\partial}{\partial t} {D_{0} (r, t) + 4 π P_{d} (r, t)},

where c is the light velocity in vacuum. In Eq. (6), P _d(r,t) stands for the virtual oscillating dipole moment as before and D ₀(r,t) denotes the electric displacement due to the regular dielectric structure of the photonic crystal. The latter is generally given by the convolution integral of the electric field E(r, t) and the dielectric response function Φ(r,t):

D_{0} (r, t) = \int_{- \infty}^{\infty} dt' Φ (r, t - t') E (r, t') .

Φ(r, t) is given by the Fourier transform of the dielectric constant ϵ(r, ω), which is now a function of the frequency as well as the spatial coordinates:

Φ (r, t) = \frac{1}{2 π} \int_{- \infty}^{\infty} dωϵ (r, ω) \exp (- iωt) .

In order to show the detail of our numerical method, let us treat the particular problem discussed by Kuzmiak et al.[7] in what follows. So, we consider a photonic crystal that contains metallic components and assume a dielectric constant of Drude type in the metal:

ϵ (r, ω) = ϵ_{\infty} [1 - \frac{ω_{p}^{2} τ}{ω (ωτ + i)}],

where ω_p is the plasma frequency, τ is the relaxation time, and ϵ _∞ is the dielectric constant at sufficiently high frequencies. In Eq. (9), we took into account the imaginary part of the dielectric constant in order to fulfill the Kramers-Kronig relation, and hence, the causality. Then Eq. (8) leads to

Φ (r, t) = ϵ_{\infty} δ (t) + ϵ_{\infty} ω_{p}^{2} τ [1 - \exp (- \frac{t}{τ})] θ (t),

where θ(t) is a step function. From Eqs. (5), (6), (7), and (10), we obtain the following wave equation for the metallic region.

\nabla \times {\nabla \times E (r, t)}

On the other hand, we assume that ϵ(r, ω) = 1 (air) outside the metal. Then for this region,

\nabla \times {\nabla \times E (r, t)} = - \frac{1}{c^{2}} [\frac{\partial^{2}}{\partial t^{2}} E (r, t) - 4 π ω^{2} μ δ (r - r_{0}) \exp (- iωt)] .

We discretized Eqs. (11) and (12) to obtain difference equations and solved them numerically with an initial condition E(r, 0) = 0 and a boundary condition

E (r + a, t) = \exp (i k \cdot a) E (r, t),

where k is a wave vector in the first Brillouin zone and a is the elementary lattice vector. The latter condition extracts the contribution to the radiated electromagnetic field from particular eigenmodes with the specified wave vector. Therefore, we can calculate the resonance frequency as a function of k, i.e., we can obtain the dispersion relations. According to Kuzmiak et al.,[7] we considered a 2D photonic crystal composed of a square array of metallic cylinders with a radius R for the numerical calculation. The following parameters were assumed: R/a = 0.472, ϵ _∞ = 1.0, and ω_p a/27πc = 1.0, where a denotes the lattice constant. Because of the boundary condition, Eq. (13), it was enough to deal with only one unit cell, and therefore, the CPU time necessary for the numerical calculation was small. In the actual calculation, the unit cell was divided into 40 × 40 parts, and one period of the oscillation was divided into 160 parts in order to discretize the wave equation. The further decrease of the size of the spatial and temporal meshes did not bring about an apparent change in the resonance frequencies.

3. Results and discussion

The numerical calculation was done for the wave vectors in the (1, 0) direction in the first Brillouin zone of the 2D square lattice. As an example, Fig. 1 shows the electromagnetic energy radiated by the oscillating dipole moment as a function of the oscillation frequency for k = 0. The abscissa represents the normalized frequency. □, ○, ◇, and • denote the accumulated electromagnetic energy after 10, 20, 50, and 100 cycles of the oscillation, respectively. A resonance at ωa/2πc = 0.745 is clearly observed in this figure.

Figure 1. The electromagnetic energy radiated by an oscillating dipole moment located at r ₀/a = (0.3,0) as a function of the oscillation frequency for k = 0. The electric field parallel to the cylinder axis (E polarization) was calculated. □, ○, ◇, and • denote the accumulated electromagnetic energy after 10, 20, 50, and 100 cycles of the oscillation, respectively. The abscissa represents the normalized frequency. According to the previous calculation by Kuzmiak et al. (Ref. [7]), the following parameters were used: R/a = 0.472, ω = 1.0, ϵ _∞ ω_p a/2πc = 1.0. We also assumed that τ^-1 = 0.01 × ω_p. A resonance at ωa/2πc = 0.745 is clearly observed.

Download Full Size | PDF

Figure 2. The dispersion relations in the (1, 0) direction (or the Γ-X direction) of the 2D photonic crystal composed of the metallic cylinders. The ordinate is the normalized frequency and the abscissa is the normalized wave vactor. The same parameters as for Fig. 1 were used.

Download Full Size | PDF

The dispersion relations thus obtained are presented in Fig. 2 where the symmetry of each eigenmode is also shown. In this figure, the number in parentheses is given in order of the ascending frequency when there are more than one modes of the same symmetry in the analysed spectral region. Note that there is no eigenmode for ωa/2πc < 0.745. We can show that this cut-off frequency is consistent with the result of the long-wavelength approximation of the Maxwell’s equations. When we compare this figure with Fig. 1(b) of Ref. [7], it is found that the overall feature is quite common to both calculations except the degeneracy of the fourth and the fifth bands at the Γ point. As will be shown later, the symmetry of the eigenfunctions of these two modes is really that of the E representation of the C _4v point group that is doubly degenerate. Therefore, we believe that our calculation is more accurate concerned with this matter.

The distribution of the electric fields of the eigenmodes at the Γ point is presented in Fig. 3. The maximum of each electric field is normalized to unity. It is very clear that each eigenmode shows its peculiar symmetry. Especially, Fig. 3(d) is a replica of Fig. 3(e) given by a 90° rotation, which implies that they are degenerate and attributed to the E representation. Finally, Fig. 4 shows the field distribution of the eigenmodes at the X point. They are consistent with the symmetry assignments given in Fig. 2.

Figure 3. The distribution of the electric fields of the (a) A ₁(1), (b) A ₁(2), (c) B ₁, and (d) (e) E modes at the Γ point that has the C _4v symmetry. The maximum of each electric field is normalized to unity. For all eigenmodes, the eigenfunctions show their peculiar symmetries. For doubly degenerate E mode, one eigenfunction is a replica of the other given by a 90° rotation.

Download Full Size | PDF

We would like to conclude this section by making three remarks. First, we have treated the frequency-dependent dielectric constant of a particular type in this paper. But it is easy to deal with other types of frequency dependence which can be understood from the formulation given in Sec. 2. Second, we can treat three-dimensional (3D) photonic crystals as well as 2D ones. In that case, the spatial mesh for the discretization is 3D and the number of the data that should be dealt with in the numerical calculation increases by a considerable amount. But we believe that the CPU time is still allowable. Thirdly, we did not calculate the decay rate γ of the eigenmodes in this paper because we were mainly interested in the evaluation of the accuracy and efficiency of the present method as a tool for the photonic band calculation. γ is inevitably non-zero when the dielectric constant depends on frequency and has the non-zero imaginary part as a result of the Kramers-Kronig relation. γ can be evaluated by analysing the width of the resonance curve, or more practically by calculating the temporal evolution of the accumulated electromagnetic energy after switching off the oscillation of the dipole.

Figure 4. The distribution of the electric fields of the (a) B ₁(1), (b) A ₁(1), (c) A ₂, (d) B ₁(2), (e) A ₁(2), and (f) B ₂ modes at the X point that has the C _2v symmetry. The maximum of each electric field is normalized to unity as before.

Download Full Size | PDF

4. Conclusion

We formulated a new and general method to calculate the dispersion relations of arbitrary photonic crystals with frequency-dependent dielectric constants based on the numerical simulation of the radiation process of a virtual oscillating dipole moment embedded in the crystals. As an example, we apllied the present method to the 2D square lattice composed of metallic cylinders with the dielectric constant of Drude type. In addition to the dispersion relations that were consistent with the symmetry of the crystal, we could calculate the wave functions of the eigenmodes.

Acknowledgments

This work was financially supported by the Sumitomo Foundation and the Research Foundation for Opto-Science and Technology.

References

1. J. D. Joannopoulos, R. D. Meade, and J. N. Winn, Photonic Crystals, (Princeton University Press, Princeton, 1995).

2. Photonic Band Gaps and Localization, edited by C. M. Soukoulis (Plenum, New York, 1993).

3. Photonic Band Gap Materials, edited by C. M. Soukoulis (Kluwer, Dordrecht, 1996).

4. S. L. McCall, P. M. Platzman, R. Dalichaouch, D. Smith, and S. Schultz, “Microwave propagation in two-dimensional dielectric lattices,” Phys. Rev. Lett. 67, 2017–2020 (1991). [CrossRef] [PubMed]

5. E. Yablonovitch, T. J. Gmitter, R. D. Meade, A. M. Rappe, K. D. Brommer, and J. D. Joannopoulos, “Donor and acceptor modes in photonic band structure,” Phys. Rev. Lett. 67, 3380–3383 (1991). [CrossRef] [PubMed]

6. E. Yablonovitch, “Photonic band-gap structures,” J. Opt. Soc. Am. B 10, 283–295 (1993). [CrossRef]

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

8. K. Sakoda and H. Shiroma, “Numerical method for localized defect modes in photonic lattices,” Phys. Rev. B 56, 4830–4835 (1997). [CrossRef]

9. K. Sakoda, T. Ueta, and K. Ohtaka, “Numerical analysis of eigenmodes localized at line defects in photonic lattices,” Phys. Rev. B 56, 14 905–14 908 (1997). [CrossRef]

10. K. Sakoda and K. Ohtaka, “Optical response of three-dimensional photonic lattices: Solutions of inhomogeneous Maxwell’s equations and their applications,” Phys. Rev. B 54, 5732–5741 (1996). [CrossRef]

Novel approach to photonic bands with frequency-dependent dielectric constants

Abstract

1. Introduction

2. Theory

3. Results and discussion

4. Conclusion

Acknowledgments

References

Cited By

Figures (4)

Equations (15)

Optics Express