Cavity-enhanced second harmonic generation via nonlinear-overlap optimization

We describe an approach based on topology optimization that enables automatic discovery of wavelength-scale photonic structures for achieving high-efficiency second-harmonic generation (SHG). A key distinction from previous formulation and designs that seek to maximize Purcell factors at individual frequencies is that our method not only aims to achieve frequency matching (across an entire octave) and large radiative lifetimes, but also optimizes the equally important nonlinear--coupling figure of merit $\bar{\beta}$, involving a complicated spatial overlap-integral between modes. We apply this method to the particular problem of optimizing micropost and grating-slab cavities (one-dimensional multilayered structures) and demonstrate that a variety of material platforms can support modes with the requisite frequencies, large lifetimes $Q>10^4$, small modal volumes $\sim (\lambda/n)^3$, and extremely large $\bar{\beta} \gtrsim 10^{-2}$, leading to orders of magnitude enhancements in SHG efficiency compared to state of the art photonic designs. Such giant $\bar{\beta}$ alleviate the need for ultra-narrow linewidths and thus pave the way for wavelength-scale SHG devices with faster operating timescales and higher tolerance to fabrication imperfections.

Introduction.-Nonlinear optical processes mediated by second-order (χ (2) ) nonlinearities play a crucial role in many photonic applications, including ultra-short pulse shaping [1,2], spectroscopy [3], generation of novel frequencies and states of light [4][5][6], and quantum information processing [7][8][9]. Because nonlinearities are generally weak in bulk media, a well-known approach for lowering the power requirements of devices is to enhance nonlinear interactions by employing optical resonators that confine light for long times (dimensionless lifetimes Q) in small volumes V [10][11][12][13][14][15][16][17][18][19]. Microcavity resonators designed for on-chip, infrared applications offer some of the smallest confinement factors available, but their implementation in practical devices has been largely hampered by the difficult task of identifying wavelength-scale (V ∼ λ 3 ) structures supporting long-lived, resonant modes at widely separated wavelengths and satisfying rigid frequency-matching and mode-overlap constraints [15,20].
In this letter, we extend a recently proposed formulation for the scalable topology optimization of microcavities, where every pixel of the geometry is a degree of freedom, to the problem of designing wavelength-scale photonic structures for second harmonic generation (SHG). We apply this approach to obtain novel micropost and grating microcavity designs supporting strongly coupled fundamental and harmonic modes at infrared and visible wavelengths with relatively large lifetimes Q 1 , Q 2 > 10 4 . In contrast to recently proposed designs based on known, linear cavity structures hand-tailored to maximize the Purcell factors or mode volumes of individual resonances, e.g. ring resonators [17,[21][22][23] and nanobeam cavities [19,24], our designs ensure frequency matching and small confinement factors while also simultaneously maximizing the SHG enhancement factor Q 2 1 Q 2 |β| 2 to yield orders of magnitude improvements in the nonlinear couplingβ described by (3) and determined by a special overlap integral between the modes. These particular optimizations of multilayer stacks illustrate the benefits of our formalism in an approachable and experimentally feasible setting, laying the framework for future topology optimization of 2D/3D slab structures that are sure to yield even further improvements.
In situations involving all-resonant conversion, where confinement and long interaction times lead to strong nonlinearities and non-negligible down conversion [12,20], the maximum achievable occurs at a critical input power [20], where χ (2) eff is the effective nonlinear susceptibility of the medium [SM], is the dimensionless quality factor (ignoring material absorption) incorporating radiative decay 1 Q rad and coupling to an input/output channel 1 Q c . The dimensionless coupling coefficientβ is given by a complicated, spatial overlap-integral involving the fundamental and harmonic modes [SM], whereǭ(r) = 1 inside the nonlinear medium and zero elsewhere. Based on the above expressions one can define the following dimensionless figures of merit, where FOM 1 represents the efficiency per power, often quoted in the so-called undepleted regime of low-power conversion [33], and FOM 2 represents limits to power enhancement. Note that for a given radiative loss rate, FOM 1 is maximized when the modes are critically coupled, Q = Q rad 2 , with the absolute maximum occurring in the absence of radiative losses, Q rad → ∞, or equivalently, when FOM 2 is maximized. From either FOM, it is clear that apart from frequency matching and lifetime engineering, the design of optimal SHG cavities rests on achieving a large nonlinear couplingβ.
Optimal designs.- Table I characterizes the FOMs of some of our newly discovered microcavity designs, involving simple micropost and gratings structures of various χ (2) materials, including GaAs, AlGaAs and LiNbO 3 . The low-index material layers of the microposts consist of alumina (Al 2 O 3 ), while gratings are embedded in either silica or air (see supplement for detailed specifications). Note that in addition to their performance characteristics, these structures are also significantly different from those obtained by conventional methods in that traditional designs often involve rings [17,18], periodic structures or tapered defects [24], which tend to ignore or sacrificeβ in favor of increased lifetimes and for which it is also difficult to obtain widely separated modes [19]. Figure 1 illustrates one of the optimized structures-a doubly-resonant rectangular micropost cavity with alternating AlGaAs/Al 2 O 3 layers-along with spatial profiles of the fundamental and harmonic modes. It differs from conventional microposts in that it does not consist of periodic bi-layers yet it supports two localized modes at precisely λ 1 = 1.5 µm and λ 2 = λ 1 /2. In addition to having large Q rad 10 5 and small V ∼ (λ 1 /n) 3 , the structure exhibits an ultra-large nonlinear couplingβ ≈ 0.018 that is almost an order of magnitude larger than the best overlap found in the literature (see Fig. 2). From an experimental point of view, the micropost system is of particular interest because it can be realized by a combination of existing fabrication techniques such as molecular beam epitaxy, atomic layer deposition, selective oxidation and electron-beam lithography [25]. Additionally, the micropost cavity can be naturally integrated with quantum dots and quantum wells for cavity QED applications [26]. Similar to other wavelength-scale structures, the operational bandwidths of these structures are limited by radiative losses in the lateral direction [10,25,27], but their ultra-large overlap factors more than compensate for the increased bandwidth, which ultimately may prove beneficial in experiments subject to fabrication imperfections and for large-bandwidth applications [1,2,6,28].
To understand the mechanism of improvement inβ, it is instructive to consider the spatial profiles of interacting modes. Figure 1b plots the y-components of the electric fields in the xz-plane against the background structure. Sinceβ is a net total of positive and negative contributions coming from the local overlap factor E 2 1 E 2 in the presence of nonlinearity, not all local contributions are useful for SHG conversion. Most notably, one observes that the positions of negative anti-nodes of E 2 (light red regions) coincide with either the nodes of E 1 or alumina layers (where χ (2) = 0), minimizing negative contributions to the integrated overlap. In other words, improvements inβ do not arise purely due to tight modal confinement but also from the constructive overlap of the modes enabled by the strategic positioning of field extrema along the structure.
Based on the tabulated FOMs (Table I), the efficiencies and power requirements of realistic devices can be directly calculated. For example, assuming χ the AlGaAs/Al 2 O 3 micropost cavity ( Fig. 1) yields an efficiency of P 2,out P 2 1 = 2.7 × 10 4 /W in the undepleted regime when the modes are critically coupled, Q = Q rad 2 . For larger operational bandwidths, e.g. Q 1 = 5000 and Q 2 = 1000, we find that P 2,out P 2 1 = 16/W. When the system is in the depleted regime and critically coupled, we find that a maximum efficiency of 25% can be achieved at P crit 1 ≈ 0.15 mW whereas assuming smaller Q 1 = 5000 and Q 2 = 1000, a maximum efficiency of 96% can be achieved at P crit Comparison against previous designs. -Table II  umes, improvements inβ are still hardly sufficient for achieving high efficiencies at low powers.
Ultra-compact nanophotonic resonators such as the recently proposed nanorings [18], 2D photonic crystal defects [16], and nanobeam cavities [19], possess even smaller mode volumes but prove challenging for design due to the difficulty of finding well-confined modes at both the fundamental and second harmonic frequencies [16]. Even when two such resonances can be found by fine-tuning a limited set of geometric parameters [18,19], the frequency-matching constraint invariably leads to sub-optimal spatial overlaps which severely limits the maximal achievableβ.
Comparing Tables I and II, one observes that for a comparable Q, the topology-optimized structures perform significantly better in both FOM 1 and FOM 2 than any conventional geometry, with the exception of the LN gratings, whose low Q rad lead to slightly lower FOM 2 . Generally, the optimized microposts and gratings perform better by virtue of a large and robustβ which, notably, is significantly larger than that of existing designs. Here, we have not included in our comparison those structures which achieve non-negligible SHG by special poling techniques and/or quasiphase matching methods [32][33][34], though their performance is still sub-optimal compared to the topology-optimized designs. Such methods are highly material-dependent and are thus not readily applicable to other material platforms; instead, ours is a purely geometrical topology optimization technique applicable to any material system.
Optimization formulation.-Optimization techniques have been regularly employed by the photonic device community, primarily for fine-tuning the characteristics of a pre-determined geometry; the majority of these techniques involve probabilistic Monte-Carlo algorithms such as particle swarms, simulated annealing and genetic algorithms [35][36][37]. While some of these gradientfree methods have been used to uncover a few unexpected results out of a limited number of degrees of freedom (DOF) [38], gradient-based topology optimization methods efficiently handle a far larger design space, typically considering every pixel or voxel as a DOF in an extensive 2D or 3D computational domain, giving rise to novel topologies and geometries that might have been difficult to conceive from conventional intuition alone. The early applications of topology optimization were primarily focused on mechanical problems [39] and only recently have they been ex-   A key realization in [41] is that instead of maximizing the LDOS at a single discrete frequency ω, a betterposed problem is that of maximizing the frequency-averaged f in the vicinity of ω, denoted by where W is a weight function defined over some specified bandwidth Γ. Using contour integration techniques, the frequency integral can be conveniently replaced by a single evaluation of f at a complex frequency ω + iΓ [41]. For a fixed Γ, the frequency average effectively shifts the algorithm in favor of minimizing V over maximizing Q; the latter can be enhanced over the course of the optimization by gradually winding down the averaging bandwidth Γ [41]. A major merit of the frequency-averaged LDOS formulation is that it features a mathematically well-posed objective as opposed to a direct maximization of the cavity Purcell factor Q V , allowing for rapid convergence of the optimization algorithm into an extremal solution.
As suggested in [41], a simple extension of the optimization problem from single-to multimode cavities maximizes the minimum of a collection of LDOS at different frequencies. Applying such an approach to the problem of SHG, the optimization objective becomes: maxǭ α min LDOS(ω 1 ), LDOS(2ω 1 ) which would require solving two separate scattering problems, M 1 E 1 = J 1 and M 2 E 2 = J 2 , for the two distinct point sources J 1 , J 2 at ω 1 and ω 2 = 2ω 1 respectively. However, as discussed before, rather than maximizing the Purcell factor at individual resonances, the key to realizing optimal SHG is to maximize the overlap integralβ between E 1 and E 2 , described by (3). Here, we suggest an elegant way to incorporateβ by coupling the two scattering problems. In particular, we consider not a point dipole but an extended source J 2 ∼ E 2 1 at ω 2 and optimize a single combined radiated power f = − Re dr J * 2 · E 2 instead of two otherwise unrelated LDOS. The advantage of this approach is that f yields precisely theβ parameter along with any resonant enhancement factors (∼ Q/V ) in E 1 and E 2 . Intuitively, J 2 can be thought of as a nonlinear polarization current induced by E 1 in the presence of the second order susceptibility tensor χ χ χ (2) , and in particular is given by J 2i =ǭ(r) jk χ (2) ijk E 1j E 1k where the indices i, j, k run over the Cartesian coordinates. In general, χ (2) ijk mixes polarizations and hence f is a sum of different contributions from various polarization-combinations. In what follows and for simplicity, we focus on the simplest case in which E 1 and E 2 have the same polarization, corresponding to a diagonal χ (2) tensor determined by a scalar χ (2) eff . Such an arrangement can be obtained for example by proper alignment of the crystal orientation axes [SM]. With this simplification, the generalization of the linear topology-optimization problem to the case of SHG becomes: and where ǫ d denotes the dielectric contrast of the nonlinear medium and ǫ m is that of the surrounding linear medium. Note thatǭ α is allowed to vary continuously between 0 and 1 whereas the intermediate values can be penalized by so-called threshold projection filters [45]. The scattering framework makes it straightforward to calculate the derivatives of f (and possible functional constraints) with respective toǭ α via the adjoint variable method [39][40][41]. The optimization problem can then be solved by any of the many powerful algorithms for convex, conservative, separable approximations, such as the well-known method of moving asymptotes [46].
For computational convenience, the optimization is carried out using a 2D computational cell (in the xz-plane), though the resulting optimized structures are given a finite transverse extension h y (along the y direction) to make realistic 3D devices (see Fig. 3). In principle, the wider the transverse dimension, the better the cavity quality factors since they are closer to their 2D limit which only consists of radiation loss in the z direction; however, as h y increases,β decreases due to increasing mode volumes. In practice, we chose h y on the order of a few vacuum wavelengths so as not to greatly compromise either Q orβ. We then analyze the 3D structures via rigorous FDTD simulations to determine the resonant lifetimes and modal overlaps. By virtue of our optimization scheme, we invariably find that frequency matching is satisfied to within the mode linewidths. We note that our optimization method seeks to maximize the intrinsic geometric parameters such as Q rad andβ of an un-loaded cavity whereas the loaded cavity lifetime Q depends on the choice of coupling mechanism, e.g. free-space, fiber, or waveguide coupling, and is therefore an external parameter that can be considered independently of the optimization. When evaluating the performance characteristics such as FOM 1 , we assume total operational lifetimes Q 1 = 5000, Q 2 = 1000. In the optimized structures, it is interesting to note the appearance of deeply sub-wavelength features ∼ 1 − 5% of λ 1 n , creating a kind of metamaterial in the optimization direction; these arise during the optimization process regardless of starting conditions due to the low-dimensionality of the problem. We find that these features are not easily removable as their absence greatly perturbs the quality factors and frequency matching.
Concluding remarks.-We have presented a formulation that allows for large-scale optimization of SHG. Applied to simple micropost and grating structures, our approach yields new classes of microcavities with stronger performance metrics over existing designs. One potentially challenging aspect for fabrication in the case of gratings is the presence of deeply sub-wavelength features, which would require difficult high-aspect-ratio etching or growth techniques. This is not an issue for the micropost cavities since each material layer can be grown/deposited to a nearly arbitrary thickness [25,26]. Another caveat about wavelength-scale cavities is that they are sensitive to structural perturbations near the cavity center, where most of the field resides. In our optimized structures, the figures of merit are robust to within ∼ ±20 nm variations (approximately one computational pixel). One possible way to constrain the optimization to ensure some minimum spatial feature and robustness is to exploit so-called regularization filters and worst-case optimization techniques [45], which we will consider in future work.
Our results provide just a glimpse of the kinds of designs that could be realized in structures with 2D and 3D variations, where we expect even better performance metrics due to the larger design space. In fact, preliminary application of our formulation to 2D systems reveals overlap factors and lifetimes at least one order of magnitude larger than those attained here. Apart from SHG optimization, our approach can be generalized to consider other nonlinear processes, even those involving more than two frequencies [SM]. Preliminary investigations reveal orders-of-magnitude improvements in the efficiency of third harmonic and sum-frequency generation processes. These findings, together with higher-dimensional applications, will be presented in future work. * zlin@seas.harvard.edu

COUPLED-MODE THEORY FOR SECOND HARMONIC GENERATION
The temporal coupled mode equations describing second harmonic generation (SHG) in a doubly-resonant cavity coupled to a channel are [1]: such that |a k | 2 is the modal energy in the cavity and |s k± | 2 is the input/output power in the waveguide, and where Q k and Q rad k denote the total and radiative quality factors corresponding to mode k. The nonlinear coupling coefficient β 1 , obtained from perturbation theory [1], is given by: with β 2 = β * 1 /2 far off from material resonances where Kleinman symmetry is valid [2], as required by conservation of energy [1]. In general, the overlap integral in the numerator is a sum of products between different E-field polarizations weighted by off-diagonal components of the nonlinear χ (2) tensor. For simplicity, however, in the main text we only consider the simple case of diagonal χ (2) involving same-polarization interactions described by an effective χ (2) eff , resulting from an appropriate orientation of the crystal axes of the nonlinear material. All of these considerations suggest a simple dimensionless normalization of β, given by: such that β 2 = 4βχ (2) eff / ǫ 0 λ 3 1 . As defined in the text,ǭ(r) = 1 for nonlinear dielectric and ǫ(r) = 0 for the surrounding linear medium.
Most SHG experiments operate in the small-signal regime of small input powers, leading to negligible down-conversion and pump depletion. Ignoring the down-conversion or β 1 term in Eq. 1, one obtain the following simple expression for the second harmonic output power: In the limit of large up-conversion and non-negligible down-conversion, solution of the coupledmode equations yields the maximum efficiency (defined as η = P out 2 /P in 1 ) and corresponding critical power [1]:

QUENCY CONVERSION PROCESS
Nonlinear frequency conversion processes can be viewed as frequency mixing schemes in which two or more constituent photons at a set of frequencies {ω n } interact to produce an output photon at frequency Ω such that Ω = n c n ω n , where the photon number coefficients {c n } can be either negative or positive, depending on whether the corresponding photons are created or destroyed in the process. In general, because the optical nonlinear response of materials is a tensor and hence the frequency conversion process mixes different polarizations [2]. However, for notational simplicity, we will describe our optimization problem only for a single component of the susceptibility tensor. If one wishes to consider the full tensorial properties, one can easily add extra optimization terms (weighted by the tensor components) by following the same approach The formulation is now given by: maxǭ f (ǭ; ω n ) = −Re J(Ω) * · E(Ω) dr , M(ǭ, ω n )E n = iω n J n , J n =ê nν δ(r − r ′ ), In practice, we maximize a frequency-averaged version of the power output f (ω) rather than f (ω) itself since the latter has poor convergence [3], i.e., we maximize f = dω ′ W(ω ′ ; ω, Γ)f (ω ′ ), where we simply choose the weighting function W to be a simple lorentzian with the desired resonance ω and a certain linewidth Γ, W(ω ′ ) = Γ/π (ω ′ − ω) 2 + Γ 2 (10) (Note that this linewidth is not necessarily the same as the intrinsic radiative linewidth of the cavity in that Γ is a computational artifice introduced to aid rapid convergence [3].) In Ref. [3], it is shown by means of contour integration that this averaging is equivalent to evaluating f at a complex frequency f (ω + iΓ), also equivalent to adding a uniform loss iΓ/ω to µ(r) and ǫ(r). In our implementation of the optimization process, we typically begin with a large Γ which affords rapid convergence to a stable geometry in a few hundred iterations. Γ is then decreased by an order of magnitude every time the optimization converges until Γ ∼ 10 −5 at which point the structure settles to within a linewidth Γ of the desired frequencies (perfect frequency matching).
Application of this formulation to the problem of second harmonic generation is straightforward and described in the main text, in which case Ω = ω 2 = 2ω 1 and J(Ω) = J 2 = ǫ(r) (ê 1j · E 1 ) 2ê 2i . Note that for the structures described in the text, we choseê 1j =ê 2i =ê y . In addition to the problem statement of Eq. 9, the optimization algorithm also benefits from gradient information of the objective function, which exploits the adjoint variable method [3][4][5]. Here, we simply quote the result for the gradient of our SHG objective function (dropping the polarization