Fabrication-constrained nanophotonic inverse design

A major difficulty in applying computational design methods to nanophotonic devices is ensuring that the resulting designs are fabricable. Here, we describe a general inverse design algorithm for nanophotonic devices that directly incorporates fabrication constraints. To demonstrate the capabilities of our method, we designed a spatial-mode demultiplexer, wavelength demultiplexer, and directional coupler. We also designed and experimentally demonstrated a compact, broadband 1 × 3 power splitter on a silicon photonics platform. The splitter has a footprint of only 3.8 × 2.5 μm, and is well within the design rules of a typical silicon photonics process, with a minimum radius of curvature of 100 nm. Averaged over the designed wavelength range of 1400–1700 nm, our splitter has a measured insertion loss of 0.642 ± 0.057 dB and power uniformity of 0.641 ± 0.054 dB.

where is the permittivity distribution, and µ 0 is the permeability of free space. For each input mode i, we then specify a set of output modes j = 1 . . . N i , whose amplitudes are bounded between α ij and β ij . If our output modes are guided modes of waveguides with modal electric fields E ij and magnetic fields H ij , this constraint can be written using a mode orthogonality relationship [2], Here,n is a unit vector pointing in the propagation direction, and r ⊥ denotes the coordinates perpendicular to the propagation direction. We can use Faraday's law ∇ × E i = −iωµ 0 H i to rewrite (S.2) purely in terms of the electric field: More generally, we can specify the output mode amplitude in terms of a linear functional L ij of the electric field E i , } is the space of all possible electric field distributions, and L ij : V → C maps electric field distributions to a complex scalar. We are thus interested in finding a permittivity distribution and electric fields E i which simultaneously satisfy (S.1) and (S.4), for all input modes i = 1 . . . M and output modes j = 1 . . . N . To ensure that the resulting device is fabricable, we will later impose additional constraints on .

Linear algebra description
Since we will solve Maxwell's equations numerically, and employ numerical optimization techniques to design our devices, it is convenient to recast the design problem in terms of linear algebra. We do this by discretizing space and making the substitutions We thus wish to find electric fields x i and a permittivity distribution z which satisfy for i = 1 . . . M and j = 1 . . . N i . Here, diag (v) refers to the diagonal matrix whose diagonal entries are given by the vector v, and u † is the conjugate transpose of u. For convenience, we further define the matrices This lets us rewrite equation (S.6) as The final problem we wish to solve is then

Parametrizing the structure
As described in the main article, we describe our structure using a two-dimensional level-set function φ(x, y) : R 2 → R , where the permittivity in the design region is given by For the purposes of numerical optimization, we discretize the level set function in space, which transforms the level set function into a two dimensional array φ ∈ R U ×V . We parametrize the permittivity distribution z with the level set φ by using a mapping function m : When the level set boundaries are not perfectly aligned with simulation grid cells, we render the structure using anti-aliasing. This allows us to continuously vary the structure, rather than being forced to make discrete pixel-by-pixel changes.

Formulating the optimization problem
We are finally in a position to construct our optimization problem. Although there are a variety of ways we could solve (S.10) and (S.11), the particular optimization problem we choose to solve is Here, we constrain the fields to satisfy Maxwell's equations, parameterize the permittivity z with the level set function φ ∈ R U ×V , and construct a penalty function for violating our field constraints from equation (S.11). The penalty f i (x i ) for each input mode is given by where I + (u) is a relaxed indicator function [3], Typically, we use q = 2 and s = max i f i (x i ).

Steepest descent optimization
We solve our optimization problem (S.14) by first ensuring that Maxwell's equations (S.10) are always satisfied. This implies that both the fields x 1 , . . . , x M and the fieldconstraint penalty F are a function of the level set φ. We then optimize the structure using a steepest descent method. Using the chain rule, the gradient of the penalty function F is given by The majority of the computational cost comes from computing the gradient dF/dz. As described in the main text, we evolve the level set function φ by advecting it with a velocity field v ∈ R U ×V . To implement gradient descent, we set the velocity field to be equal to the gradient of the penalty function:

Computing gradient of penalty function F
We now consider how to efficiently compute the gradient of the penalty function F with respect to the permittivity z, which can be written using (S.15) as Although f i is not a holomorphic function since f i : C n → R, we can compute df i /dz using the expression where we have taken the Wirtinger derivative of f i [4]. The Wirtinger derivative with respect to some complex variable w = u + iv is defined as Using this definition, the Wirtinger derivatives ∂f i /∂x i are given by Here, we have used the identity Next, we consider how to take the derivative of the electric fields x i with respect to the permittivity z. If we take the derivative of the discretized Maxwell's equations (S.6) with respect to z, we obtain where we have used our definitions of A i and B i from (S.8). Rearranging, we find the derivative to be We obtain our final expression for df i /dz by subsitututing (S.28) into (S.21): Since A i and B i are large n × n matrices, we have rearranged the expression in the final step to require only a single matrix solve rather than n solves. This method for reducing the computational cost of computing gradients is known as adjoint sensitivity analysis, and is described in detail elsewhere [5].
The cost of computing dF/dz is dominated by the cost of solving Maxwell's equations. For each input mode i = 1 . . . M , we need to solve both the forward problem x i = A −1 i b to find the electric field x i , and the adjoint problem A − † i (z) ∂f i ∂x i † from equation (S.29). Both the forward and adjoint problems can be solved by any standard Maxwell's equation solver [6,7,8]. We use a graphical processing unit (GPU) accelerated implementation of the finite-difference frequency-domain (FDFD) method [9,10].

Level set implementation 2.1 Curvature limiting
In the main text, we wrote that we implement curvature limiting by evolving the level set function φ with where d(x, y) is the Euclidean distance to the nearest element in the set Ω, We choose Ω = {(x, y)|κ(x, y) > κ 0 } to be the set of points with a local curvature greater than our threshold κ 0 . The distance function d(x, y) can be efficiently computed using the Euclidean distance transform commonly included in image processing libraries.

Numerical implementation
In our design algorithm, we apply gradient descent using the partial differential equation where v(x, y) is the local velocity, and apply curvature limiting with equation S.30. We spatially discretize equation S.34 using Godunov's scheme, and equation S.30 using central differencing, as is common practice [11]. We discretize in the time dimension using Euler's method.
To ensure that our level set equations remain well behaved, we regularly reinitialize φ to be a signed distance function [11], where |∇φ| ≈ 1. Most reinitialization schemes, however, result in subtle shifts in the interface locations, which can cause optimization to fail. We use Russo and Smerka's reinitialization scheme to avoid these issues [12].

Fabrication robustness
To obtain a better understanding of the fabrication robustness of our 1 × 3 splitter, we simulated the device for a range of over-etching and under-etching errors, which correspond to lateral growth or shrinkage of the design. We have presented the results in figure S1.

Backreflections
We also simulated the return loss for our 1 × 3 splitter, which we present in figure  S2. The backreflections into the fundamental TE 10 mode of the input waveguide are < 23 dB over the operating bandwidth of our device. Backreflections into the TE 20 mode are zero due to reflection symmetry in our device in the horizontal direction, and mode conversion to TM modes is impossible due to reflection symmetry of our structure in the vertical direction. Finally, scattering into higher order waveguide modes is negligible since the input and output waveguides are close to single-mode. Thus, backreflections into the input waveguide comprise only a small fraction of the total losses.