Inverse design of an apodized grating coupler through shape optimization

A grating coupler (GC) is a key photonic device used to couple light between optical fibers and on-chip waveguides by diffracting light at a specific angle. These devices often employ apodization, a technique that modifies the grating tooth parameters along the device length to improve efficiency and minimize reflection losses.

In this notebook, we will use gradient-based optimization to determine the optimal grating tooth distribution. We start from a uniform grating coupler and modify the tooth and gap sizes simultaneously using gradient based optimization. Feature size constraints on the parameters enable the device to remain fabricable over the course of optimization. In just a handful of iterations, we achieve a design that outperforms an apodized grating coupler designed using semi-analytical methods.

While the simulations are performed using a 2D cross section, the same approach can be modified to 3D without much further complication.

If you are interested in other inverse design examples using tidy3d, you can find many of them here.

If you are unfamiliar with inverse design, we also recommend our intro to inverse design tutorials and our primer on automatic differentiation with tidy3d.

import tidy3d as td
import tidy3d.web as web

# as we are using autograd to perform gradient calculations, we need to import numpy from autograd's wrapper
import autograd as ag
import autograd.numpy as np

# we will use optax to perform the optimization using gradient descent
import optax

import matplotlib.pyplot as plt

Basic Parameters

First, we define the basic parameters used to define a uniform grating coupler for operation primarily in the C-band.

lda0 = 1.55  # central wavelength
n_ldas = 101 # number of wavelength points
ldas = np.linspace(1.5, 1.6, n_ldas)  # wavelength range

freq0 = td.C_0 / lda0  # central frequency
freqs = td.C_0 / ldas  # frequency range

num_freqs_objective = 5
ldas_objective =  np.linspace(1.5, 1.6, num_freqs_objective)
freqs_objective = td.C_0 / ldas_objective

fwidth = 0.5 * (np.max(freqs) - np.min(freqs))  # width of the source frequency range

The material platform is the common silicon on insulator (SOI) with a 260 nm silicon thickness and 2 µm BOX thickness. The top cladding layer is 680 nm oxide. The grating will be partially etched with an etching depth of 160 nm.

# define materials from the material library
si = td.material_library["cSi"]["Palik_Lossless"]
sio2 = td.material_library["SiO2"]["Palik_Lossless"]

t_si = 0.26 # thickness of the silicon layer
etch_depth = 0.16 # etching depth
t_tox = 0.68 # top oxide layer thickness
t_box = 2 # bottom oxide layer thickness

We aim to design the GC for the standard single-mode fiber (SMF) with a mode field diameter (MFD) of 10.8 µm. The fiber is tilted at a 14.5 degree angle.

theta = np.deg2rad(14.5) # fiber tilt angle
mfd = 10.8 # mode field diameter
source_x = 4.5 # x position of the fiber

N = 16 # number of grating teeth to create
inf_eff = 1e3 # effective infinity
buffer = 1.1*lda0 # buffer spacing to pad the simulation domain

Next we compute the uniform tooth parameters (periodicity and fill fraction) to maximize efficiency.

neff_unetch = 2.9 # effective index of the slab mode of the unetched waveguide
neff_etch = 2.2 # effective index of the slab mode of the etched waveguide
n_c = 1.44 # refractive index of the cladding (SiO2)
theta_c = np.sin(theta)/n_c # incident angle in the cladding

def get_periodicity(fill_fraction: float) -> float:
    """periodicity (bragg condition) as function of fill fraction, angle, wavelength, and etching parameters."""
    return lda0/(fill_fraction*neff_unetch + (1-fill_fraction)*neff_etch - n_c*np.sin(theta_c))

f0 = 0.80
periodicity = get_periodicity(f0)
p_list = N * [periodicity]
f_list = N * [f0]

After obtaining \(p_{\text{i}}\) and \(f_{\text{i}}\), we will write a function to return the width of each air gap and each tooth. We will use these gap widths as the parameters for our design. The advantage is that we can impose bounds on these values to ensure the minimal feature size is above the fabrication constraint.

Here we aim to have a minimal feature size above 60 nm. The current design has a minimal feature size right above it.

def get_widths(p_list, f_list):
    # calculate the widths of air gaps and silicon teeth
    return np.array([item for p, f in zip(p_list, f_list) for item in [p*(1-f), p*f]])

widths = get_widths(p_list, f_list)

