## Abstract

We analyze and compare the effect of fabrication disorder on the quality factor of six well-known high-index photonic crystal cavity designs. The theoretical quality factors for the different nominal structures span more than three orders of magnitude, ranging from 5.4 × 10^{4} to 7.5 × 10^{7}, and the defect responsible for confining light is introduced in a different way for each structure. Nevertheless, among the different designs we observe similar behavior of the statistics of the disorder-induced light losses. In particular, we show that for high enough disorder, such that the quality factor is mainly determined by the disorder-induced losses, the measured quality factors differ marginally – not only on average as commonly acknowledged, but also in their full statistical distributions. This notably shows that optimizing the theoretical quality factor brings little practical improvement if its value is already much larger than what is typically measured, and if this is the case, there is no way to choose an alternative design more robust to disorder.

© 2013 Optical Society of America

## 1. Introduction

On-chip photonic technology has witnessed vigorous development in the last two decades. Photonic crystals (PhCs) in particular currently lie at the forefront of research in photonics [1–4], with a major effort directed towards fabricating high quality factor (*Q*) two-dimensional photonic crystal slab cavities with modal volumes close to the diffraction limit [5–12]. Such devices have found many applications, for example making possible one of the first observations of the strong coupling cavity quantum electrodynamics regime [13, 14] in a solid-state system [15]. Recently, there have been equally remarkable experimental breakthroughs in the all-optical field of photonics, including the realization of low-power optical switching [16–18] and optical random-access memory [19], and the demonstration of strong coupling between two distant PhC cavities [20]. The success of such experiments is due, on one hand, to the highly sophisticated fabrication technology, which allows for sub-nanometer precision in state-of-the-art silicon devices [21, 22], and on the other, to the introduction of optimized, ultra-high-*Q* cavity designs [5–12], for which the theoretical *Q* is in the range of 10^{5} − 10^{9}. However, a general feature of such designs is that the quality factor predicted for the nominal structure is much higher – sometimes by more than an order of magnitude – than the experimentally measured one, which is unanimously attributed to disorder residual in the fabrication process (most notably disorder in the size and positioning of air holes, though other contributions have also been identified [23–25]). Understanding the disorder effects has spurred some theoretical [25–27] and experimental [21, 22, 25] research, which has most notably shown – for two specific designs (the *L*3 cavity and the heterostructure cavity) – that the disorder-averaged value of the disorder-induced cavity losses scales as *σ*^{2}, where *σ* is the magnitude of the random fluctuations in one spatial direction (e.g. in hole radius and/or position). Additionally, it is by now acknowledged that said average value depends marginally on the specific cavity design. However, a better understanding of the statistics of the disorder effects (e.g. investigating the full probability distribution of *Q* values) – and how those compare among different cavity designs – is still needed.

In this paper we analyze and compare six very different, well-known PhC cavities with theoretical *Q* ranging from 5.4×10^{4} to 7.5×10^{7}, i.e. spanning more than three orders of magnitude: the *L*3 [5], *A*3 and *A*1 [8, 9], heterostructure (HS) [6], *H*1 [28, 29] and *H*0 [12, 28]. We perform a detailed statistical analysis of the disorder effects, and conclude, most notably, that whenever the disorder-induced losses are clearly dominating over the intrinsic ones, their statistical distribution depends weakly on the particular cavity design, which validates and strengthens the common belief concerning the average value of *Q*. The results are relevant for Si, GaAs, and other high-index materials, although there is no reason to believe that similar considerations would not apply to lower index contrast. In addition, in the Appendix we discuss the different numerical tools currently available for computing *Q*, and focus in particular on the applicability to high-*Q* cavities of the recently developed Bloch-mode expansion method [30], which proves to be both computationally efficient and reliable. The paper is organized as follows: in Section 2, we show the statistical analysis of disorder; Section 3 contains a discussion in view of improving the performance of PhC cavities and the importance of design optimization; in Section 4 (the Appendix), we discuss the numerical tools and the Bloch-mode expansion.

## 2. Analysis of disorder

#### 2.1. Cavity designs and disorder model

The designs of the six PhC nanocavities studied in this paper are presented in Fig. 1. All of the devices are based on a triangular lattice of circular holes of radius *R* etched in a silicon slab of thickness *d* (*R* and *d* differ among designs). The *L*3 cavity of Fig. 1(a), also studied in Refs. [22, 26], has *R* = 0.3*a* and *d* = (220/420)*a*, with *a* the lattice constant. The cavity defect is introduced by three missing holes in a row, and in addition, for *Q*-optimization, the two holes on each side of the defect (marked in red) are shifted outward by 0.16*a*, while their radius is decreased by 0.06*a*. The *A*3 and *A*1 cavities (Fig. 1(b) and Fig. 1(c)) are taken as originally defined in Ref. [8], with *R* = 0.25*a*, *d* = 0.472*a* for the *A*3 and *R* = 0.257*a*, *d* = 0.486*a* for the *A*1. The *A*3 cavity is formed by first introducing a one-dimensional waveguide defect of width 0.9*a*, and then shifting two holes (marked with red contours in panel (b) of Fig. 1) on each side of the linear defect outward by *d _{y}* = 0.0278

