## Abstract

We present a method for directly analyzing photonic nano-devices and apply it to photonic crystal cavities. Two-dimensional photonic crystals are scanned and reproduced in computer memory for Finite Difference Time Domain simuations. The results closely match experimental observations, with a fidelity far beyond that for idealized structures. This analysis allows close examination of error mechanisms and analytical error models.

© 2006 Optical Society of America

## 1. Introduction

Photonic crystals (PCs) allow exquisite control over light on small length scales, making them useful in a wide range of applications. The devices are typically designed in the following sequence: intuition-guided rough design, numerical optimization and simulation, fabrication, and testing. Discrepancies between tested and predicted device performance are usually attributed to fabrication inaccuracies. Several models have been proposed to explain these inaccuracies [1, 2, 3], but to date no direct evaluation of the effect of these deviations exists.

In this paper, we turn the problem around and simulate fabricated structures. Although the approach is general, we focus on PC cavities because their fields are highly sensitive to device errors. The PC cavities were digitized by Scanning Electron Microscopy (SEM) and their electromagnetic fields calculated by Finite Difference Time Domain (FDTD) simulations. The results predict the actual photonic properties of the PC structure far more accurately than results for idealized structures. This direct approach approximates very closely the true mode patterns inside fabricated structures and reveals the dominant loss mechanisms. In particular, we show that higher losses in our devices result primarily from global hole radius deviations that change the resonance condition. Edge roughness and material loss play only secondary roles in our structures, but they nevertheless do cause significant losses that need to be considered in device design. In addition, our analysis gives good estimates on what discritization is required in simulations to capture the over-riding device behavior. This analysis of electromagnetic fields in real structures is very general and could be of use in areas ranging from photonic crystal cavities to diffraction by butterfly wings[4].

## 2. Photonic crystal cavities

We now analyze the electromagnetic modes in fabricated planar PC cavities. Light is confined in the in-plane direction by distributed Bragg reflection (DBR), and in the vertical direction by total internal reflection (TIR). Imperfect confinement leads to losses in these directions, denoted by the rates *γ*
_{∥} and *γ*
_{⊥}, respectively. These losses are inherent to the PC cavity design for the particular mode, so we will refer to them as intrinsic. These losses remain even if the structure were fabricated perfectly. In contrast, extrinsic losses are comprised of material absorption (with rate *γ*_{M}
) and scattering from fabrication imperfections (with rate *γ*_{S}
). The total loss rate is expressed as the sum,

We can express this in terms of the quality factor *Q* = *ω* 〈*E*_{mode}
〉 / 〈*P*〉, which measures the mode confinement as the fractional energy loss per resonator cycle, where *ω* is the confined mode frequency and 〈*P*〉 and 〈*E*_{mode}
〉 are the averaged radiated power and mode energy. Eq.1 then becomes

