Resonator benchmark (COMSOL)

See on github, run on colab, or just follow along with the output below.

In this example, we reproduce the findings of Zhang et al. (2018), which is linked here.

This notebook was originally developed and written by Romil Audhkhasi (USC).

The paper investigates the resonances of silicon structures by measuring their transmission spectrum under varying geometric parameters.

The paper uses a finite element solver (COMSOL), which matches the result from Tidy3D.

diagram

(Citation: Opt. Lett. 43, 1842-1845 (2018). With permission from the Optical Society)

To do this calculation, we use a broadband pulse and frequency monitor to measure the flux on the opposite side of the structure.

[1]:
# get the most recent version of tidy3d
!pip install -q --upgrade tidy3d

# make sure notebook plots inline
%matplotlib inline

# standard python imports
import numpy as np
import matplotlib.pyplot as plt

# tidy3D import
import tidy3d as td
from tidy3d import web

Set Up Simulation

[2]:
nm = 1e-3

# define the frequencies we want to measure
Nfreq = 1000
wavelengths = nm * np.linspace(1050, 1400, Nfreq)
freqs = td.constants.C_0 / wavelengths

# define the frequency center and width of our pulse
freq0 = freqs[len(freqs)//2]
freqw = freqs[0] - freqs[-1]

# Define material properties
n_SiO2 = 1.46
n_Si = 3.52
SiO2 = td.Medium(epsilon=n_SiO2**2)
Si = td.Medium(epsilon=n_Si**2)
[3]:
# space between resonators and source
spc = 1.5

# geometric parameters
Px = Py = P = 650 * nm  # periodicity in x and y
t = 260 * nm            # thickness of silcon
g = 80 * nm             # gap size
L = 480 * nm            # length in x
w_sum = 400 * nm        # sum of lengths in y

# resolution (should be commensurate with periodicity)
dl = P / 32
dz = 20 * nm

# computes widths in y (w1 and w2) given the difference in lengths in y and the sum of lengths
def calc_ws(delta):
    """ delta is a tunable parameter used to break symmetry.
        w_sum = w1 + w2
        delta = w1 - w2
        w_sum + delta = 2 * w1
    """
    w1 = (w_sum + delta) / 2
    w2 = w_sum - w1
    return w1, w2
[4]:
# total size in z and [x,y,z]
Lz = spc + t + t + spc
sim_size = [Px, Py, Lz]

# sio2 substrate
substrate = td.Box(
    center=[0, 0, -Lz/2],
    size=[td.inf, td.inf, 2*(spc+t)],
    material=SiO2,
    name='substrate'
)

# creates a list of structures given a value of 'delta'
def geometry(delta):
    w1, w2 = calc_ws(delta)
    center_y = (w1 - w2) / 2.0

    cell1 = td.Box(
        center=[0, center_y + (g + w1)/2., t/2.],
        size=[L, w1, t],
        material=Si,
        name='cell1'
    )

    cell2 = td.Box(
        center=[0, center_y - (g + w2)/2., t/2.],
        size=[L, w2, t],
        material=Si,
        name='cell2'
    )

    return [substrate, cell1, cell2]

[5]:
# time dependence of source
gaussian = td.GaussianPulse(freq0, freqw)

# plane wave source
source = td.PlaneWave(
            source_time=gaussian,
            injection_axis='-z',
            position= Lz/2 - spc + 2*dl,
            polarization='x')

# Simulation run time.  Note you need to run a long time to calculate high Q resonances.
run_time = 7e-12
[6]:
# monitor fields on other side of structure (substrate side) at range of frequencies
monitor = td.FreqMonitor(
        center=[0., 0., -Lz/2 + spc - 2 * dl],
        size=[td.inf, td.inf, 0],
        freqs=freqs,
        store=['flux'],
        name='transmitted_fields')

Define Case Studies

Here we define the three simulations to run

  • With no resonators (normalization)

  • With symmetric (delta = 0) resonators

  • With asymmetric (delta != 0) resonators

[7]:
# normalizing run (no Si) to get baseline transmission vs freq
# can be run for shorter time as there are no resonances
sim_empty = td.Simulation(size=sim_size,
                    mesh_step=[dl, dl, dz],
                    structures=[substrate],
                    sources=[source],
                    monitors=[monitor],
                    run_time=run_time/10,
                    pml_layers=[0,0,15])

# run with delta = 0
sim_d0 = td.Simulation(size=sim_size,
                    mesh_step=[dl, dl, dz],
                    structures=geometry(0),
                    sources=[source],
                    monitors=[monitor],
                    run_time=run_time,
                    pml_layers=[0,0,15])

# run with delta = 20nm
sim_d20 = td.Simulation(size=sim_size,
                    mesh_step=[dl, dl, dz],
                    structures=geometry(20 * nm),
                    sources=[source],
                    monitors=[monitor],
                    run_time=run_time,
                    pml_layers=[0,0,15])
Initializing simulation...
Mesh step (micron): [2.03e-02, 2.03e-02, 2.00e-02].
Simulation domain in number of grid points: [32, 32, 206].
Total number of computational grid points: 2.11e+05.
Total number of time steps: 19987.
Estimated data size (GB) of monitor transmitted_fields: 0.0000.
Initializing simulation...
Mesh step (micron): [2.03e-02, 2.03e-02, 2.00e-02].
Simulation domain in number of grid points: [32, 32, 206].
Total number of computational grid points: 2.11e+05.
Total number of time steps: 199868.
Estimated data size (GB) of monitor transmitted_fields: 0.0000.
Initializing simulation...
Mesh step (micron): [2.03e-02, 2.03e-02, 2.00e-02].
Simulation domain in number of grid points: [32, 32, 206].
Total number of computational grid points: 2.11e+05.
Total number of time steps: 199868.
Estimated data size (GB) of monitor transmitted_fields: 0.0000.
[8]:
# Structure visualization in various planes

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14, 8))
sim_d0.viz_eps_2D(normal='x', position=0, ax=ax1)
sim_d0.viz_eps_2D(normal='y', position=g, ax=ax2)
sim_d0.viz_eps_2D(normal='z', position=0, ax=ax3)
plt.show()
../_images/examples_HighQ_Si_10_0.png

