1. INTRODUCTION
Nonlinear frequency conversion plays a crucial role in many photonic applications, including ultrashort pulse shaping ^{[1,2]}, spectroscopy ^{[3]}, generation of novel optical states ^{[4–6]}, and quantum information processing ^{[7–9]}. Although frequency conversion has been studied exhaustively in bulky optical systems, including large ring resonators ^{[10]} and etalon cavities ^{[11]}, it remains largely unstudied in micro and nanoscale structures where light can be confined to lengthscales of the order of or even smaller than its wavelength. By confining light over long a time and to small volumes, such highly compact devices greatly enhance light–matter interactions, enabling similar as well as new ^{[12]} functionalities compared to those available in bulky systems but at much lower power levels. Several proposals have been put forward based on the premise of observing enhanced nonlinear effects in structures capable of supporting multiple resonances at faraway frequencies ^{[13–21]}, among which are microring resonators ^{[22,23]} and photonic crystal (PhC) cavities ^{[24,25]}. However, to date, these conventional designs fall short of simultaneously meeting the many design challenges associated with resonant frequency conversion, chief among them being the need to support multiple modes with highly concentrated fields, exactly matched resonant frequencies, and strong mode overlaps ^{[26]}. Recently, we proposed to leverage powerful, largescale optimization techniques (commonly known as inverse design) to allow computeraided photonic designs that can address all of these challenges.
Our recently demonstrated optimization framework allows automatic discovery of novel cavities that support tightly localized modes at several desired wavelengths and exhibit large nonlinear mode overlaps. As a proofofconcept, we proposed doubly resonant structures, including multilayered, aperiodic micropost cavities and multitrack ring resonators, capable of realizing secondharmonic generation efficiencies exceeding ${10}^{4}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{W}}^{1}$ ^{[27,28]}. In this paper, we extend and apply this optimization approach to design extended structures, including microstructured optical fibers and PhC threedimensional metasurfaces, as shown in Fig. 1, for achieving highefficiency (second and thirdharmonic) frequency conversion. Harmonic generation, which underlies numerous applications in science, including coherent light sources ^{[29]}, optical imaging and microscopy ^{[30,31]}, and entangledphoton generation ^{[32]}, is now feasible at lower power requirements thanks to the availability of highly nonlinear ${\chi}^{(2)}$ and ${\chi}^{(3)}$ materials such as III–V semiconductor compounds ^{[33,34]} and novel types of chalcogenide glasses ^{[35]}. In combination with advances in materials synthesis, emerging fabrication technologies have also enabled demonstrations of sophisticated microstructured fibers ^{[36]} and metasurfaces ^{[37–44]}, paving the way for experimental realization of inversedesigned structures of increased geometric and fabrication complexity, which offer ordersofmagnitude enhancements in conversion efficiencies and the potential for augmented functionalities.
Figure 1.Schematic illustration of thirdharmonic generation and secondharmonic generation processes in inversedesigned microstructured fibers and metasurfaces, respectively.
Given a material system of intrinsic ${\chi}^{(2)}$ or ${\chi}^{(3)}$ nonlinear coefficient, the efficiency of any given frequencyconversion process in a resonant geometry will be determined by a few modal parameters. The possibility of confining light within small mode volumes over a long time or distance leads to significant gains in efficiency (i.e., lower power requirements), stemming from the higher intensity and cascadability of nonlinear interactions (compensating for the otherwise small bulk nonlinearities). In particular, the efficiency of such resonant processes depends on the product of mode lifetimes and a nonlinear coefficient $\beta $, given by Eqs. (6) and (8) below, which generalizes the familiar concept of quasiphasematching to situations that include wavelengthscale resonators ^{[26]}. For propagating modes, leaky or guided, the existence of a propagation phase further complicates this figure of merit, with optimal designs requiring: (i) phasematching and frequencymatching conditions, (ii) large nonlinear mode overlaps $\beta $, and (iii) large dimensionless lifetimes $Q$ (low material absorption and/or radiative losses in the case of leaky modes). The main design challenge is the difficult task of forming a doubly resonant cavity with farapart modes that simultaneously exhibit long lifetimes and large $\beta $, along with phase and frequency matching. To date, the majority of prior works on frequency conversion in fibers ^{[45–47]} and metasurfaces ^{[38–40,42,48–51]} have focused on only one of these aspects (usually phase matching) while ignoring the others. The geometries discovered by our optimization framework, in contrast, address the above criteria, revealing complex fibers and metasurfaces supporting TE or TM modes with guaranteed phase and frequency matching, long lifetimes $Q$, and enhanced overlap factors $\beta $ at any desired propagation wavevector, and resulting in ordersofmagnitude enhancements in conversion efficiencies.
Sign up for Photonics Research TOC Get the latest issue of Advanced Photonics delivered right to you！Sign up now
2. OVERVIEW OF OPTIMIZATION
The possibility of finetuning spatial features of photonic devices to realize functionalities not currently achievable by conventional optical design methodologies based on index guiding and bandgap confinement (which work exceedingly well but are otherwise limited for narrowband applications) has been a major drive behind the past several decades of interest in the topic of photonic optimization ^{[52,53]}. Among these techniques are probabilistic Monte Carlo algorithms, e.g., particle swarms, simulated annealing, and genetic algorithms ^{[54–56]}. Though sufficient for the majority of narrowband (singlemode) applications, many of these gradientfree methods are limited to typically small sets of design parameters ^{[57]} that often prove inadequate for handling wideband (multimode) problems. On the other hand, gradientbased inversedesign techniques are capable of efficiently exploring a much larger design space by making use of analytical derivative information of the specified objective and constraint functions ^{[58]}, demonstrated to be feasible for as many as ${10}^{9}$ design variables ^{[59]}. Recently, the development of versatile mathematical programming methods and the rapid growth in computational power have enabled concurrent progress in photonic inverse design, allowing theoretical (and more recently, experimental) demonstrations of complex topologies and unintuitive geometries with unprecedented functionalities that would be arguably difficult to realize through conventional intuition alone. However, to date, most applications of inverse design in photonics are confined to linear devices such as mode converters, waveguide bends, and beam splitters ^{[57,58,60–65]}. We believe that this paper along with our recent works ^{[27,28]} provides a glimpse of the potential of photonic optimization in nonlinear optics.
A typical optimization problem seeks to maximize or minimize an objective function $f$, subject to certain constraints $g$, over a set of free variables or degrees of freedom (DOFs) ^{[66]}. Generally, one can classify photonic inverse design into two different classes of optimization strategies, based primarily on the nature or choice of DOF ^{[67]}. Given a computational domain or grid, the choice of a finitedimensional parameter space not only determines the degree of complexity but also the convergence and feasibility of the solutions. One possibility is to exploit each DOF in the computational domain as an optimization parameter, known as topology optimization (TO), in which case one typically (though not always) chooses the dielectric permittivity of each pixel $\epsilon (\mathbf{r})$ as a DOF (known as a continuous relaxation parameter ^{[68]}). Another possibility, known as shape optimization, is to expand the optimization parameter space in a finite set of shapes (independent of the computational discretization), which may be freeform contours represented by socalled level sets ^{[69]} (the levelset method) or basic geometric entities with simpler parametrizations (e.g., polytopes) ^{[70]}. In the levelset method, the zeros of a levelset “function” $\mathrm{\varphi}(\mathbf{r})$ define the boundaries of “binary shapes”; the optimization then proceeds via a levelset partial differential equation characterized by a velocity field, which is, in turn, constructed from derivative information ^{[69]}. A much simpler variant (which we follow) is to choose a fixed but sufficient number of basic binary shapes whose parameters can be made to evolve by an optimization algorithm. Essentially, for such a parametrization, the mathematical representations of the shapes must yield continuous (analytic) derivatives, which is not feasible a priori due to the finite computational discretization and can instead be enforced by the use of a “smoothing kernel” (described below).
A generic TO formulation is written down as where the DOFs are the normalized dielectric permittivities ${\overline{\epsilon}}_{\alpha}\in [0,1]$ assigned to each pixel or voxel (indexed $\alpha $) in a specified volume ^{[58,60]}. The subscript $\alpha $ denotes appropriate spatial discretization $\mathbf{r}\to {(i,j,k)}_{\alpha}\mathrm{\Delta}$ with respect to Cartesian or curvilinear coordinates. Depending on the choice of background (bg) and structural materials, ${\overline{\epsilon}}_{\alpha}$ is mapped onto a positiondependent dielectric constant via ${\epsilon}_{\alpha}=(\epsilon {\epsilon}_{\text{bg}}{\overline{\epsilon}}_{\alpha}+{\epsilon}_{\text{bg}})$. The binarity of the optimized structure is enforced by penalizing the intermediate values $\overline{\epsilon}\in (0,1)$ or utilizing a variety of filter and regularization methods ^{[58]}. Starting from a random initial guess, the technique discovers complex structures automatically with the aid of powerful gradientbased algorithms such as the method of moving asymptotes (MMA) ^{[71]}. For an electromagnetic problem, $f$ and $g$ are typically functions of the electric $\mathbf{E}$ or magnetic $\mathbf{H}$ fields integrated over some region, which are in turn solutions of Maxwell’s equations under some incident current or field. In what follows, we exploit direct solution of frequencydomain Maxwell’s equationsdescribing the steadystate field $\mathbf{E}(\mathbf{r};\omega )$ in response to incident currents $\mathbf{J}(\mathbf{r},\omega )$ at frequency $\omega $. While solution of Eq. (4) is straightforward and commonplace, the key to making optimization problems tractable is to obtain a fastconverging and computationally efficient adjoint formulation of the problem. Within the scope of TO, this requires efficient calculations of the derivatives $\frac{\partial f}{\partial {\overline{\epsilon}}_{\alpha}},\frac{\partial g}{\partial {\overline{\epsilon}}_{\alpha}}$ at every pixel $\alpha $, which we perform by exploiting the adjointvariable method (AVM) ^{[58]}.
While the TO technique is quite efficient in handling the enormity of an unconstrained design space, it often leads to geometries with irregular features that are difficult to fabricate. An alternative approach that is in principle more conducive to fabrication constraints is to exploit shape optimization. In this work, we primarily focus on a simple implementation of the latter that employs a small and, hence, limited set of elementary geometric shapes, e.g., ellipses ^{[72]} and polytopes, parametrized by a few DOFs. In particular, we express the dielectric profile of the computational domain as a sum of basic shape functions with permittivities, ${\overline{\epsilon}}_{\alpha}=\sum _{\beta}{H}_{\beta}({\mathbf{r}}_{\alpha};\{{p}_{\beta}\})$, described by shape functions ${H}_{\beta}$ and a finite set of geometric parameters $\{{p}_{\beta}\}$, where $\beta $ denotes the shape index. Here, to deal with potential overlap of two or more shapes, we implement a filter function that enforces the same maximumpermittivity constraint $\overline{\epsilon}\le 1$ described above. The derivatives of a given objective function $f$ (and associated constraints) can then be obtained via the chain rule $\frac{\partial f}{\partial {p}_{i}}=\frac{\partial f}{\partial {\overline{\epsilon}}_{\alpha}}\frac{\partial {\overline{\epsilon}}_{\alpha}}{\partial {p}_{i}}$, where the smoothness of the derivatives is guaranteed by insisting that the shape functions $H$ be continuously differentiable functions. Below, we choose nonpiecewiseconstant ellipsoidal shapes with exponentially varying dielectric profiles near the boundaries, the smoothness of which is determined by a few simple parameters that can, at various points along the optimization, be slowly adjusted to realize fully binary structures upon convergence. Such a “relaxation” process ^{[70]} is analogous to the application of a binary filter in the objective function ^{[58]}.
Any nonlinear frequency conversion process can be viewed as a frequency mixing scheme in which two or more constituent photons at a set of frequencies $\{{\omega}_{n}\}$ interact to produce an output photon at frequency $\mathrm{\Omega}=\sum _{n}{c}_{n}{\omega}_{n}$, where $\{{c}_{n}\}$ can be either negative or positive, depending on whether the corresponding photons are created or destroyed in the process ^{[73]}. Given an appropriate nonlinear tensor component ${\chi}_{ijk\dots}$, with $i,j,k,\dots \in \{x,y,z\}$, mediating an interaction between the field components ${E}_{i}(\mathrm{\Omega})$ and ${E}_{1j},{E}_{2k},\dots $, we begin with a collection of point dipole currents, each at the constituent frequency ${\omega}_{n},n\in \{1,2,\dots \}$, such that ${\mathbf{J}}_{n}={\widehat{\mathbf{e}}}_{n\nu}\delta (\mathbf{r}{\mathbf{r}}^{\prime})$, where ${\widehat{\mathbf{e}}}_{n\nu}\in \{{\widehat{\mathbf{e}}}_{1j},{\widehat{\mathbf{e}}}_{2k},\dots \}$ is a polarization vector chosen so as to excite the desired electricfield polarization components ($\nu $) of the corresponding mode at an appropriate position ${\mathbf{r}}^{\prime}$. Given the choice of incident currents ${\mathbf{J}}_{n}$, we solve Maxwell’s equations to obtain the corresponding constituent electricfield response ${\mathbf{E}}_{n}$, from which one can construct a nonlinear polarization current $\mathbf{J}(\mathrm{\Omega})=\overline{\epsilon}(\mathbf{r})\prod _{n}{E}_{n\nu}^{{c}_{n}(*)}{\widehat{\mathbf{e}}}_{i}$, where ${E}_{n\nu}={\mathbf{E}}_{n}\xb7{\widehat{\mathbf{e}}}_{n\nu}$ and $\mathbf{J}(\mathrm{\Omega})$ can be generally polarized (${\widehat{\mathbf{e}}}_{i}$) in a (chosen) direction that differs from the constituent polarizations ${\widehat{\mathbf{e}}}_{n\nu}$. Here, $(*)$ denotes complex conjugation for negative ${c}_{n}$ and no conjugation otherwise. Finally, maximizing the radiated power, $\mathrm{Re}[\int \mathbf{J}{(\mathrm{\Omega})}^{*}\xb7\mathbf{E}(\mathrm{\Omega})\mathrm{d}\mathbf{r}]$, due to $\mathbf{J}(\mathrm{\Omega})$, one is immediately led to the following nonlinear optimization problem: where $\overline{\epsilon}$ is given by either the topology or shape parameterizations described above. Writing down the objective function in terms of the nonlinear polarization currents, it follows that solution of Eq. (5), obtained by employing any mathematical programming technique that makes use of gradient information, e.g., the AVM ^{[58]}, maximizes the nonlinear coefficient (mode overlap) associated with the aforementioned nonlinear optical process. The above framework can be easily extended to consider propagating modes once we take into account the appropriate Bloch boundary conditions that may arise from any desired wave vectors imposed at the requisite frequencies ^{[74]}. In the case of optical fibers or PhC metasurfaces (or, more generally, any waveguiding system), such an extension naturally guarantees perfect phase and frequency matching of the relevant modes in the optimized structure.
3. THIRDHARMONIC GENERATION IN FIBERS
Conventional microstructured fibers (e.g., Bragg and holey fibers) are typically designed based on intuitive principles like slow light ^{[47]}, index guiding, and bandgap confinement ^{[52]}, and thus often consist of periodic cross sections comprising simple shapes ^{[75,76]}. Below, we apply the aforementioned optimization techniques to propose much more complicated heterostructure fibers designed to enhance thirdharmonic generation at any desired wavelength. To achieve large thirdharmonic generation efficiencies, the fiber must support two copropagating modes of frequencies ${\omega}_{1}$ and ${\omega}_{3}=3{\omega}_{1}$ and wavenumbers that satisfy the phasematching condition ${k}_{3}=3{k}_{1}$. Furthermore, the system must exhibit low radiative/dissipative losses or, alternatively, attenuation lengths that are much longer than the corresponding interaction lengths $L$, defined as the propagation length at which 50% of the fundamental mode is upconverted. In the smallinput signal regime, the converted thirdharmonic output power ${P}_{3}\propto {P}_{1}^{2}$ and the interaction length $L=\frac{16}{3{k}_{1}{Z}_{0}{\beta}_{3}{P}_{1}}$ depend on the incident power ${P}_{1}$, vacuum impedance ${Z}_{0}$, and nonlinear overlap factor ^{[77]} which involves a complicated spatial overlap of the two modes over the crosssectional surface $S$ of the fiber. Note that the attenuation coefficient $\gamma \equiv \omega /2{v}_{g}Q$ of each mode (the inverse of their respective attenuation length) is proportional to their lifetime $Q$ and group velocity ${v}_{g}$.
We focus on fibers comprising chalcogenide/polyethersulfone (PES) composites of permittivities ${\epsilon}_{{\mathrm{As}}_{2}{\mathrm{Se}}_{3}}=5.8125$ and ${\epsilon}_{\mathrm{PES}}=2.4025$ at telecom wavelengths. Although our technique can be readily applied to design the requisite properties at any given wavenumber $k$ and for any desired polarization, we specifically focus on designs for operation at wavenumbers in the range $0.1(2\pi /\lambda )<{k}_{\mathrm{opt}}<2.3(2\pi /\lambda )$, with $\lambda $ denoting the corresponding vacuum wavelength and ${k}_{\mathrm{opt}}$ the optimized wavenumber. We consider both leaky and guided modes above and below the PES lightline $\omega =ck/\sqrt{{\epsilon}_{\mathrm{PES}}}$, respectively, along with different choices or transverse electric ${\mathrm{TE}}_{01}$ and transverse magnetic ${\mathrm{TM}}_{01}$ polarizations. ${\mathrm{TE}}_{01}$ modes are those polarized along the plane of the fiber and consist primarily of circulating ${E}_{x}$ and ${E}_{y}$ electric fields ^{[78]}, while ${\mathrm{TM}}_{01}$ modes have electric fields ${E}_{z}$ polarized mainly along the propagation direction $z$.
The top insets in Fig. 2 show an inversedesigned fiber cross section that supports phasedmatched ${\mathrm{TM}}_{01}$ fundamental and thirdharmonic modes (with profiles superimposed on the insets) at ${k}_{\mathrm{opt}}={k}_{1}=1.4(2\pi /\lambda )$. To ensure that the optimization algorithm selectively finds ${\mathrm{TM}}_{01}$ modes, we employ a magnetic current ${\mathbf{J}}_{1}\sim \nabla \times \delta (\mathbf{r})\widehat{z}$ as the source in Eq. (5), resulting in electric fields of the desired polarization. The fiber cross section is represented by a $3\lambda \times 3\lambda $ computational cell consisting of $300\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{pixel}\times 300\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{pixel}$, where the size of each pixel is $0.01\lambda \times 0.01\lambda $. From Fig. 2 (inset), it is clear that both the fundamental and thirdharmonic modes are well confined to the fiber core and exhibit substantial modal overlaps, while again, the phasematching condition is automatically satisfied by the optimization process, with ${k}_{3}=3{k}_{\mathrm{opt}}$. We find that ${{\beta}_{3}}^{2}\approx 2\times {10}^{4}({\chi}^{(3)}/{\lambda}^{4})$ is almost 4 orders of magnitude larger than what has been demonstrated in standard plain fibers, which have typical values of ${{\beta}_{3}}^{2}\lesssim 2({\chi}^{(3)}/{\lambda}^{4})$ ^{[77]}. Figure 2 shows the dispersion of the two leaky modes (solid lines), with the PES lightline represented by the gray region and their corresponding dimensionless lifetimes, around ${Q}_{1}\approx {10}^{6}$ and ${Q}_{3}\approx {10}^{5}$ at ${k}_{\mathrm{opt}}$, plotted as dashed lines. Noticeably, while the fiber is optimized to ensure phase matching at a single ${k}_{\mathrm{opt}}$, any phase mismatch remains small in the vicinity of $k\approx {k}_{\mathrm{opt}}$. In fact, even for $k\ll {k}_{\mathrm{opt}}$, the frequency difference is found to be only around 1%. Technically, the only factor limiting the lifetimes is the finite computational crosssection (imposed by the finite computational cell), with much larger lifetimes possible for larger cross sections. Away from ${k}_{\mathrm{opt}}$, the quality factors decrease while remaining relatively large over a wide range of $k$. Considering the group velocity ${v}_{g}$ around ${k}_{\mathrm{opt}}$, we find that the attenuation length of the fiber ${L}_{\mathrm{rad}}=1/\gamma \approx 2{v}_{g}Q/\omega =1.66\times {10}^{5}\lambda $. We note that while the fiber supports multiple modes around these wavelengths, the only modes near ${k}_{\mathrm{opt}}$ are those discovered by the optimization and shown in the figure.
Figure 2.Dispersion relations (solid line) and radiative lifetimes $Q$ (dashed line) versus propagation wavenumber $k$ of ${\mathrm{TM}}_{01}$ fundamental ${\omega}_{1}$ (red) and thirdharmonic ${\omega}_{3}$ (blue) modes in a chalcogenide/PES fiber optimized to achieve frequency matching ${\omega}_{3}=3{\omega}_{1}$ and large nonlinear overlaps at ${k}_{\mathrm{opt}}=1.4(2\pi /\lambda )$. The shaded area in gray indicates regions lying below the chalcogenide light cone. The top insets show the fiber cross section overlaid with corresponding power densities at ${\omega}_{1}$ (left) and ${\omega}_{3}$ (right).
Figure 3 shows the ${\beta}_{3}$ corresponding to fibers optimized for operation at different values of ${k}_{\mathrm{opt}}$ and polarizations, and obtained by application of either topology (squares or circles) or shape (triangles) optimization. The figure shows a general trend in which ${\beta}_{3}$ decreases with increasing ${k}_{\mathrm{opt}}$ for both polarizations, except that ${\mathrm{TM}}_{01}$ fibers tend to exhibit nonmonotonic behavior, with ${\beta}_{3}$ increasing sharply at an intermediate ${k}_{\mathrm{opt}}\approx 2\pi /\lambda $ below the lightline, above which it drops significantly before increasing again in the guided regime, peaking again at ${k}_{\mathrm{opt}}\approx 1.7(2\pi /\lambda )$ before plummeting once again. We suspect that this complicated behavior is not a consequence of any fundamental limitation or physical consideration, but rather stems from the optimization algorithm getting stuck in local minima. Regardless, our results provide a proofofprinciple of the existence of fiber designs with performance characteristics that can greatly surpass those of traditional, handdesigned fibers. Furthermore, Fig. 3 shows typical fiber cross sections at selective ${k}_{\mathrm{opt}}$, along with their corresponding superimposed (fundamental) mode profiles, illustrating the fabricability of the structures, in which the structure via shape optimization [Fig. 3(iv)] is easiest for fabrication.
Figure 3.Nonlinear overlap factor ${{\beta}_{3}}^{2}$ corresponding to fundamental and thirdharmonic modes in fibers that have been optimized to ensure phasematched modes (${k}_{3}=3{k}_{\mathrm{opt}}$) at various fundamentalmode propagation wavenumbers ${k}_{\mathrm{opt}}$, for both ${\mathrm{TE}}_{01}$ (blue) and ${\mathrm{TM}}_{01}$ (red) polarizations, by the application of either topology (circles or squares) or shape (triangles) optimization. The grayshaded area denotes the regime of guided modes below the chalcogenide lightline. For comparison, also shown is ${{\beta}_{3}}^{2}$ (black cross) of a standard plain fiber manually designed for operation at ${\omega}_{1}=0.914(2\pi c/\lambda )$ and ${k}_{1}=0.992(2\pi /\lambda )$ [77]. Shown as insets are fiber cross sections along with power densities of fundamental modes optimized at four different ${k}_{\mathrm{opt}}=\{0.1,1.4,1.7,2.0\}(2\pi /\lambda )$ for both ${\mathrm{TE}}_{01}$ (upper insets) and ${\mathrm{TM}}_{01}$ (lower insets), with (i)–(iii) obtained via topology optimization and (iv) via shape optimization.
Finally, we provide estimates of the power requirements associated with these fiber designs. We find that, for a ${\mathrm{TM}}_{01}$ fiber operating at ${k}_{\mathrm{opt}}=1.4(2\pi /\lambda )$ and at a wavelength of $\lambda =1\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{\mu m}$, conversion efficiencies of 50% can be attained at relatively small pump powers ${P}_{1}\approx 1.7\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mW}$ over a fiber segment $L\approx 3\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{cm}$, while the corresponding (radiative) attenuation lengths are $\approx 17\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{cm}$. For comparison, plain silica fibers ^{[77]} exhibit modeoverlap factors ${\beta}_{3}\approx 2({\chi}^{(3)}/{\lambda}^{2})$, leading to conversion efficiencies of the order of ${10}^{8}\%$ for the same input power and fiber length. (Note that typical PhC fibers rely on enhancements such as slow light effects ^{[47]}, exhibiting even poorer ${\beta}_{3}$ than that of a fiber.) Hence, the optimized structures achieve considerably ($\sim {10}^{9}$ times) higher conversion efficiencies, an improvement that is only partially due to the larger ${\chi}^{(3)}$ of chalcogenide compared to glass (approximately 440 times larger). In particular, defining the normalized interaction fiber length $L({\chi}^{(3)})$, which removes any source of material enhancement, we find that the optimized fiber leads to a factor of ${10}^{4}$ enhancement. Similarly, we find that a ${\mathrm{TE}}_{01}$ fiber operating at ${k}_{\mathrm{opt}}=0.2(2\pi /\lambda )$ results in a factor of ${10}^{3}$ enhancement compared to plain fibers.
4. SECONDHARMONIC GENERATION IN METASURFACES
Metasurfaces offer an advantageous platform for realizing complicated beam generation and wavefront shaping over extended surfaces ^{[79]} and have recently been exploited in conjunction with nonlinear materials as a means of generating and controlling light at multiple wavelengths ^{[43,48,80,81]}. A typical nonlinear metasurface can suffer from poor frequencyconversion efficiencies due to a combination of weak confinement, material absorption, and suboptimal mode overlaps. In particular, typical designs exploit plasmonic ^{[38–40,50]} or alldielectric ^{[42,49]} elements comprising simple shapes distributed over a unit cell, including split ring resonators ^{[38,40,50]}, cross bars ^{[39]}, and cylindrical posts ^{[49]}, with the main focus being that of satisfying the requisite frequency and phasematching condition ^{[82]}. Here, we show that inverse design can not only facilitate the enforcement of frequency and phasematching requirements but also allow further enhancements stemming from the intentional engineering of nonlinear modal overlaps, often neglected in typical designs.
To achieve large secondharmonic generation efficiencies, a metasurface must support two extended resonances at frequencies ${\omega}_{1}$ and ${\omega}_{2}=2{\omega}_{1}$ and wavevectors satisfying the phasematching condition ${\mathbf{k}}_{2}=2{\mathbf{k}}_{1}$. As illustrated schematically in Fig. 4(a), a typical setup consists of an incident wave of power per unit cell ${P}_{1}$ at some frequency and angle (described by wavenumber ${\mathbf{k}}_{1}$) and a corresponding output harmonic wave of power per unit cell, ${P}_{2}$. In the smallsignal regime, the output power ${P}_{2}\propto {P}_{1}^{2}$ scales quadratically with ${P}_{1}$, resulting in a conversion efficiency per unit cell of where $Q$ and ${Q}_{\mathrm{rad}}$ denote total and radiative dimensionless lifetimes and ${\beta}_{2}$ the nonlinear overlap factor:
Figure 4.(a) Schematic illustration of secondharmonic generation in a squarelattice metasurface of finite thickness $t$ and period $\mathrm{\Lambda}\times \mathrm{\Lambda}$. Shown to the right are dielectric profiles and mode profiles ${\mathbf{E}}^{2}$ corresponding to two inversedesigned metasurfaces, both over single unit cells and $z=0$ cross sections. The structures are optimized to ensure frequency and phase matching for light incident at (i) an angle $\theta =3\xb0$ or (ii) normal incidence. Dark (white) represents gallium phosphide (vacuum) regions. (b) Convergence of the objective function with respect to iteration number, leading to structure (ii).
Note that here the conversion efficiency is defined as the efficiency per unit cell for such an extended surface, hence the volume integration is performed inside a unit cell.
We now apply our optimization framework to discover new alldielectric threedimensional metasurfaces, with the permittivity of the medium ${\epsilon}_{\mathrm{GaP}}$ taken to be that of gallium phosphide near telecom wavelengths ^{[83,84]}. Note, however, that the same framework can be easily extended to design plasmonic surfaces. The metasurfaces illustrated schematically in Fig. 4 are square PhC slabs of inplane periodicity $\mathrm{\Lambda}\times \mathrm{\Lambda}$ and finite thickness $t$. To ensure fabricability, here we consider $z$invariant structures, in which the optimization parameters are taken to lie in the plane perpendicular to the $z$ axis, resulting in a structure that can be fabricated by etching. As a proofofprinciple, we consider metasurfaces suspended in air, while the same framework can be easily applied to include any substrate ^{[27]}.
Figure 4 shows cross sections of the unit cell of two GaP metasurfaces of thickness $t=612\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$ and $\mathrm{\Lambda}=480\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$, designed for operation at a fundamental frequency ${\omega}_{1}=1.57\times {10}^{15}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{rad}/\mathrm{s}$ ($\lambda =1.2\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{\mu m}$) so as to satisfy both frequency and phasematching conditions. Also shown are the corresponding fundamental and harmonic mode profiles. The structure on the left is optimized for operation at an incident angle $\theta \approx 3.6\xb0$ relative to the outofplane axis, and is found to exhibit large radiative lifetimes ${Q}_{1(2)}^{\mathrm{rad}}\approx 6(2)\times {10}^{4}$ and overlap factor ${{\beta}_{2}}^{2}=1.6\times {10}^{3}({\chi}^{(2)}/{\lambda}^{3})$. The structure on the right is instead optimized for operation at normal incidence, resulting in a slightly smaller ${{\beta}_{2}}^{2}=4\times {10}^{4}({\chi}^{(2)}/{\lambda}^{3})$. Because of the symmetry of the structure, the modes exhibit infinite lifetimes (and, hence, are technically dark modes), though in practice, fabrication imperfections necessarily lead to finite lifetimes. Furthermore, Fig. 4(b) illustrates the convergence of the TO optimization process to achieve structure (i), converged within $\sim {10}^{3}$ iterations. Table 1 compares a few of the relevant FOMs for representative metasurface designs, which include both plasmonic and dielectric structures. Although comparing ${\beta}_{2}$ appears to be impossible due to a surprising lack of relevant modal parameters in these studies ^{[38,39,49,50]}, such as the absence of radiative and dissipative quality factors, we find that the optimized designs exhibit orders of magnitude larger conversion efficiencies. While it is difficult to distinguish the relative impact of the mode lifetimes and overlap factors, arguably, the optimized structures overcome several limitations associated with previous designs. On the one hand, plasmonic structures exhibit tightly confined modes and therefore lead to large nonlinear overlaps, but absorptive losses and weak material nonlinearities imply that they suffer from small lifetimes. On the other hand, several of the proposed alldielectric metasurfaces have had negligible material losses and, hence, larger lifetimes, but have not been designed to ensure large nonlinear overlaps.
Table 1. Representative SecondHarmonic Generation FOMs for Both Hand and InverseDesigned Metasurfaces, Including ${\chi}^{(2)}$, Fundamental Wavelength ${\lambda}_{1}$, and Conversion Efficiency $\eta $ per Unit Cell^{a}
Table 1. Representative SecondHarmonic Generation FOMs for Both Hand and InverseDesigned Metasurfaces, Including ${\chi}^{(2)}$, Fundamental Wavelength ${\lambda}_{1}$, and Conversion Efficiency $\eta $ per Unit Cell^{a}
Structure  χ(2)(nm/V)  λ1(μm)  η/(χ(2))2  Gold split resonators [38]  250  10  2.1×1011  Gold split resonators [50]  1.3  3.4  3.8×1011  Gold cross bars [39]  54  8  1.4×1013  Alldielectric cylinders [49]  0.2  1.02  1.6×1017  Optimized design [Fig. 4]  0.1  1.2  9.6×1024 

5. CONCLUDING REMARKS
We have demonstrated an optimization approach for the design of nonlinear photonic fibers and metasurfaces. The optimized structures demonstrate very high leakymode lifetimes for both fundamental and harmonic modes and ordersofmagnitude larger overlap factors than traditional designs. Inverse design not only overcomes efficiency limitations of traditional index fibers and PhC metasurfaces but also greatly reduces challenges and difficulties inherent to the design process. Although in this paper we have not considered effects resulting from self or crossphase modulation, we expect no significant impact on the conversion efficiency in the smallsignal limit, since the finite bandwidth around the designated phasematched propagation wavevectors can potentially compensate for any small phase mismatch that might arise. At larger powers where these effects cannot be ignored, one could account and compensate for them through minor modifications to the optimization objective function, the subject of future work. Furthermore, we will consider extending our inversedesign framework to terahertz frequency generation and other nonlinear processes.