1. Introduction

Structural optimization has emerged as a new paradigm for the design of physical devices subject to wave and scattering physics. The typical electromagnetic wave problem consists of optimizing some quadratic field objective (e.g., the energy density at some location or the power radiated by a source) with respect to structural variations of the scattering medium. Because such design problems are non-convex, there are generally no guarantees of convergence to optimal solutions. Recently, however, a new means of probing wave equations makes possible the calculation of physical limits, or performance bounds, transforming the non-convex structural design problem into a convex and structure-agnostic field optimization problem with controlled relaxations of wave physics. Relying on the linearity of the wave equation, a bound on any quadratic wave objective is formulated as the solution of a quadratically constrained quadratic program (QCQP) by means of convex relaxations, e.g., Lagrangian duality or semi-definite programming [1–4]. This general framework has enjoyed great success in establishing new limits and scaling laws for diverse electromagnetic design objectives, from scattering cross sections to local density of states and focusing. Thus far, however, no such approach for deriving fundamental limits in the context of nonlinear photonics has been developed. Establishing limits beyond linear physics is of importance not only in guiding ongoing investigations of compact lasers, quantum controllers, and power limiters, but also in other disciplines including microfluidics and quantum mechanics.

In this paper, we propose a Lagrange duality framework [1] that is applicable to broad classes of nonlinear photonics objectives. Focusing on the problem of light incident on a $\chi ^{(2)}$ or Pockels medium, we show that under appropriate modifications, the introduction of auxiliary degrees of freedom transforms this nonlinear wave optimization problem into a QCQP susceptible to bound calculation through Lagrangian dual relaxations. Technically, proper selection of convex physical constraints is needed to ensure a feasible dual problem and non-trivial bounds. We use this approach to investigate upper limits to second harmonic generation (SHG), the maximum power that may be generated at twice the oscillating frequency of a harmonic source interacting with a Pockels medium. Apart from establishing limits for second harmonic (SH) power objectives, the bounds can predict the scalings and trends one may expect from optimal structures with respect to changes in material and device size. The method allows systematically increasing the number of physical constraints incorporated to achieve tighter limits. Finally, we propose several fabrication-ready doubly resonant slab cavities exhibiting among the highest nonlinear mode overlaps found in the literature. Extrapolating the observed performance gap in proof-of-concept two-dimensional settings to these realistic structures suggests little room for additional improvements beyond what is seen in topology-optimized cavities.

Second harmonic generation is an important phenomenon in nonlinear optics, with many practical applications including high frequency source generation [5,6], electro-optic modulators [7], as well as high resolution imaging [8] and spectroscopy [9,10], to name a few. While bulk nonlinearities are weak, the output power from a field incident on a nonlinear medium at the second harmonic can be resonantly enhanced through structuring: a cavity supporting tightly confined modes at both wavelengths leads to higher field intensities, or smaller cavity mode volumes $V=\frac {\int \varepsilon \vert{\textbf {E}}\vert ^2 \ d\textbf {r}}{\max {\left (\varepsilon \vert {\textbf {E}}\vert^2\right )}}$, and longer interaction times, or larger quality factors $Q$, along with a high nonlinear field-overlap factor $\beta = \frac {1}{4}\frac {\int \varepsilon _0\sum _{ijk}\chi ^{(2)}_{ijk} \left ( {E}_{1i}^{*}{E}_{1j}^{*}{E}_{2k}+ {E}_{1i}^{*}{E}_{2j}{E}_{1k}^{*}\right ) \ d\textbf {r} }{\left ( \int \varepsilon _0\varepsilon _1|\mathbf {E}_1|^2 \ d\textbf {r}\right ) \sqrt {\int \varepsilon _0\varepsilon _2|\mathbf {E}_2|^2\ d\textbf {r}} }$, where $\varepsilon _1$ and $\varepsilon _2$ denote the relative permittivity of the structure at the fundamental and second harmonic wavelengths, while $\textbf {E}_1$ and $\textbf {E}_2$ represent the corresponding mode profiles [11,12]; $\beta$ scales $\propto 1/\sqrt {V}$ and captures the degree of spatial confinement and mode coupling. The problem of designing resonators maximizing second harmonic generation is complicated because $Q$ and $\beta$ are often interconnected, with smaller cavities leading to larger radiative losses [13]. (Note that photonic crystals often exhibit a single narrow gap, precluding use of gap confinement at multiple wavelengths [14].) Previously, the relative ease in conceptualizing and fabricating large-etalon resonators supporting a plethora of tunable modes, proposed in the mid 1960s [15], led to devices with slower operation (larger $Q$) and smaller $\beta$. Over the past few decades, however, fabrication advances and the growing need to miniaturize devices (on a microchip) along with the promise of larger nonlinearities and greater operating speeds ushered a search for compact microcavities, from millimeter-[16] and micrometer-scale whispering gallery mode resonators [17–19] to, more recently, wavelength-scale structures designed via brute-force optimization [11]. In going toward smaller devices, the chief challenge of achieving simultaneous resonances at far-apart wavelengths poses a conceptual limitation [12,20], making quantitative analyses of expected trade-offs in temporal and spatial confinement very challenging. Recent inverse designs like the novel class of slab geometry proposed below, appear to exploit a range of physical mechanisms—from index guiding to Bragg scattering and slot confinement—to create spectrally far-apart, highly spatially confined resonances with moderately large bandwidths. This paper represents a top-down engineering perspective on the problem of enhancing frequency conversion in structured environments, setting fundamental wave constraints that complement our recent bottom-up investigations based on large-scale structural optimization [11].