l_grating = np.sum(widths) + 3 * lda0 # total length of the GC

# print the current minimal feature size
print(f"The minimal feature size is {1e3*np.min(widths):.2f} nm.")

The minimal feature size is 123.46 nm.

To ensure the minimal feature size is maintained during the optimization, we will use a design parameter in the range from -inf to inf and project it to a tanh function that is bounded by the minimal and maximal feature sizes.

min_width = 0.08 # minimal feature size we want to maintain
max_width = 1 # maximal feature size

# function to project a design parameter between -inf to inf to between min_width and max_width
def project(x):
    return 0.5*(max_width - min_width)*np.tanh(x)+0.5*(max_width + min_width)

# function to inversely project a design parameter between min_width and max_width to between -inf to inf
def inverse_project(y):
    return np.arctanh((2 * (y - 0.5 * (max_width + min_width))) / (max_width - min_width))

# project the widths to parameters between -inf to inf
params0 = inverse_project(widths)

Define Static Components of the Simulation

Next we will define the components that don’t change during the optimization. These include the cladding layer, the output waveguide, the unetched layer, the bottom oxide layer, and the silicon substrate.

mask = td.Structure(
    geometry=td.Box(size=(td.inf, td.inf, td.inf), center=(0,0,0)),
    medium=td.Medium(),
)

# create the top oxide layer
tox = td.Structure(
    geometry=td.Box.from_bounds(
        rmin=(-inf_eff, -inf_eff, 0),
        rmax=(inf_eff, inf_eff, t_tox)),
    medium=sio2
)

# create the output waveguide
slab_waveguide = td.Structure(
    geometry=td.Box.from_bounds(
        rmin=(-inf_eff, -inf_eff, 0),
        rmax=(0, inf_eff, t_si)),
    medium=si
)

# create the unetched waveguide
unetched_waveguide = td.Structure(
    geometry=td.Box.from_bounds(
        rmin=(0, -inf_eff, 0),
        rmax=(inf_eff, inf_eff, t_si-etch_depth)),
    medium=si
)

# create the bottom oxide layer
box = td.Structure(
    geometry=td.Box.from_bounds(
        rmin=(-inf_eff, -inf_eff, -t_box),
        rmax=(inf_eff, inf_eff, 0)),
    medium=sio2
)

# create the silicon substrate layer
substrate = td.Structure(
    geometry=td.Box.from_bounds(
        rmin=(-inf_eff, -inf_eff, -inf_eff),
        rmax=(inf_eff, inf_eff, -t_box)),
    medium=si
)

The source and monitor are not changed either. We use a GaussianBeam to represent the incident fiber mode. A ModeMonitor is placed at the output waveguide to measure the coupling efficiency.

gap = 1 # gap size between the source plane and the top surface of the cladding

# define a gaussian beam source
source = td.GaussianBeam(
    size=(2*mfd, td.inf, 0),
    center=[source_x, 0, t_tox+gap],
    source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
    angle_theta=theta,
    direction="-",
    waist_radius=mfd/2,
    pol_angle=np.pi/2 # 90 degree polarization angle for TE polarization
)

# define a mode monitor
mode_monitor = td.ModeMonitor(
    center=(-buffer/2, 0, t_si/2),
    size=(0, td.inf, 6 * t_si),
    freqs=freqs,
    mode_spec=td.ModeSpec(num_modes=1, target_neff=3),
    name="mode",
)

Simulate the Designed Apodized GC

To simulate the designed GC and facility future optimization, we define a function that takes in an array of design parameters and returns a `Simulation <>`__ object. Note again that each design parameter has a range from -inf to inf. It will then be projected to a range between min_width and max_width.

def make_2d_sim(design_parameters):

    # calculate the widths of air gaps and silicon teeth
    widths_si = project(design_parameters[1::2])
    widths_air = project(design_parameters[::2])

    # initialize the center and size of each silicon teeth
    center = 0
    size = 0

    # create the grating geometries from the given widths
    gratings = 0
    for width_si, width_air in zip(widths_si, widths_air):
        center += width_air+width_si/2
        size = width_si
        gratings += td.Box(center=(center, 0, t_si-etch_depth/2), size=(size, td.inf, etch_depth))
        center += width_si/2

    # create the grating structure
    gratings = td.Structure(geometry=gratings, medium=si)

    # create a box to represent the simulation domain box
    sim_box = td.Box.from_bounds(rmin=(-buffer, 0, -t_box-buffer/2),
                                rmax=(l_grating+buffer,0,t_si+buffer))

    run_time = 1e-12  # simulation run time

    # construct simulation
    sim = td.Simulation(
        center=sim_box.center,
        size=sim_box.size,
        grid_spec=td.GridSpec.auto(min_steps_per_wvl=30, wavelength=lda0), # use a fine grid to ensure the small features are well resolved
        structures=[mask, tox, gratings, unetched_waveguide, slab_waveguide, box, substrate],
        sources=[source],
        monitors=[mode_monitor],
        medium=sio2,
        run_time=run_time,
        boundary_spec=td.BoundarySpec(
            x=td.Boundary.pml(),
            y=td.Boundary.periodic(),  # set the boundary to periodic in y since it's a 2D simulation
            z=td.Boundary.pml()
        ),
    )

    return sim

Create the simulation for the designed apodized GC and visualize it.

sim0 = make_2d_sim(params0)
sim0.plot_eps(y=0, freq=freq0)
plt.show()

Also visualize the grid to ensure it resolves the small features well.

ax = sim0.plot(y=0)
sim0.plot_grid(y=0, ax=ax)
ax.set_xlim(-0.1, 1)
ax.set_ylim(0, 0.3)
plt.show()

Run the simulation.

sim_data0 = web.run(simulation=sim0, task_name="initial design")

13:09:46 EDT Created task 'initial design' with task_id
             'fdve-b1542448-f18d-4fbf-a4de-73f74814d5ab' and task_type 'FDTD'.

             View task using web UI at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-b1542448-f18
             d-4fbf-a4de-73f74814d5ab'.

13:09:48 EDT status = success

             loading simulation from simulation_data.hdf5

Result Visualization

After the simulation is done, we can extract the coupling efficiency from the mode monitor data. In addition, we also calculate two important metrics: the maximum coupling efficiency and the 1dB bandwidth. The maximal coupling efficiency is about 2.7dB and the 1dB bandwidth is about 36 nm.

# function to calculate the coupling efficiency from simulation data
def compute_ce(sim_data):
    amp = sim_data["mode"].amps.sel(mode_index=0, direction="-").values
    ce = np.abs(amp) ** 2 # transmission to the top waveguide
    return ce

def dB(ce):
    return 10*np.log10(ce)

# calculate the coupling efficiency for the designed GC
ce0 = compute_ce(sim_data0)

# plot the coupling efficiency
plt.plot(ldas, dB(ce0), c='red', linewidth=2)
plt.xlim(min(ldas), max(ldas))
plt.ylim(-10, 0)
plt.xlabel("Wavelength (µm)")
plt.ylabel("Coupling efficiency (dB)")
plt.grid()
plt.show()

# function to calculate the 1dB bandwidth from the coupling efficiency
def bandwidth(ldas, ce):
    max_ce = np.max(dB(ce))

    threshold = max_ce - 1

    within_1db = np.where(dB(ce) >= threshold)[0]

    lambda_min = ldas[within_1db[0]]
    lambda_max = ldas[within_1db[-1]]
    bandwidth = lambda_max - lambda_min
    return 1e3*bandwidth

# print the maximum coupling efficiency and 1dB bandwidth
print(f"The 1dB bandwidth is {bandwidth(ldas, ce0):.1f} nm")
print(f"The maximum coupling efficiency is {np.max(dB(ce0)):.2f} dB.")

The 1dB bandwidth is 36.0 nm
The maximum coupling efficiency is -2.74 dB.

Optimize GC with Inverse Design to Maximize the Efficiency and Bandwidth

In this part, we will apply inverse design to maximize the coupling efficiency of our coupling over the wavelength range from 1500 to 1600 nm.

To do this, we need to express our objective as a function of our parameters returning a single float value to maximize. This function will involve constructing, running, and postprocessing our simulation.

def J(design_parameters: np.ndarray) -> float:

    sim = make_2d_sim(design_parameters)
    sim_data = web.run(sim, task_name="GC_invdes", local_gradient = False, verbose=False)
    ce =  np.sum(compute_ce(sim_data)) / n_ldas

    return ce

dJ = ag.value_and_grad(J)

Before running the optimization, we can check the value and gradient of the initial design to ensure the gradient tracking is working properly.

%%time

val, grad = dJ(params0)
print(val)
print(grad)

0.3389903919270378
[-1.19602694e-02  6.20778868e-01  2.64248243e-02  5.76042908e-01
 -1.76322293e-03  1.29420074e-01 -2.18638124e-03 -2.12022992e-02
  6.51520927e-04  5.86568764e-02  1.14255341e-02  5.56881759e-02
  1.50206311e-02 -2.24280572e-02  3.85486469e-03 -7.80827528e-02
  7.60361154e-03 -7.04035656e-02  1.07565584e-03 -8.46643660e-02
  4.36607376e-03 -7.72914850e-02 -1.97607690e-03 -7.11080918e-02
  2.49066963e-04 -4.97352013e-02 -3.88170239e-04 -3.69572558e-02
 -6.04949639e-04 -2.51740880e-02 -5.48443825e-04 -1.45726593e-02]
CPU times: user 1e+03 ms, sys: 36.2 ms, total: 1.04 s
Wall time: 4.47 s

We will use the Adam optimizer to adjust the parameters for a number of iterations until we achieve a good performance.

%%time

# hyperparameters
num_steps = 25
learning_rate = 0.05

# initialize adam optimizer with starting parameters
params = params0
optimizer = optax.adam(learning_rate=learning_rate)
opt_state = optimizer.init(params)

# store history
J_history = []
params_history = []

for i in range(num_steps):
    # compute gradient and current objective function value
    value, gradient = dJ(params)

    # outputs
    print(f"step = {i + 1}")
    print(f"\tJ = {value:.3e}")
    print(f"\tgrad_norm = {np.linalg.norm(gradient):.4e}")

    # compute and apply updates to the optimizer based on gradient (-1 sign to maximize obj_fn)
    updates, opt_state = optimizer.update(-gradient, opt_state, params)
    params = optax.apply_updates(params, updates)
    params = np.array(params)

    # save history
    J_history.append(value)
    params_history.append(params)

step = 1
        J = 3.390e-01
        grad_norm = 8.8129e-01
step = 2
        J = 3.173e-01
        grad_norm = 9.1435e-01
step = 3
        J = 3.793e-01
        grad_norm = 5.7908e-01
step = 4
        J = 3.860e-01
        grad_norm = 7.1871e-01
step = 5
        J = 3.886e-01
        grad_norm = 5.8824e-01
step = 6
        J = 3.916e-01
        grad_norm = 6.4806e-01
step = 7
        J = 4.034e-01
        grad_norm = 5.9342e-01
step = 8
        J = 4.205e-01
        grad_norm = 4.2752e-01
step = 9
        J = 4.279e-01
        grad_norm = 3.6594e-01
step = 10
        J = 4.307e-01
        grad_norm = 5.1920e-01
step = 11
        J = 4.384e-01
        grad_norm = 4.4993e-01
step = 12
        J = 4.441e-01
        grad_norm = 4.7875e-01
step = 13
        J = 4.522e-01
        grad_norm = 3.3082e-01
step = 14
        J = 4.552e-01
        grad_norm = 3.7186e-01
step = 15
        J = 4.581e-01
        grad_norm = 3.1254e-01
step = 16
        J = 4.603e-01
        grad_norm = 3.4604e-01
step = 17
        J = 4.659e-01
        grad_norm = 2.4904e-01
step = 18
        J = 4.690e-01
        grad_norm = 2.0352e-01
step = 19
        J = 4.696e-01
        grad_norm = 3.1720e-01
step = 20
        J = 4.730e-01
        grad_norm = 1.9975e-01
step = 21
        J = 4.745e-01
        grad_norm = 1.8170e-01
step = 22
        J = 4.750e-01
        grad_norm = 2.7087e-01
step = 23
        J = 4.787e-01
        grad_norm = 1.5339e-01
step = 24
        J = 4.800e-01
        grad_norm = 1.7580e-01
step = 25
        J = 4.802e-01
        grad_norm = 1.5796e-01
CPU times: user 15.8 s, sys: 997 ms, total: 16.8 s
Wall time: 2min 46s

Plot the objective function (average coupling efficiency) as a function of the number of iterations. A steady increase of the average coupling efficiency is observed.

plt.plot(J_history, c='black')
plt.xlabel("Iterations")
plt.ylabel("Objective function (average coupling efficiency)")
plt.show()

Design the Apodized GC

For comparison with more established methods, we’ll simulate a device designed using the approach introduced in the reference paper. First of all, the Bragg condition is given by

\[p_{\text{i}} = \frac{\lambda_\text{c}}{(n_{\text{eff,i}} - n_{\text{c}}\sin \theta_{\text{c}})},\]

where \(p_{\text{i}}\) is the local period of the grating, \(\lambda_\text{c}\) is the central wavelength, \(n_{\text{c}}\) is the refractive index of the cladding, \(\theta_{\text{c}}\) is the incident angle in the cladding (~ 10 degrees from Snell’s law), and

\[n_{\text{eff,i}} = f_\text{i} \times n_\text{neff,nonetch} + (1-f_\text{i}) \times n_\text{neff,etch}.\]

Here \(f_\text{i}\) is the local filling fraction, \(n_\text{neff,unetch}\) is the effective index of the slab mode of the unetched waveguide (260 nm thick), and \(n_\text{neff,etch}\) is the effective index of the slab mode of the etched waveguide (100 nm thick). \(n_\text{neff,unetch}\) and \(n_\text{neff,etch}\) are determined to be about 2.9 and 2.2 by mode analysis (not shown). Furthermore, we follow a linear apodization, namely

\[f_\text{i} = f_\text{0} - R \times x,\]

where \(f_\text{0}\) is the filling fraction of the first period, \(R\) is the linear apodization factor and \(x\) is the position of each tooth from the starting point of the grating. By using an initial local filling fraction \(f_\text{0}\)=0.9 and apodization factor \(R\)=0.0025 µm\(^{-1}\) we can recursively calculate \(p_{\text{i}}\) and \(f_\text{i}\) for \(i=0, 1, ..., N-1\) using the equations above. The result is a grating with apodized pitch length and filling fraction.

For more details on the design principles, refer to the reference paper.

# and then the apodized design parameters to compare to later

R = 0.025 # linear apodization factor

p_list_apodized = [] # list to store all local periodicities
f_list_apodized = [] # list to store all local filling fractions

f0_apodized = 0.86 # filling fraction of the first teeth
f_current = f0_apodized # variable to store the current local filling fraction

# recursively calculate all local periodicities and filling fractions
for i in range(N):
    p_current = get_periodicity(f_current)
    p_list_apodized.append(p_current)
    f_list_apodized.append(f_current)
    f_current = f0_apodized - R*sum(p_list_apodized)

widths_apodized = get_widths(p_list_apodized, f_list_apodized)
params_apodized = inverse_project(widths_apodized)
sim_apodized = make_2d_sim(params_apodized)

Let’s plot all of the simulations together to see how their features compare.

params_opt = params_history[-1]
sim_opt = make_2d_sim(params_opt)

f, (ax1, ax2, ax3) = plt.subplots(3, 1, tight_layout=True, figsize=(10,6))

sim0.plot_eps(y=0, freq=freq0, ax=ax1)
sim_apodized.plot_eps(y=0, freq=freq0, ax=ax2)
sim_opt.plot_eps(y=0, freq=freq0, ax=ax3)

ax1.set_title('uniform')
ax2.set_title('apodized')
ax3.set_title('optimized')


plt.show()

Next, we’ll run each of the designs to be able to compare the coupling efficiencies directly.

sim_data_apodized = web.run(simulation=sim_apodized, task_name="apodized design")
sim_data_opt = web.run(simulation=sim_opt, task_name="optimized design")

13:12:39 EDT Created task 'apodized design' with task_id
             'fdve-5a56a3f9-a80e-442d-8c39-021abad4fd8c' and task_type 'FDTD'.

             View task using web UI at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-5a56a3f9-a80
             e-442d-8c39-021abad4fd8c'.

13:12:40 EDT status = success

13:12:41 EDT loading simulation from simulation_data.hdf5

             Created task 'optimized design' with task_id
             'fdve-43d33b82-22dc-495d-8fba-1c9b51addfe1' and task_type 'FDTD'.

             View task using web UI at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-43d33b82-22d
             c-495d-8fba-1c9b51addfe1'.

13:12:42 EDT status = success

13:12:43 EDT loading simulation from simulation_data.hdf5

ce_opt = compute_ce(sim_data_opt)
ce_apodized = compute_ce(sim_data_apodized)

plt.plot(ldas, dB(ce0), c='red', linewidth=2, label="Initial design")
plt.plot(ldas, dB(ce_apodized), c='orange', linewidth=2, label="Apodized design")
plt.plot(ldas, dB(ce_opt), c='blue', linewidth=2, label="Optimized design")
plt.xlim(min(ldas), max(ldas))
plt.ylim(-10, 0)
plt.xlabel("Wavelength (µm)")
plt.ylabel("Coupling efficiency (dB)")
plt.legend()
plt.grid()
plt.show()

We see that the optimized and apodized designs both exceed that of the uniform with -2dB efficiency. However, the optimized design far ouperforms in terms of bandwidth. We can also confirm that the minimal feature size is maintained above 85 nm.

print('Uniform:')
print(f"  The 1dB bandwidth is {bandwidth(ldas, ce0):.1f} nm")
print(f"  The maximum coupling efficiency is {np.max(dB(ce0)):.2f} dB.")
print(f"  The minimal feature size is {1e3*np.min(project(params0)):.2f} nm.")
print('')
print('Apodized:')
print(f"  The 1dB bandwidth is {bandwidth(ldas, ce_apodized):.1f} nm")
print(f"  The maximum coupling efficiency is {np.max(dB(ce_apodized)):.2f} dB.")
print(f"  The minimal feature size is {1e3*np.min(project(params_apodized)):.2f} nm.")
print('')
print('Optimized:')
print(f"  The 1dB bandwidth is {bandwidth(ldas, ce_opt):.1f} nm")
print(f"  The maximum coupling efficiency is {np.max(dB(ce_opt)):.2f} dB.")
print(f"  The minimal feature size is {1e3*np.min(project(params_opt)):.2f} nm.")

Uniform:
  The 1dB bandwidth is 36.0 nm
  The maximum coupling efficiency is -2.74 dB.
  The minimal feature size is 123.46 nm.

Apodized:
  The 1dB bandwidth is 43.0 nm
  The maximum coupling efficiency is -1.99 dB.
  The minimal feature size is 85.00 nm.

Optimized:
  The 1dB bandwidth is 50.0 nm
  The maximum coupling efficiency is -1.98 dB.
  The minimal feature size is 86.11 nm.

Take the Model Further

In this notebook, we only perform 2D simulations. A good next step would be to convert the 2D design to a 3D linear GC or focusing GC. The coupling efficiency of the 3D GC is likely slightly lower than that in 2D. Inverse design can then be applied to fine-tune it again. To minimize the footprint, one can potentially also apply shape optimization to design a compact and low-loss taper section, as demonstrated in this example.

Further visualizations

In case you would like to gain more insight into how your device has changed over optimization, the following are a few post-processing steps to make animations and visualize the grating teeth locations.

import matplotlib.animation as animation
from IPython.display import HTML

fig, ax = plt.subplots(1, 1, tight_layout=True, figsize=(10, 3))

sim0 = make_2d_sim(params_history[0])

def plot_fn(params, ax=ax):
    sim = make_2d_sim(params)
    ax = sim.plot_eps(y=0.0, ax=ax, freq=freq0)
    ax.set_aspect("equal")

def animate(i):
    plot_fn(params_history[i])

# create animation
ani = animation.FuncAnimation(fig, animate, frames=len(params_history))
plt.close()

td.config.logging_level="ERROR"

# display the animation (press "play" to start)
HTML(ani.to_jshtml())

Once Loop Reflect

<Figure size 640x480 with 0 Axes>

# uncomment to save to file
# ani.save("/users/twhughes/Desktop/apodized_grating_coupler.gif", fps=10)

<Figure size 640x480 with 0 Axes>

history_min = [[] for _ in range(len(params_history))]
history_max = [[] for _ in range(len(params_history))]

for i, params in enumerate(params_history):
    sim = make_2d_sim(params)
    for geometry in sim.structures[2].geometry.geometries:
        (xmin, _, _), (xmax, _, _) = geometry.bounds
        history_min[i].append(xmin)
        history_max[i].append(xmax)
        # print(i, j)

history_min = np.array(history_min).T
history_max = np.array(history_max).T

from matplotlib.pyplot import  fill_between

xs = np.arange(len(params_history))

_, ax = plt.subplots()

for y1, y2 in zip(history_min, history_max):
    fill_between(xs, y1, y2)

plt.xlabel('iteration number')
plt.ylabel('x position')
plt.title('tooth position and width over optimization')

plt.show()