See discussions, stats, and author profiles for this publication at: [Link]

net/publication/329391388

Adjoint Method and Inverse Design for Nonlinear Nanophotonic Devices

Article in ACS Photonics · December 2018

DOI: 10.1021/acsphotonics.8b01522

CITATIONS READS
359 743

4 authors, including:

Tyler W. Hughes Momchil Minkov

Stanford University Flexcompute Inc
57 PUBLICATIONS 3,158 CITATIONS 143 PUBLICATIONS 6,049 CITATIONS

SEE PROFILE SEE PROFILE

All content following this page was uploaded by Tyler W. Hughes on 16 May 2019.

The user has requested enhancement of the downloaded file.

 Adjoint method and inverse design for nonlinear nanophotonic devices

Tyler W. Hughes,∗ Momchil Minkov,∗ Ian A. D. Williamson, and Shanhui Fan†

Department of Electrical Engineering, and Ginzton Laboratory, Stanford University, Stanford, CA 94305, USA
(Dated: November 6, 2018)
The development of inverse design, where computational optimization techniques are used to
design devices based on certain specifications, has led to the discovery of many compact, non-
intuitive structures with superior performance. Among various methods, large-scale, gradient-based
optimization techniques have been one of the most important ways to design a structure containing
a vast number of degrees of freedom. These techniques are made possible by the adjoint method,
in which the gradient of an objective function with respect to all design degrees of freedom can
be computed using only two full-field simulations. However, this approach has so far mostly been
applied to linear photonic devices. Here, we present an extension of this method to modeling
arXiv:1811.01255v1 [[Link]] 3 Nov 2018

nonlinear devices in the frequency domain, with the nonlinear response directly included in the
gradient computation. As illustrations, we use the method to devise compact photonic switches in a
Kerr nonlinear material, in which low-power and high-power pulses are routed in different directions.
Our technique may lead to the development of novel compact nonlinear photonic devices.

In recent years, there has been significant interest in Furthermore, because in many cases the steady-state be-
using computational optimization tools to design novel havior of the system is of interest, a frequency-domain
nanophotonic devices with a wide range of applications approach is preferred as the steady state response can
[1–18]. Much of this progress [1–12] is made possible by be obtained directly, without the need for going through
the adjoint method [18, 19], a technique which allows the a large number of time steps as in a time-domain simu-
gradient of an objective function to be computed with lation. The general mathematical formalism for the ad-
respect to an arbitrarily large number of degrees of free- joint method in nonlinear systems is known in the applied
dom using only two full-field simulations. This method mathematics literature [28]. But, with the exception of a
makes large-scale gradient-based design of electromag- very recent preprint that seeks to design a nonlinear ele-
netic structures possible. When compared to brute force ment in an optical neural network [23], such a formalism
searching through the parameter space [13], and other has not been previously applied to nonlinear photonic
commonly used design methods, like stochastic global op- device optimizations.
timization algorithms [14–17], gradient-based design has In this work we outline, in detail, how the adjoint
a number of practical advantages. For example, a very method may be used to optimize the steady-state re-
large number of design parameters can be adjusted simul- sponse of a nonlinear optical device in the frequency do-
taneously, and the number of structures one is required main. We first outline a formalism for generalizing ad-
to evaluate in order to reach a high-performing structure joint problems to arbitrary nonlinear problems. Then,
can be far smaller compared with the total number of as a demonstration, we use our method to inverse-design
structures in the search space. photonic switches with Kerr nonlinearity. Our results
Up to now, in photonics the adjoint method has been may be applied more generally to other objective func-
mostly applied to gradient-based optimization of linear tions and sources of nonlinearity and provides new pos-
optical devices. The generalization of the adjoint method sibilities for designing novel nonlinear optical devices.
to nonlinear optical devices would create new possibilities
in several exciting fields such as on-chip lasers [20], fre- NONLINEAR ADJOINT METHOD
quency combs [21], spectroscopy [22], neural computing
[23], and quantum information processing [24]. To this
We first outline the formulation of the adjoint method
end, several recent works [25–27] have applied adjoint
for the inverse design of nonlinear optical devices. The
methods to engineer linear devices to display favorable
goal of inverse design is to find a set of real-valued design
properties for nonlinear optical applications, such as high
variables ϕ that maximize a real-valued objective func-
quality factors, small mode volume, or large field overlap
tion L = L(e, e∗ , ϕ), where the complex-valued vector e
between the modes of interest. However, these works do
is given by the solution to the equation
not directly optimize the nonlinear systems.
f (e, e∗ , ϕ) = 0. (1)
To solve for the adjoint sensitivity of a nonlinear sys-
tem, the standard option is to work within a time-domain For example, eq. (1) may represent the steady-state
adjoint formalism, which entails simulating an additional Maxwell’s equations with an intensity-dependent permit-
linear system with a time-varying permittivity [3]. How- tivity distribution where e is the electric field distribu-
ever, as this formalism requires the storage of the fields at tion. The solution to eq. (1) may be found with any non-
each time step, it has substantial memory requirements. linear equation solver, such as with the Newton-Raphson
 2

method [29]. We further note that the treatment of e APPLICATION TO KERR NONLINEARITY
and its complex conjugate as independent variables is
necessary for differentiation as will be shown later. We now apply the general formalism as discussed above
The aim of the optimization is to maximize the ob- to the optimization of nonlinear optical systems. Since
jective function with respect to the design variables ϕ. the formalism is applicable to linear optical systems as
For this purpose, it is essential to compute the sensitivity well, for illustration purposes here we use it to treat
of L with respect to each element of ϕ. For simplicity, both the linear and the nonlinear cases, in order to high-
we derive the derivative of the objective function with light aspects that are unique to nonlinear systems. A
respect to a single parameter ϕ, which is written schematic outlining the two cases is presented in Fig. 1.
For a linear system, Maxwell’s equations for the steady-
dL ∂L ∂L de ∂L de∗ state at a frequency ω0 may be written as
= + + ∗ . (2)
dϕ ∂ϕ ∂e dϕ ∂e dϕ
µ−1 2
0 ∇ × ∇ × E(r) − ω0 0 r (r)E(r) = −iω0 J(r), (9)
Or, in matrix form as where E(r) is the electric field, J(r) is the electric current
  source, r (r) is the relative dielectric permittivity, and we
dL ∂L   de/dϕ
= + ∂L/∂e ∂L/∂e∗ . (3) have assumed relative permeability µr = 1 everywhere.
dϕ ∂ϕ de∗ /dϕ
Compactly, and to make connection to the general for-
malism in the previous section, this can be written in
To compute de/dϕ and de∗ /dϕ, we differentiate eq. (1):
matrix form as
df ∂f ∂f de ∂f de∗ f (e, e∗ , ϕ) = A(r )e − b = 0, (10)
=0= + + ∗ . (4)
dϕ ∂ϕ ∂e dϕ ∂e dϕ
where A is a linear operator, vectors e and r now con-
Eq. (4) together with its complex conjugate then yields tain the electric fields and the relative permittivity, re-
spectively, and b is a vector proportional to the current

∂f /∂e ∂f /∂e∗

de/dϕ
 
∂f /∂ϕ
 source. The design parameters ϕ in this case is the per-
= − . (5) mittivity distribution r . Eq. (10) can be solved to ob-
∂f ∗ /∂e ∂f ∗ /∂e∗ de∗ /dϕ ∂f ∗ /∂ϕ
tain the electric fields e, as diagrammed by Fig. 1(a).
Thus, formally we can rewrite eq. (3) as We assume an objective function L that depends on
the field solution to eq. (10) and we take the linear rel-
dL ∂L ative permittivity distribution as the set of design vari-
= − (6) ables. Because ∂f /∂e = A and ∂f /∂e∗ = 0 for the linear
dϕ ∂ϕ
system, from eq. (7), the adjoint field may be written
 ∂f /∂e ∂f /∂e∗ −1 ∂f /∂ϕ
   
∂L/∂e ∂L/∂e∗ simply as the solution to the equation

.
∂f ∗ /∂e ∂f ∗ /∂e∗ ∂f ∗ /∂ϕ
T
AT (r )eaj = − (∂L/∂e) , (11)
In analogy with the linear adjoint method, we can now as shown in Fig. 1(b). For a reciprocal system, AT = A,
compute the gradient by solving an additional linear sys- thus the original and the adjoint fields are solutions to
tem. We define a complex-valued adjoint field eaj as the the same linear problem but with different source terms.
solution to Note that the source for the adjoint field, − (∂L/∂e)
T

T   "
T
# depends on both the objective function and the original
∂f /∂e ∂f /∂e∗

eaj ∂L/∂e solution.
=− T , (7)
∂f ∗ /∂e ∂f ∗ /∂e∗ e∗aj ∂L/∂e∗ Once the adjoint field is computed, the gradient of the
objective function with respect to the permittivity dis-
and the gradient of the objective function is then tribution is given, through eq. (8), by
 
dL ∂L
  dL ∂L T ∂A
T ∂f = + 2R eaj e (12)
= + 2R eaj , (8) dr ∂r ∂r
dϕ ∂ϕ ∂ϕ
∂L
− 2ω02 0 R eTaj e .


= (13)
where R denotes taking the real part. In deriving eq. ∂r
(8), we have used the fact that both L and ϕ are real. In
Having reviewed the adjoint formalism for linear opti-
the case of multiple parameters ϕ, we can simply replace
cal systems we now consider nonlinear optical systems.
∂f /∂ϕ with the matrix ∂f /∂ϕ. Since eaj only needs to
As an example, we introduce Kerr nonlinearity into the
be solved once regardless of the number of parameters,
system [30], which corresponds to an intensity-dependent
gradients may be computed with very little marginal cost
permittivity
for an arbitrary number of free parameters, making large-
2
scale, gradient-based optimization possible. ˜r (r) = r (r) + 3ω02 0 χ(3) (r) |E(r)| , (14)
 3

to determine the derivative of the objective function, is a

linear problem. The size of the adjoint problem is twice
as large as the corresponding linear problem of Eq. (10),
but it is of a similar form, with the source dependent
upon the solution e.
Once the adjoint field is computed, the gradient
dL/dr is evaluated from eq. (13) as in the linear case.
Here for simplicity we do not assume any explicit depen-
dence of the nonlinearity on the design variable. How-
ever, the formalism is straightforward to extend to that
case, as explained in the Supplementary Information.

INVERSE DESIGN OF OPTICAL SWITCHES

Figure 1. Illustration of the adjoint field computation for a We now demonstrate the use of this nonlinear adjoint
linear and a nonlinear system. (a) The linear system driven formalism to inverse design optical switches with desired
by a point source b with an objective function L given by the power-dependent performance characteristics. In Figs. 2
field intensity at a measuring point. (b) The adjoint prob-
and 3, we show the optimization procedures and perfor-
lem for the linear system: the same system driven by a point
source given by −∂L/∂e located at the measuring point. (c) mance characteristics of a 1 → 1 and 1 → 2 port device,
The nonlinear system containing a medium with Kerr nonlin- respectively. The operating frequency for both devices
earity (red). The electric fields are the solution to a nonlinear correspond to a free-space wavelength of 2µm.
equation. (d) The adjoint problem for the nonlinear system, For each device, we seek to maximize the correspond-
which is a linear system of equations for the adjoint field and ing objective function with respect to the permittivity
its complex conjugate. The Kerr medium is replaced by a lin- distribution within a fixed design region. To perform
ear region whose permittivity depends on the nonlinear fields.
the numerical optimization of the structure, we use the
finite-difference frequency-domain method (FDFD) [31],
where χ(3) (r) is the nonlinear susceptibility distribution. where the fields and operators of Eq. (9) are spatially
Other types of nonlinear terms can also be treated with discretized using a Yee lattice [32]. For simplicity, we re-
the formalism outlined above. Replacing r (r) in Eq. (9) strict our study to two-dimensional structures (i.e. struc-
with ˜r (r) in Eq. (14), our system is then described by tures with infinite extent in the third dimension), and
the equation: transverse-magnetic polarization, which has only non-
zero out-of-plane electric field components. In the opti-
f (e, e∗ , ϕ) = Anl (r , χ, e)e − b = 0, (15) mization process, we start with an initial relative permit-

 tivity in the design region. We solve the electric field dis-
where Anl ≡ A(r ) − diag χ |e|2 . Here, is
tribution in the structure by solving the nonlinear equa-
element-wise vector multiplication and diag(v) repre-
tion (eq. (15)). Then, we compute the gradient of L
sents a diagonal matrix with vector v on the main diag-
with respect to the relative permittivity distribution in
onal. The vector χ corresponds to the term 3ω02 0 χ(3) (r)
this region using eq. (13). With the gradient informa-
and |e|2 ≡ e e∗ . Again, for concreteness, the design
tion, we perform updates of the design variables using the
parameters ϕ correspond to the permittivity r . The
limited-memory BFGS [33] algorithm, although a simple
solution to this problem is diagrammed in Fig. 1(c).
gradient ascent algorithm would also suffice. This proce-
From eq. (7) we may now compute the partial deriva-
dure is repeated until convergence on a final structure.
tives of f with respect to the electric fields e, which is
We choose optimization parameters corresponding to
needed to construct the adjoint problem.
a device made from chalcogenide glass (Al2 S3 ), which ex-
∂f /∂e = A − 2diag χ |e|2 hibits a strong χ(3) response and high damage threshold


(16)
∂f ∗ /∂e = −diag (χ e∗ e∗ ) . (17) [30, 34, 35]. During the optimization, the relative per-
mittivity was constrained to lie between 1 (air) and 5.95
With this, we then express the adjoint field as a solution (Al2 S3 ). We further assume that the materials exhibit
to the linear system nonlinearity only within the design regions outlined in
T T T Fig. 2(a) and 3(a).
(∂f /∂e) eaj + (∂f ∗ /∂e) e∗aj = − (∂L/∂e) , (18) To create a more realistic final structure, the strength
which is diagrammed in Fig. 1(d). of the nonlinear susceptibility was assumed to be propor-
For a nonlinear system, to obtain the field e, one will tional to the ”density” of material, ρ, defined as
need to solve a nonlinear equation (e.g. eq. (15)). How- r (r) − 1
ever, we emphasize that the adjoint problem, as required ρ(r) = , (19)
m − 1
 4

1.0
(a) (b)

transmission
design region linear
0.5
X
nonlinear
5 m
0.0
10 4 10 2 100 102
input power (W / m)
|Ez|/Pin1/2 |Ez|/Pin1/2
(c) low power (d) high power
102 102

5 m 101 5 m 101

Figure 2. Inverse design demonstration of a 1 → 1 port switch. (a) Optical power is input into the left port (purple
arrow). The goal of optimizing the design region (blue square) is to maximize power transmission in the linear regime (blue
arrow) and minimize transmission in the nonlinear regime (red arrow). The final permittivity distribution after optimizing is
also shown. The black regions are chalcogenide with a relative permittivity of 5.95 and a χ(3) of 4.1 × 10−19 m2 /V2 . The
waveguide regions outside the design region have a width of 0.3µm. The operating wavelength is 2µm. (b) The transmission
as a function of input power, demonstrating the switching behavior at around 10−3 W/µm. The dashed black line indicates
the input power used for the high-power regime in the optimization. (c-d) The amplitude of the simulated electric field of the
final structure, in the linear (c) and nonlinear (d) regimes, respectively – with input power of 10−9 W/µm and 0.157 W/µm,
respectively. Ez corresponds to the out-of-plane electric field in the 2D simulation.

where m is the permittivity of the material. This as- its maximum value is 1. The optimization setup and the
sumption ensures that air regions do not exhibit a non- optimized structure are diagrammed in Fig. 2(a). The
linear refractive index. Eq. (19) adds an r dependence final structure resembles a resonator between two Bragg
in the nonlinear susceptibility, which is straightforwardly mirrors, effectively acting like a bistable switch [37, 38].
treated in the adjoint method, as discussed in the Supple- Fig. 2(b) shows the transmission as a function of the in-
mentary Information. Low-pass spatial filtering and pro- put power, and it clearly switches from high to low as the
jection techniques [36] were applied during optimization input power increases. This is also illustrated in panels
to create binarized (air and chalcogenide) final structures (c)-(d), where we plot the field amplitude distributions
with large, smoothed features. Additional details on this in the low-power (high-transmission) regime and in the
are described in the Supplementary Information. high-power (low-transmission) regime, respectively. The
Our first device, as shown in Fig. 2, consists of a computed power transmission coefficients for these two
waveguide-fed 1 → 1 port geometry with a central de- panels are 98.2% and 3.1%, respectively. The value of the
sign region. We optimize this design region to maximize input power used in the optimization and in panel (d) is
power transmission in the linear regime when the incident shown by a dashed line in panel (b). At this input power,
power is sufficiently low such that the nonlinear terms do the device exhibits a maximum nonlinear refractive index
not affect the transmission, and minimize transmission shift of 4.0 × 10−3 , which is below the damage thresh-
in the nonlinear regime when the incident power is at a old for Al2 S3 using sub-nanosecond pulses [39] (see Sup-
specific high value such that the nonlinearity plays a sig- plementary Information). The transmission spectrum of
nificant role. This corresponds to an objective function this structure gives a resonance peak with a full-width
of the form at half maximum of 38GHz (see Supplementary Infor-
2 2 mation). In the Supplementary Information, we also list
L(elow , ehigh ) = mT elow − mT ehigh , (20) the specific optimization parameters, and show the value
of the objective function during the optimization process.
where elow and ehigh are the simulated fields with a low
Reaching the final optimized structure shown in Fig. 2
and a high input power, respectively, m is the modal pro-
required the evaluation of 2000 structures, but a reason-
file of the electric field for the waveguide in the output
ably high-performing structure is already reached after
port, and the objective function is normalized such that
 5

1.0
(a) (b)
design 0.8
region

transmission
0.6
right port
linear bottom port
0.4

nonlinear 0.2
5 m
0.0
10 5 10 4 10 3 10 2 10 1
input power (W / m)

|Ez|/Pin1/2 |Ez|/Pin1/2
(c) low power (d) high power

102 102

5 m 101 5 m 101

Figure 3. Inverse design demonstration of a 1 → 2 port switch. (a) Optical power is input into the left port (purple
arrow). The goal of optimizing the design region (blue square) is to maximize the power transmission to the right port (blue
arrow) in the linear regime and maximize transmission to the bottom port (red arrow) in the nonlinear regime. The final
permittivity distribution after optimizing is also shown. (b) The transmission in the right (blue) and bottom (red) ports as
a function of input power, demonstrating the switching behavior at around 2 × 10−2 W/µm. The dashed black line indicates
the input power used for the high-power regime in the optimization. (c-d) The amplitude of the simulated electric field of the
final structure, in the linear (c) and nonlinear (d) regimes, respectively – with input power of 10−9 W/µm and 0.057 W/µm,
respectively.

only a few hundred iterations. maximum value of 1 for a perfect switching operation.
We also use the same technique for the inverse design The setup of the optimization problem and the final de-
of a 1 → 2 port switch where light is guided to the right sign are diagrammed in Fig. 3(a). We note that the
port in the linear regime and to the bottom port in the device displays a non-intuitive geometry while retaining
nonlinear regime. To achieve this design, we define the large features and good binarization.
objective function as In Fig. 3(b), we plot the transmission through the
2 2 right and through the bottom ports as a function of in-
L(elow , ehigh ) = mTr elow − mTr ehigh (21) put power. This clearly shows the switching of power
2 2
− mTb elow + mTb ehigh , from the right port to the bottom port in the linear and
nonlinear regimes, respectively. Specifically, in the linear
where mr and mb denote the mode profiles of the waveg- regime, the device has a power transmission of 81.8% and
uides in the right and bottom ports, respectively. These 5.9% to the right and bottom ports, respectively, while
are normalized such that the objective function has a in the nonlinear regime, at the input power marked by
 6

the dashed line in Fig. 3(b), these values are 6.1% and optimization applied to electromagnetic design,” Optics
80.8%, respectively. The electric field amplitudes for lin- Express 21, 21693–21701 (2013).
ear and nonlinear regimes are displayed in 3(c)-(d). The [2] Jiahui Wang, Yu Shi, Tyler Hughes, Zhexin Zhao,
operational bandwidth for this device is approximately and Shanhui Fan, “Adjoint-based optimization of active
nanophotonic devices,” Optics Express 26, 3236–3248
90 GHz (see Supplementary Information). The device (2018).
exhibits a maximum nonlinear refractive index shift of [3] Y. Elesin, B. S. Lazarov, J. S. Jensen, and O. Sigmund,
3.9 × 10−3 , which is also below the acceptable damage “Design of robust and efficient photonic switches using
threshold for Al2 S3 using sub-nanosecond pulses [39]. A topology optimization,” Photonics and Nanostructures -
full list of optimization parameters and a plot of the ob- Fundamentals and Applications 10, 153–165 (2012).
jective function vs. iteration number is shown in the [4] Alexander Y. Piggott, Jan Petykiewicz, Logan Su, and
Supplementary Information. Jelena Vučković, “Fabrication-constrained nanophotonic
inverse design,” Scientific Reports 7, 1786 (2017).
[5] Alexander Y. Piggott, Jesse Lu, Konstantinos G.
Lagoudakis, Jan Petykiewicz, Thomas M. Babinec, and
DISCUSSION
Jelena Vučković, “Inverse design and demonstration of
a compact and broadband on-chip wavelength demulti-
We have presented an extension to the adjoint vari- plexer,” Nature Photonics 9, 374–377 (2015).
able method applied to the optimization of an electro- [6] C. Y. Kao, S. Osher, and E. Yablonovitch, “Maximiz-
magnetic system with Kerr nonlinearity. Our approach ing band gaps in two-dimensional photonic crystals by
using level set methods,” Applied Physics B 81, 235–244
can be straightforwardly applied to other types of non-
(2005).
linearities which do not mix frequencies, such as sat- [7] Tyler Hughes, Georgios Veronis, Kent P. Wootton,
urable gain or absorption. Moreover, the methods here R. Joel England, and Shanhui Fan, “Method for com-
should be straightforwardly generalizable to treat nonlin- putationally efficient design of dielectric laser accelerator
ear problems involving frequency mixing. For example, structures,” Optics Express 25, 15414–15427 (2017).
one can imagine implementing a similar adjoint method [8] Sean Molesky, Zin Lin, Alexander Y. Piggott, Weiliang
in combination with the multi-frequency finite-difference Jin, Jelena Vučković, and Alejandro W. Rodriguez,
“Outlook for inverse design in nanophotonics,” arXiv
frequency-domain implementations for nonlinear wave in-
preprint arXiv:1801.06715 (2018), arXiv: 1801.06715.
teractions [40]. [9] O. Sigmund and J. Sondergaard Jensen, “Systematic
In addition to the design of optical switches, our for- design of phononic band-gap materials and structures
malism may prove useful for many other interesting prob- by topology optimization,” Philosophical Transactions of
lems in nonlinear photonics. For example, one could ap- the Royal Society A: Mathematical, Physical and Engi-
ply our approach to design nonlinear elements in opti- neering Sciences 361, 1001–1019 (2003).
cal neural networks [41] with specific forms of activation [10] Ren Matzen, Jakob S. Jensen, and Ole Sigmund, “Sys-
tematic design of slow-light photonic waveguides,” Jour-
functions. Another interesting application is power reg-
nal of the Optical Society of America B 28, 2374–2382
ulation in photonic networks. For example, as photonic (2011).
networks for laser-driven particle accelerators [42] must [11] Jakob S. Jensen and Ole Sigmund, “Topology optimiza-
be able to handle large input powers, it may be of interest tion of photonic crystal structures: a high-bandwidth
to use our approach to design compact optical limiters in low-loss T-junction waveguide,” Journal of the Optical
these networks. For the purposes of exploring these and Society of America B 22, 1191 (2005).
many other potential applications, we have made pub- [12] Louise F. Frellsen, Yunhong Ding, Ole Sigmund, and
Lars H. Frandsen, “Topology optimized mode multiplex-
licly available a software package that implements the
ing in silicon-on-insulator photonic wire waveguides,”
algorithms discussed here [43]. Optics Express 24, 16866–16873 (2016).
To summarize this paper, we have developed an adjoint [13] Bing Shen, Peng Wang, Randy Polson, and Ra-
method, which enables gradient optimization of nonlinear jesh Menon, “An integrated-nanophotonics polarization
photonic devices. Our work broadens the capability of beamsplitter with 2.4 × 2.4 µm2 footprint,” Nature Pho-
inverse design for producing novel nonlinear devices. tonics 9, 378–382 (2015).
This work is supported by the Gordon and Betty [14] Linfang Shen, Zhuo Ye, and Sailing He, “Design of two-
dimensional photonic crystals with large absolute band
Moore Foundation (GBMF4744); the Swiss National Sci-
gaps using a genetic algorithm,” Physical Review B 68,
ence Foundation (P300P2 177721); and the Air Force Of- 035109 (2003).
fice of Scientific Research (AFOSR) (FA9550-17-1-0002). [15] Momchil Minkov and Vincenzo Savona, “Automated op-
timization of photonic crystal slab cavities.” Scientific re-
ports 4, 5124 (2014).
[16] Momchil Minkov and Vincenzo Savona, “Wide-band slow
light in compact photonic crystal coupled-cavity waveg-
∗
These authors contributed equally. uides,” Optica 2, 631–634 (2015).
†
shanhui@[Link] [17] Yu Shi, Wei Li, Aaswath Raman, and Shanhui Fan, “Op-
[1] Christopher M. Lalau-Keraly, Samarth Bhargava, timization of Multilayer Optical Films with a Memetic
Owen D. Miller, and Eli Yablonovitch, “Adjoint shape Algorithm and Mixed Integer Programming,” ACS Pho-
 7

tonics 5, 684–691 (2018). problems involving maxwell’s equations in isotropic me-

[18] Georgios Veronis, Robert W. Dutton, and Shanhui Fan, dia,” IEEE Transactions on antennas and propagation
“Method for sensitivity analysis of photonic crystal de- 14, 302–307 (1966).
vices,” Optics Letters 29, 2288–2290 (2004). [33] Richard H Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou
[19] Michael B Giles and Niles A Pierce, “An introduction to Zhu, “A limited memory algorithm for bound constrained
the adjoint approach to design,” Flow, turbulence and optimization,” SIAM Journal on Scientific Computing
combustion 65, 393–415 (2000). 16, 1190–1208 (1995).
[20] Daiki Yamashita, Yasushi Takahashi, Takashi Asano, [34] Richard T. White and Tanya M. Monro, “Cascaded Ra-
and Susumu Noda, “Raman shift and strain effect in man shifting of high-peak-power nanosecond pulses in
high-Q photonic crystal silicon nanocavity,” Optics Ex- as2 s3 and as2 se3 optical fibers,” Optics Letters 36, 2351
press 23, 3951 (2015). (2011).
[21] Yoshitomo Okawachi, Kasturi Saha, Jacob S. Levy, [35] Michael R. Lamont, Barry Luther-Davies, Duk-Yong
Y. Henry Wen, Michal Lipson, and Alexander L. Gaeta, Choi, Steve Madden, and Benjamin J. Eggleton, “Su-
“Octave-spanning frequency comb generation in a silicon percontinuum generation in dispersion engineered highly
nitride chip,” Optics Letters 36, 3398 (2011). nonlinear (γ = 10 w/m) as2 s3 chalcogenide planar waveg-
[22] Joong Ho Moon, Jin Ho Kim, Ki-jeong Kim, Tai- uide,” Optics Express 16, 14938 (2008).
Hee Kang, Bongsoo Kim, Chan-Ho Kim, Jong Hoon [36] Mingdong Zhou, Boyan S Lazarov, Fengwen Wang, and
Hahn, and Joon Won Park, “Absolute Surface Density Ole Sigmund, “Minimum length scale in topology op-
of the Amine Group of the Aminosilylated Thin Layers: timization by geometric constraints,” Computer Meth-
Ultraviolet-Visible Spectroscopy, Second Harmonic Gen- ods in Applied Mechanics and Engineering 293, 266–282
eration, and Synchrotron-Radiation Photoelectron Spec- (2015).
troscopy Study,” Langmuir 13, 4305–4310 (1997). [37] Marin Soljačić, Mihai Ibanescu, Steven G. Johnson, Yoel
[23] Erfan Khoram, Ang Chen, Dianjing Liu, Qiqi Wang, Fink, and J. D. Joannopoulos, “Optimal bistable switch-
and Zongfu Yu, “Stochastic optimization of nonlinear ing in nonlinear photonic crystals,” Phys. Rev. E 66,
nanophotonic media for artificial neural inference,” arXiv 055601 (2002).
preprint arXiv:1810.07815 (2018). [38] Mehmet Fatih Yanik, Shanhui Fan, Marin Soljačić, and
[24] Xiang Guo, Chang-Ling Zou, Hojoong Jung, and J. D. Joannopoulos, “All-optical transistor action with
Hong X. Tang, “On-Chip Strong Coupling and Effi- bistable switching in a photonic crystal cross-waveguide
cient Frequency Conversion between Telecom and Visi- geometry,” Optics Letters 28, 2506 (2003).
ble Optical Modes,” Physical Review Letters 117 (2016), [39] Marine Chorel, Thomas Lanternier, Eric Lavastre, Nico-
10.1103/PhysRevLett.117.123902. las Bonod, Bruno Bousquet, and Jérôme Néauport, “Ro-
[25] Zin Lin, Xiangdong Liang, Marko Lončar, Steven G. bust optimization of the laser induced damage threshold
Johnson, and Alejandro W. Rodriguez, “Cavity- of dielectric mirrors for high power lasers,” Optics express
enhanced second-harmonic generation via nonlinear- 26, 11764–11774 (2018).
overlap optimization,” Optica 3, 233 (2016). [40] Yu Shi, Wonseok Shin, and Shanhui Fan, “Multi-
[26] Zin Lin, Marko Lončar, and Alejandro W. Rodriguez, frequency finite-difference frequency-domain algorithm
“Topology optimization of multi-track ring resonators for active nanophotonic device simulations,” Optica 3,
and 2d microcavities for nonlinear frequency conversion,” 1256–1259 (2016).
Optics Letters 42, 2818 (2017). [41] Yichen Shen, Nicholas C Harris, Scott Skirlo, Mihika
[27] Jorge Bravo-Abad, Alejandro Rodriguez, Peter Bermel, Prabhu, Tom Baehr-Jones, Michael Hochberg, Xin Sun,
Steven G. Johnson, John D. Joannopoulos, and Marin Shijie Zhao, Hugo Larochelle, Dirk Englund, and Marin
Soljačić, “Enhanced nonlinear optics in photonic-crystal Soljačić, “Deep learning with coherent nanophotonic cir-
microcavities,” Optics Express 15, 16161–16176 (2007). cuits,” Nature Photonics 11, 441 (2017).
[28] Gilbert Strang, Computational Science and Engineering [42] Tyler W. Hughes, Si Tan, Zhexin Zhao, Neil V. Sapra,
(Wellesley-Cambridge Press, 2007) chapter 8. Kenneth J. Leedle, Huiyang Deng, Yu Miao, Dylan S.
[29] William H Press, Saul A Teukolsky, William T Vetter- Black, Olav Solgaard, James S. Harris, Jelena Vuckovic,
ling, and Brian P Flannery, Numerical recipes 3rd edi- Robert L. Byer, Shanhui Fan, R. Joel England, Yun Jo
tion: The art of scientific computing (Cambridge univer- Lee, and Minghao Qi, “On-chip laser-power delivery sys-
sity press, 2007). tem for dielectric laser accelerators,” Physical Review
[30] Robert W. Boyd, Nonlinear Optics (Academic Press, Applied 9, 054017 (2018).
2008). [43] Tyler W Hughes, Momchil Minkov, and Ian
[31] Wonseok Shin and Shanhui Fan, “Choice of the perfectly A D Williamson, “Angler – an adjoint nonliner
matched layer boundary condition for frequency-domain gradient open-source package,” [Link]
maxwell’s equations solvers,” Journal of Computational fancompute/angler (2018).
Physics 231, 3406–3431 (2012).
[32] Kane Yee, “Numerical solution of initial boundary value
 S1

SUPPLEMENTARY INFORMATION

OPTIMIZATION DETAILS

Table S1 contains the parameters used in the inverse design demonstrations of Fig. 2 and Fig. 3 of the main text.
The values of the objective function vs. iteration are shown in Fig. S1 for both the 2-port and 3-port devices.

parameter symbol value (2-port) value (3-port) units

max relative permittivity m 5.95 5.95 -
nonlinear susceptibility χ(3) 4.1 × 10−19 4.1 × 10−19 m2 V−2
input power P0 157 57 mW/µm
free space wavelength λ0 2 2 µm
design region length L 10 6 µm
design region height H 1.6 6 µm
waveguide width w 300 300 nm
grid size g 40 40 nm
low-pass filter feature size R 160 200 nm
projection strength β 100 500 -
projection mid-point η 0.5 0.5 -

Table S1. Parameters used in the optimization study. Column ‘2-port’ refers to the device from Fig. 2. Column ‘3-port’ refers
to the device from Fig. 3
.

1.0 1.0
(a) (b)
0.8 0.8
objective function

objective function

0.6 0.6
0.4 0.4
0.2 0.2
0.0 0.0
0 500 1000 1500 2000 0 250 500 750 1000 1250 1500
iteration iteration

Figure S1. Objective function vs. iteration of the optimization for (a) 2-port device of Fig. 2 and (b) 3-port device of Fig. 3
of the main text.

PERMITTIVITY-DEPENDENT NONLINEAR SUSCEPTIBILITY

In the inverse design demonstration of the main text, we made the assumption that the nonlinear susceptibility
distribution was proportional to the density of material in the design region. Here we will derive the form of the
adjoint sensitivity with this modification.
With our assumption, the nonlinear susceptibility vector χ may be written in terms of the scalar magnitude of the
nonlinear susceptibility χ(3) and the relative permittivity vector r explicitly as
r − 1
χ = 3ω02 0 χ(3) (S1)
m − 1
where 1 is a vector of all ones and m is the maximum relative permittivity allowed in the optimization, corresponding
to the material relative permittivity.
 S2

When choosing to relative permittivity distribution as the set of design variables, ϕ = r , the nonlinear adjoint
problem requires the calculation of the partial derivatives ∂f /∂e, ∂f ∗ /∂e, and ∂f /∂r . When choosing the form of χ
from eq. (S1), the partial derivatives ∂f /∂e and ∂f ∗ /∂e are the same as derived in the main text. However, the term
∂f /∂r takes on a more complicated form given by
 
∂f ∂ r − 1
= A(r )e − 3ω02 0 χ(3) e e e∗ − b (S2)
∂r ∂r m − 1
∂A 1
e − 3ω02 0 χ(3) diag e |e|2


= (S3)
∂r m − 1
1
= diag (e) − 3ω02 0 χ(3) diag e |e|2 ,


(S4)
m − 1
 
∂A
where we make use of the fact that ∂ r
= δij δjk , where δ is the Kronecker delta.
ijk
Thus, while the adjoint field eaj will have the same form as in the main text, when computing the sensitivity as
in eq. (8) in the main text, one must insert the form of ∂f /∂r from eq. (S4). This is in contrast to the usual case
where the nonlinear susceptibility is fixed, ∂f /∂r = diag (e).

MAINTAINING MINIMUM FEATURE SIZE AND BINARIZATION

To create realistic devices with sufficiently large minimum feature sizes and binarized permittivity distributions,
we employed filtering and projection schemes during our optimization. These schemes are discussed in great detail in
[36] and related works.
Rather than updating the permittivity distribution directly, one may instead choose to update a design density ρ,
which varies between 0 and 1 within the design region. To create a structure with larger feature sizes, a low pass filter
can be applied to ρ to created a filtered density, labelled ρ̃:
P
j∈D Wij ρj
ρ̃i = P , (S5)
j∈D Wij

where D denotes the design region, and W is the spatial filter, defined for a feature size of R as
(
R − |ri − rj | if |ri − rj | ≤ R
Wij = (S6)
0, otherwise

with |ri − rj | being the distance between points i and j. This defines a low-pass spatial filter on ρ with the effect of
smoothing out features with length scale below R.
Now, for binarization of the structure, a projection scheme is used to recreate the final permittivity from the filtered
density. For this we define ρ̄ as the projected density, which is created from ρ̃ as

tanh (βη) + tanh (β [ρ̃i − η])

ρ̄i = . (S7)
tanh (βη) + tanh (β [1 − η])

Here, η is a parameter between 0 and 1 that controls the mid-point of the projection, typically 0.5, and β controls
the strength of the projection, typically around 100.
The relative permittivity can then be determined from ρ̄ as

r = (m − 1)ρ̄ + 1, (S8)

where m is the maximum permittivity.

The effect of these filtering and the binarization techniques on a sample permittivity set is illustrated in Fig. S2.
In the optimizations of the main text, these techniques were performed only within the design region and required
minimal modifications to the adjoint sensitivity. The determination of ∂r /∂ ρ̄, ∂ ρ̄/∂ ρ̃, and ∂ ρ̃/∂ρ were required to
compute the derivatives of the objective function with respect to the underlying ρ. For more details, see the software
package accompanying this work [43].
 S3

design density ( ) filtered density ( )

0 1.0 0
0.8
20 0.8 20

0.6 0.6
40 40
pixels (y)

pixels (y)
60 0.4 60 0.4

80 0.2 80 0.2

0.0 0.0
0 20 40 60 80 0 20 40 60 80
pixels (x) pixels (x)
projected density ( ) rel. permittivity ( r)
0 1.0 0 6

20 0.8 20 5

40 0.6 40 4
pixels (y)

pixels (y)
60 0.4 60 3

80 0.2 80 2

0.0 1
0 20 40 60 80 0 20 40 60 80
pixels (x) pixels (x)

Figure S2. Filtering and projection of an example design density, ρ. (top left) the original density ρ before processing. (top
right) the density after applying a low pass filter, ρ̃. (bottom left) the density ρ̄ after applying projection. (bottom right) the
final relative permittivity distribution r . The parameters used are R = 200nm, β = 100, η = 0.5.

NONLINEAR INDEX SHIFT

Here we estimate the maximum nonlinear index shift of chalcogenide (Al2 S3 ) materials. Based on [30, 34, 35], Al2 S3
has a nonlinear index (n2 ) between 3 × 10−14 and 2 × 10−13 cm2 /W. From [39], the damage threshold is 2.5 J/cm2 .
At a pulse duration of 100 ps, this damage threshold corresponds to 2.5 × 1010 W/cm2 . Together with the nonlinear
refractive index, the maximum refractive index shift sustainable by the material is approximately

∆n = n2 Idamage ≈ 7.5 × 10−4 − 5 × 10−3 (S9)

with a corresponding pulse bandwidth (for a bandwidth-limited Gaussian pulse) of ≈ 4.5 GHz. Our final structures
exhibit maximum refractive index shifts below 5 × 10−3 and their objective functions have FWHM bandwidths above
10 GHz. This suggests that they should exhibit their desired switching effects without damage using pulse durations
on the order of 100 ps and input powers on the order of 100 mW/µm.

TRANSMISSION SPECTRA

The transmission vs. frequency for the two devices in the main text is shown in Fig. S3 in the linear regime.
 S4

Figure S3. (a) Transmission spectrum through the 2-port device for the low power (linear) regime. (b) Transmission spectrum
through the right (blue) and bottom (red) ports of the 3-port device in the low power (linear) regime. The x-axis represents
the difference in frequency with respect to the design frequency.

View publication stats