Global optimization of an encapsulated Si/SiO\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2 L3 cavity with a 43 million quality factor

We optimize a silica-encapsulated silicon L3 photonic crystal cavity for ultra-high quality factor by means of a global optimization strategy, where the closest holes surrounding the cavity are varied to minimize out-of-plane losses. We find an optimal value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q_c=4.33\times 10^7$$\end{document}Qc=4.33×107, which is predicted to be in the 2 million regime in presence of structural imperfections compatible with state-of-the-art silicon fabrication tolerances.


L3 cavity optimization
We consider a silica-encapsulated silicon PC slab with a hexagonal lattice of holes of radii r = 100 nm, lattice parameter a = 390 nm and thickness d = 220 nm. Refractive indices of SiO 2 and Si are taken at telecom wavelengths, i.e., n SiO 2 = 1.44 and n Si = 3.47 , respectively. A L3 cavity is introduced by removing three holes along the ŴK direction of the lattice. In order to optimize the quality factor Q c of its fundamental mode, we adopt a global optimization approach in which only the closest holes surrounding the cavity are varied, in size r → r + dr and position (x, y) → (x + dx, y + dy) , to reduce out-of-plane losses. This technique has been extremely successful during the last few years to reach record theoretical and experimental quality factors for a wide variety  20,25,27,35,[37][38][39] . For a single optimization step, the quality factor of a given configuration must be efficiently computed. To this purpose, we use the guided mode expansion (GME) method 40 which, for PC slab structures, has proven to be as predictive as first-principle commercial Maxwell simulators like FDTD, while being several orders of magnitude faster (see methods). This allows to simulate hundreds of thousands configurations as needed by most optimization algorithms. In addition to a Maxwell solver, one needs an efficient global optimization algorithm to find the global maximum of the quality factor in the complex landscape in parameter space. For this, non-gradient-based algorithms have empirically proven to be very effective in the past. One in particular -the particle swarm (PS) algorithm -was successfully used in past PC optimization works and turned out to be very effective also in the present analysis. As an additional advantage, as discussed in what follows, the PS algorithm is embarrassingly parallel and can be run on a multicore architecture with dramatic computational advantage. We show in Fig. 1 the schematic representation of the holes to be considered in the optimization procedure, where mirror symmetry with respect to the planes x = 0 and y = 0 is assumed. Thereby, we end up with a total of 53 optimization parameters. However, after 1400 iterations of the PS algorithm we have noted that the most relevant parameters for increasing Q c are those highlighted in red in Fig. 1. This preliminary analysis allowed us to reduce the dimension of the optimization parameter space from 53 to 27, and decrease the number of function evaluations required by the algorithm to converge. We summarize in Table 1 our final results where is the linear mode volume and is the non-linear one 41 , with ǫ(r) representing the dielectric function of the system and E(r) the electric field of the cavity mode. An optimal solution of Q c = 4.33 × 10 7 (computed with FDTD 42 ) is found after 806200 function evaluations, leading to an improvement of four orders of magnitude with respect to the non-optimized cavity.
While this theoretical quality factor is around 50% of the record value achieved in glass-clad hetero-structured cavities 43 , it corresponds to the largest reported for the encapsulated L3 cavity. Moreover, the short-length mode confinement provided by the L3 geometry is more convenient for applications in ultra-compact devices. It is noteworthy that, Q c is maximized at the expense of the linear and non-linear mode volumes, in contrast to previous optimizations of the L3 cavity 25 . Nevertheless, we still get extremely large enhancement factors Q c /V l and Q 2 c /V 2 nl which are in the 10 7 and 10 13 regimes, respectively. The increase of the mode volume is clearly seen in the Fig. 2, where we plot the near-field intensity distribution of the fundamental cavity mode in the middle of   www.nature.com/scientificreports/ the slab, for the non-optimized cavity, Fig. 2(a), and the optimized one, Fig. 2(b). The holes which are actually varied are represented by magenta circles in Fig. 2(b). The optimal parameters of the cavity as well as the far-filed projection of the near-field components are reported in the Supplementary Information. The same optimization strategy can be directly applied to the air-bridge silicon L3 cavity within the same parameter space of dimension 27. For this configuration, we have obtained an FDTD quality factor Q c = 1.91 × 10 8 . Such Q c is around 20 times larger than the previous record obtained with deep neural networks 36 , where a different set of parameters were considered to optimize the objective function. While our optimization requires a much larger number of evaluations to find the maximum of the objective function, it clearly shows that there is still considerable room for further improvement of these figures of merit by better identifying the optimization parameters. Detailed results for the Si/Air L3 cavity are given in the Supplementary Information.

Disorder analysis
Realistic photonic devices are always subject to a small amount of intrinsic disorder, coming from unavoidable imperfections introduced at the fabrication stage. We model such effect by considering random Gaussian fluctuations in all hole positions and radii of our PC, where the standard deviation of the Gaussian probability distribution σ is taken as the disorder parameter [44][45][46] . Results of this analysis are shown in Fig. 3, where the averaged cavity quality factor Q c , computed over 100 independent disorder realization of the system, is plotted as a function of σ . Typical tolerances in silicon state-of-the-art fabrication techniques range between σ = 0.001a and σ = 0.002a 31,47 , which leads to an averaged Q c in the 2 million regime. Previous measured Q values on 2D PC slab cavities with low index claddings range from 0.6 × 10 6 (silica cladding) 48 to 1 × 10 6 (glass cladding) 43 , with estimated disorder magnitudes larger than σ = 0.002a . However, the confinement mechanism is based on hetero-structured PCs waveguides or effective photonic potentials which usually require a long-length cavity region, i.e., more than 20 periods of the underlying photonic lattice.

Conclusions
In conclusion, we have optimized a silica-encapsulated silicon L3 cavity by means of a global optimization strategy, where the closest holes surrounding the cavity are varied to decrease out-of-plane losses. We have found a value of Q c = 4.33 × 10 7 , corresponding to outstanding result given the short-length of the confinement region. To better relate our optimal design to realistic samples, we have also studied the effects of intrinsic disorder. When considering typical tolerances in modern fabrication techniques, the averaged quality factor of the optimized cavity remains in the 2 million regime, which is comparable to previous measurements in fully embedded heterostructured and photonic potential based designs. Apart from setting a new record for the encapsulated L3 cavity, our results open the way to a new class of optimized designs in low-index-contrast materials, such as AlN, GaN or Si 3 N 4 , holding great promise for nonlinear optical enhancement, sensing, and solid-state quantum optics.

Methods
The guided mode expansion. We employ the guided mode expansion method (GME) 40 as our main PC solver for estimating the cavity mode frequencies and quality factors at the optimization stage. The expansion was carried out in a supercell of dimensions 20a × 7 √ 3a with 8517 plane waves and one TE guided mode. During the optimization procedure we mainly focus on the growing direction of Q c within the high-dimensional space of parameters. Therefore, only one k point is needed (in the Brillouin zone of the supercell) to compute the real R{ω c } and imaginary I{ω c } parts of the mode frequencies ω c . The cavity quality factor, associated to out-of-plane losses, is thus defined as 40 Note that, while a single k point describes the topology of the optimization space with accuracy (which is enough for the PS algorithm), it does not give the actual value of Q c . This is why the GME estimation of Q c with a single k may differ significantly from the value computed with FDTD (see Supplementary information). In order to recover the absolute quality factor of the cavity with GME, Eq. (3) must be integrated over k in the Brillouin zone of the supercell.

Particles Swarm algorithm.
We have used the parallel version of the Particles Swarm (PS) algorithm available in the Global Optimization Toolbox of MATLAB. It starts from a random swarm, uniformly distributed, of 139 particles and requires 5800 iterations to converge. Every particle at each iteration corresponds to an independent evaluation of Q c with GME executed by one CPU core. The total number of cores used in the computation is 139. This results in 806200 GME simulations to find the maximum of the objective function Q c . It is worth mentioning that the number of iterations can be reduced by increasing the number of particles in the swarm, which may be advantageous if more cores are employed in the optimization.

FDTD numerical experiments.
Once the optimal parameters are found by the PS optimization, we simulate the final structure with a commercial FDTD solver 42 . The cavity mode is excited with an electric dipole oriented along the dominant mode polarization E y . This source was set to be a standard Gaussian pulse, centered at the cavity frequency with a narrow bandwidth of 0.88 THz. The total size of the PC structure was 86a × 75 √ 3a/2 with the L3 cavity localized at the center of the computational cell. PML boundary conditions and mirror symmetry with respect to the planes x = 0 , y = 0 and z = 0 were assumed. Finally, maximum mesh steps of a/40 and 0.5 √ 3a/40 were considered along x and y directions, respectively, while 40 cells per wavelength with a grading factor of 1.41421 were set along the z direction.  www.nature.com/scientificreports/