## Abstract

Two numerical techniques for analysis of defect modes in photonic crystals are presented. Based on the finite-difference time-domain method (FDTD), we use plane wave incidences and point sources for excitation and analysis. Using a total-field/scattered-field scheme, an ideal plane wave incident at different angles is implemented; defect modes are selectively excited and mode symmetries are probed. All modes can be excited by an incident plane wave along a non-symmetric direction of the crystal. Degenerate modes can also be differentiated using this method. A proper arrangement of point sources with positive and negative amplitudes in the cavity flexibly excites any chosen modes. Numerical simulations have verified these claims. Evolution of each defect mode is studied using spectral filtering. The quality factor of the defect mode is estimated based on the field decay. The far-field patterns are calculated and the Q values are shown to affect strongly the sharpness of these patterns. Animations of the near-fields of the defect modes are presented to give an intuitive image of their oscillating features.

© 2003 Optical Society of America

## 1. Introduction

It is well known that a point defect in photonic crystals can generate localized states in the band gap. The localized states can resonate in the microcavity with high Q values since they cannot escape into the surrounding bulk material. It is crucial to understand the nature of the defect modes since they are used in the design and application of many devices, such as high Q filters [1] and add/drop filters for multiplexer [2].

Numerical techniques for the study of these localized modes mainly fall into two categories: the time domain and frequency domain. The plane wave method (PWM) [3], which is the most popular frequency domain method, is able to find the eigen-frequencies and mode fields by using a supercell scheme. It assumes an infinite lattice with periodic defects, and the coupling of these defects leads to a defect band. The coupling effect between neighboring defects decays exponentially with the distance between defects, and a supercell of moderate size can give accurate information of the defect modes. However, the plane wave method cannot yield the information of their dynamics, Q values and far-field patterns, which are important for mode couplings.

Several techniques using FDTD were used to analyze the defect modes in micro-cavities. In Ref. [1], a plane wave was incident normal to the crystal and the transmission was monitored. Resonant modes were found as spectral peaks in the band gap. However, a plane wave incidence cannot excite those defect modes that are anti-symmetric to the k-vector of the plane wave [1]. Reference [4] used an initial field with low symmetry to excite all defect modes. By monitoring the field in various locations, defect modes can be found from the spectrum; however, the degenerate modes cannot be differentiated by this method. Dipole sources were proposed to excite the defect modes in Refs. [5–7], and it was implemented by adding an embedded dipole polarization field into the Maxwell’s equations. By choosing bandwidths, orientations and locations of the dipole sources, defect modes could be selectively excited.

In this paper, two approaches using the FDTD method are proposed to excite and analyze the defect modes in micro-cavities. One technique employs a plane wave incidence as the excitation. The total-field/scattered-field scheme [8] is used to implement an ideal plane wave incidence in any direction. This scheme has some advantages in photonic crystal simulations. It can reduce the computation region to contain just the area of interest; also, a plane wave incident at different angles can excite modes selectively due to the directional dependence of the defect modes. The second technique makes use of point sources at proper locations for excitation and there is no need to construct a dipole polarization field in Maxwell’s equations for implementation. Moreover, defect modes can be excited one or all at a time.

## 2. Methods of computation

A two-dimensional FDTD method for TM modes (which have only E_{z}, H_{x} and H_{y} components) based on Yee’s algorithm [9] is implemented for photonic crystal calculations, and the Perfectly Matched Layer (PML) boundary condition [10] is used to absorb outgoing waves efficiently. A total-field/scattered-field algorithm is deployed as in [8], which enables an ideal plane wave incidence at any angle from the interface of the total-field/scattered-field. In addition, ‘hard’ or ‘soft’ point sources are implemented as excitation sources. When the point sources are used, the total-field/scattered-field scheme is disabled.

Far-field patterns of defect modes are important to understand their mode coupling strengths and their directional dependences. A near-to-far field transformation [8] is implemented in the FDTD method to measure the far-field patterns for one or more frequencies in either total or scattered field region. One should be aware that the incident conditions or initial field distributions would affect the calculated far-field patterns, but the effects can be reduced by starting the transformation after the steady state is reached and excitation pulses are dissipated.

