The cavities studied here are formed in a triangular lattice with pitch a of air-holes of radius R etched in a silicon (n = 3.46) slab of thickness d. In what follows, all lengths will be expressed in units of a, as the quality factor is invariant upon a spatial rescaling of the structure. We however set parameters such that, for the typically used thickness d = 220 nm, the resonant modes lie in the telecommunication band around λ = 1.55 μm. For the simulation of a single structure, we use the GME method38 (see Methods). This method has already proven reliable for modelling high-Q cavities35,37. As a further check of its validity, the results for all final (optimized) structures are verified using the 3-D finite-difference time-domain (FDTD) method40 (see Methods). For each of the optimized designs presented here, we also analyse the probability density of Q e in the presence of fabrication imperfections (and neglecting the absorption loss contribution 1/Q a ). The disorder model is that of Gaussian distributed random fluctuations in the position and radius of each hole, with zero mean and standard deviation σ (which is a measure of the disorder magnitude). Disorder in all holes in the computational cell was included and no hole-hole correlations were taken into account. Irregular hole shapes can in general easily be included, but the resulting effect is well described by an effective radius fluctuation41. The disorder model thus captures well what is accepted as the main source of losses in silicon photonic crystal cavities12,35,36.

The first cavity design we investigate is the widely-employed L3 cavity formed by three missing holes in a row (Fig. 1(a)), with d = 0.55a and R = 0.25a. The quality factor of this cavity has already been optimized2 with respect to shifts of the positions of three neighbouring holes (marked S 1x , S 2x and S 3x in the Figure), by using a simplified approach in which, starting from the unshifted design, each of the three shifts has been varied once while keeping the two others constant. To explore the extent to which this approach is suitable, we compute a full map of the quality factor on a relevant region of the (S 1x , S 2x , S 3x )-space. The map is displayed in Fig. 1(b)–(m). There, in each panel, Q t is plotted as a function of S 2x and S 3x , while the value of S 1x increases from 0.15a to 0.37a in steps of 0.02a when going from panel (b) to panel (m). Technically, these plots already provide a global optimization of the cavity (a clear maximum of Q t can be identified), although performed in the least practical, brute-force way. If applied to the panels of Fig. 1 (though approximately, given the coarse S 1x step used in this figure), the simplified optimization procedure in Ref. 2 leads to the point marked by a white cross in panel (e) (more precisely, S 1x = 0.21a, S 2x = 0.01a, S 3x = 0.23a), i.e. far off the maximum that can be seen in panel (k) at S 1x = 0.33a, S 2x = 0.26a, S 3x = 0.10a. It is also interesting to note that within the range of the plots in Fig. 1(b–m), two maxima are visible - a local one around S 1x = 0.25a, S 2x = 0.09a, S 3x = 0.21a in panel (g) and another one which first appears at S 1x = 0.29a, S 2x = 0.18a, S 3x = −0.10a in panel (i) and then shifts to become the global maximum at S 1x = 0.327a, S 2x = 0.257a, S 3x = 0.116a. In general, an even larger number of local extrema might in principle be present, especially for larger number of parameters, thus making the search for the global extremum more difficult. This highlights the need for a global optimization procedure instead of a more conventional algorithm (e.g. the conjugate gradient) that would almost inevitably find a local rather than a global maximum.

Figure 1 (a): The design of the L3 cavity. For quality factor optimization, shifts of the positions of the five neighbouring holes in the x-direction were introduced, labeled as S 1x , S 2x , S 3x , S 4x and S 5x in the figure. (b)–(m): A parameter scan of the GME-computed quality factor values for different S 1x , S 2x and S 3x , where S 1x starts from 0.15a in panel (b) and increases in multiples of 0.02a in every consecutive panel, up to 0.37a in (m) and S 4x = S 5x = 0 in all panels. Full size image

Ideally, we would like to apply a global, stochastic (since there is no general way to come up with a “good” guess for a starting point) procedure to the problem of optimizing the cavity parameters. Thus, we choose to employ a genetic algorithm which typically relies on a range of evolution-inspired techniques to create consecutive “generations” of cavities each containing better and better “individuals”, until convergence is reached. This type of algorithm proves very efficient especially with increasing number of free parameters. We use the genetic algorithm provided in the MATLAB® Global Optimization Toolbox (see Methods), where we choose the objective function to be the GME-computed quality factor Q t . When applied to the L3 with freedom in S 1x , S 2x and S 3x , the optimal design is found for S 1x = 0.327a, S 2x = 0.257a and S 3x = 0.116a (Fig. 2(a)). This yields Q t = 2.1 × 106 (FDTD: 1.6 × 106), which is an increase by a factor of ≈ 6 as compared to the previously highest value2 of 3.3 × 105 (FDTD: 2.6 × 105), while the mode volume (see Methods) increases from 0.77(λ/n)3 to 0.94(λ/n)3. One obvious advantage of using evolutionary optimization rather than the brute-force parameter scan of Fig. 1 is the precision with which the maximum can be pinpointed; another one is that a few tens of generations with a population of 80 individuals are sufficient to reach this converged configuration (see Methods). Moreover, the design can be further improved if two more shifts (S 4x and S 5x in Fig. 1(a)) are allowed in the optimization, which is still easily handled by the genetic algorithm, although ≈ 100 generations are needed for convergence. In this case, the optimized design is found for S 1x = 0.337a, S 2x = 0.270a, S 3x = 0.088a, S 4x = 0.323a and S 5x = 0.173a (Fig. 2(b)) and has Q t = 5.1 × 106 (FDTD: 4.2 × 106) with mode volume V = 0.95(λ/n)3, i.e. an increase in Q t by one order of magnitude compared to the previous optimal values2,33, with an increase in the mode volume comparable to the three-shift case. The resonant frequency of the modes is at for both designs, which is slightly lower than the frequency of the unmodified L3 (i.e. with the same d/a and R/a but with no hole shifts).

Figure 2 (a)–(b): Electric field (E y ) profiles of the two optimized L3 designs; the holes that were displaced are marked in red. (a): S 1x = 0.327a, S 2x = 0.257a, S 3x = 0.116a. (b): S 1x = 0.337a, S 2x = 0.270a, S 3x = 0.088a, S 4x = 0.323a, S 5x = 0.173a. (c): Dependence of Q t on S 2x and S 3x for S 1x = 0.327a; the shifts in the design of panel (a) are marked by a white cross. (d): Histograms showing the probability of occurrence of different Q e -values in the design of panel (a), for two different disorder magnitudes: σ = 0.003a (red) and σ = 0.0015a (blue). The black line indicates the ideal Q t . (e): Same as (d), for the design of (b). (f): Dependence of Q t on the overall radius and the slab thickness, for the values of S 1x − S 5x corresponding to the optimal design in (b). On the right of the dashed line, the slab becomes multi-mode at the cavity frequency. Full size image

The choice of just these few free parameters is not unique, but appears optimal. Attempts were made at including other position shifts or radii variations of the holes surrounding the cavity defect, which however did not bring to a significant further increase in Q t . Having only a few hole shifts as free parameters results in a technologically friendly structure and in a more compact cavity defect, characterized by a much smaller footprint on the PhC, than waveguide-based ultrahigh-Q designs4,5,7. In addition, the present designs are as robust to fabrication imperfections as any other high-Q PhC cavity37, as can be inferred from Fig. 2(c)–(f). In panel (c) we plot, for the three-shift L3, the dependence of Q t on S 2x and S 3x as in Fig. 1, but for the value S 1x = 0.327a corresponding to our optimal design (the white cross indicates where the design lies with respect to S 2x and S 3x ). In this plot, we observe that the width of the maximum is larger than the typical uncertainty in the hole positions (smaller than 0.005a for Si42). Furthermore, in Fig. 2(d)–(e) we show, respectively for the three- and five-shift design, the computed probability of occurrence of Q e -values in presence of disorder. Each of these histograms was obtained by simulating 1000 disordered realisations of the corresponding cavity design. The blue plot in panel (d) in particular shows that for a state-of-the-art disorder magnitude σ = 0.0015a (i.e. about 0.6 nm12,35, when assuming a = 400 nm in a silicon slab), the average value lies at about Q e = 2.5 × 106, i.e. quality factors one order of magnitude larger than the previous theoretical maximum can be expected in practice, highlighting the significance of the design optimization. Finally we note that, for a given set of optimal values of the S nx parameters, the designs are also robust to small changes in the overall hole radius R and slab thickness d, that can originate from an offset in the fabrication process and/or be introduced on purpose in order to e.g. tune the resonant frequency to a desired value. To show this, in panel (f) we plot the value of Q t obtained by varying R and d while keeping the shifts S 1x − S 5x constant and set to the values obtained for the optimal design computed at d = 0.55a and R = 0.25a. We observe that Q t > 4 × 106 for a range of R and d values which is much larger than the fabrication uncertainty and which allows fine-tuning of the frequency. For certain values of d and R (to the right of the dashed line in the Figure), higher-order guided modes of the slab become non-negligible. More precisely, a TM-like band of the regular PhC becomes resonant with the cavity mode and introduces a new loss channel38. We point out that, while Q t appears to systematically increase with d in the single-mode region, it drops rapidly as soon as d increases into the multi-mode region. An analogous trend of the loss rates as a function of R/a and d/a is expected for the other cavity designs discussed in this study. In principle, one could consider d/a and/or R/a as free parameters in the optimization, but in that case setting a target wavelength, for a fixed value of d that might arise from technological requirements, becomes more difficult.

