# Inverse Design for Waveguide Dispersion with a Differentiable Mode Solver

Dodd Gray  1,2,* Gavin N. West  1 Rajeev J. Ram1  1MIT, 77 Massachusetts
Avenue, Cambridge, MA 02139, USA
2MIT Lincoln Laboratory, 244 Wood Street, Lexington, MA 02421, USA
<EMAIL_ADDRESS>
††journal: opticajournal††articletype: Research Article

Inverse design of optical components based on adjoint sensitivity analysis has
the potential to address the most challenging photonic engineering problems.
However existing inverse design tools based on finite-difference-time-domain
(FDTD) models are poorly suited for optimizing waveguide modes for adiabatic
transformation or perturbative coupling, which lies at the heart of many
important photonic devices. Among these, dispersion engineering of optical
waveguides is especially challenging in ultrafast and nonlinear optical
applications involving broad optical bandwidths and frequency-dependent
anisotropic dielectric material response. In this work we develop gradient
back-propagation through a general purpose electromagnetic eigenmode solver
and use it to demonstrate waveguide dispersion optimization for second
harmonic generation with maximized phase-matching bandwidth. This optimization
of three design parameters converges in eight steps, reducing the
computational cost of optimization by $\sim$100x compared to exhaustive search
and identifies new designs for broadband optical frequency doubling of laser
sources in the 1.3–$1.4\text{\,}\mathrm{\SIUnitSymbolMicro m}$ wavelength
range. Furthermore we demonstrate that the computational cost of gradient
back-propagation is independent of the number of parameters, as required for
optimization of complex geometries. This technique enables practical inverse
design for a broad range of previously intractable photonic devices.

##  1 Introduction

Inverse design methods based on adjoint sensitivity analysis have been
demonstrated as a promising toolkit for photonic design [1], enabling tunable
high performance devices through optimization of complex designs with large
numbers of parameters. These techniques generally rely on gradients of
electromagnetic models calculated using the adjoint method with a
computational cost that is independent of the number of optimized parameters.
A variety of mature software tools exist for inverse design of optical devices
based on differentiable finite-difference-time-domain (FDTD) [2, 3] and
finite-difference-frequency-domain (FDFD) models [4]. These have been used to
optimize light scattering between fixed input and output modes [5, 6, 7, 8]
for a wide range of applications but are poorly suited for optimization of the
modal dispersion of optical resonances and waveguide modes.

Modeling of waveguide mode fields and modal dispersion is ubiquitous and
indispensable in integrated photonics. Many important design problems,
especially those involving adiabatic or perturbative coupling of waveguide
modes require differentiable mode solving for compatibility with sensitivity
analysis and inverse design. The efficiency of nonlinear optical interactions
are particularly sensitive to dispersion [9, 10, 11] because coupling is
achieved over large optical bandwidths and long propagation lengths. Various
approaches to modal dispersion optimization have been demonstrated, including
neural networks [12, 13, 14, 15], model order reduction [16], automatic
differentiation [17, 18, 19, 20, 21], and others [22, 23, 24, 25]. None of
these simultaneously support arbitrary, non-periodic waveguide cross sections
and account for dielectric material anisotropy and frequency-dependence,
limiting their general utility for photonic design.

Surprisingly, optimization of dispersion-related objective functions with only
two or three geometry parameters can be computational expensive. For example
second harmonic generation (SHG) or three-wave mixing involve frequencies that
span at least one octave. This wide optical bandwidth, coupled with material
birefringence often requires a large number of modes to be modeled (>10) for
each set of geometric parameters and wavelengths. The underlying material
anisotropy also necessitates solving the Helmholtz Equation with tensor-valued
dielectric constants and introduces mode degeneracies which can frustrate
numerical convergence. Finally, the broadband nonlinear phenomena considered
here are extremely sensitive to small changes in geometry – a fact that has
driven the development of adaptive fabrication strategies in thin film lithium
niobate (TFLN) on insulator waveguides [26, 27] – necessitating the use of a
fine discretization mesh and a large number of simulations for brute-force
optimization of the geometry111 We consider the design problem discussed in
Section 3. On a conventional workstation (16 cores, 3 GHz clock speed, 32 GB
RAM) robustly solving the Helmholtz Equation for fundamental and second
harmonic TE00 modes in isotropic media would require $10\text{\,}\mathrm{s}$
(4 modes), however the corresponding simulations with anisotropic materials
require $60\text{\,}\mathrm{s}$ (12 modes). A coarse sweep of waveguide
geometry parameters (20 widths $\times$ 20 thicknesses $\times$ 10 etch
depths) considering SHG of a single wavelength requires 67 hours (48000 modes)
if solved in serial, or 8 hours with 8 mode solvers running concurrently. This
increases linearly with wavelengths and exponentially with parameters. .
Indeed, the existing literature on SHG bandwidth optimization only considers
two [9] or three [28, 29] waveguide geometry parameters (eg. width, thickness,
and etch depth). Additional design parameters are easily accessible
experimentally (eg. upper and lower cladding thickness and side-wall angle)
but have not been considered because of the exponential scaling of
computational effort. Clearly a more efficient approach to optimization is
required. Here, differentiable mode solvers with gradient-based optimization
are shown to make such dispersion-related optimization problems tractable.
This class of inverse design tools efficiently manages the sensitivity to
geometry by using parameter gradients to adaptively step through geometric
parameters even in the presence of strong material dispersion and anisotropy.
As we demonstrate below, these gradients are calculated efficiently for a
large number of parameters by exploiting automatic differentiation. Our
implementation enables incorporation of waveguide mode solutions as a modular
component in larger differentiable models, including concurrent evaluation of
mode solutions at multiple frequencies for time-efficient broadband dispersion
optimization.

In Section 2 we describe the implementation of gradient backpropagation
through a fully vectorial plane-wave expansion mode solver compatible with
anisotropic and frequency-dependent dielectric materials. In Section 3 we
demonstrate its application to waveguide dispersion optimization for nonlinear
optical applications. In Section 4 we describe the performance and scaling of
our differentiable mode solver. Finally in Section 5 we summarize our results
and discuss future work.

##  2 Implementation

The workflow for inverse design of waveguides with optimized modal dispersion
is depicted in Fig. 1–1. First a function is defined mapping input design
parameters to dielectric tensor values on a discretized spatial grid. An
example of such parametric waveguide design is shown in Fig. 1. The
computational steps we used to numerically model and optimize waveguide modes
for a desired dispersion are schematically depicted in Fig. 1. We use desired
modal properties illustrated in Fig. 1 such as effective index, group velocity
(or group index) and group velocity dispersion (GVD) at target optical
frequencies to determine optimal waveguide designs. To achieve this we
supplied with the value and sensitivity to design parameters of an objective
function quantifying a desired dispersion property to a generic optimization
algorithm [30]. Allowing the optimization algorithm to iteratively update
design parameters in a feedback loop rapidly improved the cost function,
typically converging on a local optimum after only a handful of iterations.
The sequence of calculations comprising a single iteration is shown in Fig.
1–1 and described below.

