Statistics of the disorder-induced losses of high-Q photonic crystal cavities

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⇥ 104 to 7.5⇥ 107, 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 OCIS codes: (230.5298) Photonic crystals; (220.4241) Nanostructure fabrication; (140.3948) Microcavity devices. 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). 3. M. Notomi, “Manipulating light with strongly modulated photonic crystals,” Rep. Prog. Phys. 73, 096501 (2010). 4. J. L. O’Brien, A. Furusawa, and J. Vuckovic, “Photonic quantum technologies,” Nat. Photonics 3, 687–695 (2009). 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). 6. B. Song, S. Noda, T. Asano, and Y. Akahane, “Ultra-high-Q photonic double-heterostructure nanocavity,” Nat. Mater. 4, 207–210 (2005). 7. D. Englund, I. Fushman, and J. Vuckovic, “General recipe for designing photonic crystal cavities,” Opt. Express 13, 5961–5975 (2005). 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). 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). 10. Y. Tanaka, T. Asano, and S. Noda, “Design of photonic crystal nanocavity with Q-factor of ⇠ 109,” J. Lightwave Technol. 26, 1532–1539 (2008). #193647 $15.00 USD Received 9 Jul 2013; revised 9 Oct 2013; accepted 31 Oct 2013; published 11 Nov 2013 (C) 2013 OSA 18 November 2013 | Vol. 21, No. 23 | DOI:10.1364/OE.21.028233 | OPTICS EXPRESS 28233 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). 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). 13. G. Khitrova, H. Gibbs, M. Kira, S. Koch, and A. Scherer, “Vacuum rabi splitting in semiconductors,” Nat. Phys. 2, 81–90 (2006). 14. H. Carmichael, Statistical Methods in Quantum Optics 2: Non-Classical Fields (Springer, 2008). 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). 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). 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). 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). 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). 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). 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). 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). 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). 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). 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). 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). 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). 28. Z. Zhang and M. Qiu, “Small-volume waveguide-section high q microcavities in 2d photonic crystal slabs,” Opt. Express 12, 3988–3995 (2004). 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). 30. V. Savona, “Electromagnetic modes of a disordered photonic crystal,” Phys. Rev. B 83, 085301 (2011). 31. M. Minkov and V. Savona, “Effect of hole-shape irregularities on photonic crystal waveguides,” Opt. Lett. 37, 3108–3110 (2012). 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). 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). 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). 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). 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). 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). 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). 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). 40. N. Descharmes, U. P. Dharanipathy, Z. Diao, M. Tonin, and R. Houdré, “Observation of backaction and selfinduced trapping in a planar hollow photonic crystal cavity,” Phys. Rev. Lett. 110, 123601 (2013). 41. M. Minkov and V. Savona, “Radiative coupling of quantum dots in photonic crystal structures,” Phys. Rev. B #193647 $15.00 USD Received 9 Jul 2013; revised 9 Oct 2013; accepted 31 Oct 2013; published 11 Nov 2013 (C) 2013 OSA 18 November 2013 | Vol. 21, No. 23 | DOI:10.1364/OE.21.028233 | OPTICS EXPRESS 28234 87, 125306 (2013). 42. V. Savona, “Erratum: Electromagnetic modes of a disordered photonic crystal [phys. rev. b 83, 085301 (2011)],” Phys. Rev. B 86, 079907 (2012).


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][2][3][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][6][7][8][9][10][11][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][17][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][6][7][8][9][10][11][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][24][25]).Understanding the disorder effects has spurred some theoretical [25][26][27] and experimental [21,22,25] research, which has most notably shown -for two specific designs (the L3 cavity and the heterostructure cavity) -that the disorder-averaged value of the disorder-induced cavity losses scales as s 2 , where s 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 L3 [5], A3 and A1 [8,9], heterostructure (HS) [6], H1 [28,29] and H0 [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.(f): Optimized H0, where the 8 holes marked in red are shifted symmetrically: the two shifts in the x-direction are 0.14a, and 0.06a, while the two shifts in the y-direction are 0.04a and 0.02a (cf.[12]).

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 L3 cavity of Fig. 1(a), also studied in Refs.[22,26], has R = 0.3a 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.16a, while their radius is decreased by 0.06a.The A3 and A1 cavities (Fig. 1(b) and Fig. 1(c)) are taken as originally defined in Ref. [8], with R = 0.25a, d = 0.472a for the A3 and R = 0.257a, d = 0.486a for the A1.The A3 cavity is formed by first introducing a one-dimensional waveguide defect of width 0.9a, 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.0278a.The A1 cavity is similar, but the width of the waveguide is 0.98a, and three types of hole shifts are introduced: by d y = 0.0214a (red holes in Fig. 1(c)), by 2d y /3 (black holes) and by d y /3 (grey holes).Next, we consider the heterostructure cavity first introduced in Ref. [6], with 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 H1 cavity formed by one missing hole, with d = 0.5a and R = 0.3a 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 H0 as defined in Ref. [12], where d = 0.6a and R = 0.26a, 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 s .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 3n 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 s = 0.001a (relevant for state-of-the-art Si structures [21,32]) to s = 0.015a, i.e. includes values relevant for different systems like GaAs, InGaAs or InGaAsP.The measured Q of a cavity can be written as

Statistics of the disorder-induced losses
where Q i is the theoretically predicted Q-factor for the particular cavity design in the absence of disorder (also called "intrinsic"), Q d is the quality factor associated to the disorder-induced losses and Q a accounts for any absorption (e.g.due to free carriers or water condensation).The absorption term Q a 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 Q i .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 Q ⌘ whUi/hPi, where w is the cavity frequency, hUi is the total electromagnetic energy of the cavity mode, and hPi is the radiation power.The physical assumption entering Eq. ( 1) (neglecting Q a ) 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 destructive.In particular, this interference can be to a good approximation neglected when one of the channels is strongly dominating, but is important when Nevertheless, in our analysis we use Eq. ( 1), bearing in mind that Q d is, in general, a purely phenomenological quality factor that accounts for the difference between Q i and the measured Q, and study the validity of interpreting Q d as a quality factor associated to disorder-induced scattering losses which do 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 d , which we show are possible for particular disorder realizations.Physically, this reflects the possibility for the disorder to cancel -due to destructive interference -some of the intrinsic radiative components, rather than introduce more losses.The possibility for negative Q d -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 A1 and the A3), 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, where G are reciprocal lattice vectors in the slab-plane, defined with respect to a supercell that contains the structure to simulate, a = 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, a) are then obtained by diagonalization of a matrix with elements defined in Ref. [33].The size of the computational supercell taken here is between 16a and 32a in the x-direction and between 14 p 3a/2 and 20 p 3a/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 max = 2 2p 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 s = 0.001a, 0.0014a, 0.002a, 0.003a, 0.004a and 0.005a, and for 400 disorder realizations for each cavity and each s .For all cavities and disorder magnitudes, the mean of the deviation from the ideal-cavity frequency w 0 , hw w 0 i, is at least one order of magnitude smaller than the standard deviation d (w) = (hw 2 i hwi 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 s 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 H1, where two modes with frequencies w 1 and w 2 exist, we define d (w) = (d (w 1 ) + d (w 2 ))/2).For some practical applications, having the smallest possible d (w) 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 h1/Q d i and standard deviation d vs. s 2 for the different cavities, where 1/Q d was computed from Eq. ( 1).The scaling h1/Q d i µ s 2 has already been numerically established for the L3 [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 s 2 for the A1 and HS in the range we consider, while for the other designs, the dependence is linear for all but the smallest values of s .
An intuitive model for the disorder-induced scattering, which matches to some extent the computed results, is to assume that 1/Q , where d d is a random variable drawn from a Gaussian distribution of zero mean and standard deviation s .Thus, h1/Q d i = A c s 2 with A c a design-specific constant, as observed in Fig. 2 . the standard deviation also scales linearly with s 2 , as observed in Fig. 2(c).However, there are two notable deviations from this simplified expectation.One is the fast increase of d (1/Q d ) for the A3, L3, H0 and H1 for low s .This, however, occurs in the s range for which h1/Q d i 1 >⇡ Q i ; in this situation, as discussed earlier, Q d 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 d ) > h1/Q d i, meaning that negative values for 1/Q d 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 d (1/Q d )/h1/Q d i = p 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 L3, A3, A1, HS, H1 and H0, respectively.This means that d d is actually not Gaussian-distributed, as hd 4 d i < 3hd 2 d i 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 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 d distribution.
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 c for the different cavities are 0.159, 0.121, 0.107, 0.101, 0.273 and 0.230 for the L3, A3, A1, HS, H1 and H0, respectively.With the exception of the HS, these decrease with increasing cavity mode volume; the HS has a volume slightly lower than the A3 and A1 combined with the lowest A c .Additionally, a cavity with a lower A c 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/Q d , plotted for s = 0.0014a, s = 0.005a and s = 0.015a in Fig. 3(a)-(c), we find that the difference is marginal whenever h1/Q d i 1 < Q i .The latter holds for all three values of s for the A1 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 A3 plot in Fig. 2(a) is visibly higher than of the A1 and HS, already for s = 0.005a (panel (b) of Fig. 3) it can be seen that the three distributions become very close, illustrating that the difference in h1/Q d i is actually minor.This reflects the fact that if disorder-induced losses are dominating, which is already the case for the A3 with s = 0.005a, for which Q i = 1.5⇥10 6 and h1/Q d i 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 i is more than two orders of magnitude smaller than the A1 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 s , when Q i is important.However, when Q d 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 h1/Q d i µ s 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 L3, A1, A3 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 = 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 s as in Fig. 3, the distribution of P is always peaked around 1, and all values are above 0.999 for s = 0.0014a, 0.99 for s = 0.005a, and 0.9 for s = 0.015a.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, A1 and A3 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 d in one cavity vs. 1/Q d 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 s ) is indeed found, with a correlation coefficient as high as 0.98 between the A1 and the HS for s = 0.015a.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 d = 1/Q des + 1/Q wm , where 1/Q des 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/Q wm reflects the fact that the waveguide modes themselves become lossy.Thus, 1/Q des is design-dependent, while 1/Q wm is largely ubiquitous, especially since the near-field profiles of Fig. 1 are very similar; the strong correlations with increasing s then show that the losses become increasingly waveguide-mode dominated.Finally, we note that the A1-HS correlation is very high even for s = 0.0014a, which highlights the fact that those two defect cavities are extremely similar.

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 i ), the statistics of 1/Q d are to a good approximation design-independent.Thus, optimizing Q i is important only if the measured Q (or, better, the hQi measured over many samples) is not much lower than Q i .If that is the case, there are just two ways to significantly increase the measured Q.One is to decrease the magnitude of disorder, which for all studied designs is related to the extrinsic losses via h1/Q d i µ s 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 d 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 Q = 7.4 ⇥ 10 6 -the average value with s = 0.001a -even within the s = 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 i 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.
Our results are in line with the idea -commonly acknowledged within the communitythat 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 i 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][37][38][39][40].

The Bloch-mode expansion
Since four of the cavities studied here can be considered as a perturbation over a PhC waveguide (including the L3 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.

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 finitedifference 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 32a (24a for the L3 cavity) in the x-direction (see Fig. 1), and length 16 p 3/2a (10 p 3/2a for the L3) 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 firstorder 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).

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 (A1 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 s = 0.0014a 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.

Fig. 1 .
Fig. 1.Designs and profiles (computed with the guided-mode expansion) of the ycomponent of the electric field, E y (r), for the different cavities studied here.(a) Improved L3 (see text); (b) A3, where the red holes are shifted outwards; (c) A1 with the red, black, and grey holes correspondingly shifted outwards by d y , 2d y /3 and d y /3; (d) Heterostructure, where the lattice constant is (420/410)a in the region between the dashed lines and (415/410)a in the two regions between a dashed and a dashed-dotted line; (e) Optimized H1: the holes marked in red are shifted outward by 0.12a and their radius is decreased by 0.06a, the ones marked black have a radius decreased by 0.03a, and the ones marked grey are shifted inward by 0.26a.Only the y-polarized of the two degenerate modes is plotted;(f): Optimized H0, where the 8 holes marked in red are shifted symmetrically: the two shifts in the x-direction are 0.14a, and 0.06a, while the two shifts in the y-direction are 0.04a and 0.02a (cf.[12]).

Fig. 2 .
Fig. 2. For the six different cavities, (a): Dependence of the standard deviation in the resonant mode frequency w on the disorder magnitude s ; (b) and (c): Dependence of the mean and standard deviation of 1/Q d on s 2 .In (a) and (b), the solid lines are determined from a linear interpolation, while in (c) the dashed lines only connect the individual points and serve as a guide for the eye.Inset in (c): zoom-in over the low-s region.

Fig. 3 .
Fig. 3. (a)-(c): Histograms showing the occurrence of 1/Q d for the six different cavities, computed using GME for 1000 disorder realizations with (a): s = 0.0014a, (b): s = 0.005a, and (c): s = 0.015a.(d)-(f): Histograms showing the occurrence of Q with the same values of s as in (a)-(c).The insets in (d) and (e) are a zoom-in over the low-Q range, where the L3, H0 and H1 values are sitting.

#Fig. 4 .
Fig.4.For the three waveguide-based cavities (A1, A3 and HS), each panel shows 1/Q d in one cavity vs. 1/Q d in another for the same underlying disorder configuration, with 1000 different configurations.First row: HS vs. A1; second row: A1 vs. A3; third row: HS vs. A3.Across columns, s changes from 0.0014a to 0.005a to 0.015a.The correlation coefficient r is indicated in the top right corner of each panel.Note that the range of the x-axis is twice larger than that of the y-axis due to the larger A3 spread.

Fig. 5 .
Fig.5.Convergence of the quality factor with the number of bands included in the BME for four different cavities; thick blue line shows the disorder-less case, while the remaining lines represent 20 different disorder realizations with s = 0.0014a, a value relevant to stateof-the-art silicon based systems.The FDTD, FEM and GME results for the disorder-less case are also given for reference.

Fig. 6 .
Fig. 6.Histograms showing the occurrence of Q for the four cavities, computed for 1000 disorder realizations with s = 0.0014a.On the scale of the main plots, the L3 histogram would appear very narrow, and thus it is shown in an inset.(a) two-band BME results; (b) 100-band BME results; (c) GME results, which correspond to all-band BME.