## Abstract

We present a first-principles method to compute radiation properties of ultra-high quality factor photonic crystal cavities. Our Frequency-domain Approach for Radiation (FAR) can compute the far-field radiation pattern and quality factor of cavity modes ∼ 100 times more rapidly than conventional finite-difference time domain calculations. We explain how the radiation pattern depends on the perturbation used to create the cavity and on the Bloch modes of the photonic crystal.

© 2012 Optical Society of America

The high quality factors (*Q*) and small modal volumes of photonic crystal (PC) cavities make them ideally suited for applications requiring strong optical field enhancement, such as low-energy optical switching [1], strongly coupled cavity quantum electrodynamics (QED) [2] and harmonic generation [3]. Recent interest has focused on high-*Q* cavities where the far-field radiation pattern is engineered to emit vertically, enabling free-space mode excitation [4–6] in cavity QED [5] and harmonic generation [7] experiments.

Photonic crystal cavity design uses established theoretical ideas [8–12] to maximize quality factors, in conjunction with finite difference time domain (FDTD) calculations to compute the cavity mode. Even with improvements in speed and accuracy [13,14], time domain calculations are by nature computationally intensive for the long life-times of ultra-high *Q* cavities, taking hours or days per design on a supercomputer. Optimizing both the quality factor and the radiation pattern can require the exploration of a large parameter space [5, 6], further increasing the computation effort. These severe computational demands, and the inability of FDTD to provide insight into the underlying physics, point to the need for an alternative method. Here we provide such a method. Our first-principles Frequency-domain Approach for Radiation (FAR) is ∼100 times more efficient than FDTD calculations since we do not compute radiative modes directly. It consists of two parts: we initially approximate the cavity mode as a bound mode after which the radiation is obtained using perturbation theory, thus avoiding the most time-consuming part of FDTD calculations. The FAR also provides a design strategy for achieving cavities with specified radiation patterns without requiring exhaustive simulations.

We apply the FAR to double heterostructure photonic crystal cavities [10], which are formed by perturbing a photonic crystal waveguide (PCW) in a slab geometry. The two geometries are shown in Fig. 1(a)–(b); in the *photosensitive cavity* (Fig. 1(a)), the refractive index of a strip around the PCW (yellow shading) is uniformly increased by Δ*n _{p}*, as can be achieved in chalcogenide glass [15, 16]. In the

*fluid infiltrated cavity*(Fig. 1(b)), the refractive index of the holes is increased by Δ

*n*in a strip-like region (red shading), typically by fluid infiltration [17, 18]. These two cavities are therefore complementary,

_{i}*i.e.*, in one only the background is perturbed, while in the other only the holes. We have found that considerable qualitative insight into the radiation pattern of the cavity mode can be obtained by examining a single term in the equation that governs the radiation from the cavity. This term has the form

*Ã*(

**r**)

**D**

*(*

^{a}**r**), where

*Ã*(

**r**) is associated with the perturbation that creates the cavity and

**D**

*(*

^{a}**r**) is the bound approximation for the cavity mode. This term is shown for a

*z*= 0 slice through the PC slab in Figs. 1(c)–1(d). The perturbation term

*Ã*(

**r**) is only non-zero in the background for the photosensitive cavity (Fig. 1(c)) and only non-zero in the holes for the fluid infiltrated cavity (Fig. 1(d)). The Fourier components in the light cone of this product are peaked near the edge of the light cone for the photosensitive cavity (Fig. 1(e)) corresponding to radiation being directed towards the horizon. However, for the fluid infiltrated cavity, the Fourier transform is strongest near

*k*=

_{x}*k*= 0 (Fig. 1(f)), and thus it predominantly radiates vertically.

_{y}Our theory uses a Hamiltonian formulation [19] to construct cavity modes by superposing a basis of bound PCW modes expressed in terms of the **B**(**r**) and **D**(**r**) fields, so any superposition is divergence-free. The Hamiltonian for a dielectric PC cavity with relative permittivity *ε*(**r**) is

*ε*(

**r**) =

*ε*̄(

**r**) +

*ε*̃(

**r**), where

*ε*̄(

**r**) is the permittivity of the PCW, and

*ε*̃(

**r**) the small permittivity change that creates the cavity. We then expand the cavity mode using the normalized PCW modes [19] below the light cone

*E*(

_{y}**r**) is even in

*y*. We can include more modes, but ultra-high

*Q*cavity modes are typically gently confined and thus different bands couple weakly.

Substituting (2) into (1) we obtain an approximation for the Hamiltonian of the PC cavity

*γ*(

**r**) = 1/(2

*ε*

_{0})[1/

*ε*(

**r**) − 1/

*ε*̄(

**r**)], and we dropped non-rotating wave terms involving ${a}_{k}^{\u2020}{a}_{{k}^{\prime}}^{\u2020}$ and

*a*

_{k}a_{k}_{′}. Diagonalizing ℋ

_{1}determines an eigenvalue, the energy

*h*̄

*ω*

_{0}of a photon in the cavity mode, while its eigenfunction

*v*

_{0}(

*k*) gives the cavity mode in the basis of PCW modes:

*Q*factor cavities of interest here is small, and we have found that

**D**

*(*

^{a}**r**) is a good approximation for the shape of the cavity mode. Similarly, the eigenvalues of (3) approximate the real part of the frequency of the cavity mode well. We thus use

**D**

*(*

^{a}**r**) to find a first approximation for the polarization field

**P**(

**r**) within the light cone.

The polarization field **P**(**r**) = *ε*_{0} [*ε*(**r**) − 1]**E**(**r**) of a mode with frequency *ω* satisfying the macroscopic Maxwell equations is also a solution to the integral equation

**r**′ due to an oscillating polarization source at

**r**. We use the formalism for layered media [20], in which we deal with a sheet of polarization. Since we need to compute the out-of-plane (

*z*-direction) radiation of a PC cavity, this formalism is particularly useful as it separates propagating modes in the

*z*-direction, with ${\left|\mathit{\kappa}\right|}^{2}\equiv {k}_{x}^{2}+{k}_{y}^{2}\le {k}_{0}^{2}$, from evanescent modes with ${\left|\mathit{\kappa}\right|}^{2}>{k}_{0}^{2}$, where

*k*

_{0}=

*ω*

_{0}/

*c*.

Defining Γ(**r**) = [*ε*(**r**) − 1]*/ε*(**r**), the polarization field of the cavity is approximated by **P*** ^{a}*(

**r**) = Γ(

**r**)

**D**

*(*

^{a}**r**) ≡ (Γ̄(

**r**) + Γ̃(

**r**))

**D**

*(*

^{a}**r**), where again the over-bar denotes a quantity for the PCW and the tilde denotes the perturbation creating the cavity. Since

**D**

*(*

^{a}**r**) has no Fourier components within the light cone, neither does Γ̄(

**r**)

**D**

*(*

^{a}**r**), Γ̄(

**r**) being periodic with the lattice. However, Γ̃(

**r**)

**D**

*(*

^{a}**r**) does have components within the light cone, providing a starting point for calculating the radiative polarization. We relate the actual polarization field of the cavity mode

**P**(

**r**) to

**P**

*(*

^{a}**r**) by writing

**P**(

**r**) =

**P**

*(*

^{a}**r**) +

**P**

*(*

_{c}**r**), whose radiative components are ${\mathbf{P}}^{\text{rad}}(\mathbf{r})=\tilde{\mathrm{\Gamma}}(\mathbf{r}){\mathbf{D}}^{a}(\mathbf{r})+{\mathbf{P}}_{c}^{\text{rad}}(\mathbf{r})$, where

**P**

*(*

_{c}**r**) is the correction to the polarization field, while

*rad*refers only to Fourier components in the light cone. We write the complex cavity mode frequency

*ω*as

*ω*=

*ω*

_{0}+

*ω*̃. We next perform a Taylor expansion about

*ω*

_{0}of the Green function, substitute into (5), and use the fact that ${\mathbf{P}}_{c}^{\text{rad}}(\mathbf{r})$ and the variables with tildes are small. After some manipulation and keeping only terms with Fourier components in the light cone, we obtain a first order expression for

**P**

^{rad}

*Ã*(

**r**)

**D**

*(*

^{a}**r**) gives good qualitative insight into the far-field radiation. In general though, Eq. (6) is a Fredholm integral equation of the second kind, in which the Green function ensures a self-consistent interaction between the dipoles.

By solving (6), we obtain the full quantitative radiative polarization components of the cavity mode, from which the far-field radiation can be determined using the Green function in (5). In the far-field we write the electric field as **E**_{far}(**r**) = **e**_{±}(** κ**̄)

*e*

^{ik0r}/

*r*, where

**̄ ≡**

*κ**k*

_{0}

**r̂**· (

**x̂x̂**+

**ŷŷ**), with, above (+) and below (−) the slab,

*s*polarization, and with a similar expression for

*p*polarization. Here

**R**= (

*x*,

*y*) and $w={\left({k}_{0}^{2}-{\left|\mathit{\kappa}\right|}^{2}\right)}^{1/2}$. Equation (7) is thus a planar (

*x*and

*y*) Fourier Transform, integrated over the thickness of the slab (

*z*) with appropriate phases. Each (

*k*,

_{x}*k*) of the polarization field inside the light cone corresponds to a unique far-field direction. The far-field electric field gives the Poynting vector, and therefore the quality factor of the cavity mode can be computed.

_{y}To obtain numerical solutions to (6), we further assume that Fourier components inside the light cone do not couple to those outside the light cone. Since inside the light cone *k _{x,y}* are small, this lets us use a coarse discretization in

*x*and

*y*, reducing the size of the problem. By using an efficient iterative bi-conjugate gradient method [21, 22], Eq. (6) can be solved to within a tolerance of 10

^{−5}in 20–100 iterations, each of which take less than 10 seconds. Our MATLAB code typically solves (3) and (6) in under 15 minutes on a work station. In contrast, the FDTD calculations for each point in Fig. 2 took tens of hours on a 32 core cluster.

In our simulations for the photosensitive cavity (Fig. 1(a)), we take a *W*1 PCW with background index of *n _{b}* = 2.7, slab thickness,

*t*= 0.7

*d*and hole radius

*a*= 0.3

*d*, where

*d*is the period and Δ

*n*= 0.02, 0.04. For the fluid infiltrated cavity (Fig. 1(b)) we use a

_{p}*W*0.98 silicon PCW (background index

*n*= 3.46), with slab thickness,

_{b}*t*= 0.49

*d*, hole radius

*a*= 0.26

*d*, Δ

*n*= 0.2, 0.4, 0.6. In Fig. 2 we show the

_{i}*Q*-factor versus cavity length calculated using the FAR method (red) and using FDTD (blue). In Fig. 2(a), which is for photosensitive cavities, the efficiency of our theory allows us to vary the cavity length continuously. This is impractical for FDTD calculations, so we only have results at even integer values of the cavity length and at some intervening points. The agreement between the results is excellent: the

*Q*-factors agree to within 30% (or their logarithms by 2%), making them suitable for examining trends in

*Q*. The strong oscillations in

*Q*correspond to a factor of 8. In Fig. 2(b), which is for fluid infiltrated cavities, we only calculated

*Q*for even integer cavity lengths. The agreement for these cavities is good: the results have the same trends and never differ by more than a factor two.

Having demonstrated the reliability of the FAR, we now exploit its semi-analytic nature to explain the far-field radiation pattern in terms of the physical properties of the cavity, and to investigate how the pattern is changed for the different types of perturbation. Figure 3 shows good agreement between the far-field radiation patterns computed using the FAR (left) and FDTD (right), for photosensitive cavities (Figs. 3(a), 3(b)) and fluid infiltrated cavities (Figs. 3(c), 3(d)) of different lengths *L*. Note that (i) the number of lobes in the radiation pattern increases as the cavity gets longer; and that (ii) as discussed earlier, photosensitive cavities radiate predominantly at large declination angles (*θ*), while the fluid infiltrated cavities radiate mostly vertically. Both features can be explained by examining *Ã*(**r**)**D*** ^{a}*(

**r**). The effect of

*Ã*(

**r**) on

**D**

*(*

^{a}**r**) is to introduce nodes and anti-nodes due to Fabry-Perot effects in the cavity. Point (ii) is more subtle: returning to Fig. 1, since the cavity modes are dielectric modes, for the photosensitive cavity the product

*Ã*(

**r**)

**D**

*(*

^{a}**r**) (Fig. 1(c)) introduces sidelobes in the Fourier transform of

**D**

*(*

^{a}**r**). The overlap of these sidelobes with the light cone (Fig. 1(e)), is dominated by (

*k*,

_{x}*k*) values at the edge of the light cone, maximizing the

_{y}*Q*factor and leading to radiation at large declination angles.

For fluid infiltrated cavities, *Ã*(**r**)**D*** ^{a}*(

**r**) (Fig. 1(d)), is nonzero only inside holes. Its Fourier transform within the light cone (Fig. 1(f)) peaks at the origin, because the cavity length is such that

*Ã*(

**r**)

**D**

*(*

^{a}**r**) has a strong non-zero

*DC*Fourier component. This is clear from the fields in the holes in Fig. 1(d): four holes have strong positive fields and only two have strong negative fields because the cavity mode is dominated by the Bloch mode at

*kd*=

*π*, which changes sign each period. This cavity therefore radiates mostly vertically. The qualitative behaviour of the cavity far-field is determined by the term

*Ã*(

**r**)

**D**

*(*

^{a}**r**). This suggests a simple design rule for engineering the radiation pattern of a cavity: construct a perturbation such that

*Ã*(

**r**)

**D**

*(*

^{a}**r**) has a Fourier transform which peaks at (

*k*,

_{x}*k*) values corresponding to the desired direction. A full investigation of such a design strategy is left to future work.

_{y}Similar arguments can be used to explain the variations in *Q* observed in Fig. 2(a): the cavity mode is a superposition of Bloch functions centred about the Brillouin zone edge (*kd* = *π*). It is thus not surprising that the period of the oscillations in Fig. 2(a) corresponds to the period of the central Bloch function. The details of the oscillations in *Q* depend on the superposition of the Bloch modes in *Ã*(**r**)**D*** ^{a}*(

**r**) overlapping with the light cone.

The generalization of our method to non-waveguide cavity types, such as the L3-type cavity [9] is also possible. In such cavities, PC slab modes form the basis for the mode expansion; the perturbation for writing the cavity is larger, hence the perturbation terms Γ̃(**r**) and *ε*̃(**r**) are no longer small, and the integral equation (6) contains additional inhomogeneous terms at leading order. These factors make solving the integral equation more computationally intensive, with the convergence ultimately determined by how well expansion (4) approximates the cavity mode.

We have presented a frequency-domain approach for radiation (FAR) that allows the efficient calculation of the radiative properties of ultra-high *Q* PC cavities. Both *Q*-factors and the radiation patterns are in good agreement with fully numerical FDTD calculations. The orders-of-magnitude improvement in computation speed will enable the application of powerful optimization algorithms, potentially transforming PC cavity design. The FAR lets us directly predict the radiation pattern through its link to the cavity’s refractive index. Although we applied the theory to cavities created by refractive index changes, extensions allow the treatment of other cavity types, created, for example, by shifting inclusions [23] or stretching the lattice [10].

## Acknowledgments

The authors thank A. Rahmani, M.J. Steel and C. Husko for discussions. This work was produced with the assistance of the Australian Research Council (ARC) under the ARC Centres of Excellence program, and was supported by the Flagship Scheme of the National Computational Infrastructure of Australia, and by the Natural Sciences and Engineering Research Council of Canada (NSERC).

## References and links

**1. **K. Nozaki, T. Tanabe, A. Shinya, S. Matsuo, T. Sato, H. Taniyama, and M. Notomi, “Sub-femtojoule all-optical switching using a photonic-crystal nanocavity,” Nat. Photonics **4**, 477–483 (2010). [CrossRef]

**2. **A. Yoshie, J. Hendrickson, G. Khitrova, H. M. Gibbs, G. Rupper, C. Ell, O. B. Shchekin, and D. G. Deppe, “Vacuum Rabi splitting with a single quantum dot in a photonic crystal nanocavity,” Nature **432**, 200–203 (2004). [CrossRef] [PubMed]

**3. **M. W. McCutcheon, J. F. Young, G. W. Rieger, D. Dalacu, S. Frédérick, P. J. Poole, and R. L. Williams “Experimental demonstration of second-order processes in photonic crystal microcavities at submilliwatt excitation powers,” Phys. Rev. B **76**, 245104 (2007). [CrossRef]

**4. **N. Tran, S. Combrié, P. Colman, A. De Rossi, and T. Mei, “Vertical high emission in photonic crystal nanocavities by band-folding design,” Phys. Rev. B **82**, 075120 (2010). [CrossRef]

**5. **D. Englund, A. Majumdar, A. Faraon, M. Toishi, N. Stoltz, P. Petroff, and J. Vučković, “Resonant excitation of a quantum dot strongly coupled to a photonic crystal nanocavity,” Phys. Rev. Lett. **104**, 073904 (2010). [CrossRef] [PubMed]

**6. **S.L. Portalupi, M. Galli, C. Reardon, T. F. Krauss, L. O’Faolain, L. C. Andreani, and D. Gerace, “Planar photonic crystal cavities with far-field optimization for high coupling efficiency and quality factor,” Opt. Express **18**, 16064–16073 (2010). [CrossRef] [PubMed]

