Controlling the Minimal Feature Sizes in Adjoint Optimization of Nanophotonic Devices Using B-spline Surfaces

TL;DR

This paper introduces a level-set method based on B-spline surfaces for adjoint optimization in nanophotonic device design, enabling precise control over minimal feature sizes to improve manufacturability.

Contribution

It presents a novel B-spline surface-based level-set approach that allows explicit control of feature sizes during the inverse design of nanophotonic devices.

Findings

Controlled feature sizes in device designs

Successful design of wavelength demultiplexer

Implementation of neural computing structure

Abstract

Adjoint optimization is an effective method in the inverse design of nanophotonic devices. In order to ensure the manufacturability, one would like to have control over the minimal feature sizes. Here we propose utilizing a level-set method based on b-spline surfaces in order to control the feature sizes. This approach is first used to design a wavelength demultiplexer. It is also used to implement a nanophotonic structure for artificial neural computing. In both cases, we show that the minimal feature sizes can be easily parameterized and controlled.

Equations

$$\begin{gathered}N_{(u)}^{i,0}=\left\{\begin{array}[]{ll}1\quad u_{i}<u<u_{i+1}\\ 0\quad otherwise\end{array}\right.\\ N_{(u)}^{i,k}=\frac{u-u_{i}}{u_{i+k}-u_{i}}N_{(u)}^{i,k-1}+\frac{u_{i+k+1}-u}{u_{i+k+1}-u_{i+1}}N_{(u)}^{i+1,k-1}\end{gathered}$$

$$C_{(x)}=\sum _{i=0}^{n_x}{N}_{(x)}^{i,k}{P}_{(i)}$$

$$\begin{array}{c}{\varphi }_{(x,y)}=\sum _{i=0}^{n_x}\sum _{j=0}^{n_y}{N}_{(x)}^{i,k_x}{N}_{(y)}^{j,k_y}{P}_{(i,j)}\end{array}$$

$$\varepsilon(r)=\left\{\begin{array}[]{ll}\varepsilon_{1}\quad\phi(r)<0\\ \varepsilon_{2}\quad\phi(r)>0\end{array}\right.$$

$$\begin{array}{c}\frac{\partial J}{\partial P_{(i,j)}}={\int }_s(\frac{\partial J}{\partial {\epsilon }_{({\varphi }_{(x,y)})}})(\frac{\partial {\epsilon }_{({\varphi }_{(x,y)})}}{\partial {\varphi }_{(x,y)}})(\frac{\partial {\varphi }_{(x,y)}}{{\partial }_{(i,j)}})\\ \Rightarrow \frac{\partial J}{\partial P_{(i,j)}}={\int }_s(\frac{\partial J}{\partial {\epsilon }_{({\varphi }_{(x,y)})}})({\delta }_{{\varphi }_{(x,y)}})(N_{(x)}^{i,k_x}{N}_{(y)}^{j,k_y})\end{array}$$

$$\begin{gathered}P^{t+1}=P^{t}-\gamma\frac{\partial J}{\partial P}\Big{\rvert}_{P=P^{t}}\\ P^{t+1}=\frac{P^{t+1}}{||P^{t+1}||_{\infty}}\end{gathered}$$

$${\partial }_t\varphi +v(x,y)|\nabla \varphi |=0$$

Full text

††preprint: APS/123-QED

INTRODUCTION

Adjoint method has become a very useful tool in the inverse design of photonicsLalau-Keraly et al. (2013); Callewaert et al. (2018); Piggott et al. (2017, 2015); Frei et al. (2007); Su et al. (2018); Pestourie et al. (2018); Frellsen et al. (2016); Jensen and Sigmund (2011); Michaels and Yablonovitch (2018); Hughes et al. (2018); Frandsen and Sigmund (2016); Frei et al. (2008); Molesky et al. (2018); Vercruysse et al. (2019); Miller and Yablonovitch (2015) and among other optimization methods such as genetic algorithmsOffrein et al. (1998), swarm optimizationMak et al. (2016), or nonlinear optimizationShen et al. (2015), adjoint method is particularly efficient as it computes the gradient of a cost function w.r.t. a large number of design parameters.

The inverse design process mostly involves the optimization of the dielectric constant distribution $\varepsilon_{(r)}$. In a gradient-based approach, $\varepsilon$ is treated as a continuous variable. However, in practice, the permittivity can only take discrete values determined by the types of materials used. To bridge this gap, one can apply the method of level-set or density-based topology optimization. However, even with level-set, or density-based topology optimization to address the issue of discrete material parameters, we still end up with structural features that might be too small for today’s nanofabrication. Additional constraints need to be applied to control the size of minimum features. For example RefShen et al. (2015) defines the medium as a collection of pixels with a fixed size larger than the critical dimension of fabrication; while Hughes et al. (2018) uses the method proposed in Zhou et al. (2015): a spatial low pass filter which controls the feature sizes. Ref Su et al. (2017) uses the convolution between the gray-scale medium and a disk that defines the minimum feature size. Refs Frei et al. (2007, 2008) use radial basis functions to control the feature size. In this work, we are going to use b-splines for topology optimization, which has been a highly successful and proven approach for size control in structural shape optimization in the field of mechanics Qian (2013).

B-SPLINES

Splines are functions generally used for interpolating data and designing curves and surfaces in computer aided designs and computer graphics. These functions are generated by a linear combination of piece-wise polynomial functions. Their strength is in the fact that they are capable of accurate interpolation of data with low degree polynomial components, that would otherwise require a high degree non-local polynomial. The degree of the spline is the degree of its highest order polynomial.

One way of representing a spline is by defining a set of basis functions, called basis-spline functions or b-spline functions. In order to define any one of these basis functions, the domain over which the spline is going to be defined, has to be segmented with a set of knots $U=[u_{0},...,u_{m}]$. In this work we only work with knots that appear only once in this set, and are equidistant from one another. Another property that is required to define a b-spline function is its degree (k). With both knots and degree defined, basis-spline functions can be expressed using the Cox-de Boor recursion formula De Boor (1972); Cox (1972).

[TABLE]

In the equation above $N_{(x)}^{i,k}$ shows a b-spline of degree $k$. As it can be seen, a b-spline of degree zero is only non-zero on one of the segments of the domain, while higher degrees span more segments of the domain. As a matter of fact, a b-spline of degree $k$ defined over a domain segmented by $m+1$ knots spans $k+1$ one of these intervals ($m+1-k$ such b-splines are defined over the domain). This can be seen in Fig.1 where a domain has been segmented by 4 knots.

With the basis spline functions, a curve can be generated using a linear combination of these basis functions. In Eq.2, the $P_{(i)}$’s are the set of coefficients that determine the contribution of each b-spline component to the curve.

[TABLE]

We now discuss the control of feature sizes. The degree of the basis function(k) plays a critical role in controlling the minimal feature size. This can be clearly seen in Fig.1: as the degree of the basis function is increased, the minimum area that it affects also increases and so the minimum size of the features increases. Another way to control the minimal feature size is to tune the spatial density of the knots. A similar effect happens when the knots are distributed more sparsely over the domain. This effect can be seen in Fig.2: the denser the knots are packed together, the smaller the area each of these functions can affect becomes. Consequently, more nuances can be generated in a curve. This is the approach we have taken in this work to control the critical dimensions of our designs.

FORMULATION

B-splines can be used for a density based approach or a level-set one. Next we use 2 dimensional (2D) photonic design to illustrate this method in the context of the level-set method. In 2D, we use basis spline functions to construct a surface $\phi_{(x,y)}$ using the following equation:

[TABLE]

Similar to Eq.2, $N_{(x)}^{i,k_{x}}$ and $N_{(y)}^{j,k_{y}}$ are basis-spline functions of degree $k_{x}$ and $k_{y}$ respectively, and $P_{(i,j)}$ are the set of coefficients. If we were to limit the value of the coefficients between 0 and 1, the generated surface could be used in a density-based topology optimization. Here, we focus on the level-set method.

In order to use the level-set method for topology optimization, we first need to define a surface $\phi_{(x,y)}$ over the area that the optimized device will occupy (The value of each point on this surface can vary across a limited range around zero). Then the medium can be shaped according to Eq.1 based on this surface.

[TABLE]

Now, optimizing the medium with a level-set method consists of evolving the boundaries between the two materials (all the points where $\phi_{(r)}=0$) in a manner that improves the performance of the device. With a b-spline surface for this task, we have direct control on this surface, and consequently on the boundary between the constituent materials. Having defined the level-set surface as a b-spline surface, we can now directly optimize a cost function J with respect to the control parameters of the b-spline surface.

[TABLE]

In Eq.5, the first term can be calculated with the adjoint method, the second term(a Dirac delta function which is nonzero only where $\phi_{(x,y)}$ passes zero),can be computed with marching squares algorithmKroon (2011), and the final term is a simple matrix multiplication.

Once the gradient of the cost function is calculated, gradient descent can be used to update the b-spline coefficients. However, in order to ensure the stability of the process, in each update iteration the coefficients are normalized by the highest value among themBernard et al. (2009). The mentioned steps can be summarized as shown in Eq.6, where $P^{t}$ is the matrix of all the b-spline coefficients at iteration $t$ and $||P^{t}||_{\infty}$ is the maximum absolute value in this matrix.

[TABLE]

APPLICATION

In order to demonstrate the application of this approach, we utilized it to design two nanophotonic devices: a wavelength demultiplexer Piggott et al. (2015) and a nanophotonic structure for artificial neural computing Khoram et al. (2019).

The wavelength demultiplexer is a three port device that guides the light coming from the input port to one of the two output ports based on its wavelength. In our design, this device is of size $4\,\mu m\times 2\,\mu m$ designed to demultiplex $1500\,nm$ and $1550\,nm$ wavelengths. We optimize three versions of this device with the intention of showcasing how the critical dimensions of the device can be tuned by the basis spline functions.

The b-spline surface is defined on the same grid as the one used for the electromagnetic simulation and over the area to be optimized. We can use either the distribution of the knots on the grid, or the degree of the b-splines to control the critical dimension. In this work we opt for the former. We define our basis functions as cubic b-splines, and use three different knot distributions to realize three version of the devices with different feature sizes. The first design has the knots placed at every second point on the simulation grid, the second design at every fifth point, and the third design at every ninth point. This effectively provides us with 4559, 629, and 152 b-spline coefficients respectively. Once the b-spline surface is set up, we use Eq.4 to define the medium, where $\varepsilon_{1}$ is set to the permittivity of $SiO_{2}$ and $\varepsilon_{2}$ is set to that of $Si$. Finally, a cost function is necessary to optimize the structure. The cost function we define here tries to maximize the the transmission of the time-averaged energy flux to the output port corresponding to the input wavelength.

The devices were optimized with respect to the cost function mentioned above using Eq.4 and Eq.5. The results for the optimization of these three instances have been depicted in Fig.3. As it was mentioned earlier, since placing the knots more sparsely results in an increase in the area that each basis-spline function affects, the critical dimension of the device increases. This can clearly be seen in Fig.3. As we move from Fig.3(a) to Fig.3(c) the feature sizes distinctly expand. However, even for the case of the sparsely distributed knots, the occasional nuances in the b-spline surface can cause the appearance of small features in the device. This problem is mitigated by applying erosion and dilation processes on the device in each update iteration to eliminate these occasional smaller features. We should point out that this control over the feature sizes of the device comes at a price: there is a trade off between the performance of the device and the critical dimensions of the structure. As we use larger critical dimensions, the degree of freedom, i.e. the number of b-Spline functions used, decreases, and the optimization is further constrained, consequently the performance generally degrades. This is also evident in Fig.3, where the transmission drops moving from the first device to the last one.

Nanophotonic artificial neural computing This device can accomplish the same task of neural inference as in Khoram et al. (2019) . The objective is to recognize the value of a handwritten digit 0-9. The input light from the handwritten digit is focused by a nanophotonic structure to different locations according to the value of the digit. The medium then has to map all the different writing styles of the same digit to the specific location corresponding to that digit; furthermore, the medium must perform well on handwritten digits that it has not yet seen.

The setup for this task is as follows: the image of the handwritten digit is first vectorized and mapped to an array of light sources placed on the left side of the nanophotonic structure. Then a set of ten receivers corresponding to the each of the ten classes 0-9 are placed on the right side of the medium. Inside the device, the input light scatters and reflects off of the numerous material interfaces, and then focuses on a spatial location on the output side. On the output side, the receivers measure the light intensity at their respective locations and the class with the maximum measured value is selected as the correct class for the input image.

The device is made possible by a much more extensive optimization process than that for the wavelength demultiplexer. This optimization process is equivalent to the training process in digital neural networks. We define a cost function and use adjoint method to calculate the gradient. At the time of training, we measure the intensity of light at the location of all the receivers. Then we use a cross-entropy cost function between the normalized output array and a $1\times 10$ one-hot vector(a vector that is all zeros except at the index of the correct class), to calculate the cost value for that specific instance. Details of the working principle can be found in RefKhoram et al. (2019). In this work we implement the level-set function based on a b-spline surface to control the minimal critical dimension. In our previous work Khoram et al. (2019), the level-set method was implemented based on a surface set as a signed distance function. Here, the initial medium is created by generating a random set of values for the b-spline coefficients, and then the medium is created based on the generated surface with the surface evolving afterwards by updating the b-spline coefficients. On the other hand, in RefKhoram et al. (2019) the medium was generated by randomly distributing inclusions made up of a second material inside the host medium, and then a surface $\phi_{(r)}$ is generated based on this initial medium (similar to the previous section, the surface relates to the medium with Eq.4). The structure then evolves by updating this signed-distance surface at each training iteration using a modified version of the following equation.

[TABLE]

In this equation $v(x,y)$ is the velocity with which each point on the zero crossing curves of the level-set function moves normal to each of those curves. This velocity is set equal to the gradient acquired with the adjoint method. This implementation of the level-set method for the image classification task turned out to be quite dependant on the initial distribution of the inclusions, and therefore we had to initialize the medium with many small inclusions to get a good performance from the device. Doing so resulted in small feature sizes in the device which make the fabrication of the device rather difficult. However, with the b-spline approach to level-set function this problem does not happen as the optimization process is less dependant on the initial shape of the medium.

We optimize two versions of this device, one with the signed distance function approach, and the other with the b-spline approach. Both devices were optimized for the wavelength $1\,\mu m$ with the size $42\,\mu m\times 30\,\mu m$ and made up of $SiO_{2}$ as the host material and $air$ as the inclusions. For the signed-distance implementation 2000 inclusions of size $5\,\mu m\times 5\,\mu m$ are spread throughout the device for initialization, and for the b-spline implementation the knots are distributed at every tenth point on the simulation grid in both axial dimensions. The accuracy on the test set for the two devices turned out to be $76.8\%$ and $77.2\%$. The performance of the b-spline approach is slightly better in this case which is not of great consequence as the gap between the two values can be reduced with better initialization of the signed-distance approach. However, the important point that is apparent in Fig.4 is the difference between the shape of the two structures. As it can be seen, the first device shown in Fig.4(a) has much smaller features whereas the second device depicted in Fig.4(b) has a much larger critical dimension. This makes the second device much easier to fabricate, and thus the b-spline approach much more desirable. Finally, Fig4.(c) and Fig.4(d) show the device in action. As it can be seen, light from different class of images is focused at distinct locations corresponding to that class. Moreover, different image samples of the same digit produce different field distributions; however, these different field distributions result in the same output class.

CONCLUSION

In this work, we adopt a proven technique to size control that is used in the field of mechanics to the inverse design of photonics. It use B-spline basis to describe the distribution of the dielectric constant. The minimal feature size is controlled by the knots distribution and the degree. This method is very easy to implement and interpret. We demonstrated the proposed method for a wavelength demultiplexer, and a nanophotonic structure for artificial neural computing. The first class of nanophotonic devices we implemented was used to demonstrate how to handle the trade off between the complexity and their performance. For the second class, two variations of the level-set methods were used to achieve the same goal. The comparison between the two cases shows how having a well-defined set of control parameters can help us in designing a nanophotonic device in one optimization stage, without the need to have a good initial guess of the variables we are trying to optimize.

I Acknowledgements

The authors thank W. Shin for his help on improving the computational speed of the implementation of this method. This work was partially supported by the National Science Foundation under Grant No. 1561917. It was also partially supported by the Defense Advanced Research Projects Agency (DARPA), under agreement HR00111820046. The views, opinions and/or findings expressed are those of the author and should not be interpreted as representing the official views or policies of the Department of Defense or the U.S. Government.