Often, obtaining the highest possible theoretical Q t is not the main goal of optimization. In fact, when Q t gets above a limit set by the material and the fabrication process (currently ≈ 5 × 106 in silicon12), the experimentally measured Q e is always dominated by losses due to disorder and/or absorption and so weakly affected by further increase in Q t . The potential of an automated optimization procedure is therefore best exploited when applied to other attractive properties. One example consists in maximizing Q t while having the smallest possible mode volume, so that Q t /V is as high as possible, since the latter is a figure of merit for applications in both cavity QED25,27,28 and non-linear optics26. With this in mind, the second design we focus on is the H0 cavity6,9 (sometimes also named “zero-cell” or “point-shift”), namely the simple defect cavity with the smallest known mode volume.

The design of the H0 is shown in 3(a); the defining defect is the shift of two holes away from each other (S 1x ), the thickness of the slab taken here is d = 0.5a, while the hole radius is R = 0.25a. For the optimization, we also use the consecutive shifts S 2x − S 5x , as well as two shifts in the vertical direction, S 1y and S 2y . Using S 1x , S 2x , S 3x , S 1y and S 2y , the cavity has already been optimized6 (following the same approach already discussed for the L32) to a quality factor of Q t = 2.8 × 105, with a corresponding mode volume V = 0.23(λ/n)3. Here, we improve on this result by on one hand using the genetic optimization and on the other by including S 4x and S 5x . It should be mentioned that in the optimization, an allowed range of variation for each parameter is set. For the H0, we find that the maximum allowed S 1x is a very important parameter, increasing which produces several different optimized designs. All of those are interesting, as increasing S 1x increases the mode volume but also the Q t of the cavity. More precisely, the designs in 3(b)–(d) were obtained by imposing the following restrictions in the genetic algorithm: S 1x ≤ 0.25a, S 1x ≤ 0.3a and S 1x ≤ 0.4a, respectively. The ensuing optimal parameters [S 1x , S 2x , S 3x , S 4x , S 5x , S 1y , S 2y ] are as follows: [0.216a, 0.103a, 0.123a, 0.004a, 0.194a, −0.017a, 0.067a] (panel (b)); [0.280a, 0.193a, 0.194a, 0.162a, 0.113a, −0.016a, 0.134a] (panel(c)); and [0.385a, 0.342a, 0.301a, 0.229a, 0.116a, −0.033a, 0.093a] (panel(d)). The corresponding quality factors are Q t = 1.04 × 106 (FDTD: 1.04 × 106), Q t = 1.88 × 106 (FDTD: 1.66 × 106), and, remarkably, Q t = 8.89 × 106 (FDTD: 8.29 × 106), while the respective mode volumes are V = 0.25(λ/n)3, V = 0.34(λ/n)3 and V = 0.64(λ/n)3. The first among these three designs (panel (b)) has a mode volume only slightly larger than the previous most optimal design6, combined to a quality factor almost four times larger. The last of the three designs (panel (d)) instead shows a more significant increase of the mode volume, but associated to an almost 30-fold increase of Q t with respect to the value obtained in Ref. 6. The resonance frequencies of the three designs decrease with the increase of V and are , and , respectively, while the original cavity with S 1x = 0.14 (and no other shifts) of Ref. 9 has .