**7. **M. Galli, D. Gerace, K. Welna, T. F. Krauss, L. O’Faolain, and G. Guizzetti, “Low-power continuous-wave generation of visible harmonics in silicon photonic crystal nanocavities,” Opt. Express **18**, 26613–26624 (2010). [CrossRef] [PubMed]

**8. **S. G. Johnson, S. Fan, A. Mekis, and J. D. Joannopoulos, “Multipole-cancellation mechanism for high-Q cavities in the absence of a complete photonic band gap,” Appl. Phys. Lett. **78**, 3388–3390 (2001). [CrossRef]

**9. **Y. Akahane, T. Asano, B.S. Song, and S. Noda,” Nature **425**, 944–947 (2003). [CrossRef] [PubMed]

**10. **B. S. Song, S. Noda, T. Asano, and Y. Akahane, “Ultra-high-Q photonic double-heterostructure nanocavity,” Nat. Mater. **4**, 207–210 (2005). [CrossRef]

**11. **C. Sauvan, G. Lecamp, P. Lalanne, and J. P. Hugonin, “Modal-reflectivity enhancement by geometry tuning in Photonic Crystal microcavities,” Opt. Express **13**, 245–255 (2005). [CrossRef] [PubMed]

**12. **D. Englund, I. Fushman, and J. Vučković, “General recipe for designing photonic crystal cavities,” Opt. Express **13**, 5961–5975 (2005). [CrossRef] [PubMed]

**13. **V. A. Mandelshtam and H.S. Taylor, “Harmonic inversion of time signals and its applications,” J. Chem Phys . **107**, 6756–6769 (1997). [CrossRef]

**14. **A. Farjadpour, D. Roundy, A. Rodriguez, M. Ibanescu, P. Bermel, J. D. Joannopoulos, S. G. Johnson, and G. W. Burr, “Improving accuracy by subpixel smoothing in the finite-difference time domain,” Opt. Lett. **31**, 2972–2974 (2006). [CrossRef] [PubMed]

**15. **S. Tomljenovic Hanic, M. J. Steel, C. M. de Sterke, and D. J. Moss, “High-Q cavities in photosensitive photonic crystals,” Opt. Lett. **32**, 542–544 (2007). [CrossRef] [PubMed]

**16. **M. W. Lee, C. Grillet, S. Tomljenovic-Hanic, E. C. Mägi, D. J. Moss, B. J. Eggleton, X. Gai, S. Madden, D. Y. Choi, D. A. P. Bulla, and B. Luther-Davies, “Photowritten high-Q cavities in two-dimensional chalcogenide glass photonic crystals,” Opt. Lett. **34**, 3671–3673 (2009). [CrossRef] [PubMed]

**17. **S. Tomljenovic-Hanic, C. M. de Sterke, and M. J. Steel, “Design of high-Q cavities in photonic crystal slab heterostructures by air-holes infiltration,” Opt. Express **14**, 12451–12456 (2006). [CrossRef] [PubMed]

**18. **U. Bog, C. L. C. Smith, M. W. Lee, S. Tomljenovic-Hanic, C. Grillet, C. Monat, L. O’Faolain, C. Karnutsch, T. F. Krauss, R. C. McPhedran, and B. J. Eggleton, “High-*Q* microfluidic cavities in silicon-based two-dimensional photonic crystal structures,” Opt. Lett. **33**, 2206–2208 (2008). [CrossRef] [PubMed]

**19. **P. Chak, R. Iyer, J. S. Aitchison, and J. E. Sipe, “Hamiltonian formulation of coupled-mode theory in waveguiding structures,” Phys. Rev. E **75**, 016608 (2007). [CrossRef]

**20. **J. E. Sipe, “New Green-function formalism for surface optics,” J. Opt. Soc. Am. B **4**, 481–489 (1987). [CrossRef]

**21. **S. L. Zhang, “GPBi-CG: Generalized product-type methods based on Bi-CG for solving nonsymmetric linear systems,” SIAM Journal on Scientific Computing **18**, 537–551 (1997). [CrossRef]

**22. **P.C. Chaumet and A. Rahmani, “Efficient iterative solution of the discrete dipole approximation for magnetodi-electric scatterers,” Opt. Lett. **34**, 917–919 (2009). [CrossRef] [PubMed]

**23. **E. Kuramochi, M. Notomi, S. Mitsugi, A. Shinya, T. Tanabe, and T. Watanabe, “Ultrahigh-Q photonic crystal nanocavities realized by the local width modulation of a line defect,” Appl. Phys. Lett. **88**, 041112 (2006). [CrossRef]