Figure 1:  1 Rib waveguide design parameterized by width, thickness and etch
depth. 1 Block diagram showing computational modeling of waveguide modes and a
cost function dependent on their 1 dispersion, then back-propagation of cost
function gradient components to the input design parameters. One optimization
step consists of the following steps: 1 instantiate the parametric shapes and
assign bulk materials. Next 1 generate a smooth dielectric function. 1
Eigenmodes are solved using a plane-wave expansion-based solver. The modal
indices and field distributions are used to calculate 1 the modal dispersion
(here the effective group index) of the specified waveguide In the final step
of the forward calculation we compute a dispersion and/or mode-field dependent
objective function quantifying waveguide performance. Gradients then propagate
backwards from the objective function via automatic differentiation to 1
gradients with respect to the computed waveguide modes and their dispersion.
To back-propagate gradient components through mode solutions we compute 1
corresponding adjoint fields for each mode impacting the objective function.
From the adjoint fields we calculate 1 Gradients with respect to the smoothed
dielectric tensor at each grid point. Finally we use automatic differentiation
to calculate 1 gradients with respect to the input design parameters, which
are shape parameters in this example. The initial parameters, objective
function value, and parameter gradients are passed to a gradient-based
optimizer which computes updated parameter values for the next iteration.

###  2.1 Parametric Geometry and Materials

To model resonant modes of a dielectric waveguide we write a function mapping
a set of shape and material parameters to dielectric tensor values at each
spatial grid point in our simulation. Figure 1 illustrates a geometry
instantiated through a constructive geometry library [31] which provides
parametric shape primitives to programatically generate waveguide geometries
as lists of shapes with assigned bulk materials.

We implement an interface for arbitrary models of material dielectric tensors
as a function of frequency and other parameters, such as temperature, alloy
concentration, and orientation. The first and second derivatives w.r.t
frequency of all material bulk dielectric tensors, $d\varepsilon/d\omega$ and
$d^{2}\varepsilon/d\omega^{2}$, are calculated using a computer algebra
package [32]. These derivatives are later used in the calculation of group
velocity and group velocity dispersion. Analytic pre-computation is
inexpensive, improves numerical precision, and reduces the burden on later
auto-differentiation steps.

###  2.2 Dielectric Smoothing

Once a geometry is specified and material dielectric tensors are computed, the
dielectric distribution must be smoothed across geometric/material boundaries,
seen in Fig. 1. This ensures that the resonant frequencies and mode fields of
electromagnetic eigenmode solutions vary smoothly with design parameters [33],
and corrects for spatial discretization errors to first order. We implement
the anisotropic dielectric smoothing algorithm developed by Kottke, et al.
[34] to calculate values of the smoothed dielectric tensor, here denoted
$\tilde{\varepsilon}$, and extend the algorithm to calculate first and second
frequency derivatives of the smoothed tensor. Near a dielectric interface with
surface normal $\hat{n}$ between materials with dielectric tensors
$\varepsilon_{1}$ and $\varepsilon_{2}$,

$\tilde{\varepsilon}=\mathbf{R}_{\hat{n}}\left[\tau^{-1}\left(r\tau(\mathbf{R}_{\hat{n}}^{T}\varepsilon_{1}\mathbf{R}_{\hat{n}})+(1-r)\tau(\mathbf{R}_{\hat{n}}^{T}\varepsilon_{2}\mathbf{R}_{\hat{n}})\right)\right]\mathbf{R}_{\hat{n}}^{T}$
(1)

where $\mathbf{R}_{\hat{n}}$ is a rotation matrix mapping between the grid
axes and a coordinate frame with the first coordinate along the local surface
normal of the material interface $\hat{n}$, $r$ is the fill fraction of
material 1 in the grid pixel or voxel, and $\tau(\varepsilon)$ and
$\tau^{-1}(\tau)$ are 3$\times$3 matrix-valued functions of 3$\times$3
matrices derived in [34]. The operator $\mathbf{R}_{\hat{n}}$ is constructed
as a function of $\hat{n}$ by ortho-normalization of successive cross products

$\mathbf{R}_{\hat{n}}(\hat{n})=\begin{bmatrix}\hat{n}&\frac{\hat{n}\times\hat{z}}{|\hat{n}\times\hat{z}|}&\frac{\hat{n}\times\frac{\hat{n}\times\hat{z}}{|\hat{n}\times\hat{z}|}}{\left|\hat{n}\times\frac{\hat{n}\times\hat{z}}{|\hat{n}\times\hat{z}|}\right|}\end{bmatrix}$
(2)

which is a differentiable function of $\hat{n}$ as long as $\hat{n}$ is not
parallel to $\hat{z}$. For 3D models a generalized form Eq. 2 can be used to
be continuous for all $\hat{n}$ [35]. Equations 1-2 define a differentiable
map from $r$, $\hat{n}$, $\varepsilon_{1}$ and $\varepsilon_{2}$ to
$\tilde{\varepsilon}$,

$\mathrm{vec}(\tilde{\varepsilon})=F_{\tilde{\varepsilon}}(\mathrm{vec}(r,\hat{n},\varepsilon_{1},\varepsilon_{2}))$
(3)

where $\mathrm{vec}(\dots)$ denotes flattening into a vector. We use a
computer algebra system [32] to analytically compute and pre-compile functions
for the Jacobian and Hessian matrices $J_{\tilde{\varepsilon}}$ and
$H_{\tilde{\varepsilon}}$ of the function $F_{\tilde{\varepsilon}}$ in Eqn. 3.
With these we calculate the first two frequency derivatives of the smoothed
dielectric tensor $d\tilde{\varepsilon}/d\omega$ and
$d^{2}\tilde{\varepsilon}/d\omega^{2}$ at each spatial grid point as functions
of the bulk dielectric tensors $\varepsilon_{1}$ and $\varepsilon_{2}$ and
their first two frequency derivatives

$\displaystyle d\tilde{\varepsilon}/d\omega=$ $\displaystyle
J_{\tilde{\varepsilon}}\vec{v}_{1}$ (4) $\displaystyle
d^{2}\tilde{\varepsilon}/d\omega^{2}=$
$\displaystyle\vec{v}_{1}^{T}H_{\tilde{\varepsilon}}\vec{v}_{1}+J_{\tilde{\varepsilon}}\vec{v}_{2}$
(5)

where $\vec{v}_{1}$ and $\vec{v}_{2}$ are vectors containing the first and
second derivatives of the material dielectric tensors w.r.t frequency and
zeros corresponding to $r$ and $\hat{n}$ which are frequency independent,

$\displaystyle\vec{v}_{1}=$
$\displaystyle\mathrm{vec}(0,[0,0,0],d\varepsilon_{1}/d\omega,d\varepsilon_{2}/d\omega))$
(6) $\displaystyle\vec{v}_{2}=$
$\displaystyle\mathrm{vec}(0,[0,0,0],d^{2}\varepsilon_{1}/d\omega^{2},d^{2}\varepsilon_{2}/d\omega^{2}))$
(7)

