Resonator benchmark (Lumerical)

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

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

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

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

The paper uses a finite-difference time-domain (Lumerical), which matches the result from Tidy3D.

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]:
Nfreq = 1000
wavelengths = np.linspace(8, 12, Nfreq)
freqs = td.constants.C_0 / wavelengths
freq0 = freqs[len(freqs)//2]
freqw = freqs[0] - freqs[-1]

# Define material properties
n_BaF2 = 1.45
n_Ge = 4
BaF2 = td.Medium(epsilon=n_BaF2**2)
Ge = td.Medium(epsilon=n_Ge**2)
[3]:
# space between resonators and source
spc = 8

# geometric parameters
Px = Py = P = 4.2
h = 2.53
L1 = 3.036
L2 = 2.024
w1 = w2 = w = 1.265

# resolution (should be commensurate with periodicity)
dl = P / 32
[4]:
# total size in z and [x,y,z]
Lz = spc + h + h + spc
sim_size = [Px, Py, Lz]

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

# Define structure

cell1 = td.Box(
        center=[(L1/2)-L2, -w1/2, h/2],
        size=[L1, w1, h],
        material=Ge,
        name='cell1'
    )

cell2 = td.Box(
        center=[-L2/2, w2/2, h/2],
        size=[L2, w2, h],
        material=Ge,
        name='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 = 3e-11
[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 two simulations to run

  • With no resonator (normalization)

  • With Ge resonator

[7]:
# normalizing run (no Ge) 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, dl],
                    structures=[substrate],
                    sources=[source],
                    monitors=[monitor],
                    run_time=run_time/10,
                    pml_layers=[0,0,15])

# run with Ge nanorod
sim_actual = td.Simulation(size=sim_size,
                    mesh_step=[dl, dl, dl],
                    structures=[substrate,cell1,cell2],
                    sources=[source],
                    monitors=[monitor],
                    run_time=run_time,
                    pml_layers=[0,0,15])
Initializing simulation...
Mesh step (micron): [1.31e-01, 1.31e-01, 1.31e-01].
Simulation domain in number of grid points: [32, 32, 190].
Total number of computational grid points: 1.95e+05.
Total number of time steps: 13188.
Estimated data size (GB) of monitor transmitted_fields: 0.0000.
Initializing simulation...
Mesh step (micron): [1.31e-01, 1.31e-01, 1.31e-01].
Simulation domain in number of grid points: [32, 32, 190].
Total number of computational grid points: 1.95e+05.
Total number of time steps: 131875.
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))
td.viz.viz_eps_2D(sim_actual, normal='x', position=0, ax=ax1)
td.viz.viz_eps_2D(sim_actual, normal='y', position=-0.1, ax=ax2)
td.viz.viz_eps_2D(sim_actual, normal='z', position=0.1, ax=ax3)
plt.show()
../_images/examples_HighQ_Ge_10_0.png

Run Simulations

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

def get_flux(sim, task_name='Ge_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='Ge_normalization')
flux_actual = get_flux(sim_actual, task_name='Ge_resonator')
Uploading the json file...
Project 'Ge_normalization-16b6e20c5ec6fc90' status: success...

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

Compute source modes time (s):     0.0518
Compute monitor modes time (s):    0.0606

Rest of setup time (s):            0.0942

Starting solver...
- Time step    280 / time 6.37e-14s (  2 % done), field decay: 1.00e+00
- Time step    527 / time 1.20e-13s (  4 % done), field decay: 6.73e-01
- Time step   1055 / time 2.40e-13s (  8 % done), field decay: 4.19e-12
Field decay smaller than shutoff factor, exiting solver.

Solver time (s):                   1.5500
Post-processing time (s):          0.0045

Uploading the json file...
Project 'Ge_resonator-16b6e219232b5480' status: success...

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

Compute source modes time (s):     0.0560
Compute monitor modes time (s):    0.0797

Rest of setup time (s):            0.2238

Starting solver...
- Time step    280 / time 6.37e-14s (  0 % done), field decay: 1.00e+00
- Time step   5274 / time 1.20e-12s (  4 % done), field decay: 2.42e-02
- Time step  10549 / time 2.40e-12s (  8 % done), field decay: 7.43e-03
- Time step  15824 / time 3.60e-12s ( 12 % done), field decay: 1.34e-03
- Time step  21099 / time 4.80e-12s ( 16 % done), field decay: 1.58e-03
- Time step  26374 / time 6.00e-12s ( 20 % done), field decay: 1.26e-03
- Time step  31649 / time 7.20e-12s ( 24 % done), field decay: 1.21e-03
- Time step  36924 / time 8.40e-12s ( 28 % done), field decay: 9.62e-04
- Time step  42199 / time 9.60e-12s ( 32 % done), field decay: 7.86e-04
- Time step  47474 / time 1.08e-11s ( 36 % done), field decay: 5.50e-04
- Time step  52749 / time 1.20e-11s ( 40 % done), field decay: 4.15e-04
- Time step  58024 / time 1.32e-11s ( 44 % done), field decay: 2.84e-04
- Time step  63299 / time 1.44e-11s ( 48 % done), field decay: 1.52e-04
- Time step  68574 / time 1.56e-11s ( 52 % done), field decay: 1.32e-04
- Time step  73849 / time 1.68e-11s ( 56 % done), field decay: 1.68e-04
- Time step  79124 / time 1.80e-11s ( 60 % done), field decay: 1.44e-04
- Time step  84399 / time 1.92e-11s ( 64 % done), field decay: 9.44e-05
- Time step  89674 / time 2.04e-11s ( 68 % done), field decay: 8.29e-05
- Time step  94949 / time 2.16e-11s ( 72 % done), field decay: 1.15e-04
- Time step 100224 / time 2.28e-11s ( 76 % done), field decay: 1.07e-04
- Time step 105499 / time 2.40e-11s ( 80 % done), field decay: 6.79e-05
- Time step 110774 / time 2.52e-11s ( 84 % done), field decay: 2.04e-05
- Time step 116049 / time 2.64e-11s ( 88 % done), field decay: 3.08e-05
- Time step 121324 / time 2.76e-11s ( 92 % done), field decay: 7.56e-05
- Time step 126599 / time 2.88e-11s ( 96 % done), field decay: 8.18e-05
- Time step 131874 / time 3.00e-11s (100 % done), field decay: 4.11e-05

Solver time (s):                   33.3560
Post-processing time (s):          0.0051

[11]:
# plot normalizing run
plt.plot(wavelengths, np.abs(flux_empty), color='k')
plt.ylim([0,1])
plt.title('normalizing run (no Ge)')
plt.xlabel('wavelength ($\mu m$)')
plt.ylabel('Transmission')
plt.show()
../_images/examples_HighQ_Ge_14_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, 1 - flux_actual / flux_empty, 'k', label='R')
plt.plot(wavelengths, flux_actual / flux_empty, 'r--', label='T')
plt.xlabel('wavelength ($\mu m$)')
plt.ylabel('Magnitude')
plt.xlim([8.8, 12])
plt.ylim([0.0, 1.0])
plt.legend()
plt.show()
../_images/examples_HighQ_Ge_16_0.png
[ ]: