Inverse Design of Grating Coupler (2D)

In this example, we will use the Inverse Design toolbox (lumopt) to design a silicon-on-insulator (SOI) grating coupler. Compared to other optimization methods such as particle swarm optimization (PSO), this optimization algorithm enables obtaining the best solution in just a few iterations. The optimum design will be exported into a GDS file for further simulation and/or fabrication.

download example

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 TE silicon on insulator (SOI) grating coupler design with maximized coupling 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:

• Photonic Inverse Design Overview - Python API
• Getting Started with lumopt - python API
• Optimizable Geometry - Python API

Step 1: Define base simulation in MODE

The goal for this initial step is to find a good starting point for the optimization. Using MODE FDE solver, we determine the initial linearly apodized grating based on the optimization parameters (central wavelength, etch depth, fiber angle).

Step 2: Optimization of the linearly apodized grating using 2D FDTD

The second step is to optimize the grating using a limited set of parameters. We use the apodized grating defined in step 1 as a starting point. The number of gratings is fixed at 25 and the etched depth at 80nm. A set of 4 parameters obtained from step 1, defining the grating apodization function, is used to extract the initial set of parameters to maximize transmission/coupling.

From the initial condition found in the previous step, we now run an optimization using the x coordinate of each wall position of the grating as a free parameter (2x 25 parameters).

Include manufacturing constraints

We run the 2D optimization enabling the minimum feature size constraint (set to 100nm) to make sure the device will be manufacturable.

Step 3: Extract 3D design and GDS file

The optimized grating coupler component shape with curved gratings will be exported into GDS II format and into a 3D simulation file which can be used for further simulation and/or fabrication (mask design). We used circular gratings that cover the source area to make sure that all the fiber light is collected. We also used a short linear taper to connect to a 450nm wide waveguide. The geometry of the gratings and taper can be further optimized by the designer.

Photonic Inverse Design: Grating Coupler (3D)

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 the base simulation file based on your needs.

For a given center wavelength \( \lambda_c \) and a given etch depth \( e \), we can follow ref. [1] to design an optimized apodized grating based on physical arguments. Here, we only briefly sketch the derivation, for more details and precise definitions of the variables, see the reference.

The main idea is to use the local effective index of an apodized grating. Specifically, for each unit cell, the effective index is approximated using a linear interpolation as:

$$ n_{eff} = F n_0 + (1-F) n_E $$

Where \( F \) is the filling fraction, \( n_O \) and \( n_E \) denote the effective indices of a mode propagating in a waveguide with the full height and the height of the partially etched structure, respectively.

To determine those effective indices, we use the MODE FDE solver. The calculation is performed for a silicon waveguide with height \( h = 220nm \), using a refractive index of \( 3.47668 \) for silicon, and \( 1.44401 \) for the background material (SiO2).

We can write the Bragg condition for a grating with periodicity \( \Lambda \) as

$$ \Lambda = \frac{\lambda_c}{n_{eff} -n_{bg} \sin \theta} $$

Where \( \theta \) is the angle of the incoming beam.

For a linearly apodized grating, the filling factor varies with the position along the waveguide

$$ F ( x ) = F_0 - R \cdot x $$

\( R \) is the apodization factor and \( F_0 \) the initial filling fraction. In the ideal case, it should be 1 as we start with a solid waveguide, but we will use F0=0.95 to avoid extremely small trenches.

Using the linear apodization in the interpolated effective index:

$$ n_{eff} ( x ) = F ( x ) n_0 + ( 1 - F ( x ) ) n_E = n_E + F ( x ) \Delta n $$

We get the spatially varying period of the grating:

$$ \Lambda ( x ) = \frac{ \lambda_c }{ (n_E - n_{bg} \sin \theta) + F(x) \Delta n } $$

To find a good initial condition, we discretize the filling fraction \( F_i \) and the pitch \( \Lambda_i \) for each period.

So far, we only have 1 unknown, the factor \( R \). For a slightly more efficient optimization, we can include 4 optimization parameters:

$$ p = [ x_0, R, a, b ] $$

With

$$ F_i = F_0 - R ( x_{i-1} - x_0) $$

$$ \Lambda_i = \frac{\lambda_c}{a + F_i b} $$

And

$$ a = n_E - n_{bg} \sin \theta $$

$$ b = \Delta n $$

Note \( x_0 \) and \( R \) are to be defined separately. We will set \( x_0 \) to \( –2.5\mu m \) and we will use \( R = 0.03\mu m^{-1} \) (see ref. [1]). If you change the materials, geometries, or polarization, you should perform sweeps of these values to find a suitable initial point (grating).

The script pid_grating_coupler_preliminary_grating_design.lsf will calculate the initial parameters and save them in the pid_grating_coupler_initial_params.json file to be used in the next step.

Step 2: Optimization of the linearly apodized grating using 2D FDTD

We first need to find the optimum position for the fiber.

1. [[step]]Open the python script file pid_grating_coupler_sweep_2D.py in FDTD script editor
2. If needed, modify the parameters
3. Run the script

This script will load the initial design from step 1 and run a parameter sweep fiber to determine the fiber position for optimal coupling.

The optimum position is found at \( 4.2 \mu m \), with a transmission over 60%.

4. [[step]]Open the python script pid_grating_coupler_2D_apodized.py from FDTD script editor
5. If needed, update the parameters, then run the script
6. Open the script pid_grating_coupler_2D_1etch.py
7. Run the script

In this step, we will use the parameters obtained from 1 to 3 and then extract the x coordinate of each wall position as free parameters. More specifically, the parameters will be the starting position of the grating, the width of the etch, and the width of the tooth of each grating.

The script pid_grating_coupler_2D_apodized.py will save to the file pid_optim_1.json these parameters, which can be used as a starting point.

We use the minimum feature size constraint in this step by setting min_feature_size=0.1, which will ensure that no trench is less than 100nm wide.

The final transmission is around 60% and now the grating coupler will be easily manufacturable. The optimized grating parameters are saved in pid_optim_final.json .

Step 3: Extract 3D design and GDS

1. [[step]]Open the script pid_grating_coupler_3D.py
2. If needed, modify the settings
3. Run the script

This script will load the results of the optimization, set the corresponding 3D simulation, and extract the GDS from the structure.

The 3D simulation can be used for further analysis and/or optimization. Note we should first determine the optimum fiber position with this 3D structure. The simulation file already includes a parameter sweep for that purpose.

Important model settings

Description of important objects and settings used in this model

The main parameters of the optimization are set at the beginning of the scripts, and in the base simulation file pid_grating_coupler_2D_TE_base.fsp .

# Optimization global parameters
lambda_c = 1.55e-6 
bandwidth_in_nm = 0 #< Only optimize for center frequency of 1550nm
F0 = 0.95
height = 220e-9
etch_depth = 80e-9
y0 = 0
x_begin = -5.1e-6
x_end = 22e-6
n_grates = 25

indexSi = 3.47668
indexSiO2 = 1.44401

Polarization : The chosen refractive index values are representative of an SOI chip manufacturing process. Because of the high refractive index contrast between silicon and silicon oxide, there is a large difference between effective indices of the two fundamental modes – TE and TM – of the integrated waveguide. For this reason, SOI grating couplers are strongly polarization selective. The presented design excites a TE mode since this is the most common choice, however, it is also possible to set up an optimization targeting the TM mode. To change the polarization to TM, or higher order mode, one needs to change the mode number in pid_grating_coupler_preliminary_design.lsf and in the ModeMatch class of the pid_grating_coupler python files. Also one should change the polarization angle in base fsp files. For TM optimization the mode number is 2 and the polarization angle is 0 degrees in this configuration.

Tilting angle : The coupling efficiency depends heavily on how the fiber meets the top silicon oxide cladding layer. In this example, we use a Gaussian beam as an approximation of the output of the beam from the fiber. The end of the fiber is assumed to be polished at a small angle so that the fiber tilts as it is mounted on the top cladding. Such tilting prevents reflections into the fiber.

Etch depth : The coupling efficiency is highly sensitive to the grating’s pitch, duty cycle, and etch depth. For simplicity, a fixed etch depth is employed here, however, it can also be varied if the available fabrication process provides that degree of freedom.

Substrate : If a silicon substrate is present in the fabricated device, it should be included in the simulations. The substrate will have a noticeable impact on how light is coupled and cannot be omitted as is often done in other device designs.

Optimization figure of merit (FOM) : The purpose of the design is to have the best possible coupling at a desired wavelength, the optimization figure of merit is chosen to be the transmission through the grating coupler at 1550nm, and the optimization algorithm will try to maximize this value. Lumopt uses mode power coupling, and when sweeping the source position we use the power transmission.

Definition of optimizable geometry : The optimizable geometry is defined as a polygon that has fixed values for the y coordinates of its boundary points and the points’ x coordinates can be modified to obtain the optimal geometry.

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 the 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 grating coupler with appropriate dimensions.

Height of the geometry objects : The height (depth) of the geometry objects is chosen based on the fabrication/foundry process. This is especially important for the 3D simulation to ensure correct geometry setup. This is done by adjusting the “height” and “etch_depth” in the scripts. These values are also used in the script exporting the structure to a GDS file.

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 FunctionDefinedPolygonI.

Optimization spectral range : The spectral range of the optimization can be specified as the array “wavelengths” as well as the number of frequency points considered in optimization. Note that choosing a large number of frequency points will make the optimization much slower and can cause issues when transferring a large amount of data between FDTD and the Python environment. It is recommended to keep this number as low as possible especially for 3D simulations. For a single frequency optimization, simply choose the same start and stop values with points set to 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.

Target FOM : The target values for FOM at different wavelengths can be specified with the "target_T_fwd" input argument when defining the FOM. This should be an array with the same length as wavelength data and contain target FOM values (forward transmission) for each wavelength which can be less than 1. The default is 1 for all wavelengths. This gives users the option to aim for a desired value of transmission within the design spectral range.

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 grating coupler, including SOI device layer thickness (grating height) and etch depth, corresponding changes should be made in the base simulation file and/or optimization setup Python script. This might need changes in object span and location in the base script, the polygon’s “etch_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.

In this example, we use the function FunctionDefinedPolygon to set the grating coupler geometry as a polygon. Alternatively, we can use ParametrizedGeometry instead, giving more flexibility in how the geometry is set. The scripts offer the possibility to switch from one to the other by commenting/uncommenting the corresponding lines in the function runGratingOptimization :

def runGratingOptimization(bandwidth_in_nm, etch_depth, n_grates, params, working_dir):

bounds = [(-4,3), #< Starting position (in um)
 (0,0.05), #< Scaling parameter R
 (1.5,3), #< Parameter a
 (0,2)] #< Parameter b

 # geometry = ParameterizedGeometry(func = etched_grating, 
 # initial_params = params, 
 # bounds = bounds, 
 # dx = 1e-5)
 
 geometry = FunctionDefinedPolygon(func = grating_params_pos, 
 initial_params = params, 
 bounds = bounds, 
 z = 0.0, 
 depth = 110e-9, 
 eps_out = indexSiO2 ** 2, 
 eps_in = indexSi ** 2, 
 edge_precision = 5, 
 dx = 1e-5)

Materials : The permittivity of the materials included in the simulation (the material making the optimizable geometry (core) and the material surrounding it (cladding)) for the desired simulation wavelength should be passed to the optimizer when defining the geometry by calling the function FunctionDefinedPolygon (here they are presented as refractive index (n) squared). The refractive index (not permittivity) of the waveguides (which is the same as the optimizable geometry) should also be defined in the base simulation setup script.

Grating geometry : We used circular shape gratings for 3D simulations and GDS export. One can use more sophisticated grating shapes similar to focusing grating couplers .

Taking the model further

Information and tips for users that want to further customize the model

Taper optimization : The 3-D coupler model uses a linear taper section to connect to the integrated waveguide at the start of the grating. The coupling efficiency can also be improved by optimizing the taper shape (see SOI taper design).

High efficiency grating couplers : Couplers with an efficiency higher than 90% over a large bandwidth have been designed with FDTD using more sophisticated gratings and a mixed 2-D/3-D optimization strategy (see references).

Additional resources

Additional documentation, examples, and training material

Related publications

1. R. Marchetti, C. Lacava, A. Khokhar, X. Chen, I. Cristiani, D. J. Richardson, G. T. Reed, P. Petropoulos and P. Minzioni, “High-efficiency grating-couplers: demonstration of a new design strategy,” Scientific Reports 7, Article number: 16670, 2017. ( https://www.nature.com/articles/s41598-017-16505-z );
2. Neil V. Sapra et.al, “Inverse design and demonstration of broadband grating couplers” IEEE Journal of Selected Topics in Quantum Electronics ( Volume: 25, Issue: 3, May-June 2019 )
3. D. Taillaert, F. Van Laere, M. Ayre, W. Bogaerts, D. Van Thourhout, P. Bienstman and R. Baets, “Grating Couplers for Coupling between Optical Fibers and Nanophotonic Waveguides,” Japanese Journal of Applied Physics, vol. 45, no. 8a, pp. 6071-6077, 2006.
4. T. Watanabe, M. Ayata, U. Koch, Y. Fedoryshyn and J. Leuthold, “Perpendicular Grating Coupler Based on a Blazed Antiback-Reflection Structure,” Journal of Lightwave Technology, vol. 35, no. 21, pp. 4663- 4669, 2017.

See also

• Photonic Inverse Design Overview - Python API
• Getting Started with lumopt - python API
• Optimizable Geometry - Python API
• Grating coupler
• Broadband grating coupler
• MATLAB-driven grating coupler optimization
• SOI taper design
• Spot size converter
• Python API
• S-parameter simulator (SPS)
• S-parameter matrix sweep feature
• S-parameter file formats

Related Ansys Innovation 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 of 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 modeling 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 optimization module can be used to perform a fully automated optimization.

S-parameter extraction from component level simulations of passive photonic devices

The S-parameter matrix formalism is a common approach to build compact models of photonic devices to be used in circuit-level simulations . Assuming the response of the device to optical signals is linear, it can be modeled by a network (black box) with multiple network ports, where each of them receives an incoming signal and scatters or reflects an outgoing signal.

We assume the device has N physical ports (which typically are input/output waveguide channels) and each of these supports a certain number of modes of interest. Each element of the S-parameter matrix is the ratio between the complex envelopes of the optical modes in two physical ports. For example, to model the behavior of the fundamental TE and TM modes at the input/output waveguide channels of a 1x2 MMI we need a 6x6 S-parameter matrix, given that there are two modes at each of the three physical ports. Each column of the S-parameter matrix corresponds to a fixed input physical port and mode and each row corresponds to a fixed output physical port and mode.

The approach to extract the S-parameter matrix of a photonic device depends on the solver being used. In both FDTD and EME solvers, port objects can be set up at each of the physical ports. The main differences between the solvers for the S-parameter extraction are:

• FDTD injects one mode in one input port object per simulation, while EME can solve for all the elements of the S-matrix in one simulation;
• FDTD is broadband so many frequency points can be calculated in one simulation, while results from EME are for a single frequency per simulation.

Therefore, in both solvers, it is necessary to run sweeps to obtain the full S-parameter matrix at multiple frequency values, either by sweeping over the input physical ports and modes or over the frequency points.

Use the examples in the Application Gallery to determine the appropriate solver for your component. This section describes the recommended approach for S-parameter extraction and export to INTERCONNECT once you have determined the appropriate solver.

In this section we describe how to:

How to extract S-parameters (FDTD)

How to extract S-parameters (EME)

S-parameter extraction for group delay option in INTERCONNECT

Some photonic circuits rely on controlling the phase of the optical signal traveling through the different components; therefore, the phase of the S-parameters can be very important. There are a couple of challenges to capture the phase correctly in the Optical N-port S-parameter primitive in INTERCONNECT:

• A spectral sampling of the phase in the S-parameter data: Typically, the variation of the phase as a function of frequency is mostly linear with a slope that is proportional to the group delay (the delay of the optical signal between two ports). Therefore, in long devices, the phase can change very quickly with frequency, and a fine sampling in frequency would be necessary to avoid phase changes greater than \( 2 \pi \) so that the phase has the correct slope after unwrapping it. This is important for both frequency- and time-domain simulations in INTERCONNECT.
• Capturing the correct group delay and frequency dependence of the S-parameters in INTERCONNECT time-domain simulations: For time-domain simulations, the most common type of filter used in the Optical N-port S-parameter primitive is the finite impulse response (FIR) one, which has multiple options for the estimation of taps. The choice depends on what is the most significant aspect of the element response to be captured in the circuit simulations. In elements where the phase is important, typically it is necessary to capture both the group delay and the frequency dependence of the S-parameters correctly.

The “group delay” option for the “number of taps estimation” in the digital filter settings of the Optical N-port S-parameter primitive provides solutions to both challenges:

• The group delay is captured correctly without requiring a fine spectral sampling of the S-parameters by using group delay estimates provided by the user for each S-parameter matrix element at the center frequency/wavelength.
• In most cases, the FIR filter in time-domain simulations can capture the group delay and frequency dependence of the element response more accurately with the "group delay" option, compared to other options for tap estimation. With the "group delay" option, the number of taps available to the filter depends on the sample rate of the INTERCONNECT simulation and the group delay estimate provided by the user: for a given sample rate, there are more taps available for larger group delays. Since most photonic devices are long compared to the wavelength of the signal, typically the number of available taps is not a concern. Only in cases where the device is relatively short, it might be necessary to increase the sample rate.

Compared to the regular file format for S-parameter files in INTERCONNECT, the group delay option requires a slightly different format , which includes the group delay estimate in the header of each S-parameter data block; furthermore, the phase provided in the file needs to be adjusted to remove the linear contribution associated with the group delay estimate. The previous section describes how this type of file can be generated when using FDTD or EME for the S-parameter extraction. For more information on the group delay option see the discussion here .