General Recipe for Designing Photonic Crystal Cavities

We describe a general recipe for designing high-quality factor (Q) photonic crystal cavities with small mode volumes. We first derive a simple expression for out-of-plane losses in terms of the k-space distribution of the cavity mode. Using this, we select a field that will result in a high Q. We then derive an analytical relation between the cavity field and the dielectric constant along a high symmetry direction, and use it to confine our desired mode. By employing this inverse problem approach, we are able to design photonic crystal cavities with Q>4*10^6 and mode volumes V ~ (\lambda/n)^3. Our approach completely eliminates parameter space searches in photonic crystal cavity design, and allows rapid optimization of a large range of photonic crystal cavities. Finally, we study the limit of the out-of-plane cavity Q and mode volume ratio.


Introduction
One of the most interesting applications of photonic crystals (PhCs) is the localization of light to very small mode volumes -below a cubic optical wavelength. The principal confinement mechanism is Distributed Bragg Reflection (DBR), in contrast to, e.g., microspheres or microdisks, which rely solely on Total Internal Reflection (TIR). However, since three-dimensional (3D) PhCs (employing DBR confinement in all 3D in space) have not been perfected yet, the mainstream of the PhC research has addressed 2D PhCs of finite depth, employing DBR in 2D and TIR in the remaining 1D. Such structures are also more compatible with the present microfabrication and planar integration techniques. These cavities can still have small mode volumes, but the absence of full 3D confinement by DBR makes the problem of the high-quality factor (high-Q) cavity design much more challenging. The main problem is out-of-plane loss by imperfect TIR, which becomes particularly severe in the smallest volume cavities.
Over the past few years, many approaches have been proposed to address this issue, but they all focused on the optimization of a particular cavity geometry and a particular mode supported by it [1,2,3,4,5,6,7,8,9,10,11]. Some of the proposed cavities seem promising and have already been proven useful as components of lasers or add/drop filters. However, a general recipe for designing optimized nanocavities is missing, and it is also not known if cavities better than the ones presently known are possible; these are exactly the issues that we will address in this article. In Section 2, we will first derive a simple set of equations for calculation of cavity Q and mode volume. In Section 3 we estimate the optimum k-space distribution of the cavity mode field. In Section 4, we will finally address the question of finding PhC cavities with maximum possible figures of merit for various applications. In this process, we start from the optimum k-space distribution of the cavity field; then we derive an approximate analytical relation between the cavity mode and the dielectric constant along a direction of high symmetry, and use it to create a cavity that supports the selected high-Q mode in a single step. Thus, we eliminate the need for trial and error or other parameter search processes, that are typically used in PhC cavity designs. Furthermore, we study the limit of out-of-plane Q factor for a given V of the cavity mode with a particular field pattern.

Simplified relation between Q of a cavity mode and its k-space Distribution
In order to simplify PhC cavity optimization, we derive an analytical relation between the nearfield pattern of the cavity mode and its quality factor in this section. Q measures how well the cavity confines light and is defined as where ω is the angular frequency of the confined mode. The mode energy is The difficulty lies in calculating P, the far-field radiation intensity. Following our prior work [2], we consider the in-plane and out-of-plane mode loss mechanisms in two-dimensional photonic crystals of finite depth separately: In-plane confinement occurs through DBR. For frequencies well inside the photonic band gap, this confinement can be made arbitrarily high (i.e., Q || arbitrarily large) by addition of PhC layers. On the other hand, out-of-plane confinement, which dictates Q ⊥ , depends on the modal k-distribution that is not confined by TIR. This distribution is highly sensitive to the exact mode pattern and must be optimized by careful tuning of the PhC defect. Assuming that the cavity mode is well inside the photonic band gap, Q ⊥ gives the upper limit for the total Q-factor of the cavity mode. Given a PhC cavity or waveguide, we can compute the near-field using Finite Difference Time Domain (FDTD) analysis. As described in Reference [3], the near-field pattern at a surface S above the PhC slab contains the complete information about the out-of-plane radiation losses of the mode, and thus about Q ⊥ (Fig. 1). The total time-averaged power radiated into the half-space above the surface S is: where K(θ , φ ) is the radiated power per unit solid angle. In the appendix, we derive a very simple form for K in terms of 2D Fourier Transforms (FTs) of H z and E z at the surface S, after expressing the angles θ , φ in terms of k x and k y : Here, η ≡ µ o ε o , λ is the mode wavelength in air, k = 2π/λ , and k || = (k x , k y ) = k(sin θ cos φ , sin θ sin φ ) and k z = k cos(θ ) denote the in-plane and out-of-plane k-components, respectively. In Cartesian coordinates, the radiated power (5) can thus be re-written as the integral over the light cone, k || < k. Substituting (6) into (5) gives This is the simplified expression we were seeking. It gives the out-of-plane radiation loss as the light cone integral of the simple radiation term (6), evaluated above the PhC slab. Substituting Eq. (2) and Eq. (7) into Eq. (1) thus yields a straightforward calculation of the Q for a given mode. In the following sections, when considering the qualitative behavior of (7), we will restrict ourselves to TE-like modes, described at the slab center by the triad (E x , E y , H z ), that have H z even in at least one dimension x or y. For such modes, the term |FT 2 (H z )| 2 in (7) just above the slab is dominant, and |FT 2 (E z )| 2 can be neglected in predicting the general trend of Q.
The figure of merit for cavity design depends on the application: for spontaneous emission rate enhancement through the Purcell effect, one desires maximal Q/V ; for nonlinear optical effects Q 2 /V ; while for the strong coupling regime of cavity QED, maximizing ratios g/κ ∼ Q/ √ V and g/γ ∼ 1/ √ V is important. In these expressions, V is the cavity mode volume: , g is the emitter-cavity field coupling, and κ and γ are the cavity field and emitter dipole decay rates, respectively [2].
Thus, for a given mode pattern, we have derived a simple set of equations that allow easy evaluations of the relevant figures of merit. In the next section, we address the problem of finding the field pattern that optimizes these figures of merit and later derive the necessary PhC to support it.

Optimum k-space field distribution of the cavity mode field
The photonic crystals considered here are made by periodically modulating the refractive index of a thin semiconductor slab (see Fig. 2). The periodic modulation introduces an energy bandstructure for light in two dimensions r = (x, y). The forbidden energy bands are the source of DBR confinement.
The periodic refractive index ε( r) = ε( r + R) can be expanded in a Fourier sum over spatial frequency components in the periodic plane: Here G are the reciprocal lattice vectors in the (k x , k y ) plane and are defined by G · R = 2πm for integer m. The real and reciprocal lattice vectors for the square and hexagonal lattices with periodicity a are: Square Lattice: Hexagonal Lattice: where m, j, q and l are integers. The electromagnetic field corresponding to a particular wave vector k inside such a periodic medium can be expressed as [12]: By introducing linear defects into a PhC lattice, waveguides can be formed. Waveguide modes have the discrete translational symmetry of the particular direction in the PhC lattice and can therefore be expanded in Bloch states E k . In Fig. 3, we plot the dispersion of a waveguide in the ΓJ direction of a hexagonal lattice PhC.
Cavity modes can then be formed by closing a portion of a waveguide, i.e., by introducing mirrors to confine a portion of the waveguide mode. In case of perfect mirrors, the standing wave is described by H = a k 0 H k 0 + a −k 0 H −k 0 . (Here we focus on TE-like PhC modes, and discuss primarily H z , although similar relations can be written for all other field components). Imperfect mirrors introduce a phase shift upon reflection; moreover, the reduction of the distance between the mirrors (shortening of the cavity) broadens the distribution of k vectors in the mode to some width ∆k. The H z field component of the cavity (i.e., a closed waveguide) mode can then be approximated as: A similar expansion of H z can be made at the surface S directly above the PhC slab (see Fig. 1), which is relevant for calculation of radiation losses. The Fourier transform of the above To reduce radiative losses, the mapping of components into the light cone should be minimized [3]. Therefore, the center of the mode distribution k 0 should be positioned at the edge of the first Brillouin zone, which is the region in k-space that cannot be mapped into the light cone by any reciprocal lattice vector G (see Fig. 2). For example, this region corresponds to k 0 = ± j x π ak x ± j y π ak y for the square lattice, where j x , j y ∈ 0, 1. Clearly, | j x | = | j y | = 1 is a better choice for k 0 , since it defines the edge point of the 1st Brillouin zone which is farthest from the light cone. Thus, modes centered at this point and k-space broadened due to confinement, will radiate the least. Similarly, the optimum k 0 for the cavities resonating in the ΓJ direction of the hexagonal lattice is k 0 = ± π ak x (as it is for the cavity from Ref. [4], and for the ΓX direction is k 0 = ± 2π a √ 3k y (as it is for the cavity from Refs. [2] and [3]). Assuming that the optimum choice of k 0 at the edge of the first Brillouin zone has been made, the summation over G can be neglected in Eq. (12), because it only gives additional Fourier components which are even further away from the light cone and do not contribute to the calculation of the radiation losses. In that case, we can express the field distribution as: where A(k x , k y ) is the Fourier space envelope of the mode, which is some function centered at k 0 = (±k 0x , ±k 0y ) with the full-width half-maximum (FWHM) determined by ∆ k = (∆k x , ∆k y ) in the k x and k y directions respectively. Eq. (13) implies that H z (x, y) and A(k x , k y ) are related by 2D Fourier transforms. For example, if A(k x , k y ) can be approximated by a Gaussian centered at k 0 = (k 0x , k 0y ) and with the FWHM of (∆k x , ∆k y ), the real space field distribution H z (x, y) is a function periodic in the x and y directions with the spatial frequencies of k 0x and k 0y , respectively, and modulated by Gaussian envelope with the widths ∆x ∼ 1/∆k x and ∆y ∼ 1/∆k y . Therefore, the properties of the Fourier transforms imply that the extent of the mode in the Fourier space ∆k is inversely proportional to the mode extent in real space (i.e., the cavity length), making the problem of Q maximization even more challenging when V needs to be simultaneously minimized. This has already been attempted in the past for a dipole cavity [2,3] and a linear defect [4,9,11], by using extensive parameters space search. In the following sections, we will design high Q cavities by completely eliminating the need for parameter space searches and iterative trial and error approaches. There are two main applications of Eq. (7). First, this formulation of the cavity Q factor allows us to investigate the theoretical limits of this parameter and its relation to the mode volume of the cavity. Second, it allows us to quantify the effect of our perturbation on the optimization of Q using only one or two layers of the computational field and almost negligible computational time compared to standard numerical methods. We applied Eq. (7) to cavities obtained from an iterative parameter space search. These cavities were previously studied in [3] and [4]. The results for Q using Eq. (7) at S, as well as full first-principle FDTD simulations are shown in Fig. 4. A good match is observed. Therefore, our expression (7) is a valid measure of the radiative properties of the cavity and can be used to theoretically approach the design problem; we can also use this form to speed up the optimization of the cavity parameters. The discrepancy between Eq. (7) and FDTD is primarily due to discretization errors. Fig. 4. Comparison of Q factors derived from Eq. (7) (squares) to those calculated with FDTD (circles). Left: cavity made by removing three holes along the ΓJ direction confining the B oe mode. The Q factor is tuned by shifting the holes closest to the defect as shown by the red arrow. The x-axis gives the shift as a fraction of the periodicity a. Right: the X dipole cavity described in [3]. The Q factor is tuned by stretching the center line of holes in the ΓX direction, as shown by the arrow. The x-axis gives the dislocation in terms of the periodicity a.

Inverse problem approach to designing PhC cavities
In the inverse approach, we begin with a desired in-plane Fourier decomposition of the resonant mode, FT 2 ( H( r)), chosen again to minimize radiation losses given by Eq. (7). The difficulty lies with designing a structure that supports the field which is approximately equal to the target field, H( r).
In this section, we first estimate the general behavior of Q/V for structures of varying mode volume. Then we present two approaches for analytically estimating the PhC structure ε( r) from the desired k-space distribution FT 2 ( H( r)). As mentioned in Sec.2, we restrict the analysis to TE-like modes B eo , B oe , and B ee (Fig.3(b)) for which we can approximate the trend of the radiation (7) by considering only H z at the surface S just above the PhC slab. Moreover, to make a rough estimate of the cavity dielectric constant distribution from the desired H z field on S, we approximate that H z at S is close to H z at the slab center.

General trend of Q/V
The simplification described above allows us to study the general behavior of Q/V for a cavity with varying mode volume. Here, we assume that a structure has been found to support the desired field H z .
We again start from the expression for radiated power, Eq. (7), and calculateQ using Eq. (1). All that is required of the cavity field is that its FT at the surface S above the slab be distributed around the four points k x0 = ±π/a, k y0 = ± 2π √ 3a , to minimize the components inside the light cone. As an example, we choose a field with a Gaussian envelope in Fourier space, as described in Sec. 3. For now, let us consider mode symmetry B oe . The Fourier Transform of the H z field is then given by where σ x and σ y denote the modal width in real space. The mode and its FT are shown in Fig.5 (d-e). We use Eq.7 without the E z terms to estimate the trend in Q, as described above. As the mode volume grows, the radiation inside the light cone shrinks exponentially. This results in an exponential increase in Q. This relationship is shown in Fig.5(g) for field B oe at frequency a/λ = 0.248. At the same time, the mode volume grows linearly with σ x . The growth of Q/V is therefore dominantly exponental , and we can find the optimal Q for a particular choice of mode volume (i.e. σ x ) of the Gaussian mode cavity. According to Fig.5(g), very large Qs can be reached with large mode volumes and there does not seem to be an upper bound on Q ⊥ . As the mode volume of the Gaussian cavity increases, the radiative Fourier components vanish exponentially, but are never zero. A complete lack of Fourier components in the light cone should result in the highest possible Q. As an example of such a field, we propose a mode with a sinc envelope in x and a Gaussian one in y. The FT of this mode in Fig.5(b) is described by where Rect(k x , ∆k x ) is a rectangular function of width ∆k x and centered at k x . The Fourier-transform implies that the cavity mode is described by a sinc function in x whose width is inversely proportional to the width of the rectangle Rect(k x , ∆k x ). To our knowledge, this target field has not been previously considered in PhC cavity design. This field is shown in Fig.5(a-c). Though it has no out-of-plane losses, this field drops off as 1 r and therefore requires a larger structure than the Gaussian field for confinement.
Over the past years, many new designs with ever-higher theoretical quality factors have been suggested [11]. In light of our result that Q/V increases exponentially with mode size, these large Qs are not surprising. On the other hand, direct measurements on fabricated PhC cavities indicate that Qs are bounded to currently ∼ 10 4 by material absorption and surface roughness [13,14], while indirect measurements using waveguide coupling indicate values up to 6 · 10 5 [9,11].
It is interesting to note that Eq. (6) also allows one to calculate the field required to radiate with a desired radiation distribution. For example, many applications require radiation with a strong vertical component; waveguide modes with even H z can be confined for this purpose so

Estimating Photonic Crystal Design from k-space Field Distribution
Now we introduce two analytical ways of estimating the dielectric structure ε( r) that supports a cavity field that is approximately equal to the desired field H c . These methods directly calculate the dielectric profile from the desired field distribution, without any dynamic tuning of PhC parameters, and are thus computationally fast. We focus on TE-like modes, since they see a large bandgap and exhibit electric-field maxima at the slab center. For TE-like modes, H c = H cẑ at the center of the slab, and H c ≈ H cẑ at the surface. First, we relate H c to one of the allowed waveguide fields H w . The fields H c and H w at the center of the PhC slab (z = 0) are solutions to the homogeneous wave equation with the corresponding refractive indices ε c and ε w , respectively.
C is a positive constant of integration, and H c = H w H e , where H w is the known waveguide field and H e is the desired field envelope. C can be chosen by fixing the value of ε c at some x, leading to a particular solution for ε c . In our cavity designs we chose C such that the value of ε c is close to ε w at the cavity center. To implement this design in a practical structure, we need to approximate this continuous ε c by means of a binary function with low and high-index materials ε l and ε h , respectively. We do this by finding, in every period j, the air hole radius r j that gives the same field-weighted averaged index on the x-axis: where E c is estimated from a linear superposition of waveguide modes as E c ∝ ∇ × H c . We assume that the holes are centered at the positions of the unperturbed hexagonal lattice PhC holes. The radii r j thus give the required index profile along the x symmetry axis. The exact shape of the holes in 3D is secondary -we choose cylindrical holes for convenience. Furthermore, we are free to preserve the original hexagonal crystal structure of the PhC far away from the cavity where the field is vanishing.
To illustrate the power of this inverse approach, we now design PhC cavities that support the Gaussian and sinc-type modes of Eq. (14), (15). In each case, we start with the waveguide field B oe of Fig.3(b) confined in a line-defect of a hexagonal PhC. The calculated dielectric structures and FDTD simulated fields inside them are shown in Fig. 6. The FT fields on S also show a close match and very little power radiated inside the light cone ( Fig.6(c,f)). This results in very large Q values, estimated from Q ⊥ to limit computational constraints. These estimates were done in two ways, using first principles FDTD simulations [2], and direct integration of lossy components by Eq. (7). The results are listed in Table 1 and show an improvement of roughly three orders of magnitude over the unmodified structure of Fig. 4 with as small increase in mode volume. Furthermore, a fit of the resulting field pattern to a Gaussian envelope multiplied by a Sine, yielded a value of σ x /a ≈ 1.6, which, according to the plot in Fig. 5 g., puts us at the attainable limit of Q ⊥ at this mode volume.
In our FDTD simulations, we verified that Q ⊥ correctly estimates Q by noting that Q ⊥ did not change appreciably as the number of PhC periods in the x− and y− directions, N x and N y , was increased: for the Gaussian-type (sinc-type) mode, increasing the simulation size from N x = 13, N y = 13 (N x = 21, N y = 9) PhC periods to N x = 25, N y = 13 (N x = 33, N y = 13) changed quality factors from Q || = 22 · 10 3 , Q ⊥ = 1.4 · 10 6 (Q || = 17 · 10 3 , Q ⊥ = 4.2 · 10 6 ) to Q || = 180 · 10 3 , Q ⊥ = 1.48 · 10 6 (Q || = 260 · 10 3 , Q ⊥ = 4.0 · 10 6 ). (The number of PhC periods in the x-direction in which the holes are modulated to introduce a cavity is 9 and 29 for Gaussian and sinc cavity, respectively, while both cavities consist of only one line of defect holes in the y-direction.) Thus, with enough periods, the quality factors would be limited to Q ⊥ , as summarized in the table. In the calculation of Q, the vertically emitted power P || was estimated from the fields a distance ∼ 0.25 · λ above the PhC surface. Note that the frequencies a/λ closely match those of the original waveguide field B oe (a/λ cav =0.251), validating the assumption in the derivation. Table 1. Q values of structures derived with inverse-approach 1 a/λ cav Q cav (freq. filter) Q cav (Eq. (7)) V mode ( λ n ) 3 Gaussian 0.248 1.4 · 10 6 1.6 · 10 6 0.85 Sinc 0.247 4.2 · 10 6 4.3 · 10 6 1.43 Unmodified 3-hole defect 0.251 6.6 · 10 3 6.4 · 10 3 0.63 We now derive a closed-form expression for ε c (x, y) that is valid in the whole PhC plane (instead of the center line only). Again, begin with the cavity field H( r) =ẑH c consisting of the product of the waveguide field and a slowly varying envelope, H c = H w H e , and treat the cavity dielectric constant as: 1 ε c = 1 ε pert + 1 ε w . In the PhC plane, Eq. (16, 17) for a TE-like mode can be rewritten as Multiplying the last equation by H e , subtracting from the first, and recalling that ω c ∼ ω w yields