Similarly to what we have done above for the L3 cavity, in Fig. 3(e)–(g) we present the probability of occurrence of Q e values, computed using 1000 random disorder realizations for each design and each disorder magnitude. From these histograms it appears clearly that, even though design 3 has the highest theoretical Q t /V, it might not be the best choice in practice. According to Eq. (1) in fact, depending on the amount of disorder, the maximum value of the actual ratio Q e /V will in general be achieved for a design having an intermediate value of Q t . For example, in the case with σ = 0.003a (red curves in panels (e)–(g)), the average values of Q e (neglecting absorption) computed from the simulations are 3.97 × 105, 5.23 × 105 and 6.49 × 105, respectively, meaning that the highest Q e /V would in practice be achieved by design 1. On the other hand, for the smaller disorder σ = 0.0015a (blue curves in panels (e)–(g)), the corresponding average values of Q e are 7.22 × 105, 1.12 × 106 and 2.02 × 106 and the highest average Q e /V is achieved by design 2. For both values of σ, the expected Q e /V for the five-shift L3 cavity is lower than that for any of the three H0 designs, thus the latter should be the cavity of choice for applications where Q/V is the most important figure of merit. We expect that design 3 will dominate for yet smaller values of σ, which however appear to be currently not achievable in practice. Additionally, design 3 could be the preferred choice for purely manufacturing reasons if quantum dots are to be inserted in the cavity, as the field maxima are further away from the hole edges than in the other two cavities. An important final remark is that cavities 1 and 2 have already been fabricated39 and experimental quality factors consistent with the σ = 0.003a disorder histograms have been measured (unloaded Q e = 485, 000 for design 2), resulting in optical bistability with ultra-low threshold power in silicon.

Figure 3 (a): The design of the H0 cavity. For quality factor optimization, shifts of the positions of the five neighboring holes in the x-direction and the two neighboring holes in the y-direction were introduced, labeled as S 1x , S 2x , S 3x , S 4x , S 5x , S 1y and S 2y in the Figure. (b)–(d): Electric field (E y ) profiles of three optimized designs with increasing Q and V. (b): S 1x = 0.216a, S 2x = 0.103a, S 3x = 0.123a, S 4x = 0.004a, S 5x = 0.194a, S 1y = −0.017a, S 2y = 0.067a. (c): S 1x = 0.280a, S 2x = 0.193a, S 3x = 0.194a, S 4x = 0.162a, S 5x = 0.113a, S 1y = −0.016a, S 2y = 0.134a. (d): S 1x = 0.385a, S 2x = 0.342a, S 3x = 0.301a, S 4x = 0.229a, S 5x = 0.116a, S 1y = −0.033a, S 2y = 0.093a. (e): Histograms showing the probability of occurrence of different Q e -values in the design of panel (b), for two different disorder magnitudes: σ = 0.003a (red) and σ = 0.0015a (blue). The black line indicates the value of Q t . (f)–(g): Same as (e), for the designs of (c)–(d). The σ = 0 line in panel (g) is not visible as Q t occurs beyond the axis boundary. Full size image