The results in this paper are obtained using a standard 2^{nd}-order FDTD scheme. A 4^{th} order FDTD scheme was also used to reduce accumulated numerical dispersions, but we found the 2^{nd}-order algorithm was accurate enough in this case.

As an example, a two-dimensional 5×5 supercell containing 25 GaAs rods in air is considered, and the setup is shown in Fig. 1. The dielectric constant of the GaAs rods is 11.56, and the radii of the regular rods and defect rod are 0.20*a* and 0.60*a*, respectively, where *a* is the lattice constant and *a*=1.0 µm. Each unit cell is divided into a 40×40 mesh and the whole FDTD region is 200×200 with Δx=Δy=0.025 µm. 10 PML layers are used outside the structure to absorb all energy flowing into the boundary.

Gaussian pulses covering the spectrum of the whole band gap are used: $s\left(t\right)=A\phantom{\rule{.2em}{0ex}}\mathrm{exp}\left\{-{\left(\frac{t-{t}_{0}}{{t}_{w}}\right)}^{2}\right\}\mathrm{sin}\left\{2\pi {f}_{0}\left(t-{t}_{0}\right)\right\}$

where *A* is the amplitude, t_{0}=1000Δt, t_{w}=250Δt, Δt=5.89e-11µs, f_{0}=0.35c/a and *c* is the speed of light in vacuum. The total number of time steps is 50,000. All pulsed sources will be removed after 5000 time steps.

The crystal has a band gap for TM modes from 0.29c/a to 0.42c/a and the micro-cavity supports four resonant modes as given in Refs. [1, 5] or they can be calculated using plane wave method as discussed in Refs. [3, 11]. Their frequencies and near-field symmetries are listed in Table 1. The symmetry property of the mode is described as ‘symmetry to x -symmetry to y’, thus ‘odd-even’ indicates that the mode is odd to the x-plane and even to y-plane. Special attention should be paid to the degenerate modes, which are the doubly degenerate hexapoles Mode 4(1) and Mode 4(2). Their mode patterns obtained using the frequency domain plane wave method (PWM) are odd-even and even-odd to ±45° [1], and they can also exist as their linear combination in the forms of constructive and destructive interferences. Therefore, the symmetries for these two modes can also be odd-even and even-odd as listed in Table 1.

Since the defect modes are well defined, detectors should be properly placed to obtain the information of all defect modes. In the calculation, detector 1 is placed at the center of the defect, detector 2 is placed in the direction of 45°, and detector 3 is placed on the x-axis outside the micro-cavity. According to the symmetries of these defect modes, detector 1 will only detect Mode 3 since other modes have no field at this point; detector 2 will not detect Mode 2 and detector 3 will not detect Mode 1 since they are located on the nodal plane of these modes.

## 3. Plane wave excitation

A plane wave incident at any non-symmetric direction of the crystal can have non-zero coupling to all defect modes, therefore it is able to excite all of them. On the contrary, a plane wave incident at some special symmetric directions of the crystal can have zero coupling with specific modes, therefore it will not excite them. This property can be used to reveal all defect modes and differentiate the degenerate modes in a cavity.

The space of interest is covered by the scattered-field region and is indicated as the dashed line in Fig. 1. A near-to-far-field transformation is performed along the dotted lines to measure the far-field pattern of the scattered field after the steady state is reached; in our case, it is specifically after 10,000 time steps. The plane wave pulse with a Gaussian shape is incident from the left at the total-field/scattered-field interface with an angle from 0° to 90° to the x-axis as shown in Fig. 1.

Four numerical simulations were carried out for incident angles 0°, 45°, 60° and 90°. Since the crystal has three symmetry directions along 0°, 45° and 90°, which are designated as Γ-x, x-M, M-Γ directions in its 1^{st} Brillouin zone, plane wave incidence along these directions will not excite specific modes. Plane wave along other directions, such as incidence at 60° in the example excites all defect modes. Since the structure is symmetric for 90° rotation, 0° and 90° plane wave incidence will yield identical spectral information. However, 0° plane wave incidence will not excite Mode 4(2) and 90° plane wave incidence will not excite Mode 4(1), hence, degenerate modes can be differentiated.

