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_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)–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_y = 0 (Fig. 1(f)), and thus it predominantly radiates vertically.

 

Fig. 1 Schematic of (a) photosensitive cavity with locally increased background index and of (b) fluid infiltrated cavity with increased hole index. Ã(r)D^a(r) (arb. units) for (c) photosensitive cavity (d) fluid infiltrated cavity, both with a length of L = 4d. (e)–(f) Modulus of the Fourier transform of (c)–(d) respectively, with non-radiating components removed.

Download Full Size | PDF

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

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

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

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 ω = ω₀ +ω̃. We next perform a Taylor expansion about ω₀ 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

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^ik₀r/r, where κ̄ ≡ k₀r̂ · (x̂x̂ + ŷŷ), with, above (+) and below (−) the slab,

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⁻⁵ 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.

 

Fig. 2 Quality factor versus cavity length for (a) the photosensitive cavity (Fig. 1(a)); (b), the fluid infiltrated cavity (Fig. 1(b)). Red symbols are computed using the FAR method, while blue ones are computed by FDTD.

Download Full Size | PDF

In our simulations for the photosensitive cavity (Fig. 1(a)), we take a W1 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 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_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.

Download Full Size | PDF

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_y) values corresponding to the desired direction. A full investigation of such a design strategy is left to future work.

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].