Moving forward we will refer to the smoothed dielectric tensor simply as
$\varepsilon$ to simplify notation.

###  2.3 Electromagnetic Eigenmode Solver

Our differentiable electromagnetic mode solver is adapted from the plane-wave
expansion algorithm in the open-source MIT Photonic Bands (MPB) [36] software
package. The algorithm leverages discrete Fourier transforms (DFTs) [37] to
efficiently apply curl and dielectric tensor operations to the electromagnetic
mode field in reciprocal and real space, respectively. We used the Julia
programming language [38] to implement a matrix-free discretized Helmholtz
operator and used generic iterative solvers based on Krylov subspace methods
[39, 40, 41] for partial eigen-decomposition to solve for a specified number
of approximate eigenvector-eigenvalue pairs. Our primary motivation for using
Julia was compatibility with automatic differentiation, discussed below in
Sec. 2.5. In this work we focus on laterally confined waveguide modes, but the
following analysis does not require this assumption. For reasons discussed
below it is most convenient to solve Eq. 9 with periodic boundary conditions,
despite the fact that we are interested in dielectric waveguide geometries
which are not periodic in the dimensions transverse to the propagation axis
(here the $x-y$ plane).

We consider a Helmholtz Equation for the magnetic field $H_{\rm n}$ of an
electromagnetic waveguide mode propagating in the $+\hat{z}$ direction (along
the waveguide axis), parameterized by dielectric tensor values
$\varepsilon(x,y,z)$ at each point on a regular spatial grid $\left\\{x_{\rm
i},y_{\rm i},z_{\rm i}\right\\}$ and a wave vector $k\hat{z}$ with magnitude
$k$

$\mathbf{M}(k,\varepsilon)H_{\rm n}=\omega_{\rm n}^{2}H_{\rm n}$ (8)

The Helmholtz operator $\mathbf{M}$ consists of two curl operations
sandwiching multiplication by the inverse dielectric tensor
$\varepsilon^{-1}$.

$\mathbf{M}(k,\varepsilon)=\left[\nabla\times\right]\left[\varepsilon^{-1}\cdot\right]\left[\nabla\times\right]$
(9)

Expressing the mode field and curl operators in a plane wave basis we have

$\mathbf{M}(k,\varepsilon)\tilde{h}_{\rm
n}=\left[\mathbf{C}^{\dagger}\varepsilon^{-1}\mathbf{C}\right]\tilde{h}_{\rm
n}=\omega^{2}_{\rm n}\tilde{h}_{\rm n}$ (10)

where $\tilde{h}_{\rm n}$ is the projection of $H_{\rm n}$ onto transversely
polarized plane wave components at each reciprocal lattice point and
$\mathbf{C}$ represents the fused DFT, Cartesian-to-transverse projection, and
curl operators. The electromagnetic Helmholtz equation as expressed in Eq. 10
enables computationally efficient approximate numerical eigenmode solutions
[36]. In many dispersion engineering applications it is important to account
for the dependence of the dielectric tensor $\varepsilon$ on optical frequency
$\omega$. In this case the Helmholtz equation becomes a nonlinear eigenproblem
with an eigenvalue-dependent operator.

$\mathbf{M}\left(k,\varepsilon(\omega_{\rm n})\right)H_{\rm n}=\omega_{\rm
n}^{2}H_{\rm n}.$ (11)

Solutions to Eq. 11 are found by iteratively solving Eq. 10 using a Newton
method to update $k$. We will refer to Eq. 10 as the “fixed-$k$” eigenproblem
and Eq. 11 as the “fixed-$\omega$” eigenproblem. Supporting derivations and
the definition of $H_{\rm n}$, $\tilde{h}_{\rm n}$, and $\mathbf{C}$ are
available in Supplement 1. For mode sorting and subsequent calculations, we
transform $\tilde{h}_{\rm n}$ to a real-space electric field mode $E_{n}$ as
shown in Fig. 1.

###  2.4 Modal Group Velocities, GVDs and the Objective Function

To calculate frequency dispersion relationships such as shown in Fig. 1 we
derive AD-compatible expressions for the frequency derivatives of the modal
index. The modal group index is calculated as the ratio between Poynting flux
and energy density in mode volume [42]. We express the modal group index in
terms of the plane-wave basis fields as

$n_{\rm g,n}=\frac{\omega_{\rm n}^{2}+\frac{1}{2}\omega_{\rm
n}\sum\tilde{h}_{\rm
n}^{\dagger}\left[\mathbf{C}^{\dagger}\varepsilon^{-1}\frac{d\varepsilon}{d\omega}\varepsilon^{-1}\mathbf{C}\right]\tilde{h}_{\rm
n}}{\omega_{\rm n}\sum\tilde{h}_{\rm
n}^{\dagger}\frac{\partial\mathbf{M}}{\partial k}\tilde{h}_{\rm n}}$ (12)

where $\frac{d\varepsilon}{d\omega}$ is the first frequency derivative of the
smoothed dielectric tensor defined in Eq. 4.

The modal group velocity dispersion (GVD)
$\text{GVD}_{n}=\frac{\text{d}n_{g,n}}{\text{d}\omega}$ can also be computed
from a single mode solution using the adjoint field corresponding to group
index and the second frequency derivative of the smoothed dielectric tensor,
$\frac{d^{2}\varepsilon}{d\omega^{2}}$, defined in Eq. 5. Expressions for the
GVD are derived in Supplement 1. Group index and GVD values calculated
directly from a single mode solution avoid numerical error from finite
difference approximation [43] and can be optimized using fewer total eigen-
and adjoint problem solutions.

In the final step of the forward calculation we compute an objective function
quantifying the waveguide performance. The objective function can depend in an
arbitrary way on the calculated mode fields and dispersion quantities as well
as the input parameters and intermediate quantities such as the dielectric
function.

###  2.5 Automatic Differentiation

Efficient calculation of the objective function value and gradient is highly
desirable when optimizing of objective functions in high-dimensional parameter
spaces [43]. The standard Julia automatic differentiation framework provides
extensive compatibility with native Julia functions [44], enabling convenient
compiled gradient computations of arbitrary user-defined functions.

We defined “pullback” functions for the fixed-$k$ and fixed-$\omega$ Helmholtz
eigen-problems using the adjoint method to map partial derivatives with
respect to mode fields and wave vector magnitudes back to partial derivatives
w.r.t the optical frequency, dielectric tensor elements, and material
dispersion at each spatial grid point (the Helmholtz operator parameters),
described in Sec. 2.6. These custom AD rules enable differentiation of cost
functions with arbitrary dependence on mode solutions. They do not require
adjustment or re-definition when the mode simulation or cost function is
modified. We note that an AD-compatible constructive geometry package is
available in Julia [31], enabling propagation of gradients w.r.t the
dielectric function (Fig. 1) to gradients in the shape parameters (Fig. 1).