Fourier analysis of the detected fields can disclose the resonant frequencies. As an example, the spectra obtained by the three detectors using a plane wave incident at 60° are shown in Fig. 2. The resonant frequencies are in excellent agreement with the frequencies given in Table 1. Just as predicted above, not all modes are detected by a single detector.

Mode symmetries can be discerned from the four simulations. The results are listed in Table 2. As predicted above, some modes cannot be excited by a plane wave along a special direction. From Table 2, mode symmetry can be deduced. Taking Mode 1 as an example, we can see that it has zero coupling to plane waves along the directions of 0° and 90°, so we can conclude easily that Mode 1 is an odd-odd mode.

## 4. Point sources

Point sources are more flexible than the plane wave sources to selectively excite the defect modes since their k-vectors are two-dimensional. Proper arrangements of point sources with equal positive and negative amplitudes can form dipoles, quadrupoles and even more complex patterns. This method is able to excite all, some or any individual defect mode in the cavity. Table 3 shows the arrangement of the point sources to achieve excitation of all or any individual mode.

As in Table 3, a single point source in any non-symmetric direction can excite all modes. To excite a single mode only, point sources are arranged according to its mode symmetry, which is assumed to be known first. Since the defect modes in a cavity are orthogonal to each other [12], such an arrangement of point sources will produce an excitation with the same symmetry as the defect mode and orthogonal to other modes; therefore it will excite only the selected mode.

Five numerical simulations as in Table 3 are conducted and the resulting spectra are shown in Fig. 3. The results show that a simple arrangement of point sources with positive and negative amplitudes can selectively excite any combination of the defect modes. Note that due to the weak coupling, the peak for Mode 2 in the figure for excitation of all modes is small.

The spectrum of a defect mode has a Lorentzian shape. The non-Lorentzian oscillations in the spectra are produced by the transient effect, which depends on the excitation condition. Those frequency components which do not meet resonance condition will decay in a short time. Under some excitation conditions, the decay time could be longer, and it is reflected in the spectrum as oscillations. To reduce or eliminate these effects in the spectra, we can take Fourier transform after a certain time. For example, in the graph at the bottom right of Fig. 3, the Fourier transform is taken after 10,000 time steps; all the oscillations disappear, confirming that they are indeed due to the transient effect.

## 5. Far-field patterns

The far-field patterns can disclose their diffraction strength at each direction, which is related to the directional dependence of the coupling strength. The far-field amplitudes of the defect modes are calculated and shown in Fig. 4 using the plane wave incidence scheme. The point source scheme also gives similar results, but is not shown here. Note that the far-field pattern for Mode 4 using a plane wave incident at 45° will be the mix of the two degenerate modes and is actually the eigen-mode calculated using PWM method. Note also that their amplitudes along symmetric angles are not exactly equal because the calculated far-field patterns are still affected a little by the incidence condition.

The far-field is dependent upon the Q values, a higher Q values will have a more directional mode. As an example, a 1^{st} order monopole mode as in [1] produced by a microcavity with the center rod removed is excited by a point source at the center and shown in Fig. 5, with 3×3, 5×5 and 7×7 lattices, respectively. The Q values increase almost exponentially with the supercell size. When the Q value increases, the far-field pattern becomes more directional and the coupling to other modes gets weaker. The far-field pattern for 1^{st} order and 2^{nd} order monopole are similar since the far-field pattern only measures its outgoing diffraction strength.

## 6. Mode dynamics and their Q values

The defect mode evolution can be obtained from the field variation on the detectors by using a narrow bandwidth Butterworth digital filter. The results are shown in Fig. 6 for the plane wave incident at an angle of 60°. The evolution for mode 4(1), 4(2) and the mixed Mode 4 are similar and only one is calculated. The field in the defect is enhanced at the eigen-frequencies due to the coherent interference in the microcavity. After the steady state is reached and no external source is providing energy, the resonant mode decays exponentially with time. We can evaluate the decay rate and the Q values of the micro-cavity based on the field evolution.