Run Simulations

[9]:
# define a function that runs the simulations on our server and computes the flux

def get_flux(sim, task_name='Si_resonator'):
    project = web.new_project(sim.export(), task_name=task_name)
    web.monitor_project(project['taskId'])
    print('Downloading results')
    web.download_results(project['taskId'], target_folder='output')
    sim.load_results('output/monitor_data.hdf5')
    with open('output/tidy3d.log', 'r') as f:
        print(f.read())
    return np.squeeze(sim.data(monitor)['flux'])

[10]:
# run all simulations, take about 2-3 minutes each with some download time
flux_empty = get_flux(sim_empty, task_name='normalization')
flux_delta0 = get_flux(sim_d0, task_name='Si-resonator-delta-0')
flux_delta20 = get_flux(sim_d20, task_name='Si-resonator-delta-20')
Uploading the json file...
Project 'normalization-16b6e1d68c55d7e8' status: success...

Downloading results
Applying source normalization to all frequency monitors using source index 0.
Simulation domain Nx, Ny, Nz: [32, 32, 206]
Applied symmetries: [0, 0, 0]
Number of computational grid points: 2.1094e+05.
Using subpixel averaging: True
Number of time steps: 19987
Automatic shutoff factor: 1.00e-05
Time step (s): 3.5023e-17

Compute source modes time (s):     0.0517
Compute monitor modes time (s):    0.0546

Rest of setup time (s):            0.0932

Starting solver...
- Time step    318 / time 1.11e-14s (  1 % done), field decay: 1.00e+00
- Time step    799 / time 2.80e-14s (  4 % done), field decay: 1.18e-05
- Time step   1598 / time 5.60e-14s (  8 % done), field decay: 4.34e-13
Field decay smaller than shutoff factor, exiting solver.

Solver time (s):                   1.6589
Post-processing time (s):          0.0047

Uploading the json file...
Project 'Si-resonator-delta-0-16b6e1e153d4b898' status: success...

Downloading results
Applying source normalization to all frequency monitors using source index 0.
Simulation domain Nx, Ny, Nz: [32, 32, 206]
Applied symmetries: [0, 0, 0]
Number of computational grid points: 2.1094e+05.
Using subpixel averaging: True
Number of time steps: 199868
Automatic shutoff factor: 1.00e-05
Time step (s): 3.5023e-17

Compute source modes time (s):     0.0600
Compute monitor modes time (s):    0.0772

Rest of setup time (s):            0.3107

Starting solver...
- Time step    318 / time 1.11e-14s (  0 % done), field decay: 1.00e+00
- Time step   7994 / time 2.80e-13s (  4 % done), field decay: 6.99e-03
- Time step  15989 / time 5.60e-13s (  8 % done), field decay: 5.25e-04
- Time step  23984 / time 8.40e-13s ( 12 % done), field decay: 5.74e-04
- Time step  31978 / time 1.12e-12s ( 16 % done), field decay: 3.27e-05
- Time step  39973 / time 1.40e-12s ( 20 % done), field decay: 1.33e-04
- Time step  47968 / time 1.68e-12s ( 24 % done), field decay: 1.53e-05
- Time step  55963 / time 1.96e-12s ( 28 % done), field decay: 2.24e-05
- Time step  63957 / time 2.24e-12s ( 32 % done), field decay: 9.00e-06
Field decay smaller than shutoff factor, exiting solver.