###  2.6 Adjoint Sensitivity Analysis

To provide an unbroken chain of derivatives to map sensitivities between the
objective function and input parameters, we apply the adjoint method to
solutions of Eq. 11. We begin by following the approach of [45] to derive the
adjoint problem for the fixed-$k$ Helmholtz equation, then extend this
derivation to the fixed-$\omega$ Helmholtz equation. The sensitivities of the
objective function with respect to the wave vector magnitude $k$ and smoothed
permittivity $\varepsilon$ for the linear problem are given by

$\displaystyle\frac{\partial g}{\partial k}=$ $\displaystyle\frac{\partial
g}{\partial\omega_{\rm n}^{2}}\frac{\partial\omega^{2}_{\rm n}}{\partial
k}+\sum_{j}\Tr\left(\matrixquantity[-(v_{j\rm z}\hat{n}_{j}+n_{j\rm
z}\vec{v}_{j})^{T}\\\ (v_{j\rm z}\hat{m}_{j}+m_{j\rm
z}\vec{v}_{j})^{T}]\left[\tilde{E}_{\rm n}\tilde{h}_{\lambda\rm
n}^{\dagger}+\tilde{E}_{\rm\lambda n}\tilde{h}_{\rm
n}^{\dagger}\right]\right)$ (13) $\displaystyle\frac{\partial
g}{\partial\varepsilon}=$ $\displaystyle E_{\rm\lambda n}E_{\rm n}^{\dagger}$
(14)

where $\tilde{h}_{\lambda\rm n}$ is the adjoint mode field shown in Fig. 1 and
where we have used the simplifying transformations

$\displaystyle\frac{\partial\omega_{n}^{2}}{\partial k}=$ $\displaystyle
2\frac{\omega_{n}}{n_{g,n}}$ $\displaystyle\vec{v}_{j}=$ $\displaystyle
k\hat{z}+\vec{g}_{j}$

where each $\vec{g}_{\rm j}$ is a reciprocal lattice point in the plane wave
basis. We introduce the nonlinearity of Eq. 11 through an equality constraint
(see Supplement 1) and differentiate to find

$\displaystyle\frac{\partial k}{\partial\varepsilon}=$
$\displaystyle-\frac{\partial k}{\partial\omega^{2}}E_{\rm n}E_{\rm
n}^{\dagger}$ (15) $\displaystyle\frac{\partial g}{\partial\omega^{2}}=$
$\displaystyle\left[\frac{\partial g}{\partial k}+\frac{\partial
g}{\partial\tilde{h}_{\rm n}}\frac{\partial\tilde{h}_{\rm n}}{\partial
k}\right]\frac{\partial k}{\partial\omega^{2}}$ (16)
$\displaystyle\frac{\partial g}{\partial\varepsilon}=$
$\displaystyle\frac{\partial
g}{\partial\omega^{2}}\frac{\partial\omega^{2}}{\partial\varepsilon}$ (17)

where in Eq. 16

$\frac{\partial g}{\partial\tilde{h}_{\rm n}}\frac{\partial\tilde{h}_{\rm
n}}{\partial k}=\sum_{j}\Tr\left(\matrixquantity[-(v_{j\rm
z}\hat{n}_{j}+n_{j\rm z}\vec{v}_{j})^{T}\\\ (v_{j\rm z}\hat{m}_{j}+m_{j\rm
z}\vec{v}_{j})^{T}]\left[\tilde{E}_{\rm n}\tilde{h}_{\lambda\rm
n}^{\dagger}+\tilde{E}_{\rm\lambda n}\tilde{h}_{\rm
n}^{\dagger}\right]\right)$ (18)

as in Eq. 13 and requires calculation of the adjoint field. Using Eqs. 15-17
we define a $\left\\{\frac{\partial g}{\partial k},\frac{\partial
g}{\partial\tilde{h}_{\rm n}}\right\\}\mapsto\left\\{\frac{\partial
g}{\partial\omega_{\rm n}},\frac{\partial g}{\partial\varepsilon}\right\\}$
pullback map for gradient back-propagation through the fixed-$\omega$
eigenproblem.

##  3 Optimization of Second Harmonic Conversion Bandwidth

To demonstrate the utility of this design approach we explore the optimization
of an etched TFLN waveguide for maximized second harmonic generation (SHG)
phase-matching bandwidth. We chose this design goal because broadband
simultaneous phase matching of optical three-wave mixing in TFLN waveguides is
of great interest for a wide range of applications[9, 46, 28, 47, 48, 49, 50,
51, 52, 53, 54, 55]. In this example the rib waveguide geometry is
parameterized by core width, core thickness and normalized etch fraction. The
phase mismatch between fundamental and second harmonic quasi-TE00 modes is
assumed to be compensated by periodic poling of the TFLN [56]. We perform the
optimization purely based on the desired group indices and calculate the
required poling period from the corresponding phase indices. This is practical
because inverting the TFLN domains has no impact on the linear dielectric
polarizability of LiNbO3 and hence the optimized modal dispersion. To lowest
order the bandwidth of the phase-matched SHG process is inversely proportional
to the modal group velocity mismatch between first and second harmonic
frequencies.

$\displaystyle\Delta\omega\propto$
$\displaystyle\left|n_{g,2\omega}-n_{g,\omega}\right|^{-1}$ (19)

###  3.1 Single Frequency Optimization

We first attempted to maximize SHG bandwidth around target optical frequencies
by minimizing the magnitude of the modal group index mismatch in Eq. 19 at the
target frequency.

We defined an objective function equal to the square of the group index
difference to avoid discontinuities in the objective function derivative which
would be introduced by the absolute value

$\displaystyle g=$ $\displaystyle(n_{g,2\omega}-n_{g,\omega})^{2}.$ (20)

This objective function is chosen for simplicity and empirically works well.
We note that $\Delta\omega$ in Eq. 19 is limited by higher order dispersion
terms when $\left|n_{g,2\omega}-n_{g,\omega}\right|^{-1}$ approaches zero.
Direct optimization over a large bandwidth may be achieved by choosing a
different objective function, for example the sum of square group index
differences at multiple wavelengths or a combination of group index mismatch
and group velocity dispersion mismatch.

Figure 2:  Evolution of the cross section, fundamental wavelength mode, and
effective conversion envelope for 8 iterations at a fundamental wavelength of
2 $1.25\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and 2
$1.95\text{\,}\mathrm{\SIUnitSymbolMicro m}$. The waveguide was parameterized
by its 2 top width (having a $\sim$ $67\text{\,}\deg$ sidewall angle), 2
thickness, and 2 partial etch fraction. 0 is un-etched and 1 is fully-etched.
2 The associated objective function values $(n_{g,2\omega}-n_{g,\omega})^{2}$.
2 Required poling period to compensate the mismatch in effective index and 2
the resulting FWHM over which the SHG process is quasi phase matched for a
$1\text{\,}\mathrm{c}\mathrm{m}$ length.

The frequency dependence and anisotropy of lithium niobate’s dielectric
susceptibility are crucial factors in this design problem, highlighting the
novel capability of our implementation. Figure 22 and 22 show the evolution of
the waveguide cross section optimization for fundamental wavelengths of
$1.25\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and
$1.95\text{\,}\mathrm{\SIUnitSymbolMicro m}$. The corresponding evolution of
parameter and objective function values is shown in Fig. 22-2 and derived
values in Fig. 22-2. Beginning with randomly chosen values for the waveguide
geometry, in only 8 iterations the optimizer increases the SHG conversion
bandwidth by a factor of 3 (20) around
$1.25\text{\,}\mathrm{\SIUnitSymbolMicro m}$
($1.95\text{\,}\mathrm{\SIUnitSymbolMicro m}$) compared to bulk dispersion. We
estimate that comparable optimizations based on exhaustive search require
$\sim$100$\times$ as many eigenmode calculations for similar results[9, 28,
29], even in this low-dimensional parameter space. We expect this performance
gap to grow exponentially with number of design parameters, as discussed in
Sec. 4.

We performed similar optimizations of SHG bandwidth and modal GVD for many
target vacuum wavelengths spanning $1.25\text{\,}\mathrm{\SIUnitSymbolMicro
m}$–$2.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$ to establish the robustness
and generality of our approach for optimization of waveguide dispersion. These
results are presented in Supplement 1. We note that our SHG bandwidth
optimizations centered at $1.95\text{\,}\mathrm{\SIUnitSymbolMicro m}$ (Fig.
2(b)) and $2.03\text{\,}\mathrm{\SIUnitSymbolMicro m}$ (Fig. S1(g)) reach
$3\text{\,}\mathrm{dB}$ bandwidths of $78\text{\,}\mathrm{nm}$ and
$70\text{\,}\mathrm{nm}$ respectively for a $1\text{\,}\mathrm{cm}$
interaction length, which is comparable to the $110\text{\,}\mathrm{nm}$
bandwidth recently experimentally demonstrated in a $6\text{\,}\mathrm{mm}$
exhaustive-search-optimized TFLN waveguide [9] ($\sim$$67\text{\,}\mathrm{nm}$
for $1\text{\,}\mathrm{cm}$ length).

###  3.2 Broadband Optimization

Figure 3:  3 Normalized waveguide SHG transfer functions resulting from multi-
frequency waveguide geometry optimization for maximal quasi-phase matching
bandwidth around $1325\text{\,}\mathrm{nm}$ (orange),
$1350\text{\,}\mathrm{nm}$ (green), $1375\text{\,}\mathrm{nm}$ (pink) and
$1400\text{\,}\mathrm{nm}$ (blue) fundamental wavelengths. In each case the
modeled device length is $1\text{\,}\mathrm{cm}$. The minimax approach was
used to minimize the fundamental/second harmonic group velocity mismatch at
several optical frequencies as described in the text. The
$1350\text{\,}\mathrm{nm}$ SHG QPM spectrum in periodically poled bulk LiNbO3
is plotted in black to highlight the $\approx$100$\times$ bandwidth
improvement enabled by waveguide dispersion optimization. The TE00 modal
intensity distributions for fundamental and second harmonic frequencies are
shown in the inset along with a refractive index map showing the optimized
waveguide cross-section geometry. 3 SHG phase-mismatch $\Delta k$ spectra
corresponding to the optimized SHG transfer function spectra
$\propto\textrm{sinc}^{2}(\Delta kL)$ plotted in 3. 3 Parameter traces showing
the evolution of the optimized three waveguide design parameters over 12
optimization steps.

Single-frequency SHG bandwidth optimizations showed limited perforamce at
design wavelengths shorter than $1.6\text{\,}\mathrm{\SIUnitSymbolMicro m}$
(see Fig. S1 in Supplement 1). In this wavelength range higher order modal
dispersion terms tend to be large in magnitude and limit SHG bandwidths to
smaller values even when the fundamental and second-harmonic modes have
matched group velocities. To achieve larger SHG QPM bandwidths at shorter
wavelengths we re-cast the design problem as a minimax optimization[57, 58].
The minimax SHG bandwidth optimization attempts at each step to minimize the
largest magnitude group velocity mismatch computed at a discrete set of
optical frequencies $\omega_{\rm i}$. This approach has been demonstrated to
improve optimization results for other types of multi-frequency photonic
design problems [3, 59]. An epigraph formulation is used to achieve this with
a differentiable cost function. In this form the a dummy variable $t$ is
introduced and a set of differentiable inequality constraints $c_{\rm i}(x,t)$
are enforced as the cost function $g(x,t)$ is optimized w.r.t the design
parameters $x$.

$\begin{gathered}g(x,t)=t\\\ s.t.\>\>\>c_{\rm i}(x,t)=\left|n_{g,2\omega_{\rm
i}}-n_{g,\omega_{\rm i}}\right|-t\leq 0\end{gathered}$ (21)

We ran separate minimax SHG bandwidth optimizations using sets of three
regularly spaced fundamental optical frequencies with center fundamental
wavelengths of $1325\text{\,}\mathrm{nm}$, $1350\text{\,}\mathrm{nm}$ ,
$1375\text{\,}\mathrm{nm}$ and $1400\text{\,}\mathrm{nm}$. The results of
these minimax dispersion optimizations are shown in Fig. 3. In each case a
broadband phase matching condition is found in the allowed parameter range.
The optimized designs have SHG QPM full-width-half-max bandwidths ranging
approximately 40–$150\text{\,}\mathrm{nm}$ (6–$25\text{\,}\mathrm{THz}$) for
$1\text{\,}\mathrm{cm}$ interaction lengths, 30–100$\times$ larger than the
corresponding values for bulk PPLN.

Our optimized designs have SHG QPM bandwidths similar to or exceeding
previously published results [9, 28, 29] (6–$25\text{\,}\mathrm{THz}$ for
$1\text{\,}\mathrm{cm}$ interaction length) at shorter wavelengths where
material dispersion makes the design problem more challenging.

The benefit of the multi-frequency approach for broadband dispersion
engineering is highlighted in Fig. 22. In single-frequency optimizations shown
in Fig. 2(b) and Fig. S1 the GVM (the slope of $\Delta k$) at the target
wavelength was routinely driven to negligible levels near the precision limit
of the simulation ($|n_{g,2\omega}-n_{g,\omega}|\leq 10^{-5}$), but the
curvature of $\Delta k$ and higher-order dispersion mismatch terms were
uncontrolled. In contrast the minimax-optimized designs reach more moderate
GVM magnitudes at the target wavelength ($|n_{g,2\omega}-n_{g,\omega}|\approx
10^{-3}$) but the curvature and higher order dispersion of $\Delta k$ are
simultaneously minimized around the target wavelength, further enhancing the
simultaneous phase matching bandwidth. To our knowledge these designs
represent the shortest wavelength broadband SHG phase matching conditions
predicted in a TFLN platform and the first for the
1.3–$1.4\text{\,}\mathrm{\SIUnitSymbolMicro m}$ wavelength range, where the
commercial availability of tunable and swept laser sources [60, 61] as well as
fiber amplification [62, 63] make efficient broadband SHG highly desirable.

##  4 Performance and Accuracy

We performed a series of numerical experiments to study the performance and
accuracy of our differentiable mode solver. We defined a parametric waveguide
ridge geometry with a variable number of parameters $N_{\rm p}$ determining a
smooth surface morphology, as shown in the inset of Fig. 4. We then calculated
the group index of the highest effective-index waveguide mode and its gradient
w.r.t parameters with varied spatial grid resolutions (square grids
32$\times$32–512$\times$512 with fixed spatial extent) and number of geometry
parameters $3\leq N_{\rm p}\leq 100$. We used 10 randomly generated sets
parameter values for each combination of each grid size and number of
parameters and separately recorded the calculation times of each primal and
gradient calculation. To assess the gradient accuracy we compared gradients
calculated using the adjoint method to nominally equivalent finite difference
gradient approximations calculated using a 5-point central difference scheme
[64]. This comparison is shown in Fig. 4, where each point represents one
component of a computed gradient vector and points lying along $y=x$ indicate
agreement between adjoint-based and finite-difference gradients. Some
discrepancies are observed but the vast majority of the $\sim$10000 gradient
elements compared show good agreement, especially those with the largest
magnitude.

Figure 4:  4 Values of gradient magnitudes as calculated by finite difference
and via automatic differentiation. 4 In the present implementation the
gradient calculation dominates total simulation time. Scaling for the
underlying calculation is driven by DFT computation time for our plane wave
expansion-based solver. Run time is given as “wall time” on a personal
workstation with 32 GB RAM and a 4-core 3.2 GHz Intel processor. 4 Ratio of
time required to compute the gradients vs. the underlying variable (group
index or GVD).

The group index and gradient calculation run times (as wall times) are plotted
as a function of grid size in Fig. 4. The run time for the primal group index
calculation is dominated by the underlying eigenmode solve. The observed run-
time scaling with respect to grid size appears to be consistent with the
$\mathcal{O}(N\mathrm{log}(N))$ scaling (where $N$ is the total number of grid
points) seen with the software package [36] from which our solver is adapted.
In our implementation the gradient calculation is a factor of $\sim$50–100
slower than the primal for all spatial grid sizes and is roughly independent
of the number of design parameters $N_{\rm p}$. This is not fundamental to our
approach and we expect that improvements of the implementation will yield
group index/gradient run time ratios of $\sim$1.

Table 1. Comparison of Waveguide Dispersion Engineering Approaches
---
Work |  Optimization Approach | Mode Solves Per Iter. $[N]$ | Materials |  Geometry Definition
Anisotropy | Dispersion
Talenti, _et al_. [16] |  $k\cdot p$ reduced model | $\mathcal{O}(1)^{a}$ | No | No |  1-parameter PhC
Castelló, _et al_. [23] | Hellmann-Feynman | $\mathcal{O}(1)$ | No | No |  Limited Param.
Minkov, _et al_. [17] | Auto-Diff. | $\mathcal{O}(1)$ | No | No |  Limited Param.
Vial and Hao [21] | Auto-Diff. | $\mathcal{O}(1)$ |  Partialb | No |  Binary Topology
Hameed, _et al_. [24] | Trust Region Algo. | $\mathcal{O}(N^{3})$ | Yes | No |  Limited Param.
This Work | Auto-Diff. | $\mathcal{O}(1)$ | Yes | Yes |  Arbitrary Param.
Table 1:  Abbreviations: Automatic differentiation (“Auto-Diff.”), Trust-
Region Algorithm (“Trust Region Algo.”), Limited parameterization (“Limited
Param.”), Arbitrary Parameterization (“Arbitrary Param.”), Photonic Crystal
(“PhC”), Hellmann-Feynman Theorem (“Hellmann-Feynman”). “Limited
parameterization” geometries are, for example, the position and size of holes
in a PhC or the width and height of a rectangular waveguide. In contrast,
“arbitrary parametization” allows for any geometrical shape which can be
described by a finite number of input parameters. We describe “limited” and
“arbitrary” parametrizations based on what is reported in each work; these may
not be fundamental to the approach. Additionally, we consider inclusion of
“dispersion” to be explicit use of material information to determine spectral
derivatives without finite differences in frequency.
aOnly geometries parameterized in 1 variable are supported, then
$\mathcal{O}(1)$ 3D FDTD solves are required; additionally $\sim 100$s of
(fast) reduced model solves are performed.
bAnisotropic materials supported for uniaxial materials with z-oriented
anisotropy.

The ratios between corresponding group index and gradient calculation run
times are plotted as a function of the number of design parameters $N_{\rm p}$
in Fig. 4. We observe that the ratio between gradient and primal calculation
times approximately doubles for a 30-fold increase in $N_{\rm p}$. This
indicates that the gradient calculation run-time does not substantially depend
with $N_{\rm p}$, as expected for gradients computed using the adjoint method.

We computed independent eigen- and adjoint problem solutions in parallel for
objective functions depending on multiple waveguide mode solutions sharing a
common geometry (for example modes at several distinct optical frequencies).
Thus using a parallel computing environment we were able to scale the spectral
bandwidth or resolution of multi-frequency objective functions without
appreciably increasing simulation time.

Table 1 compares this work with several previous demonstrations of inverse
design for dispersion engineering. While many are based on similar frameworks,
none support the dispersive and anisotropic materials, limiting their
relevance to key nonlinear optical applications. Furthermore, our approach
back-propagates gradients w.r.t both the eigenvalues _and_ eigenvectors (mode
fields), a limitation of [23].

##  5 Conclusion

In summary we have demonstrated a differentiable electromagnetic eigenmode
solver suitable for a wide variety of modeling and design applications and
applied it to optimize waveguide geometries for maximal SHG phase matching
bandwidth. We demonstrated optimized designs for broadband SHG in TFLN
waveguides in the 1.3–$1.4\text{\,}\mathrm{\SIUnitSymbolMicro m}$ wavelength
range, enabling widely tunable visible light sources based on commercially
available infrared lasers. Our results for a practical waveguide dispersion
optimization with three parameters demonstrate a $\sim$100-fold reduction in
the number of required mode solutions compared to commonly used exhaustive
search approaches. To enable this we have derived the adjoint problem and
back-propagation equations for waveguide mode solutions accounting for
frequency-dependent and anisotropic dielectric materials and sub-pixel
smoothing. We believe this is the first mode solver capable of back-
propagating gradients with respect to mode fields and effective indices, group
indices and GVDs calculated from single eigenmode solutions. In this work we
focused on 2D waveguide cross-section models, however the solver and the
adjoint method described in our paper also applies to 3D eigenmode models.

This type of differentiable mode solver can be modularly combined with FDTD,
eigenmode expansion (EME) and multi-physics models to enable optimizable end-
to-end models of integrated photonic and optoelectronic components and
systems. Our work enables application of inverse design methods to a wide
variety of integrated photonic design problems based on waveguide mode models.

Funding This research was supported in part by an appointment to the
Intelligence Community Postdoctoral Research Fellowship Program at the
Massachusetts Institute of Technology administered by Oak Ridge Institute for
Science and Education (ORISE) through an interagency agreement between the
U.S. Department of Energy and the Office of the Director of National
Intelligence (ODNI), and in part by the DARPA DSO Robust Optical Clock Netowrk
(ROCkN) Program.

Acknowledgments The authors would like to thank Steven G. Johnson for helpful
discussions. The authors acknowledge the MIT SuperCloud and Lincoln Laboratory
Supercomputing Center for providing HPC resources that have contributed to the
research results reported within this paper.

Disclosures The authors declare no conflicts of interest.

Data Availability Statement The presented data is available upon reasonable
request.

Supplemental document See Supplement 1 for supporting content.

##  6 References

## References

  * [1] S. Molesky, Z. Lin, A. Y. Piggott, _et al._ , “Inverse design in nanophotonics,” Nature Photonics 12, 659–670 (2018).
  * [2] A. F. Oskooi, D. Roundy, M. Ibanescu, _et al._ , “MEEP: A flexible free-software package for electromagnetic simulations by the FDTD method,” Computer Physics Communications 181, 687–702 (2010).
  * [3] A. M. Hammond, A. F. Oskooi, M. Chen, _et al._ , “High-performance hybrid time/frequency-domain topology optimization for large-scale photonics inverse design,” Optics Express 30, 4467–4491 (2022).
  * [4] A. Y. Piggott, J. Petykiewicz, L. Su, and J. Vučković, “Fabrication-constrained nanophotonic inverse design,” Scientific reports 7, 1786 (2017).
  * [5] A. Y. Piggott, J. Lu, K. G. Lagoudakis, _et al._ , “Inverse design and demonstration of a compact and broadband on-chip wavelength demultiplexer,” Nature Photonics 9, 374–377 (2015).
  * [6] A. Michaels and E. Yablonovitch, “Inverse design of near unity efficiency perfectly vertical grating couplers,” Optics express 26, 4766–4779 (2018).
  * [7] H. Chung and O. D. Miller, “High-NA achromatic metalenses by inverse design,” Optics Express 28, 6945–6965 (2020).
  * [8] R. E. Christiansen, J. Michon, M. Benzaouia, _et al._ , “Inverse design of nanoparticles for enhanced raman scattering,” Optics Express 28, 4444–4462 (2020).
  * [9] M. Jankowski, C. Langrock, B. Desiatov, _et al._ , “Ultrabroadband nonlinear optics in nanophotonic periodically poled lithium niobate waveguides,” Optica 7, 40–46 (2020).
  * [10] Z. Jia, W. Qarony, J. Park, _et al._ , “Interpretable inverse-designed cavity for on-chip nonlinear photon pair generation,” Optica 10, 1529–1534 (2023).
  * [11] L. Ledezma, R. Sekine, Q. Guo, _et al._ , “Intense optical parametric amplification in dispersion-engineered nanophotonic lithium niobate waveguides,” Optica 9, 303–308 (2022).
  * [12] W. Ma, Z. Liu, Z. A. Kudyshev, _et al._ , “Deep learning for the design of photonic structures,” Nature Photonics 15, 77–90 (2021).
  * [13] G. Alagappan and C. Png, “Group refractive index via auto-differentiation and neural networks,” Scientific Reports 13, 4450 (2023).
  * [14] D. D. El-Mosalmy, M. Hameed, N. F. Areed, and S. Obayya, “Novel neural network based optimization approach for photonic devices,” Optical and Quantum Electronics 46, 439–453 (2014).
  * [15] G. N. Malheiros-Silveira and H. E. Hernandez-Figueroa, “Prediction of dispersion relation and pbgs in 2-d pcs by using artificial neural networks,” IEEE Photonics Technology Letters 24, 1799–1801 (2012).
  * [16] F. R. Talenti, S. Wabnitz, I. Ghorbel, _et al._ , “Fast dispersion tailoring of multimode photonic crystal resonators,” Physical Review A 106, 023505 (2022).
  * [17] M. Minkov, I. A. Williamson, L. C. Andreani, _et al._ , “Inverse design of photonic crystals through automatic differentiation,” Acs Photonics 7, 1729–1741 (2020).
  * [18] I. Nasidi, R. Hao, S. Jin, and E. Li, “Inverse design of a photonic moiré lattice waveguide towards improved slow light performances,” Applied Optics 62, 2651–2655 (2023).
  * [19] D. Vercruysse, N. V. Sapra, K. Y. Yang, and J. Vuckovic, “Inverse-designed photonic crystal circuits for optical beam steering,” ACS Photonics 8, 3085–3093 (2021).
  * [20] E. Nussbaum, E. Sauer, and S. Hughes, “Inverse design of broadband and lossless topological photonic crystal waveguide modes,” Optics Letters 46, 1732–1735 (2021).
  * [21] B. Vial and Y. Hao, “Open-source computational photonics with auto differentiable topology optimization,” Mathematics 10, 3912 (2022).
  * [22] H. Men, K. Y. Lee, R. M. Freund, _et al._ , “Robust topology optimization of three-dimensional photonic-crystal band-gap structures,” Optics express 22, 22632–22648 (2014).
  * [23] D. Castelló-Lurbe, E. Silvestre _et al._ , “Inverse dispersion engineering in silicon waveguides,” JOSA B 31, 1829–1835 (2014).
  * [24] M. F. O. Hameed, A.-k. S. Hassan, A. E. Elqenawy, and S. S. Obayya, “Modified trust region algorithm for dispersion optimization of photonic crystal fibers,” Journal of Lightwave Technology 35, 3810–3818 (2017).
  * [25] R. Stainko and O. Sigmund, “Tailoring dispersion properties of photonic crystal waveguides by topology optimization,” Waves in Random and Complex Media 17, 477–489 (2007).
  * [26] C. Xin, S. Lu, J. Yang, _et al._ , “Wavelength-accurate and wafer-scale process for nonlinear frequency mixers in thin-film lithium niobate,” arXiv preprint arXiv:2404.12381 (2024).
  * [27] P.-K. Chen, I. Briggs, C. Cui, _et al._ , “Adapted poling to break the nonlinear efficiency limit in nanophotonic lithium niobate waveguides,” Nature Nanotechnology 19, 44–50 (2024).
  * [28] H. R. Fergestad and K. Gallo, “Ultrabroadband second harmonic generation at telecom-wavelengths in lithium niobate waveguides,” in _Integrated Photonics Research, Silicon and Nanophotonics,_ (Optica Publishing Group, 2020), pp. ITu4A–13.
  * [29] C. Zhang, X. Tian, H. Liu, _et al._ , “Ultrabroadband second-harmonic generation based on a dispersion-engineered lithium niobate on insulator waveguide,” JOSA B 39, 3304–3310 (2022).
  * [30] S. G. Johnson, “The NLopt nonlinear-optimization package,” (2014).
  * [31] S. G. Johnson and contributors, “GeometryPrimitives.jl,” https://github.com/stevengj/GeometryPrimitives.jl (2022).
  * [32] S. Gowda, Y. Ma, A. Cheli, _et al._ , “High-performance symbolic-numerics via multiple dispatch,” arXiv preprint arXiv:2105.03949 (2021).
  * [33] A. F. Oskooi, C. Kottke, and S. G. Johnson, “Accurate finite-difference time-domain simulation of anisotropic media by subpixel smoothing,” Opt. Lett. 34, 2778–2780 (2009).
  * [34] C. Kottke, A. Farjadpour, and S. G. Johnson, “Perturbation theory for anisotropic dielectric interfaces, and application to subpixel smoothing of discretized numerical methods,” Physical Review E 77, 036611 (2008).
  * [35] Y. Zhou, C. Barnes, J. Lu, _et al._ , “On the continuity of rotation representations in neural networks,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR),_ (2019).
  * [36] S. G. Johnson and J. D. Joannopoulos, “Block-iterative frequency-domain methods for maxwell’s equations in a planewave basis,” Optics express 8, 173–190 (2001).
  * [37] M. Frigo and S. G. Johnson, “The design and implementation of FFTW3,” Proceedings of the IEEE 93, 216–231 (2005).
  * [38] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, “Julia: A fresh approach to numerical computing,” SIAM review 59, 65–98 (2017).
  * [39] “IterativeSolvers.jl,” https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl (2021).
  * [40] A. Montoison, D. Orban, and contributors, “Krylov.jl: A Julia basket of hand-picked Krylov methods,” https://github.com/JuliaSmoothOptimizers/Krylov.jl (2020).
  * [41] M. F. Herbst, A. Levitt, and E. Cancès, “DFTK: A Julian approach for simulating electrons in solids,” Proc. JuliaCon Conf. 3, 69 (2021).
  * [42] J. D. Jackson, _Classical electrodynamics_ (Wiley, 1998), 3rd ed.
  * [43] J. Nocedal and S. J. Wright, _Numerical optimization_ (Springer, 2006), 2nd ed.
  * [44] “ChainRules.jl,” https://github.com/JuliaDiff/ChainRules.jl (2021).
  * [45] S. G. Johnson, “Notes on adjoint methods for 18.335,” https://math.mit.edu/~stevenj/18.336/adjoint.pdf (2021).
  * [46] B. S. Elkus, K. Abdelsalam, A. Rao, _et al._ , “Generation of broadband correlated photon-pairs in short thin-film lithium-niobate waveguides,” Optics express 27, 38521–38531 (2019).
  * [47] J. Mishra, M. Jankowski, A. Y. Hwang, _et al._ , “Ultra-broadband mid-infrared generation in dispersion-engineered thin-film lithium niobate,” Optics Express 30, 32752–32760 (2022).
  * [48] H. Du, X. Zhang, L. Wang, and F. Chen, “Highly efficient, modal phase-matched second harmonic generation in a double-layered thin film lithium niobate waveguide,” Optics Express 31, 9713–9726 (2023).
  * [49] L. Ledezma, A. Roy, L. Costa, _et al._ , “Octave-spanning tunable infrared parametric oscillators in nanophotonics,” Science Advances 9, eadf9711 (2023).
  * [50] P. Aashna, H.-L. Lin, E. H. Dogheche, _et al._ , “Theoretical design for broadband parametric amplification in thin film lithium niobate,” in _2023 IEEE Photonics Conference (IPC),_ (IEEE, 2023), pp. 1–2.
  * [51] R. Sekine, R. M. Gray, L. Ledezma, _et al._ , “Multi-octave frequency comb from an ultra-low-threshold nanophotonic parametric oscillator,” arXiv preprint arXiv:2309.04545 (2023).
  * [52] A. Roy, L. Ledezma, L. Costa, _et al._ , “Visible-to-mid-ir tunable frequency comb in nanophotonics,” Nature Communications 14, 6549 (2023).
  * [53] J. Williams, R. Nehra, E. Sendonaris, _et al._ , “Ultra-short pulse biphoton source in lithium niobate nanophotonics at 2$\backslash$textmu m,” arXiv preprint arXiv:2402.05163 (2024).
  * [54] T.-H. Wu, L. Ledezma, C. Fredrick, _et al._ , “Visible-to-ultraviolet frequency comb generation in lithium niobate nanophotonic waveguides,” Nature Photonics pp. 1–6 (2024).
  * [55] C. Huang, Y. Luo, Y. Zhao, _et al._ , “A conformal mapping approach to broadband nonlinear optics on chip,” Nature Photonics pp. 1–7 (2024).
  * [56] M. M. Fejer, G. Magel, D. H. Jundt, and R. L. Byer, “Quasi-phase-matched second harmonic generation: tuning and tolerances,” IEEE Journal of quantum electronics 28, 2631–2654 (1992).
  * [57] K. Svanberg, “A class of globally convergent optimization methods based on conservative convex separable approximations,” SIAM Journal on Optimization 12, 555–573 (2002).
  * [58] S. P. Boyd and L. Vandenberghe, _Convex optimization_ (Cambridge university press, 2004).
  * [59] M. Chen, R. E. Christiansen, J. A. Fan, _et al._ , “Validation and characterization of algorithms and software for photonics inverse design,” JOSA B 41, A161–A176 (2024).
  * [60] “Santec TSL-570 tunable laser datasheet,” https://www.santec.com/dcms_media/other/instruments_TSL-570-C-E.pdf (2022).
  * [61] V. Jayaraman, G. Cole, M. Robertson, _et al._ , “High-sweep-rate 1310 nm MEMS-VCSEL with 150 nm continuous tuning range,” Electronics letters 48, 867–869 (2012).
  * [62] N. K. Thipparapu, A. Umnikov, P. Barua, and J. Sahu, “Bi-doped fiber amplifier with a flat gain of 25 dB operating in the wavelength band 1320–1360 nm,” Optics letters 41, 1518–1521 (2016).
  * [63] Y. Wang, N. K. Thipparapu, D. J. Richardson, and J. K. Sahu, “Ultra-broadband bismuth-doped fiber amplifier covering a 115-nm bandwidth in the O and E bands,” Journal of Lightwave Technology 39, 795–800 (2021).
  * [64] “FiniteDifferences.jl,” https://github.com/JuliaDiff/FiniteDifferences.jl (2022).