*a*. The

*A*1 cavity is similar, but the width of the waveguide is 0.98

*a*, and three types of hole shifts are introduced: by

*d*= 0.0214

_{y}*a*(red holes in Fig. 1(c)), by 2

*d*/3 (black holes) and by

_{y}*d*/3 (grey holes). Next, we consider the heterostructure cavity first introduced in Ref. [6], with

_{y}*R*= (110/410)

*a*and

*d*= (220/410)

*a*. It is again based on a linear waveguide (of width

*a*), but the cavity defect is introduced by changing the lattice constant to (420/410)

*a*in the central region and to (415/410)

*a*in a region on each side (Fig. 1(d)). Fig. 1(e) shows the currently best-optimized [29] design for an

*H*1 cavity formed by one missing hole, with

*d*= 0.5

*a*and

*R*= 0.3

*a*and nearby-hole shifts introduced as explained in the figure caption. This cavity in fact has two degenerate modes, and in the figure only the y-polarized one is plotted. The last cavity is the optimized

*H*0 as defined in Ref. [12], where

*d*= 0.6

*a*and

*R*= 0.26

*a*, and the defining defect is outward shifts of neighboring holes. In all simulations of disorder, random fluctuations in the positions and radii of the holes were considered, with an underlying Gaussian distribution with zero mean and standard deviation

*σ*. The magnitude of the fluctuations was assumed to be spatially uniform, and no hole-hole correlations were considered. Spatial correlations between hole size and positions, as well as spatially non-uniform distribution of the disorder parameters could easily be included in our simulations. However, since there is no specific experimental characterization of the magnitude of such distribution properties and their influence on the PhC, such an analysis is beyond the scope of the present work. The disorder model can also be readily extended to include irregular hole shapes [31] (i.e. side-wall roughness), but previous results show quite unambiguously that the disorder effects are to a large extend determined by the magnitude of the

*area*fluctuations rather than the particular hole shapes [31, 32], thus captured very well by the simpler model used here. Everywhere below, by “disorder realization” we understand one particular set of 3

*n*independent pseudo-random numbers drawn from a Gaussian distribution, giving the shifts in the

*x*and

*y*positions and in the radius

*R*of each hole (for

*n*holes included in the system). The magnitude of disorder included here ranges from

*σ*= 0.001

*a*(relevant for state-of-the-art Si structures [21,32]) to

*σ*= 0.015

*a*, i.e. includes values relevant for different systems like GaAs, InGaAs or InGaAsP.

#### 2.2. Statistics of the disorder-induced losses

The measured *Q* of a cavity can be written as

*Q*is the theoretically predicted

_{i}*Q*-factor for the particular cavity design in the absence of disorder (also called “intrinsic”),

*Q*is the quality factor associated to the disorder-induced losses and

_{d}*Q*accounts for any absorption (e.g. due to free carriers or water condensation). The absorption term

_{a}*Q*could come from non-linear absorption and thus depend on e.g. the mode volume or field intensity in the cavity, but in any case it does not correspond to a stochastic effect, and for a particular cavity, it is expected to be constant for the different disorder realizations. Since the main interest of our analysis is the statistics of the disorder effects, this additional constant can then be incorporated into the constant

_{a}*Q*. The same consideration applies to any other systematic effect that results in losses which are independent of the disorder realization, e.g. tilt of the hole walls. The expression of Eq. (1) is commonly employed [7,21,22,25,27] and stems from the definition of

_{i}*Q*≡

*ω*〈

*U*〉/〈

*P*〉, where

*ω*is the cavity frequency, 〈

*U*〉 is the total electromagnetic energy of the cavity mode, and 〈

*P*〉 is the radiation power. The physical assumption entering Eq. (1) (neglecting

*Q*) is then that, upon small perturbations, the frequency and total energy of the cavity mode are approximately constant, while the total radiated power is obtained by summing the contributions of two radiation channels: one due to the intrinsically radiative components of a cavity mode and the other due to disorder-induced scattering into modes radiating above the surface (i.e. modes with frequency inside the light-cone). This summation is in fact not necessarily correct: the interference between the fields emitted through the two different channels has to be taken into account, which could potentially be

_{a}*destructive*. In particular, this interference can be to a good approximation neglected when one of the channels is strongly dominating, but is important when

*Q*≈

_{i}*Q*. Nevertheless, in our analysis we use Eq. (1), bearing in mind that

_{d}*Q*is, in general, a purely phenomenological quality factor that accounts for the difference between

_{d}*Q*and the measured

_{i}*Q*, and study the validity of interpreting

*Q*as a quality factor associated to disorder-induced scattering losses which do

_{d}*not*interfere with the intrinsic cavity losses. One example for which this interpretation is obviously not true can be found in the negative values of

*Q*, which we show are possible for particular disorder realizations. Physically, this reflects the possibility for the disorder to

_{d}*cancel*– due to destructive interference – some of the intrinsic radiative components, rather than introduce more losses. The possibility for negative