The quality factor Q can be roughly estimated using the data from Fig. 6. According to Ref. [1], the energy in the micro-cavity decays in the form of
${e}^{-\left[\frac{{\omega}_{0}}{Q}\right]t}$
, so the E field amplitude will decay in a form of
${e}^{-\left[\frac{{\omega}_{0}}{2Q}\right]t}$
. Therefore, Q can be evaluated using
$Q=-\left[\frac{{\omega}_{0}\left({N}_{1}-{N}_{0}\right)\Delta t}{2\phantom{\rule{.2em}{0ex}}\mathrm{ln}\frac{{E}_{1}}{{E}_{0}}}\right]$
if we know the amplitudes *E*
_{1} and *E*
_{0} at time step *N*
_{1} and *N*
_{0}. Taking Mode 1 as an example, it reaches an amplitude *E*
_{0}=1.0 at *N*
_{0}=14000 time steps, and its amplitude decays to *E*
_{1}=0.40 at time step *N*
_{1}=50,000. This gives a Q value of 648. The Q values for Mode 2, 3 and 4 are estimated to be 276, 466 and 2936, respectively. The two degenerate modes and the mixed Mode 4 have similar Q values.

The near-field pattern obtained by PWM method is the field snapshot at a time spot and generally the phase is fixed at 0°. Time domain method can show the oscillations of the mode in real time. Animations for each defect mode are shown in Fig. 7. Each animation consists of 50 sequential field images with a time interval of 50 time steps; at one time spot each image records the Ez field components at all points in the whole space of study (200×200). In this case, the point source scheme is used.

## 7. Conclusions

In conclusion, we have presented two numerical techniques for excitation and analysis of the defect modes in two-dimensional cavities using plane wave incidence and point sources. The two schemes do not require incorporation of dipole sources in Maxwell’s equations. By using a plane wave incidence along different directions or arranging point sources with positive or negative amplitudes in a proper way, we have selectively excited the defect modes. The plane wave scheme cannot excite an arbitrarily chosen mode, however, it can be used to probe the symmetry of defect modes in the cavity. A point source scheme requires the knowledge of mode symmetry first, but it can excite any chosen modes. The dynamics of the defect modes inside the cavity have been observed using spectral filtering. Q values have been obtained from the decaying field. Far-field patterns have also been calculated simultaneously and higher Q values make the far-field patterns more direction dependent.

## Acknowledgement

This research is supported by NASA Langley Research Center through NASA-University Photonics Education and Research Consortium (NUPERC). S. Guo thanks his colleagues, Qin Hu, Feng Wu and Bing Xiao for many fruitful discussions.

## References and links

**1. **P R Villeneuve et al, “Microcavities in photonic crystals: mode symmetry, tunability, and coupling efficiency,” Phy. Rev. B **54**, 7837 (1996) [CrossRef]

**2. **Shanhui Fan et al, “Channel drop filters in photonic crystals,” Opt. Express, **3**, 4 (1998), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-3-1-4 [CrossRef]

**3. **K. M. Ho et al, “Existence of a photonic gap in periodic dielectric structures,” Phy. Rev. Lett. **65**, 3152 (1990) [CrossRef]

**4. **Min Qiu et al, “Numerical method for computing defect modes in two-dimensional photonic crystals with dielectric or metallic inclusions,” Phy. Rev. B **61**, 12871 (2000) [CrossRef]

**5. **Kazuaki Sakoda et al, “Optical response of three-dimensional photonic lattices: solution of inhomgeneous Maxwell’s equations and their applications,” Phy. Rev. B **54**, 5732 (1996) [CrossRef]

**6. **Vladimir Kuzmiak et al, “Localized defect modes in a two-dimensional triangular photonic crystal,” Phy. Rev. B **57**, 15242 (1998) [CrossRef]

**7. **Kazuaki Sakoda et al, “Numerical method for localized defect modes in photonic lattices,” Phy. Rev. B **56**, 4830 (1997) [CrossRef]

**8. **Allen Taflove, *Computational electrodynamics, the finite difference time domain method* (Artech House1995)

**9. **K. S Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antenna and Propagation **14**, 302 (1966) [CrossRef]

**10. **J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Computational Physics, 185, (1994) [CrossRef]

**11. **S. Guo and S. Albin, “Simple plane wave implementation for photonic crystal calculations,” Opt. Express **11**, 167 (2003), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-11-2-167 [CrossRef] [PubMed]

**12. **J. D. Joannopoulos et al, *Photonic crystals - Molding the flow of light* (Princeton University Press1995)