We note that while there have been prior attempts at establishing bounds on Raman scattering of single molecules [21], those bounds exploit special simplifying characteristics of that problem—the nonlinear process occurs at a single fixed point in space—to factor the nonlinear photonics objective into a product of two linear objectives; such a factorization is not possible in more general settings such as problems involving extended nonlinear media. Independent linear photonics bounds were then evaluated using solely the constraint of passivity [22] and multiplied together to obtain a bound on the Raman scattering power; this neglects potential trade-offs in performance when simultaneously optimizing both linear objectives [23], and passivity alone does not adequately capture multiple scattering effects in photonic structures [13,24]. The framework presented below allows for a more complete account of wave physics as embodied by Maxwell’s equations [1] and consideration of spatially distributed nonlinearities subject to the typical restrictions of finite resonance strength and non-depletion of the pump; hence, it is generalizable to many other wave objectives and nonlinear processes.

2. Problem formulation

Consider a structure, as shown schematically in Fig. 1, represented by spatially varying local linear susceptibilities $\overline {\chi }_j(\textbf {r}) = \chi _j \overline {\rho }(\textbf {r})$ for $j=\{1,2\}$, at both fundamental $\omega _1$ and second-harmonic $\omega _2 = 2\omega _1$ frequencies, and by a local nonlinear susceptibility $\overline {\chi }^{(2)}(\textbf {r}) = \chi ^{(2)} \overline {\rho }(\textbf {r})$. Here, $\overline {\rho }(\textbf {r})$ denotes the distribution of material and is a scalar field that is allowed to vary continuously from $0$ to $1$ within some design volume $\Omega$, while $\chi _1$, $\chi _2$, and $\chi ^{(2)}$ are tensors representing the bulk susceptibilities of the medium comprising the structure. The design domain is illuminated by an incident electric field $\textbf {E}_1^{inc}(\textbf {r})$ created by a current source $\textbf {J}^{inc}(\textbf {r})$ oscillating at frequency $\omega _1$. For brevity, we refer to $\chi$ as the bulk material susceptibility and $\overline {\chi }$ as its spatial distribution, and omit spatial dependencies of both fields and distributions for the remainder of the manuscript. Generally, this nonlinear problem is challenging to solve, leading for instance to an infinite sequence of "up-converted" and "down-converted" waves at multiples of the incident frequency. However, the weak nature of nonlinearities allows for several simplifications. First, the scarcity and challenge of designing perfectly matched resonances at several, far-apart wavelengths means that one can focus exclusively on energy transfer between the fundamental ($\textbf {E}_1$) and second harmonic ($\textbf {E}_2$) fields. Second, under the undepleted pump approximation [25], one may neglect the down-converted field at the fundamental frequency. Finally, to lowest order in the nonlinearity, one may treat the induced nonlinear polarization at the second harmonic as a free source generating fields at $\omega _2$. These simplifying assumptions made it recently possible to frame the problem of maximizing radiative power at second harmonic as a sequence of coupled scattering problems [11], described by the following structural optimization program:

Here, $\varepsilon _0$ and $\mu _0$ are the vacuum permittivity and permeability respectively, the superscript $( ^*)$ denotes conjugated vector fields, and $\overline {\chi }^{(2)}\textbf {E}_1 \textbf {E}_2$ is the action of the rank-$3$ tensor field $\overline {\chi }^{(2)}(\textbf {r})$ on the vector fields $\textbf {E}_1(\textbf {r})$ and $\textbf {E}_2(\textbf {r})$. The $e^{-i\omega t}$ time convention is used to represent harmonic fields, with the vacuum Maxwell operator given by $\mathbb {M}_{l} = \left (-\mathbb {I}+ \frac {1}{\omega _l^2 \varepsilon _0}\boldsymbol{\nabla \times} \frac {1}{\mu _0}\boldsymbol{\nabla \times} \right ), \ l={1,2}$; where $\mathbb {I}$ is the identity operator. Finally, we define $\textbf {J}_1$ and $\textbf {J}_2$ as the induced currents resulting from nonlinear down-conversion and up-conversion, respectively.

 

Fig. 1. Schematic illustration of the second-harmonic generation problem under investigation: a device (representative design obtained through structural optimization) contained within some design region efficiently converts light incident at $\omega _1$ to the second harmonic frequency $\omega _2=2\omega _1$ by addressing the multi-prong challenge of supporting two spectrally far-apart and long-lived but tightly confined and spatially overlapping modes [11]. The proposed optimization framework establishes shape-independent limits on the largest achievable conversion efficiencies that may be realized for a given material and design region.

Download Full Size | PDF

In line with the arguments presented in Ref. [11], $\textbf {J}_2$ may be considered here as a free source at $\omega _2$ while $\textbf {J}_1$, included in (1) for technical reasons discussed further below, may be omitted under the undepleted pump approximation. Ignoring down-conversion, (1) amounts to the maximization of a quadratic objective satisfying a set of coupled partial differential equations, the solution of which maximizes the radiative SHG. Effectively, by Poynting’s theorem, $-\frac {1}{2} \operatorname {Re}\left [\int _\Omega \textbf {J}_2^* \cdot \textbf {E}_2 \ d\textbf {r}\right ]$ is the power extracted from the polarization field at the second harmonic and the absorptive loss at SH is given by $-\frac {1}{2}\omega _2 \varepsilon _0 \operatorname {Im}\left [\int _\Omega \textbf {E}_2^* \cdot \overline {\chi }_2^*\textbf {E}_2 \ d\textbf {r}\right ]$ while the difference between the two denotes the radiative SH power, $\Phi _2$. Even under the above stated simplifications, several difficulties must be addressed before one can arrive at a corresponding bound optimization problem, chief among them being the need to derive an analogous QCQP with dual feasibility. As detailed in [1] and further below, the transformation from a structural to a bound optimization problem requires shifting from structural $\overline {\rho }$ to polarization degrees of freedom, which leads to an objective that is evidently not quadratic. We solve this problem by promoting $\textbf {J}_2$ from a free source to an optimization degree of freedom, enforcing the nonlinear relation $\textbf {J}_2 = -i\omega _2\varepsilon _0 \overline {\chi }^{(2)}\textbf {E}_1 \textbf {E}_1$ instead as an auxiliary constraint. Additional relaxations and modifications are needed to ensure that this auxiliary constraint leads to tractable dual solutions. Finally, we show that including the down-conversion term in (1) leads to a power conservation constraint that contains negative semi-definite quadratic forms for all polarization fields, which makes finding feasible points for the dual straightforward. Further details follow below.

We begin by defining the following optimization degrees of freedom:

These volumetric identities are complex generalizations of the familiar Poynting’s theorem, establishing energy conservation across $\Omega _k$ and written here explicitly in terms of the bulk susceptibility ${\chi }_1$, and ${\chi }_2$, and vacuum electric Green’s function $\mathbb {G}_l = \mathbb {M}_l^{-1}$, $l=1,2$. $\Omega _k$ may be any sub-domain within the design domain, with the introduction of constraints over smaller sub-domains (down to the pixel level in a computational grid) forcing solutions to resolve wave physics at increasingly smaller length scales.

Apart from (5) and (6) which enforce energy conservation at each frequency separately, one may also impose similar constraints across fundamental and second harmonic frequencies that reflect the fact that the shape of any structure is the same for both frequencies [23,28]:

Notably, both sets of constraints are quadratic functions of the fields. Further, as noted above, the objective of (1) is also manifestly quadratic if $\textbf {P}^{(2)}$ is considered as an optimization degree of freedom and its nonlinear relation to the fields,

The relaxation in (13) relies on the assumption of a discretized basis representation for the vector fields; consequently (13) includes an explicit resolution factor $\mathcal {R}$ that will generally depend on the discretization scheme used to carry out the integrals: for example, given a finite difference grid with voxel volume $\Delta v$, $\mathcal {R}=1/\Delta v$. We refer to the Supplement 1 for discussions on normed relaxations and $\mathcal {R}$. Effectively, the result in (13) follows by taking the 2-norm ($\lVert {.}\rVert $) on both sides of (12) and taking advantage of the triangle inequality for norms, along with the Cauchy-Schwarz inequality, and the fact that the $l^4$ norm of a finite-dimensional vector is no greater than its $l^2$ norm [29,30]. We deduce the required norm constraint on the components of $\textbf {P}^{(2)}$; $\lVert {P_\alpha ^{(2)}}\rVert \leq \sum _{\beta \gamma } |\chi ^{(2)}_{\alpha \beta \gamma }| \lVert {E_{1_\beta } \odot E_{1_\gamma }}\rVert \leq \sum _{\beta \gamma } |\chi ^{(2)}_{\alpha \beta \gamma }| \left ( \sum _n |E_{1_\beta,n}|^4 \right )^{1/4} \left ( \sum _n |E_{1_\gamma,n}|^4 \right )^{1/4} \sqrt {\Delta v} \leq \sum _{\beta \gamma } |\chi ^{(2)}_{\alpha \beta \gamma }| \lVert {E_{1_\beta }}\rVert \lVert {E_{1_\gamma }}\rVert \frac {1}{\sqrt {\Delta v}}$, where the $\sum _n$ is a summation over all voxels in the sub-region $\Omega _j$. Using $E_{1_\alpha } = \chi _{1_{\alpha \alpha }}^{-1} P_{1_\alpha }$ we arrive at (13), where the $\int _{\Omega _j} P_{1,\beta }^* \cdot P_{1,\beta } \ d\textbf {r}$ terms can now be further bounded by the following auxiliary optimization

Finally, denoting $B_{j,\beta }$ as the Lagrangian dual bound on (14) yields the quadratic constraint

Note the complete freedom in choosing the set of sub-regions $\Omega _j$, $\Omega _k$, $\Omega _l$ associated with the different constraints.

For computational simplicity, we shall present results where $\chi _1$ is isotropic, and the only non-trivial component of $\chi ^{(2)}$ is $\chi ^{(2)}_{zzz}$ (this is the case, for example, in lithium niobate [17,25], or even some recently developed polaritonic metasurfaces [31]), and consider 2D problems with the fields polarized in the out-of-plane direction $z$. In this case (15) reduces to the following simpler form:

Putting together the aforementioned relaxations and modifications of the structural optimization problem (1), we arrive at the following transformed QCQP problem for maximizing radiative SHG:

Since any physically realizable set of polarization fields $\textbf {P}_1$, $\textbf {P}_2$, $\textbf {P}^{(2)}$ must satisfy all constraints in (17), its global optimum is an upper bound of the structural optimization problem (1). While it is still difficult to exactly determine the global optimum of (17), convex relaxations such as Lagrangian duality can be used to further bound the global optimum; see Supplement 1 for more details on formulating and evaluating the Lagrangian dual.

The convex nature of the Lagrangian dual problem implies that so long as one can find a feasible point for the dual, application of gradient based algorithms guarantees convergence to the global optimum. To find an initial feasible point of the Lagrangian dual problem, it is convenient to have a constraint or a set of constraints in the primal problem that collectively contain negative definite quadratic forms for all primal optimization degrees of freedom, $\textbf {P}_1$, $\textbf {P}_2$, and $\textbf {P}^{(2)}$. Assuming the nonlinear susceptibility $\chi ^{(2)}$ to be real and to have full permutation symmetry [25] and including the down-conversion term, albeit small, an elegant approach to construct such a constraint is by combining Maxwell’s equations at both the fundamental and the second harmonic frequencies and deriving a corresponding power-conservation constraint (see section 1 in Supplement 1 for details):

Here, the superscript $(^a)$ denotes the asymmetric part of an operator, $\mathbb {O}^a = (\mathbb {O}-\mathbb {O}^\dagger )/2i$, and $( ^\dagger )$ stands for conjugated transpose of an operator. The semi-definiteness of the asymmetric part of Green’s function [32] combined with the definiteness of the imaginary part of the linear responses of the material makes all quadratic forms negative definite when shifted to the left-hand side of (18). Physically, (18) is precisely the statement that power drawn from the incident field is equal to the sum of radiative and absorptive loss mechanisms at both frequencies. We note that near the dual optimum, this definite constraint, and thus, down conversion can actually be dropped without leaving the feasible domain or affecting the optimum value. While the central derivations leading to (17) pertain to second harmonic generation, similar procedures can be applied to formulate bounds on higher-order nonlinear processes. The reduction of higher than quadratic field dependencies to quadratic forms via the introduction of auxiliary degrees of freedom and relaxed nonlinear constraints should, in principle, be straightforward to carry out in other nonlinear settings, including Raman scattering [21,33] and Kerr media [34,35].

3. Numerical results

As proof of concept, we compute bounds for two representative problems, with results shown in Fig. 2: namely, maximizing the radiative second-harmonic power produced by either (a) a normally incident planewave or (b) a dipolar source in the vicinity of a structure contained within a square design region of size $d \times d$. For computational expedience, the calculations are carried out in 2D with sources and fields polarized in the out-of-plane direction (TM polarization). We assume an isotropic, lossy medium with a weak nonlinearity and normalize resulting power quantities by either (a) the net radiation of the dipole source in vacuum or (b) the power carried by the planewave incident on the largest of the design regions.

 

Fig. 2. Upper bounds on the maximum radiative second harmonic power possible from either an (a) incident planewave or (b) dipolar current source $0.1\lambda$ away from the center of the edge of a square design region of varying size $d \times d$. The non-dispersive medium under consideration has isotropic linear susceptibility $\chi _1=\chi _2 = 3+0.1i$ and Pockels susceptibility $\chi ^{(2)}_{zzz} = 398$ pm/V; the vacuum wavelength at $\omega _1$ is set at $\lambda = 2\pi c/\omega _1 = 1.5$ $\mu$m. The radiated power at the second harmonic frequency $\omega _2=2\omega _1$, denoted by $\Phi _2 = \frac {1}{2} \omega _2 \varepsilon _0 \operatorname {Im}\left [\int _\Omega \textbf {P}_2^* \cdot \mathbb {G}_2 \textbf {P}_2 \ d\textbf {r}\right ]$, is normalized by either the (a) power from flux $\Phi _{1} = \frac {1}{2}\operatorname {Re}\left [\int d\textbf {S}\cdot \left ( \textbf {E}_1^{inc} \times \textbf {H}_1^{inc*} \right )\right ] = 0.064$ mW/$\mu$m carried by the planewave incident on the largest design size considered ($d=2\lambda$), or (b) the net power emitted from the dipolar source in vacuum, $\Phi _1 = -\frac {1}{2} \operatorname {Re}\left [\int \textbf {J}^{inc*} \cdot \textbf {E}_1^{inc} \ d\textbf {r}\right ] = 0.197$ mW/$\mu$m. See Supplement 1 for more discussion on the valuse of incident field and $\chi ^{(2)}_{zzz}$ used. All calculations were restricted to 2D fields resulting from TM-polarized (out of the plane) sources. The solid lines only incorporate global constraints across the design region $\Omega$, consistent with the choice $\{\Omega _k\}=\{\Omega _j\}=\{\Omega \}$ in (17); filled circles demonstrate how the imposition of additional norm constraints (15) (larger sets of $\{\Omega _j\}$) can lead to tighter bounds, here only shown for a single representative $d=0.5\lambda$. Also shown are values of $\Phi _2/\Phi _1$ achieved by optimized structures (open circles). Representative dielectric profiles and the corresponding modes at both wavelengths are depicted in (c) and (d), with black/white denoting dielectric/vacuum regions and red/blue/white denoting positive/negative/zero values of the real part of the modal electric fields. The cavity mode-overlap factor $\beta$ is given in units of $\frac {\chi ^{(2)}}{\sqrt {\varepsilon _0\lambda ^2}}$.

Download Full Size | PDF

The magnitude of the incident fields and $\chi ^{(2)}$ were chosen to be within range of values seen in practical applications; we note that within the undepleted pump regime for a fixed structure, the SHG power scales quadratically with both the incident power at the fundamental frequency $\Phi _1$ and $|\chi ^{(2)}_{zzz}|$, so bounds for one set of $(\Phi _1, \chi ^{(2)}_{zzz})$ can be obtained through appropriate scaling from bounds for a different set of $(\Phi _1, \chi ^{(2)}_{zzz})$.

Bound calculations were performed by enforcing only the loosest of constraints, consistent with a choice of global domains $\{\Omega _k\}=\{\Omega _j\}=\{\Omega \}$ (solid line), with the addition of more localized constraints (filled circles) enforced by introducing smaller sub-domains leading to further tightening of the bounds (illustrated only for a single representative design region). See also Fig. S1 in Supplement 1 for a comparison with bounds obtained by imposing only the norm constraint in (13) for a global domain, $\{\Omega _j\}=\{\Omega \}$, and passivity as in [21]; these are found to be many orders of magnitude looser than the bounds in Fig. 2 even for materials with considerable loss. The bounds are compared to the performances of inverse designs (open circles) produced by solving the structural optimization problem of enhancing radiative SHG via gradient-based algorithms, including the method of moving asymptotes [36].

As shown in Fig. 2, bound calculations not only place relatively tight constraints on additional improvements that may be gained from further fine-tuning of structural optimizations (coming within one or at most two orders of magnitude of inverse designs) but also anticipate several noteworthy trends seen in optimized designs with increasing system size. Both bounds and the inverse designs seem to grow polynomially for subwavelength devices with resonances appearing in optimized structures around $d=\lambda$ (vacuum wavelength), with larger devices leading to higher quality factors $Q_1$, $Q_2$. The quality factors remain relatively small due to the high material loss, which explains the performance saturation with increasing system size (in more realistic settings, such as in 3D structured slabs, radiative losses are expected to play a similar role as material loss in this 2D example, leading to decreased performance or potentially saturation, a subject for future investigations). The mode profiles of optimized structures exhibit significant nonlinear overlaps in either scenario: in the particular case of a design region of size $2\lambda$ x $2\lambda$, $\beta \approx 0.14 \frac {\chi ^{(2)}}{\sqrt {\varepsilon _0 \lambda ^2}}$ for the planewave source and $\beta \approx 0.07 \frac {\chi ^{(2)}}{\sqrt {\varepsilon _0 \lambda ^2}}$ for the corresponding dipole source. Inspection of inverse structures also reveals a lack of consistent structural features in designs optimized for planewave compared to dipole sources: essentially, while there are many ways to "interfere" with a propagating wave comprising a narrow range of spatial wavelengths, dipolar fields include fast-decaying evanescent components (near fields) that require more sharply varying polarization profiles, restricting possible structures. While the bound calculations above focused on illustrative 2D examples, similar investigations may be carried out in more realistic settings. For example, Fig. 3 depicts optimized three-dimensional structures obtained by addressing the more challenging problem of maximizing second-harmonic power for on-chip light routed by a waveguide into a design domain. Taking advantage of the same simplifying assumptions of weak nonlinearity, undepleted pump, and bimodal approximations in (1), we consider maximization of the second harmonic Poynting flux, $\frac {1}{2}\operatorname {Re}\left [\int _{\mathrm {wvg}} \mathrm {d}\mathbf {S}\cdot \left [ \mathbf {E}_2 \times \mathbf {H}_2^{*} \right ]\right ]$, transmitted into a prescribed output waveguide. Schematics showing lateral cross sections of the geometries are depicted in Fig. 3, which consist of GaP on a thermal SiO$_2$ substrate as an integrated photonic platform to realize high conversion efficiency in robust, compact, and wide-bandwidth cavities [18]. A single port system is considered, with either 3(a) a single or 3(b) separate waveguides guiding light at the input and output wavelengths. For zincblende crystal phase of GaP, $\chi ^{(2)}$ is nonzero only for mutually perpendicular field components; consequently, the fundamental mode must be TE-like (E-field is mostly in-plane) and the second-harmonic mode TM-like (E-field is mostly out-of-plane). Figure 3 also depicts lateral cross-sections of the TE-like mode profiles at the fundamental wavelength $\lambda _1=1550$ nm and TM-like mode profiles at the second-harmonic wavelength $\lambda _2=775$ nm for both the structures. For ease of fabrication, the structures are designed to be invariant along the thickness dimension and to have minimal in-plane feature sizes $\gtrsim 60$ nm, leading to a vertical-to-lateral etch aspect ratio well within current experimental capabilities. Both cavities support resonances with moderate radiative quality factors $\sim 10^3$ and tightly confined modes, leading to large spatial overlaps, $\beta \approx 0.01 \frac {\chi ^{(2)}}{\sqrt {\varepsilon _0 \lambda ^3}}$, while ensuring near-critical coupling to nearby input/output waveguides.

 

