First-principles method for high-$Q$ photonic crystal cavity mode calculations

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 $\sim 100$ times more rapidly than conventional finite-difference time domain calculations. It also provides a simple rule for engineering the cavity's far-field radiation pattern.

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,5,6] in cavity QED [5] and harmonic generation [7] experiments.
Photonic crystal cavity design uses established theoretical ideas [8,9,10,11,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 [6,5], 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 i in a strip-like region (red shading), typically by fluid infiltration [17,18]. These two cavities are therefore complementary, 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)-(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 y = 0 ( Fig.  1(f)), and thus it predominantly radiates vertically. We later return to this insight and use it to provide a general design rule for engineering the radiation pattern of a cavity mode.
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 Since we use PCW modes as a basis, it is convenient to define ε(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 where we only include modes of the even PCW band, i.e. those for which 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 where γ(r) = 1/(2ε 0 ) [1/ε(r) − 1/ε(r)], and we dropped non-rotating wave terms involving a † k a † k and a k a k . Diagonalizing H 1 determines an eigenvalue, the energyhω 0 of a photon in the cavity mode, while its eigenfunction v 0 (k) gives the cavity mode in the basis of PCW modes: We now have an approximate expression for the cavity mode in terms of a basis with Fourier components outside the light cone. The Fourier content within the light cone of the ultra-high 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 where the Green tensor expresses the electric field at 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 |κ κ κ| 2 ≡ k 2 x + k 2 y ≤ k 2 0 , from evanescent modes with |κ κ κ| 2 > k 2 0 , where k 0 = ω 0 /c. Defining Γ(r) = [ε(r) − 1] /ε(r), the polarisation 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 P rad (r) =Γ(r)D a (r) + P rad c (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 P rad c (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 P rad 1 (r) − ε 0 (ε(r) − 1) dr G(r − r ; ω 0 )P rad 1 (r ) =Ã(r)D a (r) ≡ Γ (r) +ε where the driving term, which contains information about the cavity via the parameters with a tilde, couples to Fourier components inside the light cone. As discussed earlier in this paper, we have found thatÃ(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.
for s polarization, and with a similar expression for p polarization. Here R = (x, y) and w = (k 2 0 − |κ κ κ| 2 ) 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 y ) 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.
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.7d and hole radius a = 0.3d, where d is the period and ∆n p = 0.02, 0.04. For the fluid infiltrated cavity ( Fig. 1(b)) we use a W 0.98 silicon PCW (background index n b = 3.46), with slab thickness, t = 0.49d, hole radius a = 0.26d, ∆n i = 0.2, 0.4, 0.6. In Fig. 2 we show the 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 design desirable far-field radiation properties. Figure 3 shows good agreement between the farfield radiation patterns computed using the FAR (left) and FDTD (right), for photosensitive cavities (Figs. 3(a),(b)) and fluid infiltrated cavities (Figs. 3(c),(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)) has the effect of merely introducing 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 y ) values at the edge of the light cone, maximizing the Q factor and leading to radiation at large declination angles. Fig. 3. Symmetric quadrants of far-field radiation (S r ) for (a),(b) cavities in Fig. 1(a) with ∆n p = 0.02 and (c),(d) those in Fig. 1(b) with ∆n i = 0.2. Left frames are computed using the FAR method while right frames are computed using FDTD. Colors are as in Fig. 1(d).
Angles φ and θ are azimuthal and declination angles respectively.
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. Since examiningÃ(r)D a (r) is sufficient for qualitative insight into the far-field, the requirement for steering the radiation of cavity modes is thus simple: construct a perturbation such thatÃ(r)D a (r) has a Fourier transform which peaks at (k x , k y ) values corresponding to the desired direction.
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.
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 ordersof-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].
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).