where we used *γ* = 2*ω*/*Q* and defined *Q*_{M}
= 4*πn*/*λα*_{M}
, denoting the commonly used optical absorption coefficient *α*_{M}
= *γ*_{M}
/(*c*/*n*), with *c*/*n* the speed of light in the material with refractive index *n*.

To analyze these losses by FDTD in 3D, consider first the structure of Fig. 1(a, insent). This is a single-defect cavity with small mode volume *V*_{mode}
= 0.5(*λ*/*n*)^{3}, modified to maximize the quality factor *Q* for the *x*-dipole mode (Fig. 1(d)). At the design stage, only intrinsic losses were considered. The optimal structure was obtained by changing the six holes neighboring the defect, yielding *Q*
_{∥} = 58,000 and *Q*
_{⊥} = 45,000 for structures with seven hole rings around the defect [5]. Since *Q*
_{∥} can always be increased arbitrarily high by addition of hole rings, the total *Q* is limited by *Q*
_{⊥} — in this case, to *Q* ~ 45,000.

This cavity design was fabricated on a GaAs substrate, for a target wavelength of *λ*
_{0} = 960 nm and lattice periodicity *a* = 256nm. The design was written with ten hole rings, a number sufficient to guarantee *Q*
_{⊥} ≪ *Q*
_{∥} and *Q* ~ *Q*
_{⊥}. The fabrication consisted of electron-beam lithography (step size 1.5 nm) on Poly-methyl-methacrylate (PMMA) resist, development of the resist, and pattern transfer by dry etching. Subsequently, a sacrificial Al_{0.94}Ga_{0.06}As layer was removed (by wet etching) to create the PC membrane. The PC membrane also contained a central layer of self-assembled InAs quantum dots (QDs) with density 200/*μ*m^{2}, which served as a convenient broad-band internal illumination source. The absorption due to the QD layer was negligible in comparison to other loss mechanisms, as we will confirm later.

The cavity resonances were measured using a confocal microscope setup, as described in [6]. The pump laser excited the QDs at 750 nm (above the GaAs bandgap) with a Ti-Saph laser. When the QDs were pumped at high intensity (~ 1.0 kW/cm^{2}), the QD emission broadened across a wide spectrum (900-950 nm), revealing the Lorentzian lineshape of the cavity as QD emission sees both higher coupling to PC leaky modes, and is emitted faster, as given by the Purcell effect (see Appendix).

In Fig. 1(c), we show the measured cavity resonance of structure S1 (Fig. 1(a)) based on the design discussed above. The polarization dependence confirms the dipole-nature of this resonance. However, the cavity’s quality factor, *Q*
_{s1} ~ 2512, falls well short of the predicted value of 45,000, and the resonance at frequency *a*/*λ*
_{s1} = 0.290 differs significantly from the target wavelength of *a*/*λ*
_{0} = 0.2847, where *a* is the lattice periodicity and *λ* the wavelength in air.

To better understand these discrepancies between prediction and experiment, we simulated the fabricated structure. This was done by converting the SEM scan (Fig. 1(a)) into a three-dimensional dielectric structure in computer memory, applying a threshold that distinguishes between high-dielectric material (‘bright’), and vacuum (‘dark’). The placement of the threshold must be calibrated for a set of structures which were imaged under identical conditions in the SEM. We will discuss this issue in more depth later. The computational structure is symmetric normal to the PC plane, matching the SEM profile, and topped by an air layer with thickness 0.75a. Then Maxwell’s equations are propagated on the space using FDTD. The result is shown in the inset of Fig. 1(d) as the ∣*B*_{z}
∣ component of the confined field at the slab center. From the energy loss relation *Q* = ω〈*E*_{mode}
〉/〈*P*〉, the quality factor is estimated as ${Q}_{S1}^{\mathit{\text{sim}}}$ = 2861, closely matching the experimentally measured value of 2512. We were also interested in the cavity resonance spectrum. As mentioned earlier, the large simulation space prohibits a detailed spectrum for a 3D simulation. A 2D simulation matches the spectrum qualitatively, although the cavity linewidths do not reliably predict the *Q* factor as vertical losses are not accounted for. In this calculation, the 3D structure is approximated by an “effective” 2D structure by calculating an effective modal index *n*_{eff}
[7]. The result is shown in Fig. 1(d). The positions of the resonances match the experimental spectrum. The simulation also replicates other low-*Q* resonances (e.g., see small unpolarized peak just right of the dipole peak). The *Q* values for these resonances were obtained from 3D simulations to account for the dominant out-of-plane loss.

The fabrication error of the photonic crystal, shown in Fig. 1(b) as the difference between ideal and fabricated structures, elucidates the reason for the resonance frequency shift of the experimental structure with respect to the design. The fabricated holes are too large (by ~ 11% over the average design diameter of 154 nm) and also deviate from the ideal circle near the cavity due to proximity effects in electron-beam exposure. The mode now overlaps more with the air holes, resulting in a resonance blue-shift and diminished vertical confinement.

We can use this example to test the applicability of first-order perturbation theory in such error analysis. If *E*⃗ denotes the electric field of the resonant mode at frequency *ω* of the ideal PC structure *ε*, then the small change Δ*ω* is given by the first-order perturbation

This expression is well-defined since in the simulations, *ε* is averaged over neighboring cells, so that *E*⃗ can be approximated as continuously varying[8]. We calculate Δ*λ*/*λ*
_{0} ~ -0.024, or Δ*λ* ~ -23 nm, close to the actual difference of *λ*_{meas}
- *λ*_{design}
= 943 - 960 nm= -17 nm.

From a practical standpoint, ensuring that all holes meet specification is difficult. One of the main obstacles results from the proximity effect in electron beam lithography: scattered electrons expose nearby areas, so the patterns in the proximity of a region affect its electron dose. In principle, this proximity effect can be corrected by local modification of the exposure dose, but we found this difficult for these non-periodic photonic crystal structures. The SEM-based simulation suggests an easier alternative: modify the design to have more uniformly sized holes.

## 2.1. Simulating Electromagnetic Fields in Fabricated Structures

Following this idea, the best design is shown in Fig. 2(a) and has *Q* = 46,000 and mode volume *V*_{mode}
= 0.5(*λ*/*n*)^{3}. This design was optimized using a parametric search on the positions of the cavity’s six nearest-neighbor holes, constraining the hole diameters to be within 10% of the unchanged values. Over 300 simulation parameter points were calculated using an automated, iterative optimization procedure; in Fig. 2(b), we show just two slices in parameter space where the four holes in the ΓJ directions were shifted by Δ*x*, Δ*y*, and their radii varied. The design that results in fabricated cavities with highest *Q* is one that has reasonably high theoretical *Q* and is tolerant to errors in the hole radii and positions. For the dipole cavity mode, we identify the best design as the one shown in Fig. 2(a).

The fabrication-friendly design of the single-defect cavity led to structures with moderately high *Q* factors that are easy to fabricate. We will now analyze a set of seven instances of this cavity design. In Fig. 3(b), we show one such fabricated cavity (S2) with its measured PL spectrum (a). We analyzed this and the other structures in the following way.

First, the SEM picture was converted into the dielectric structure shown in Fig. 3(b). The scanning electron micrographs were acquired with a resolution of 12 nm and recorded at 50 pixels per lattice parameter *a* = 256 nm. In 2D - FDTD simulations, the image was directly converted to a dielectric with 50 program points per *a*; in 3D simuations, the SEM images were averaged down to 25 pts/*a*. In both cases, the native resolution limit in the SEM scans requires a calibration for SEM images for differentiating between air and dieletric. This calibration is found from the pixel intensity histogram (Fig. 4(b)) of structure S2 (scaled from 0:black to 1:white), where a threshold value in intensity distinguishes between GaAs and vacuum. This threshold is varied until simulated and observed resonances match (see Fig. 4(a)). For the calibration structure S2, best agreement is found with the threshold at 0.55. Note that the *Q* values are roughly unchanged and are near the actual value of 2100, except for the lowest threshold value where a lower *Q* results because the cavity resonance cames close to the bandgap edge. In that case, in-plane confinement suffers and *Q* drops.

Another important consideration in digitizing the structures is the discretization, which should be low for computational ease, while still guaranteeing to capture the important features for *Q* and *λ*
_{0}. Simulations of S2 showed that discretizing periods into more than 25 points (or ~ 10 nm) did not considerably change the simulation outcomes.

Having digitized cavity S2 with these conditions for structure conversion, we simulated the fields in two ways. We used a 2D-FDTD simulation with a large number of time steps (2^{18} ~ 260,000) so that the simulated spectrum reliably reproduced the measured one (Fig. 3(c)). To obtain cavity loss rates and filter accurate mode patterns, the simulation was repeated in 3D with fewer time steps (2^{12}). This yielded *Q* = 2434 for the dipole mode (Fig. 3(d)).

We repeated this analysis for six other structures on the same sample, using the calibration data from structure S2. The structures followed the same design as S2, except that the electron-beam exposure was incremented in steps of 2%, resulting in larger holes and hence higher resonance frequencies. The resonances and *Q* values are summarized in Fig. 4(b,c). All experimental resonances matched the directly simulated ones within 2%, far better than the idealized structure simulations (Fig. 4(b)). For *Q* values, the SEM-based simulations show this improvement even more dramatically (Fig. 4(c)).

## 3. Analyzing loss in photonic crystal cavities

Having established the close agreement between experiment and simulation, we can now use the simulated fields to analyze mode losses described by Eq.1.

#### 3.1. Material Losses

The simulations of scanned structures discussed above take no material losses into account, i.e., the full loss term in the simulation is given by *γ*_{sim}
= *γ*
_{∥} + *γ*
_{⊥} + *γ*_{S}
or ${Q}_{\mathit{\text{sim}}}^{-1}$ = ${Q}_{\parallel}^{-1}$ + ${Q}_{\perp}^{-1}$ + ${Q}_{S}^{-1}$. Thus, the material losses are proportional to ${Q}_{\mathit{\text{measured}}}^{-1}$ - ${Q}_{\mathit{\text{sim}}}^{-1}$; by the simulation results, these losses are insignificant. In passing, we note also that this justifies our earlier assumption that the QDs do not influence the total *Q*.

Even though material losses play an insignificant role in our cavities, it is interesting to approximate how they limit the best possible cavity *Q* and *Q*/*V*. To answer this question, we return to an earlier analysis that considered the *Q* of a near-optimal Gaussian cavity field pattern [9]. A basic result from was that a larger structure allows lower uncertainty in the *k*-space distribution of the cavity mode, which results in lower overlap with the lossy light-cone where TIR does not confine light in the vertical direction. Without extrinsic losses, it was found that *Q* and *Q*/*V* increased roughly exponentially with cavity extent *σ*_{x}
, as shown in Fig. 5(a). If losses are taken into account according to Eq.2, then *Q* asymptotically approaches *Q*_{M}
. Interestingly, this implies that *Q*/*V* has an absolute maximum. Fig. 5(b) shows plots of *Q* and *Q*/*V* vs. cavity extent *σ*_{x}
for pure GaAs at 900 nm (below bandgap) at room temperature, where *α* = 10/cm [10]. *Q*/*V* reaches a maximum at a cavity length of ~ 1.6*a*, where *a* is the PC’s lattice periodicity. Thus, in designing a PC structure to exploit cavity quantum electrodynamic effects, the optimal structure should be only a few lattice periods in size (the optimal length has a roughly logarithmic dependence on the extrinsic loss and is therefore rather insensitive to the exact value of the loss).

#### 3.2. Extrinsic Scattering Losses

Extrinsic scattering losses result from disorder-induced scattering, due mostly to surface irregularities of the hole sidewalls. Since we have the exact field pattern, we can analyze these losses directly. In Fig. 6, we analyze the scattering loss for scanned structure S3. This structure has a sharp defect near the cavity center. The first column of Fig. 6 shows this structure, together with computed field (row 1), its Fourier transform (row 2), and the vertical component of the Poynting vector, *S*_{z}
, just above the membrane surface. We compare these results against a simulation where the defect was removed in the image (i.e., the hole is rounded out), as shown in the second column. The differences between true and ‘corrected’ structures are displayed in the third column. From the field pattern, we see that the defect shifts the cavity mode slightly (top row). The defect also causes considerable scattering, as seen from the field’s Fourier transform in the second column. Clearly visible are regions of high scatter where waves are re-directed upward and lost out of the cavity (bottom row, right). From the field patterns, we calculate the difference in the loss *S*_{z}
, as shown in Fig. 6(bottom right). By integrating this loss, we find that removing the defect causes a 4.3% decrease in losses, corresponding to an increase of *Q* of ~ 4.5%. This calculation is supported by the simulation of the ‘polished’ structure, which predicts an increase in *Q* of 5.7%. This increase in *Q* is accompanied by a frequency increase of Δ(*a*/*λ*) ~ 4.4 · 10^{-4}.

The structures analyzed here have low surface irregularities, typically far below those of structure S3. Since even in S3, surface roughness only makes up < 10% of losses, we believe that extrinsic scattering losses are not responsible for the drop in *Q* compared to ideal simulated structures. The largest drop in *Q* arises from global deviations from the design: the structures considered here have systematic radius deviations on the order of 10% and hole center deviations on the order of 5%. This deviation is enough to change cavity *Q* by up to five times, according to the parametric search (Fig. 2(b)). These deviations cause increased vertical losses due to imperfect TIR [9].

#### 3.3. Intrinsic Scattering Losses

Based on our observations on intrinsic scattering losses, we would like to briefly examine if smoothing the dielectric structure can improve *Q* values. Intrinsic scattering losses are unavoidable because the confinement in a 2D PC is not perfect: the incomplete TIR causes out-of-plane losses in direction *k*⃗ = (*k*_{x}
, *k*_{y}
, *k*_{z}
) that are described by [9]

where $\eta =\sqrt{\frac{{\mu}_{0}}{{\epsilon}_{0}}}$. Integrated over the light cone, this expression gives the total vertical loss. To reduce intrinsic losses, it it therefore important to limit the *k*- components of *B*_{z}
and *E*_{z}
, evaluated in a plane just above the PC slab. One way to reduce these losses is to reduce scattering of high *k* components back into the light cone. This can be seen by first considering the infinite PC crystal. The Fourier-transform of the wave-equation is [11]

for all *k*⃗, where we expanded the dielectric as

This infinite set of equations for all *k*⃗ can be divided into many subsets, for each wave vector *K*⃗, and containing terms of form *A*⃗ (*K*⃗) and *A*⃗(*K*⃗ - *G*⃗) with all reciprocal lattice vectors *G*⃗. Now note that if *ε*(*r*⃗) consists of only low Fourier components, then high Fourier- components of the field will be decoupled and eliminated. For example, if *ε*(*r*⃗) = cos(*xπ*/*a*), then each subset of equations for a particular *K*⃗ consists of only three equations, for *K* = 0, *K* = ±*π*/*a*.

This picture suggests that smoothing of the PC structure should result in less scattering of high-frequency *k*- components into the light cone, decreasing vertical losses. We tested this with FDTD. The smoothing was achieved by applying a filtering function *F*_{S}
(*l*), where *l* is the FWHM of a Gaussian convolving envelope. In Fig. 7, we compare the field patterns of a single defect cavity with and without smoothing, where smoothing is applied with *l* = 0.2*a*. The field is shown as ∣*FT*
_{2}(*H*_{z}
)∣ + ∣*FT*
_{2}(*E*_{z}
)∣ (in logarithmic scaling). It is clear that for the smoothed structure, high *k* components are sharply diminished.

The smoothed structure results in higher air overlap of the confined donor-type cavity resonance, shifting it to higher frequency. This changes the resonance condition and *Q* value. Thus, to show that smoothing truly results in better confinement, the smoothed structure must be re-optimized by changing hole radii and positions. Though we do not pursue this here, we note that this should lead to a decrease of intrinsic losses in PC cavities.

Smooth variation of the dielectric index could be experimentally obtained in several ways. For one, digital chemical etching[12] could be used to round the hole-walls, causing a more gradual transition in the effective dielectric index. A truly sinusoidal transition could be obtained using, e.g., porous silicon, where the position-dependent index of refraction can be controlled by a sinusoidal application of current[13]. Such a structure would seem particularly interesting for parametric up-and down-conversion, where coupling between only few frequencies is desired.

## 4. Conclusion

In conclusion, we have demonstrated a novel method for analyzing discrepancies between designed and fabricated photonic structures. The method makes few idealizations; the planar profile is derived directly from SEM images of the structure, and is translated in the third dimension by the known depth of the structure. With more refined imaging techniques, the third dimension may be scanned separately as well. The calibration from SEM image to dielectric structure in the computer is straightforward and repeatable, and we showed that discretization on the order of 10 nm gives good results for the PC structures considered here. The analysis of real structures allows direct examination of fabrication errors and their effects on the electromagnetic modes. In addition, by linking the experimentally observed results with simulations of the same structure, this method allows precise evaluation of analytical error models.

The authors would like to thank Stanford undergraduate Tim Chang for help with the image processing and analysis. This work was supported by the MURI Center for photonic quantum information systems (ARO/ARDA Program DAAD19-03-1-0199). Dirk Englund was also supported by the NDSEG fellowship.

#### A. Measuring the cavity resonance from Purcell-modified spontaneous emission

Cavity-coupled emission is enhanced over decoupled emission in two ways. Firstly, by the increased emission rate of dipoles *μ*⃗ (at position *r*⃗ and emitting at wavelength *λ*) by the Purcell factor

with ${F}_{\mathit{cav}0}\equiv \frac{3}{4{\pi}^{2}}\frac{{\lambda}_{\mathit{cav}}^{3}}{{n}^{3}}\frac{Q}{{V}_{\mathit{mode}}}$. Secondly, it is increased by enhanced out-coupling efficiency *η*_{cav}
to the objective lens, due to the cavity mode’s coupling to photonic crystal leaky modes[6]. In addition to the background signal from uncoupled QD emission, there is therefore a signal proportional to *η*_{cav}
${\int}_{0}^{2\pi}$
*dθ* ∫_{A}
*d*
^{2}
*r*⃗*ρ* (*r*⃗,*λ*,*μ*⃗) *F*_{cav}
(*λ*,*r*⃗,*μ*⃗), integrated over the observation region A and dipole angle *θ* with respect to the cavity, where *ρ*(*r*⃗,*λ*,*μ*⃗) describes the density of QD emission over position, wavelength, and dipole orientation. In this expression, we assumed a continuous laser pump with intensity far above the saturation intensity of discrete lines, as mentioned in the text. At high pump power, *ρ* broadens to *ρ*(*r*⃗,*λ*,*μ*⃗) = *ρ*
_{1}(*r*⃗)*ρ*
_{2}(*λ*)*ρ*
_{3}(*μ*) ∝ = *ρ*
_{2}(*λ*) as the other components are roughly constant over position and polarization. Thus, the PL is proportional to *η*_{cav}*ρ*(*λ*)*F*_{cav}
(*λ*), and provided *ρ*(*λ*) is smooth on the order of the cavity linewidth, the PL follows a Lorentzian lineshape that gives the cavity *Q* directly.

## References and links

**1. **S. G. Johnson, M. I. Povinelli, M. Soljacic, A. Karalis, S. Jacobs, and J. D. Joannopoulos, “Roughness losses and volume-current methods in photonic-crystal waveguides,” Appl. Phys. B: Lasers and Optics **81**, 283 – 293 (2005). [CrossRef]

**2. **S. Hughes, L. Ramunno, J. F. Young, and J. E. Sipe, “Extrinsic optical scattering loss in photonic Crystal waveguides: role of fabrication disorder and photon Group velocity.” Phys. Rev. Lett. **94**, 033,903 – 4 (2005). [CrossRef]

**3. **A. F. Koenderink, A. Lagendijk, and W. L. Vos, “Optical extinction due to intrinsic structural variations of photonic crystals,” Phys. Rev. B. **72**, 153,102 – (2005). [CrossRef]

**4. **P. Vukusic and J. R. Sambles, “Photonic structures in biology,” Nature **424**, 852–55 (2003). [CrossRef] [PubMed]

**5. **J. Vuckovic and Y. Yamamoto, “Photonic crystal microcavities for cavity quantum electrodynamics with a single quantum dot.” Appl. Phys. Lett. **82**, 2374 – 6 (2003). [CrossRef]

**6. **D. Englund, D. Fattal, E. Waks, G. Solomon, B. Zhang, T. Nakaoka, Y. Arakawa, Y. Yamamoto, and J. J. VučkoviĆ “Controlling the Spontaneous Emission Rate of Single Quantum Dots in a Two-Dimensional Photonic Crystal,” Phys. Rev. Lett. **95**(013904) (2005). [CrossRef] [PubMed]

**7. **K. Kawano and T. Kitoh, *Introduction to Optical Waveguide Analysis: Solving Maxwells Equation and the Schrödinger Equation* (Wiley-Interscience Publications, New York, 2001). [PubMed]

**8. **S. G. Johnson, M. Ibanescu, M. A. Skorobogatiy, O. Weisberg, J. D. Joannopoulos, and Y. Fink, “Perturbation theory for Maxwells equations with shifting material boundaries,” Physical Review E **65**(066611) (2002). [CrossRef]

**9. **D. Englund, I. Fushman, and J. VučkoviĆ, “General Recipe for Designing Photonic Crystal Cavities,” Opt. Express **12**, 5961–75 (2005). [CrossRef]

**10. **J. H.C. Casey, D. D. Sell, and K. W. Wecht, “,” Journal of Applied Physics **46**, 250 (1975). [CrossRef]

**11. **A. Yariv and P. Yeh, *Optical Waves in Crystals* (Wiley Interscience, 2003).

**12. **G. C. D.*et al*., “Wet Chemical Digital Etching of GaAs at Room Temperature,” J. Electrochemical Soc. **143**, 3652–56 (1996). [CrossRef]

**13. **J. Schilling, F. Muller, S. Matthias, R. B. Wehrspohn, U. Gosele, and K. Busch, “Three-dimensional photonic crystals based on macroporous silicon with modulated pore diameter,” Appl. Phys. Lett. **78**, 1180–2 (2001). [CrossRef]