Fig. 3. Lateral cross-sections of 3D optimized doubly resonant cavities—gallium phosphate (GaP) structures sitting on top of silica substrates—designed for efficient second harmonic generation under critical coupling with either (a) a single or (b) separate input/output waveguide(s). Black/white represent GaP/vacuum regions. Device (a) has lateral size $4.65\times 6.2$ $\mu$m$^2$ with a vertical thickness of 250 nm and an input waveguide of the same vertical thickness and width 457 nm; device (b) has lateral size $6.229\times 6.229$ $\mu$m$^2$ with a vertical thickness of 200 nm and input/output waveguides of the same vertical thickness and width 465 nm. Both devices support resonances at the fundamental $\lambda _1=1550$ nm and second-harmonic $\lambda _2=775$ nm wavelengths with radiative quality factors (a) $Q_1 \approx 3000$ and $Q_2 \approx 1000$, and (b) $Q_1 \approx 1100$ and $Q_2 \approx 760$, respectively. Also shown are $E_x$ and $E_z$ field profiles of the fundamental ($\omega _1$) and second harmonic ($\omega _2$) modes, respectively, with red/blue/white denoting positive/negative/zero values of the real part of the modal electric fields.

Download Full Size | PDF

4. Concluding remarks

We presented an approach to establish fundamental limits on nonlinear photonics objectives and applied it to investigate second harmonic generation in wavelength-scale structures. The limits not only quantify potential room for improvements but may be used to study scaling characteristics of optimal devices with respect to elemental design criteria, such as material choice and device size, without reference to particular geometric or physical enhancement mechanisms. A thorough and more focused exploration of limits on SHG in more realistic settings, such as the slab geometry of Fig. 3, and consideration of polaritonic or plasmonic metasurfaces in experimental setups [31,37], remains a challenging but important problem for future work. In addition to overcoming computational challenges needed to incorporate greater numbers of more finely resolved spatial constraints to obtain tighter predictions, there are several interesting directions worth pursuing. First, it should be possible to exploit symmetries to enforce that bounds respect fabrication constraints, including minimum-feature sizes and etchable geometries. Second, extension of the proposed framework to incorporate finite-bandwidth objectives should enable analysis of nonlinear space-bandwidth limitations, including trade-offs in achieving maximum spatial confinement and operating speeds at both wavelengths. Finally, detailed comparisons with established geometries like ring resonators [38–40] and inverse designs, including performance comparisons of dielectric, metallic, or even heterogeneous [41] media, should provide a more comprehensive view of the optimal design and performance landscape.