*Q*-values is intuitive and by no means an artifact of our simulations: in the same way that small, controlled structural changes can dramatically increase the quality factor of a cavity (consider for example the difference between the

_{d}*A*1 and the

*A*3), particular configurations of the random disorder can do just the same.

In this section, we use the guided-mode expansion [33] (GME) to compute the *Q* of a cavity with disorder, a method that has already proven reliable for similar studies [22] (a comparison of the accuracy and efficiency between this and other popular methods is given in the Appendix). Essentially, the method consists of expanding the PhC modes on the basis of the guided modes of the spatially uniform dielectric slab,

**G**are reciprocal lattice vectors in the slab-plane, defined with respect to a supercell that contains the structure to simulate,

*α*= 1, 2,... is a guided-mode index, and

**k**is a Bloch vector associated to the modes (thus for each

**k**a quality factor is computed, and the final

*Q*is the average over

**k**, see Ref. [33]). The coefficients

*c*(

**k**+

**G**,

*α*) are then obtained by diagonalization of a matrix with elements defined in Ref. [33]. The size of the computational supercell taken here is between 16

*a*and 32

*a*in the x-direction and between $14\sqrt{3}a/2$ and $20\sqrt{3}a/2$ in the y-direction, depending on the spread of the electric field of the different cavities, and the reciprocal space was truncated to ${G}_{\mathit{max}}=2\frac{2\pi}{a}$, resulting in 4080 to 9120 reciprocal lattice points, depending on the supercell size. Convergence was checked versus all the parameters. It is important to remark that GME automatically imposes periodic boundary conditions along the x and y directions as opposed to first-principle methods, where absorbing boundaries are used (further discussion of the different numerical tools can be found in the Appendix).

Figure 2 shows the influence of random disorder on the resonant frequencies and quality factors for the six different designs, computed for *σ* = 0.001*a*, 0.0014*a*, 0.002*a*, 0.003*a*, 0.004*a* and 0.005*a*, and for 400 disorder realizations for each cavity and each *σ*. For all cavities and disorder magnitudes, the mean of the deviation from the ideal-cavity frequency *ω*_{0}, 〈*ω* − *ω*_{0}〉, is at least one order of magnitude smaller than the standard deviation *δ* (*ω*) = (〈*ω*^{2}〉 − 〈*ω*〉^{2})^{1/2}, thus, to a good approximation, the frequency distribution can be considered centered around the ideal-cavity frequency. In contrast, the standard deviation increases linearly with increasing *σ* as shown in Fig. 2(a). This dependence matches previous results for the heterostructure [27], and is here observed for all six cavities (for the *H*1, where two modes with frequencies *ω*_{1} and *ω*_{2} exist, we define *δ* (*ω*) = (*δ* (*ω*_{1}) + *δ* (*ω*_{2}))/2). For some practical applications, having the smallest possible *δ* (*ω*) is essential, but the main focus of our work are the corresponding quality factors. For a given design, the statistics of the measured *Q* are determined by the statistics of *Q _{d}*, which is why in Fig. 2(b) and Fig. 2(c) we plot the mean 〈1/

*Q*〉 and standard deviation

_{d}*δ*(1/

*Q*) = (〈(1/

_{d}*Q*)

_{d}^{2}〉 − 〈1/

*Q*〉

_{d}^{2})

^{1/2}vs.

*σ*

^{2}for the different cavities, where 1/

*Q*was computed from Eq. (1). The scaling 〈1/

_{d}*Q*〉 ∝

_{d}*σ*

^{2}has already been numerically established for the

*L*3 [26] and the HS [27], and is confirmed by our results (Fig. 2(b)), for those two as well as the other four cavities. The straight lines in the plot show linear interpolation based on the data. Additionally we find that the standard deviation (Fig. 2(c)) also scales linearly with

*σ*

^{2}for the

*A*1 and HS in the range we consider, while for the other designs, the dependence is linear for all but the smallest values of

*σ*.

An intuitive model for the disorder-induced scattering, which matches to some extent the computed results, is to assume that
$1/{Q}_{d}={A}_{c}{\delta}_{d}^{2}$, where *δ _{d}* is a random variable drawn from a Gaussian distribution of zero mean and standard deviation

*σ*. Thus, 〈1/

*Q*〉 =

_{d}*A*

_{c}σ^{2}with

*A*a design-specific constant, as observed in Fig. 2(b), while $\u3008{\left(1/{Q}_{d}\right)}^{2}\u3009={A}_{c}^{2}\left({\delta}_{d}^{4}\right)={A}_{c}^{2}\times 3{\sigma}^{4}$, which means $\delta \left(1/{Q}_{d}\right)=\sqrt{2}{A}_{c}{\sigma}^{2}$, i.e. the standard deviation also scales linearly with

_{c}*σ*

^{2}, as observed in Fig. 2(c). However, there are two notable deviations from this simplified expectation. One is the fast increase of

*δ*(1/

*Q*) for the

_{d}*A*3,

*L*3,

*H*0 and

*H*1 for low

*σ*. This, however, occurs in the

*σ*range for which 〈1/

*Q*〉

_{d}^{−1}> ≈