Solver time (s):                   16.4182
Post-processing time (s):          0.0058

Uploading the json file...
Project 'Si-resonator-delta-20-16b6e1f5210d2e18' status: success...

Downloading results
Applying source normalization to all frequency monitors using source index 0.
Simulation domain Nx, Ny, Nz: [32, 32, 206]
Applied symmetries: [0, 0, 0]
Number of computational grid points: 2.1094e+05.
Using subpixel averaging: True
Number of time steps: 199868
Automatic shutoff factor: 1.00e-05
Time step (s): 3.5023e-17

Compute source modes time (s):     0.0577
Compute monitor modes time (s):    0.0684

Rest of setup time (s):            0.3133

Starting solver...
- Time step    318 / time 1.11e-14s (  0 % done), field decay: 1.00e+00
- Time step   7994 / time 2.80e-13s (  4 % done), field decay: 7.52e-03
- Time step  15989 / time 5.60e-13s (  8 % done), field decay: 1.05e-02
- Time step  23984 / time 8.40e-13s ( 12 % done), field decay: 6.16e-03
- Time step  31978 / time 1.12e-12s ( 16 % done), field decay: 5.09e-03
- Time step  39973 / time 1.40e-12s ( 20 % done), field decay: 6.53e-03
- Time step  47968 / time 1.68e-12s ( 24 % done), field decay: 2.04e-03
- Time step  55963 / time 1.96e-12s ( 28 % done), field decay: 1.41e-03
- Time step  63957 / time 2.24e-12s ( 32 % done), field decay: 4.56e-03
- Time step  71952 / time 2.52e-12s ( 36 % done), field decay: 3.37e-03
- Time step  79947 / time 2.80e-12s ( 40 % done), field decay: 1.25e-03
- Time step  87941 / time 3.08e-12s ( 44 % done), field decay: 2.30e-03
- Time step  95936 / time 3.36e-12s ( 48 % done), field decay: 3.46e-03
- Time step 103931 / time 3.64e-12s ( 52 % done), field decay: 2.51e-03
- Time step 111926 / time 3.92e-12s ( 56 % done), field decay: 9.81e-04
- Time step 119920 / time 4.20e-12s ( 60 % done), field decay: 7.87e-04
- Time step 127915 / time 4.48e-12s ( 64 % done), field decay: 1.58e-03
- Time step 135910 / time 4.76e-12s ( 68 % done), field decay: 1.44e-03
- Time step 143904 / time 5.04e-12s ( 72 % done), field decay: 5.84e-04
- Time step 151899 / time 5.32e-12s ( 76 % done), field decay: 8.73e-04
- Time step 159894 / time 5.60e-12s ( 80 % done), field decay: 2.06e-03
- Time step 167889 / time 5.88e-12s ( 84 % done), field decay: 2.05e-03
- Time step 175883 / time 6.16e-12s ( 88 % done), field decay: 1.12e-03
- Time step 183878 / time 6.44e-12s ( 92 % done), field decay: 7.67e-04
- Time step 191873 / time 6.72e-12s ( 96 % done), field decay: 1.17e-03
- Time step 199867 / time 7.00e-12s (100 % done), field decay: 1.03e-03

Solver time (s):                   49.5821
Post-processing time (s):          0.0057

Get Results and Plot

[11]:
# plot normalizing run (transmission through substrate alone)
plt.plot(wavelengths / nm, np.abs(flux_empty), color='k')
plt.ylim([0,1])
plt.title('normalizing run (no Si)')
plt.xlabel('wavelength ($nm$)')
plt.ylabel('Transmission')
plt.show()
../_images/examples_HighQ_Si_15_0.png

The normalizing run computes the transmitted flux for an air -> SiO2 interface, which is just below unity due to some reflection.

While not technically necessary for this example, since this transmission can be computed analytically, it is often a good idea to run a normalizing run so you can accurately measure the change in output when the structure is added. For example, for multilayer structures, the normalizing run displays frequency dependence, which would make it prudent to include in the calculation.

[12]:
# plot transmission, compare to paper results, look similar
fig, ax = plt.subplots(1, 1, figsize=(6, 4.5))
plt.plot(wavelengths / nm, flux_delta0 / flux_empty, color='red', label='$\delta=0$')
plt.plot(wavelengths / nm, flux_delta20 / flux_empty, color='blue', label='$\delta=20~nm$')
plt.xlabel('wavelength ($nm$)')
plt.ylabel('Transmission')
plt.xlim([1050, 1400])
plt.ylim([0, 1])
plt.legend()
plt.show()
../_images/examples_HighQ_Si_17_0.png

Results Comparison

Compare this plot to published results computed using COMSOL (FEM):

diagram

(Citation: Opt. Lett. 43, 1842-1845 (2018). With permission from the Optical Society)

[ ]: