Using parametric optimization (PO), and python interoperability design a silicon-on-insulator waveguide crossing in varFDTD 2.5D and 3D FDTD to rapidly converge on a solution. This optimization algorithm is based on the adjoint method, and yields optimal solution in just a few iterations compared to other schemes such as particle swarm optimization (PSO). The final design will then be exported into a GDS file for further simulation and/or fabrication.

## Overview

Understand the simulation workflow and key results

Lumerical's inverse design capability provides unparalleled optimization performance by combining the power of gradient based optimization routine with efficiencies found in fundamental properties of Maxwell equations.

This example will demonstrate how to use the inverse design method to generate a crossing shape that delivers near perfect transmission efficiency. Moreover, we'll demonstrate how to modify the example with your parameters so you can reuse this approach for your own design.

This example draws extensively from the LumOpt framework:

### Step 1: Define base simulation parameters

The goal for this initial step is to define the basic parameters of the FDTD simulation that will be used for the adjoint optimization. Here, key simulation settings for the crossing waveguides, monitors for field gradient and figure of merit (FOM) calculations, source and mesh are parametrized for ease of running the optimization routine. It is possible to provide these necessary parameters with either

- MODE project file
- lumerical script (.lsf format)
- Python function handle.

In this example we will use a lumerical script file to set up the simulation that works with the LumOpt optimization routines.

### Step 2: Define the optimizable geometry

The second step is to define a parametrized polygon geometry that will represent the body of the crossing connecting the input and output waveguides (the optimizable geometry). This is defined as a python function that accepts a set of optimization parameters and generates corresponding polygon geometry appended to the base FDTD simulation geometry defined in step 1. Initial parameter values and their limits will also get defined in this step.

### Step 3: Run the parametric optimization in 2.5D

Run the python script with the base simulation script and parametrized polygon object as an input. This will be a 2.5D simulation for rapid delivery of an optimal design which can be used as the initial guess for the 3D optimization in the next step. To find the best set of parameters that define the optimal crossing shape, the optimization routine will use the specified field monitors in the base simulation to calculate:

- The field gradients: the gradients of the field due to perturbation of permittivity resulting from slight changes in geometry of the optimizable region
- The gradient of the figure of merit (FOM): the gradient of the mode overlap to the fundamental TE mode of the waveguides as a result of the shape changes

After optimization, the optimized crossing component shape will be exported into GDS II format which can be used for further simulation and/or fabrication (mask design).

### Step 4: Run the parametric optimization in 3D [optional]

As an optional step, the best solution found by the 2.5D optimization can be used as an initial guess for a 3D optimization to verify the 2.5D simulation results and to further improve the performance of the design if possible. If one does not have a varFDTD license it is possible to start with the 3D optimization; however, without an initial solution the optimization may progress slowly. We have included symmetric boundary conditions to reduce the sim time by a factor of four.

## Run and results

Instructions for running the model and discussion of key results

### Step 1: Define base simulation parameters

This step is not necessary if you want to run the example as is, however it is required for modification of base simulation parameters based on your needs.

- Open the base simulation script file (varFDTD_crossing.lsf) in FDTD
- Modify the desired base simulation parameters
- Save and run the script file to make sure the base simulation is setup as intended and the script runs successfully. This is just a sanity check and the generated simulation is not required to be saved for the next step

### Base simulation setup

The resulting base simulation which includes the FDTD simulation region, source, monitors, mesh override for optimizable geometry, and input and output waveguides is shown in the figure below.

### Step 2: Define the optimizable geometry

This step is not necessary if you want to run the example as is, however it is required for modification of base simulation parameters based on your needs.

- Open the parametric optimization python script file (crossing_opt_varFDTD.py) in FDTD script editor
- Define the initial values of the input parameters of the optimization (the Y coordinates of the nodes of the optimizable polygon) through the array initial_params
- Define the bounds of the optimization parameters within the array bounds
- Make any other necessary modifications to the optimization parameters. See “Important model settings” for more information about these parameters
- Save the script

### Optimization python script

Below is a screenshot of the optimization python script within the script editor window in FDTD

### Step 3: Run the parametric optimization in 2D

- Before running the optimization, make sure the current working directory for FDTD is set to the location where the example files are located
- Run the optimization script file from the script editor in FDTD. Alternatively, you can run the optimization script from the command line using the Python executable provided with the Lumerical installation .

### 2D optimization results

The geometry of the crossing at the start of the optimization with default parameters is shown in the figure below:

A command prompt (in windows operating system) or terminal (in macOS or Linux) window will open showing the progress of the optimization. A screenshot of the command prompt window at the end of the optimization is shown below which contains information about the optimization progress

During the optimization, a new session of varFDTD will be opened and closed by the optimization routine to perform simulations on the guesses made by the optimization algorithm. Another window will also report the optimization progress providing various plots for evolution of the optimization geometry, figure of merit, parameters and other related information. These CAD windows can be hidden by specifying 'hide_fdtd_cad=true' in the optimizer. Another window will also report the optimization progress providing various plots for evolution of the optimization geometry, figure of merit, parameters and other related information. The figure below shows the optimization report window at the end of optimization process. The geometry of the final solution can be seen in the plot and it is shown that the FOM will have little improvement after about 25 iterations which suggests the number of iterations can be further reduced for future runs. The field gradients (due to the changes in geometry) and forward direction field distributions within the device’s structure can also be seen in the plot. In addition, it is clear that the changes in the optimized parameter values become less intense after a few iterations since the FOM gets closer to the optimum value. The changes in FOM gradient (based on changes in the parameter values) also decrease on average as seen in gradient evolution plot as the iteration continues.

From the plot, it can be seen that the optimization started with an initial design providing an FOM of ~0.75 and ended up with a geometry with a FOM of about 0.96, resulting in about 27 % improvement through the optimization process in just a few iterations which highlights the advantage of this method compared to other optimization algorithms such as particle swarm which need many iterations to achieve an optimum solution. Also, note that FOM will have little improvement after about 5 iterations which suggests higher number of iterations probably won't achieve a significant improvement. The field gradients (due to the changes in geometry) and forward direction field distributions within the device’s structure can also be seen in the plot. In addition, the changes in the optimizable parameter values become less noticeable after a few iterations as the FOM gets closer to the optimum value. The changes in FOM gradient (based on changes in the parameter values) also decrease on average as seen in gradient evolution plot as a function of the iterations.

At the end of optimization, the final parameter values will be saved in a text file (2D_parameters.txt) in the same folder as the example files which will be used in the 3D optimization for the next step. In addition, copies of the simulation files for each iteration, a text file containing the history of parameters and FOM values for each iteration and still images of the optimization report window for each iteration will be saved in a folder (with a name starting with “opts”) in the same folder as example files for future reference. Each optimization run creates a new folder with the same name but a numbered suffix to avoid overwriting previous runs data.

Here is the final geometry optimized in 2D:

The optimization python script file will also export the optimal design to a GDS file at the end of optimization. The GDS file will be saved in the same folder as the example files. The example files come with an encrypted script files with .lsfx extension that is needed for the GDS export operation. They should be saved in the same location as the python script file.

### Step 4: Run the parametric optimization in 3D (optional)

- Before running the optimization, make sure the current working directory for FDTD is set to the location where the example files are located
- Make sure the base simulation script file (FDTD_crossing.lsf) is setup correctly
- Open and run the optimization python script file (crossing_opt_3D.py) in FDTD script editor. Alternatively, you can run the optimization script from the command line using the Python executable provided with the Lumerical installation . The script will read the final parameters from the 2D optimization (stored in 2D_parameters.txt) if available. Otherwise, it will start with default values set in the script.

### 3D optimization results

Note that the 3D optimization might take a long time to complete due to the large size of the simulation.

Looking at the results of the 3D optimization will reveal that the crossing’s transmission has been initially predicted to be about 82% by 3D simulation (which is more realistic than the 2D case) and has been improved to 91% at the end of optimization. It can be seen that the FOM does not improve significantly after 6 iterations suggesting that running a longer optimization may not result in much better performance.

## Important model settings

Description of important objects and settings used in this model

** Definition of optimizable geometry ** : The optimizable geometry is defined as a polygon which has fixed values for the x coordinates of its boundary points and the points’ y coordinates can be modified to obtain the optimal geometry. These y coordinates are defined as the nodes of a cubic spline (as shown below), which makes it possible to perform cubic spline interpolation on the points in order to obtain a smooth polygon geometry. For the waveguide crossing, this optimizable polygon needs to be mirrored rotationally in order to form a complete crossing, meaning that only 1/8 of the optimizable geometry is independently parametrized. It is important to make sure that the dimensions of the optimizable geometry and the waveguides set in the base simulation match so that the crossing geometry is properly formed.

** Optimization fields monitor ** : The optimization fields DFT monitor (opt_field) is used to collect the field data within the optimizable geometry which is then used to calculate field gradients used in optimization algorithm. As such, the location of this monitor is of great importance and should cover the entire optimizable geometry.

** Figure of merit fields monitor ** : Since this monitor is used to calculate the figure of merit required for optimization (mode overlap to the fundamental TE mode of the output waveguide), it should be located within the output waveguide of the crossing with appropriate dimensions.

** Height of the geometry objects ** : The height (depth) of the geometry objects should match in both base simulation (waveguides) and the optimizable geometry (polygon). For the waveguides, this is set through the base simulation setup script (by setting their z span) and for the polygon, it is adjusted as the input argument “depth” when calling the cross function to form the optimizable geometry. This value should also be passed to the GDS export portion of the script as depth:

** Optimization parameters bounds and initial values ** : The range over which the optimization algorithm is allowed to vary the parameters and also the initial values for those parameters as a starting point are defined by initial_params and bounds arrays as input arguments when calling the function FunctionDefinedPolygon

** Simulation bandwidth ** : The inverse design PO method supports broadband simulations. The bandwidth is defined using start and stop wavelengths complemented by the number of frequency points in the given interval. The wavelength specification is passed to the optimization script when defining the FOM. Therefore, the wavelength settings defined in the base simulation script will be overridden.

Additionally, each frequency point can have their target FOM specified independently. This means that the optimization routine can aim to match custom FOM profile as a function of wavelength rather than simply maximizing the FOM across the full bandwidth.

Here is a code section from this Y-branch example that specifies the bandwidth from 1300 to 1800nm with 21 frequency points:

wavelengths = Wavelengths(start = 1300e-9, stop = 1800e-9, points = 21)

The bandwidth settings is then passed to the optimization routine during FOM definition. Notice that we specified FOM target of 1 for all frequency points to maximize the transmission across the desired bandwidth:

fom = ModeMatch(monitor_name = 'fom',

mode_number = 'fundamental mode',

direction = 'Forward',

target_T_fwd = lambda wl: np.ones(wl.size),

norm_p = 1)

** Maximum number of iterations ** : While the algorithm has the capability to stop the optimization once the gradient of the figure of merit falls below a certain threshold, the “max_iter” variable used for defining the optimization algorithm can be used to limit the number of iterations the algorithm is able to perform.

** Symmetry ** : As of 2019a r4 the inverse design PO method supports symmetric boundaries. This can be particularly useful for 3D simulations where it can decrease the simulation time by factor of 4.

## Updating the model with your parameters

Instructions for updating the model based on your device parameters

** Geometry ** : If you need to have your own geometry defined for the crossing, including SOI device layer thickness (waveguide height) and waveguide width, corresponding changes should be made in base simulation setup script and/or optimization setup python script. This might need changes in objects span and location in the base script, the polygon’s “depth” parameter in python script which its value is also passed to the GDS export portion of the script. The source, simulation region, mesh override, and field gradient and fom monitors dimensions should also be adjusted accordingly to make sure they cover the entire structure properly. The coordinates of the initial values and bounds of the optimizable geometry at points where it is connected to the waveguide might need to be adjusted to ensure proper connection between the optimizable region and the waveguides.

def cross(params): y_end = params[-1] x_end = 0 - y_end points_x = np.concatenate(([-2.01e-6], np.linspace(-2e-6, x_end, 10))) points_y = np.concatenate(([0.25e-6], params)) n_interpolation_points = 50 polygon_points_x = np.linspace(min(points_x), max(points_x), n_interpolation_points) interpolator = sp.interpolate.interp1d(points_x, points_y, kind = 'cubic') polygon_points_y = [max(min(point, 1e-6), -1e-6) for point in interpolator(polygon_points_x)] pplu = [(x, y) for x, y in zip(polygon_points_x, polygon_points_y)] ppld = [(x, -y) for x, y in zip(polygon_points_x, polygon_points_y)] ppdl = [(-y, x) for x, y in zip(polygon_points_x, polygon_points_y)] ppdr = [(y, x) for x, y in zip(polygon_points_x, polygon_points_y)] pprd = [(-x, -y) for x, y in zip(polygon_points_x, polygon_points_y)] ppru = [(-x, y) for x, y in zip(polygon_points_x, polygon_points_y)] ppur = [(y, -x) for x, y in zip(polygon_points_x, polygon_points_y)] ppul = [(-y, -x) for x, y in zip(polygon_points_x, polygon_points_y)] polygon_points = np.array(pplu[::-1] + ppld[:-1] + ppdl[::-1] + ppdr[:-1] + pprd[::-1] + ppru[:-1] + ppur[::-1] + ppul[:-1]) return polygon_points

It is recommended to set the coordinates such that there is a slight overlap between the optimizable region and the edges of the waveguides to ensure they are connected in the simulation.

** Materials ** : The materials included in the simulation (the material making the optimizable geometry and the material surrounding it) are defined in Step 1 via lsf, project file or python function. In this example we have taken the index of silicon and silicon dioxide at the desired simulation wavelength of 1550nm. One should make sure that the material defined (or given index) in the initialization is the same as that passed to the optimizer when defining the geometry as shown below;

eps_in = Material(name = 'Si: non-dispersive', mesh_order = 2)

eps_out = Material(name = 'SiO2: non-dispersive', mesh_order = 3)

depth = 220.0e-9

polygon = FunctionDefinedPolygon(func = splitter,

initial_params = initial_points_y,

bounds = bounds,

z = 0.0,

depth = depth,

eps_out = eps_out,

eps_in = eps_in,

edge_precision = 5,

dx = 1.0e-13)

With varFDTD it is now possible to use the material library to define the material properties; however, the optimization technique has not been tested with dispersive materials; thus, the effective index is solved at the center frequency and passed to the optimizer.

## Additional resources

Additional documentation, examples and training material

### See also

- Inverse design of Y-branch (PO)
- Y Branch PSO
- Y Branch S parameter extraction
- KX posts with tips and related content
- Python API

### Related Lumerical University courses

## Appendix

Additional background information and theory

### Parametric Optimization (Adjoint method)

In a typical device optimization, the goal is to minimize or maximize a figure or merit \( F(p) \) that depends on a set of \( N \) design parameters \( p=(p_1,...,p_n) \). The figure of merit can be any quantity used to benchmark the device’s performance and the design parameters can be any quantity that alters the device’s response. To evaluate \( F(p) \) for a given set of values of \( p \), a physics based modeling tool - such as FDTD - is used to map values of \( p \) to values of \( F \). The model usually contains a description of the device’s geometry, materials and all relevant environmental conditions. In a parameterized shape optimization, the goal is to allow some parts of the device’s geometry to vary to minimize or maximize the figure of merit. For this purpose, a shape function \( f(p) \) is typically defined to specify part of the device’s geometry; an example is shown in the following figure:

Y-splitter with a parameterized junction boundary \(f(p)\) defined using connected splines.

The shape function is typically user defined and must be fed to the modelling tool to construct the device’s geometry. Once the figure of merit \( F(f(p)) \) can be evaluated for a given shape function, a minimizer can make multiple evaluations of the figure of merit to find the optimal shape. There is a very large body of minimization tools and techniques that can be used to try to find the optimal shape parameters. When selecting a minimization technique, it is important to consider that running a device simulation to evaluate the figure of merit can be a slow operation. So, the best suited minimization techniques are those that limit the number of evaluations of the figure of merit. With a good starting point, gradient based methods can quickly find the optimal parameters with a relatively small number of figure of merit evaluations provided that it is also possible to evaluate the gradient of the figure of merit:

$$∇_pF=(\frac{\partial F}{\partial p_1},⋯,\frac{\partial F}{\partial p_n})$$

To evaluate the gradient without having to perform one simulation for each partial derivative, a technique called the adjoint method can be employed. This technique exploits the properties of linear partial differential equations to evaluate the gradient using only two simulations independently of the number of optimization parameters. The Python based shape optimization module in FDTD (lumopt) allows users to define a shape function and run the corresponding device simulation to evaluate the figure of merit. The module also implements the adjoint method to provide the gradient. In this way, any of the gradient based minimizers in SciPy’s optimize module can be used to perform a fully automated optimization.