Another interesting cavity with potential applications for polarization-entangled photon generation43,44 and quantum dot spin readout45 is the H19,13,34 – also known as “single point-defect cavity” – formed by one missing hole in the lattice (Fig. 4(a)). The modes of this cavity preserve the underlying hexagonal symmetry and for a wide range of parameters, the fundamental (lowest-frequency) resonance is given by two degenerate dipole modes. Based on the electric field polarization in the center of the cavity, those are usually referred to as the y-polarized (Fig. 4 (b)) and x-polarized (Fig. 4 (c)) modes, although in fact the electric field of each of those modes also has a non-vanishing component oriented in the orthogonal direction; however, since the near-field in the very center of the cavity as well as the far-field emission in the z-direction (perpendicular to the slab plane) are both truly y (correspondingly x) polarized, this labeling is in many cases appropriate, in particular for applications in which a quantum dot is placed in the center of the cavity. Here, we optimize the quality factor of the dipole mode. We take d = 0.55a and R = 0.23a and the parameters used for design optimization, labeled S 1 , S 2 and S 3 in Fig. 4(a), are an increase in the side-length of the three consecutive hexagonal “rings” around the cavity (which is also equivalent to the increase of the distance from a vertex of a hexagon to the cavity center). The previously most-optimal design was achieved using only the holes at the vertices of the hexagons (but including variations of the hole radii) and has a moderate Q t = 6.2 × 104 and a mode volume V = 0.47(λ/n)38. Here, using the genetic algorithm with the shifts as outlined in Fig. 4(a), we find an optimal design at S 1 = 0.213a, S 2 = 0.070a and S 3 = −0.009a, with Q t = 1.05 × 106 (FDTD: 0.97 × 106) and V = 0.62(λ/n)3, i.e. we find a 19-fold increase in Q t coupled to an increase in V by 32%. For this cavity as well, the disorder analysis (panel (d)) suggests that Q e -values close to a million can be expected experimentally in state-of-the-art silicon systems, i.e. more than an order of magnitude larger than the previous theoretical values. The modes lie at a frequency , which is, as expected, slightly lower than that of the unmodified cavity, .

Figure 4 (a): The design of the H1 cavity. The size of the first three hexagonal “rings” of holes is varied for quality factor optimization, with the increase of the distance from a hexagon vertex to the center of the cavity given by S 1 , S 2 and S 3 (marked). (b): Electric field (E y ) profile of the y-polarized and (c): electric field (E x ) profile of the x-polarized mode, for the optimal design S 1 = 0.213a, S 2 = 0.070a, S 3 = −0.009a. (d): Histograms showing the probability of occurrence of different Q e -values for two different disorder magnitudes: σ = 0.003a (red) and σ = 0.0015a (blue). The black line indicates the ideal Q t . (e): Histograms showing the probability of occurrence of the wavelength splitting Δλ between the two modes, which are degenerate in the disorder-less cavity. Full size image

We note that while the degeneracy of the two dipole modes is an attractive feature of the H1, it is lifted by disorder. This is why, in panel (e) of Fig. 4, we study the probability of occurrence of the splitting between the modes, based on the 1000 disorder realizations that were used for the disorder analysis in panel (d). It is important to note that there is no absolute way to define an x − y reference frame, as three equivalent frames (rotated 60° from one another) exist due to the hexagonal symmetry of the cavity. This symmetry is broken if the cavity presents preferential orientations of the axes (e.g. introduced by lithography). In the case where only random disorder is considered, x- and y-polarized modes can turn out to be oriented along either of the three x − y reference frames. Thus, what we plot in panel (e) is only the difference between the resonant wavelengths of the higher- and the lower-frequency modes, without any reference to their polarization. What is important to note is that the splitting is of the order of hundreds of picometers, i.e. two orders of magnitude larger than the linewidth corresponding to the typical Q e -s (see top x-axis of panel (d)). This implies that for applications for which overlap in frequency between the two modes is needed, some form of post-fabrication tuning is required, the possibility for which has already been demonstrated43,46. This can also be combined with an additional modulation of the neighboring holes which increases the emission in the vertical direction47, which would decrease the Q e and make the tuning easier while simultaneously increasing the intensity of the collected radiation, which is beneficial for entangled-photon generation44,48.

We finally address the hexapole mode of the H1 cavity, that typically lies at higher frequency than the dipolar modes studied above. This mode has previously been optimized to Q t = 1.6 × 10611,13 by varying only S 1 , with S 2 = S 3 = 0 in the sketch of Fig. 4(a). Here, we run a global optimization of Q t by varying the three shifts, with d = 0.55a and R = 0.22a. We could improve the previous value to Q t = 3.2 × 106 (FDTD: 3.1 × 106), obtained for the optimal values S 1 = 0.271a, S 2 = 0.039a and S 3 = 0.018a. The resonance frequency of this mode is , while in the unmodified cavity of the same d and R, this mode is not present.

The figures of merit of the seven designs that were optimized above are summarized in Table I. Notice that for the higher disorder σ 1 , the values of 〈Q e 〉 do not differ dramatically among the different designs. This reflects the fact that, as has been shown in a previous analysis37, in the disorder-dominated regime (when ), the losses are nearly design-independent.