*Q*; in this situation, as discussed earlier,

_{i}*Q*is not determined in a straightforward way by the absolute value of the disorder-induced scattering losses, but instead also contains interference effects. Indeed, in this range we also observe

_{d}*δ*(1/

*Q*) > 〈1/

_{d}*Q*〉, meaning that negative values for 1/

_{d}*Q*are likely, thus the statistics lie beyond the simple, uncorrelated assumption. The second deviation is the fact that within the simple Gaussian-squared model, we obtain $\delta \left(1/{Q}_{d}\right)/\left(1/{Q}_{d}\right)=\sqrt{2}$ for all cavities. In our simulations, this value is much lower: 0.48, 0.43, 0.36, 0.40, 0.56 and 0.67 for the

_{d}*L*3,

*A*3,

*A*1, HS,

*H*1 and

*H*0, respectively. This means that

*δ*is actually not Gaussian-distributed, as $\u3008{\delta}_{d}^{4}\u3009<3{\u3008{\delta}_{d}^{2}\u3009}^{2}$. A likely explanation is that while it is reasonable to expect that the statistics of the losses due to the fluctuation of a single hole follow a Gaussian model, the statistics of the total losses (due to the fluctuations of

_{d}*all*holes) are

*not*given by a direct sum over single-hole contributions. Our results show that this cooperative effect systematically results in a decrease of the variance of the 1/

*Q*distribution.

_{d}The constant *A _{c}* (the slope of the lines in Fig. 2(b)) can in principle be interpreted as a measure of the robustness of a cavity to disorder effects. The values of

*A*for the different cavities are 0.159, 0.121, 0.107, 0.101, 0.273 and 0.230 for the

_{c}*L*3,

*A*3,

*A*1, HS,

*H*1 and

*H*0, respectively. With the exception of the HS, these decrease with increasing cavity mode volume; the HS has a volume slightly lower than the

*A*3 and

*A*1 combined with the lowest

*A*. Additionally, a cavity with a lower

_{c}*A*than the HS has been demonstrated [34], but at the cost of an even higher mode volume. However, by investigating the full probability distribution of 1/

_{c}*Q*, plotted for

_{d}*σ*= 0.0014

*a*,

*σ*= 0.005

*a*and

*σ*= 0.015

*a*in Fig. 3(a)–(c), we find that the difference is marginal whenever 〈1/

*Q*〉

_{d}^{−1}<

*Q*. The latter holds for all three values of

_{i}*σ*for the

*A*1 and HS, and indeed the distributions for those two cavities are almost identical in all the panels. Despite the fact that the slope of the

*A*3 plot in Fig. 2(a) is visibly higher than of the

*A*1 and HS, already for

*σ*= 0.005

*a*(panel (b) of Fig. 3) it can be seen that the three distributions become very close, illustrating that the difference in 〈1/

*Q*〉 is actually minor. This reflects the fact that if disorder-induced losses are dominating, which is already the case for the

_{d}*A*3 with

*σ*= 0.005

*a*, for which

*Q*= 1.5×10

_{i}^{6}and 〈1/

*Q*〉

_{d}^{−1}= 3.4×10

^{5}, then the full statistical distribution of those losses – and not only the average

*Q*– tends towards being

*design-independent*. This conclusion is further supported by the data for the other cavities, despite the fact that their

*Q*is more than two orders of magnitude smaller than the

_{i}*A*1 or HS values. Of course, the distributions of the quality factors (panels (d)–(f) in Fig. 3) differ vastly among the different designs for small

*σ*, when

*Q*is important. However, when

_{i}*Q*is strongly dominant (Fig. 3(f)), the distributions come close together, and it can be expected that they will be even closer for yet higher disorder, but this is a case which is, first of all, irrelevant to state-of-the-art devices, second, in which the scaling 〈1/

_{d}*Q*〉 ∝

_{d}*σ*

^{2}is expected to break down [26], and third, in which the disorder becomes comparable to the hole displacements defining the cavities.

The GME provides not only the resonant frequency and quality factor of a cavity, but also the full (near-field) mode profile. A parameter which is for example interesting for cavity-QED applications is the degree of polarization in the center of a cavity (where a quantum dot would be placed), which for the *L*3, *A*1, *A*3 and HS lies entirely along the y-direction in the ideal case. If we define the degree of polarization in the center of a cavity as
$P=\frac{{\left|{E}_{y}(0)\right|}^{2}-{\left|{E}_{x}(0)\right|}^{2}}{{\left|{E}_{y}(0)\right|}^{2}+{\left|{E}_{x}(0)\right|}^{2}}$, so that *P* = 1 for fully y-polarized field (ideal-cavity case) and *P* = −1 for fully x-polarized field, we obtain that the influence of disorder is not very profound: for the 1000 realizations with each of the three values of *σ* as in Fig. 3, the distribution of *P* is always peaked around 1, and all values are above 0.999 for *σ* = 0.0014*a*, 0.99 for *σ* = 0.005*a*, and 0.9 for *σ* = 0.015*a*. Note that the effect of disorder on the far-field is expected to be less trivial [27], but since our GME implementation is currently not suited for computing the far-field profile, this analysis is left for future work.

The fact that the 1/*Q _{d}* distribution for the HS,

*A*1 and

*A*3 are almost identical, and that all three cavities are based on a local defect over a linear waveguide suggests the question of whether the same disorder realization affects the cavities in a similar way. To investigate this, in Fig. 4 we plot, for the same three values of disorder as in Fig. 3, 1/

*Q*in one cavity vs. 1/

_{d}*Q*in another for the same random configuration of hole position and size shifts (taken separately from the defining cavity defect). Strong correlation (which gets stronger with increasing

_{d}*σ*) is indeed found, with a correlation coefficient as high as 0.98 between the

*A*1 and the HS for

*σ*= 0.015

*a*. To understand why correlation is expected, we note that the modes of these cavities can be expanded on the basis of Bloch modes of the waveguide (see also the Appendix). In the ideal case, the main contribution comes from the modes close to the edge of the main guided band, which are perfectly lossless, and the lossy contributions come from small components of modes above the light line [7, 11]. When disorder is introduced in the waveguide (without a cavity), the ideally-lossless modes acquire a finite quality factor [26, 30, 33]. We can thus further split the disorder losses of a cavity into 1/

*Q*= 1/

_{d}*Q*+ 1/

_{des}*Q*, where 1/

_{wm}*Q*is due to the way the disorder in the defining defect increases the coupling to modes above the light line (i.e. disorder makes the design sub-optimal), while 1/

_{des}*Q*reflects the fact that the waveguide modes themselves become lossy. Thus, 1/

_{wm}*Q*is design-dependent, while 1/

_{des}*Q*is largely ubiquitous, especially since the near-field profiles of Fig. 1 are very similar; the strong correlations with increasing

_{wm}*σ*then show that the losses become increasingly waveguide-mode dominated. Finally, we note that the

*A*1-HS correlation is very high even for

*σ*= 0.0014

*a*, which highlights the fact that those two defect cavities are extremely similar.

## 3. Discussion and conclusion

Our study of the disorder-induced losses has an important implication on future PhC cavity optimization. It is evident from Fig. 3 that, despite the fact that the designs of the cavities we study and their theoretical *Q _{i}*-s differ vastly, in the regime in which the disorder-induced losses are dominating (

*Q*<

_{d}*Q*), the statistics of 1/

_{i}*Q*are to a good approximation design-independent. Thus, optimizing

_{d}*Q*is important only if the measured

_{i}*Q*(or, better, the 〈

*Q*〉 measured over many samples) is not much lower than

*Q*. If that is the case, there are just two ways to significantly increase the measured

_{i}*Q*. One is to decrease the magnitude of disorder, which for all studied designs is related to the extrinsic losses via 〈1/

*Q*〉 ∝

_{d}*σ*

^{2}. The second (albeit not reproducible) way is to explore the broad high-

*Q*tail of the probability distribution, which is visible in Fig. 3(d)–(f). The thickness of the tail is related to the fact that in the disorder-dominated regime, the distribution of 1/

*Q*is close to a squared Gaussian, as illustrated in Fig. 3(a)–(c). This results, for example, in a 5.7% chance to find an HS cavity with

_{d}*Q*= 7.4 × 10

^{6}– the average value with

*σ*= 0.001

*a*– even within the

*σ*= 0.0014 case. Of course, it should be kept in mind that for such high values of the nominal

*Q*, the measured value could be limited not only by disorder, but also by linear absorption, or two-photon absorption and other non-linear effects [35], which all introduce an additional effective loss rate.

Since the parameter space of hole positions and sizes in a typical PhC structure is not only large but also inseparable (i.e. the effect of hole fluctuations on *Q* is highly correlated among different holes), and presents many local maxima of *Q _{i}*, the ideal tool for optimizing a design would be a

*global*optimization algorithm – e.g. a genetic optimization. In the Appendix (Section 4), we propose a fast and qualitatively reliable Bloch-mode expansion method for the computation of

*Q*and the field profiles of waveguide-based cavities. We believe that this is the first numerical tool which is fast enough to make a global optimization algorithm effective, which is an objective of our future work. Fine-tuning can in principle still be carried out using a first-principle method, but as we discuss in the Appendix, this will probably bring little improvement.

_{i}Our results are in line with the idea – commonly acknowledged within the community – that disorder severely limits the measured *Q*, but to our knowledge represent the first study rigorously comparing the statistics of the disorder effects among different designs. We show that once the *Q _{i}* of a cavity has reached a certain value (≈ 10

^{7}for current silicon-based devices, and lower for other materials with higher magnitude of the fabrication imperfections), the best improvement strategy is not to push

*Q*higher still, but rather to improve the fabrication process, or the features of the design other than the quality factor, e.g. by focusing on cavities with much smaller mode volume [12], or with designs allowing for more sophisticated practical applications [36–40].

_{i}## 4. Appendix: Discussion of PhC modeling methods

## 4.1. The Bloch-mode expansion

Since four of the cavities studied here can be considered as a perturbation over a PhC waveguide (including the *L*3 cavity, although in this case the perturbation is rather strong, consisting in introducing additional holes along the waveguide), we explore the possibility to use the recently developed Bloch-mode expansion [30] (BME) method for computationally fast analysis of high-*Q* structures. The BME relies on the fact that the Bloch modes of a disorder-less waveguide are fast to compute (here this is again done using the GME), and the mode of the cavity, potentially including disorder, can be expanded on the basis of a set of the Bloch bands. In the limit when all bands are included, the BME computation here is thus formally equivalent to a GME one, but the potential advantage of BME lies in the possibility to truncate the set of bands to the ones closest in frequency to the spectral range of interest. Indeed, a few-band computation is very fast (of the order of minutes on a personal computer), and has already been successfully applied (and proven reliable) to various waveguide-based systems [30, 31, 41]. Here, we test the reliability of a few-band computation for estimating the quality factor of a high-*Q* cavity, which requires great precision, since the long photon lifetime in these devices is due to extremely delicate destructive interference of radiating modes.

## 4.2. Comparison between numerical methods

In Fig. 5 we show the convergence of the BME-computed *Q* with the number of bands included in the computation. The results for the disorder-less designs are shown as a thick blue line, and compared vs. the corresponding GME, 3D finite-element method (FEM), and 3D finite-difference time-domain (FDTD) results. The FEM values were obtained using a 3D-FEM solver to compute the eigenvalues of the Maxwell equations for a computational volume including the PhC cavity and 2*μ*m of air above and below the membrane. Fine and adaptive meshes were applied to the holes of the device and first-order absorbing boundary conditions were used at the truncated device edges. Convergence was checked versus hole mesh density, device width and device length for attaining the final *Q*-factor of each cavity. The FDTD results were taken from the corresponding references [6, 8]. The BME and GME results were computed using a supercell of length 32*a* (24*a* for the *L*3 cavity) in the *x*-direction (see Fig. 1), and length
$16\sqrt{3}/2a$ (
$10\sqrt{3}/2a$ for the *L*3) in the *y*-direction. As expected, when all bands are included, BME and GME yield the same value, which is also in good agreement with the FEM and FDTD results, even though a general trend appears to be that the predicted *Q*-s are always, in increasing order, given by FEM, FDTD, and GME/BME. The reason for this is likely to be found in the different boundary conditions imposed: in the FEM computation we used first-order absorbing boundary conditions; the imperfect absorption at the boundaries together with the fact that the design of the HS cavity is a global variation of hole positions (i.e. affecting all holes as opposed to the other three designs) is most likely the reason why the FEM result for that cavity is too low. The FDTD simulations, on the other hand, commonly apply perfectly matched layer conditions, which absorb any incident radiation. Finally, since the GME and BME are both mode-expansion methods, periodic boundary conditions are imposed. Regarding computational efficiency, obtaining a converged *Q*-value using 3D-FEM, when running on a single core of a last-generation CPU, takes of the order of tens of hours and 100GB of memory. The same computation with a common 3D-FDTD software would take ≈ 5 − 15 hours and several GB of memory. A similar amount of memory is needed for a GME computation, which, however, is much faster – taking less than an hour. Finally, a few-band BME computation takes of the order of several minutes and requires less than 1GB of memory, illustrating its potential applicability. In addition, contrary to the other methods, the memory required by BME is moderate even for very large structures, allowing for a full simulation of long waveguides and/or various integrated structures (as opposed to splitting the individual components).

## 4.3. Potential applications of the BME

It should be noted that the irregular variation of the BME-computed *Q* with number of bands, which is visible in some of the plots, is expected: adding more bands introduces new radiative components that can interfere either constructively or destructively with the components due to the lower bands. Admittedly, the convergence is slow especially for the highest-*Q* cavities (*A*1 and HS), for which a high number of bands is needed to obtain the quantitatively correct quality factor. However we infer from our simulations that an expansion on a low number of bands predicts relative changes of *Q* upon small perturbations very reliably, even if the prediction of the absolute *Q* value is not converged. This is illustrated by plotting in Fig. 5 the BME-computed *Q* values for 20 disorder realizations, as a function of the number of bands. The curves lie essentially parallel to each other. This shows that, when changing from one disorder realization to another, the relative change in *Q* is the same independently of the number of bands included in the simulation. This conclusion is further supported by the data in Fig. 6, where we plot histograms of the probability of occurrence of *Q*-values, computed over 1000 different disorder realizations for each cavity, when 2 (Fig. 6(a)), 100 (Fig. 6(b)), or all 285 bands (Fig. 6(c)) were included. While the statistical distributions shift to higher-*Q* values in every consecutive panel, it is evident that their shape, and the relative variations between different designs, are already captured by the two-band computation. This brings to one of the main conclusions of the present work: the BME method provides a fast and reliable tool for cavity comparison and optimization, as well as for analysis of robustness vs. disorder. In particular, a few-band computation can be used to explore the vast parameter space of the holes’ positions and sizes and find an optimal design, after which the converged-BME (or any of the other methods) can be used for fine-tuning and obtaining the correct *Q* value. An important remark is that the minimum number of bands taken should not be one but two, since it pays to include *both* guided bands of the waveguide. This is because, despite the fact that the modes of the second guided band are odd with respect to a reflection along the *x–z* plane and thus orthogonal to the modes of the ideal cavities, disorder introduces mixing that can have a pronounced effect on the losses [42]. It is also important to note that the magnitude of disorder *σ* = 0.0014*a* used in this section is close to the one quoted in Ref. [21], and that the HS histogram of Fig. 6(c) is in *quantitatively* good agreement with the experimental results of that work, validating once more the use of the GME in Section 2.

## Acknowledgments

The authors would like to thank Dario Gerace for fruitful discussions, and for kindly providing the FDTD *Q*-factor of the *L*3 cavity. The authors acknowledge the financial support from the
Swiss National Centre of Competence in Research Quantum Photonics and the
Swiss National Science Foundation Projects No.
200020_140406 and
200020_132407.

## References and links

**1. **J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, *Photonic Crystals: Molding the Flow of Light* (Princeton University, 2008).

**2. **S. Noda, M. Fujita, and T. Asano, “Spontaneous-emission control by photonic crystals and nanocavities,” Nat. Photonics **1**, 449–458 (2007). [CrossRef]

**3. **M. Notomi, “Manipulating light with strongly modulated photonic crystals,” Rep. Prog. Phys. **73**, 096501 (2010). [CrossRef]

**4. **J. L. O’Brien, A. Furusawa, and J. Vuckovic, “Photonic quantum technologies,” Nat. Photonics **3**, 687–695 (2009). [CrossRef]

**5. **Y. Akahane, T. Asano, B. Song, and S. Noda, “High-Q photonic nanocavity in a two-dimensional photonic crystal,”Nature **425**, 944–947 (2003). [CrossRef]

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

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

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

**9. **T. Tanabe, M. Notomi, E. Kuramochi, A. Shinya, and H. Taniyama, “Trapping and delaying photons for one nanosecond in an ultrasmall high-Q photonic-crystal nanocavity,” Nat. Photonics **1**, 49–52 (2007). [CrossRef]

**10. **Y. Tanaka, T. Asano, and S. Noda, “Design of photonic crystal nanocavity with Q-factor of ∼ 10^{9},” J. Lightwave Technol. **26**, 1532–1539 (2008). [CrossRef]

**11. **M. Felici, K. A. Atlasov, A. Surrente, and E. Kapon, “Semianalytical approach to the design of photonic crystal cavities,” Phys. Rev. B **82**, 115118 (2010). [CrossRef]

**12. **M. Nomura, K. Tanabe, S. Iwamoto, and Y. Arakawa, “High-q design of semiconductor-based ultrasmall photonic crystal nanocavity,” Opt. Express **18**, 8144–8150 (2010). [CrossRef] [PubMed]

**13. **G. Khitrova, H. Gibbs, M. Kira, S. Koch, and A. Scherer, “Vacuum rabi splitting in semiconductors,” Nat. Phys. **2**, 81–90 (2006). [CrossRef]

**14. **H. Carmichael, *Statistical Methods in Quantum Optics 2: Non-Classical Fields* (Springer, 2008). [CrossRef]

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

**16. **C. Husko, A. D. Rossi, S. Combrie, Q. V. Tran, F. Raineri, and C. W. Wong, “Ultrafast all-optical modulation in GaAs photonic crystal cavities,” Appl. Phys. Lett. **94**, 021111 (2009). [CrossRef]

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

**18. **T. Volz, A. Reinhard, M. Winger, A. Badolato, K. J. Hennessy, E. L. Hu, and A. Imamoglu, “Ultrafast all-optical switching by single photons,” Nat. Photonics **6**, 605–609 (2012). [CrossRef]

**19. **K. Nozaki, A. Shinya, S. Matsuo, Y. Suzaki, T. Segawa, T. Sato, Y. Kawaguchi, R. Takahashi, and M. Notomi, “Ultralow-power all-optical RAM based on nanocavities,” Nat. Photonics **6**, 248–252 (2012). [CrossRef]

**20. **Y. Sato, Y. Tanaka, J. Upham, Y. Takahashi, T. Asano, and S. Noda, “Strong coupling between distant photonic nanocavities and its dynamic control,” Nat. Photonics **6**, 56–61 (2012). [CrossRef]

**21. **Y. Taguchi, Y. Takahashi, Y. Sato, T. Asano, and S. Noda, “Statistical studies of photonic heterostructure nanocavities with an average q factor of three million,” Opt. Express **19**, 11916–11921 (2011). [CrossRef] [PubMed]

**22. **S. L. Portalupi, M. Galli, M. Belotti, L. C. Andreani, T. F. Krauss, and L. O’Faolain, “Deliberate versus intrinsic disorder in photonic crystal nanocavities investigated by resonant light scattering,” Phys. Rev. B **84**, 045423 (2011). [CrossRef]

**23. **Y. Tanaka, T. Asano, Y. Akahane, B.-S. Song, and S. Noda, “Theoretical investigation of a two-dimensional photonic crystal slab with truncated cone air holes,” Appl. Phys. Lett. **82**, 1661–1663 (2003). [CrossRef]

**24. **U. K. Khankhoje, S.-H. Kim, B. C. Richards, J. Hendrickson, J. Sweet, J. D. Olitzky, G. Khitrova, H. M. Gibbs, and A. Scherer, “Modelling and fabrication of gaas photonic-crystal cavities for cavity quantum electrodynamics,” Nanotechnology **21**, 065202 (2010). [CrossRef] [PubMed]

**25. **T. Asano, B.-S. Song, and S. Noda, “Analysis of the experimental Q factors ( 1 million) of photonic crystal nanocavities,” Opt. Express **14**, 1996–2002 (2006). [CrossRef] [PubMed]

**26. **D. Gerace and L. C. Andreani, “Effects of disorder on propagation losses and cavity q-factors in photonic crystal slabs,”Photonics
Nanostruct. Fundam. Appl. **3**, 120–128 (2005). [CrossRef]

**27. **H. Hagino, Y. Takahashi, Y. Tanaka, T. Asano, and S. Noda, “Effects of fluctuation in air hole radii and positions on optical characteristics in photonic crystal heterostructure nanocavities,” Phys. Rev. B **79**, 085112 (2009). [CrossRef]

**28. **Z. Zhang and M. Qiu, “Small-volume waveguide-section high q microcavities in 2d photonic crystal slabs,” Opt. Express **12**, 3988–3995 (2004). [CrossRef] [PubMed]

**29. **H. Takagi, Y. Ota, N. Kumagai, S. Ishida, S. Iwamoto, and Y. Arakawa, “High q h1 photonic crystal nanocavities with efficient vertical emission,” Opt. Express **20**, 28292–28300 (2012). [CrossRef] [PubMed]

**30. **V. Savona, “Electromagnetic modes of a disordered photonic crystal,” Phys. Rev. B **83**, 085301 (2011). [CrossRef]

**31. **M. Minkov and V. Savona, “Effect of hole-shape irregularities on photonic crystal waveguides,” Opt. Lett. **37**, 3108–3110 (2012). [CrossRef] [PubMed]

**32. **N. L. Thomas, Z. Diao, H. Zhang, and R. Houdre, “Statistical analysis of subnanometer residual disorder in photonic crystal waveguides: Correlation between slow light properties and structural properties,” J. Vac. Sci. Technol. B **29**, 051601 (2011). [CrossRef]

**33. **L. C. Andreani and D. Gerace, “Photonic-crystal slabs with a triangular lattice of triangular holes investigated using a guied-mode expansion method,” Phys. Rev. B **73**, 235114 (2006). [CrossRef]

**34. **K. Welna, S. Portalupi, M. Galli, L. O’Faolain, and T. Krauss, “Novel dispersion-adapted photonic crystal cavity with improved disorder stability,” IEEE J. Quantum Electron. **48**, 1177–1183 (2012). [CrossRef]

**35. **T. Uesugi, B.-S. Song, T. Asano, and S. Noda, “Investigation of optical nonlinearities in an ultra-high-q si nanocavity in a two-dimensional photonic crystal slab,” Opt. Express **14**, 377–386 (2006). [CrossRef] [PubMed]

**36. **K. Hennessy, C. Hogerle, E. Hu, A. Badolato, and A. Imamoglu, “Tuning photonic nanocavities by atomic force microscope nano-oxidation,” Appl. Phys. Lett. **89**, 041118 (2006). [CrossRef]

**37. **D. Dorfner, T. Zabel, T. Hrlimann, N. Hauke, L. Frandsen, U. Rant, G. Abstreiter, and J. Finley, “Photonic crystal nanostructures for optical biosensing applications,” Biosens. Bioelectron. **24**, 3688–3692 (2009). [CrossRef] [PubMed]

**38. **J. Jágerská, H. Zhang, Z. Diao, N. L. Thomas, and R. Houdré, “Refractive index sensing with an air-slot photonic crystal nanocavity,” Opt. Lett. **35**, 2523–2525 (2010). [CrossRef] [PubMed]

**39. **S. Vignolini, F. Riboli, D. S. Wiersma, L. Balet, L. H. Li, M. Francardi, A. Gerardino, A. Fiore, M. Gurioli, and F. Intonti, “Nanofluidic control of coupled photonic crystal resonators,” Appl. Phys. Lett. **96**, 141114 (2010). [CrossRef]

**40. **N. Descharmes, U. P. Dharanipathy, Z. Diao, M. Tonin, and R. Houdré, “Observation of backaction and self-induced trapping in a planar hollow photonic crystal cavity,” Phys. Rev. Lett. **110**, 123601 (2013). [CrossRef]

**41. **M. Minkov and V. Savona, “Radiative coupling of quantum dots in photonic crystal structures,” Phys. Rev. B **87**, 125306 (2013). [CrossRef]

**42. **V. Savona, “Erratum: Electromagnetic modes of a disordered photonic crystal [phys. rev. b 83, 085301 (2011)],” Phys. Rev. B **86**, 079907 (2012). [CrossRef]