2D ring resonator
Contents
2D ring resonator#
View this project in Tidy3D Web App.
This is a simple example of using Tidy3D to perform a 2D simulation of a ring resonator side coupled to a dielectric waveguide.
With a center wavelength of 500 nm and 10 nm resolution, this is a challenging FDTD problem because of the large simulation size and the large number of time steps needed to capture the resonance of the ring. However, it can be easily handled with Tidy3D.
[1]:
# standard python imports
import numpy as np
from numpy import random
import matplotlib.pyplot as plt
# tidy3D import
import tidy3d.web as web
import tidy3d as td
[23:36:34] WARNING This version of Tidy3D was pip installed from the 'tidy3d-beta' repository on __init__.py:103 PyPI. Future releases will be uploaded to the 'tidy3d' repository. From now on, please use 'pip install tidy3d' instead.
INFO Using client version: 1.9.0rc1 __init__.py:121
Initial setup#
Our ring resonator will include a ring centered at (0,0) with a waveguide just above the ring spanning the x direction.
[2]:
# define geometry
wg_width = 0.25
couple_width = 0.05
ring_radius = 3.5
ring_wg_width = 0.25
wg_spacing = 2.0
buffer = 2.0
# compute quantities based on geometry parameters
x_span = 2*wg_spacing + 2*ring_radius + 2*buffer
y_span = 2*ring_radius + 2*ring_wg_width + wg_width + couple_width + 2*buffer
wg_insert_x = ring_radius + wg_spacing
wg_center_y = ring_radius + ring_wg_width/2. + couple_width + wg_width/2.
Define frequency parameters for sources and monitors based on the wavelength range of interest. Since in this case study we are investigating the system’s response over a wide range of frequencies we will use a broadband mode source that takes into account frequency-dependent character of waveguide modes. The Tidy3D’s broadband feature is designed in such a way that it produces the most accurate results in the range (freq0 - 1.5 * fwidth, freq0 + 1.5 * fwidth), thus we define parameters
freq0 and fwidth such that the frequency range of interest is exactly covered by (freq0 - 1.5 * fwidth, freq0 + 1.5 * fwidth).
[3]:
# wavelength range of interest
lambda_beg = 0.4
lambda_end = 0.6
# define pulse parameters
freq_beg = td.C_0 / lambda_end
freq_end = td.C_0 / lambda_beg
freq0 = (freq_beg + freq_end) / 2
fwidth = (freq_end - freq0) / 1.5
min_steps_per_wvl = 30
run_time = 1e-11
Define materials. (docs)
[4]:
n_bg = 1.0
n_solid = 1.5
background = td.Medium(permittivity=n_bg**2)
solid = td.Medium(permittivity=n_solid**2)
Define structures. (docs)
[5]:
# background of entire domain (set explicitly as a box)
background_box = td.Structure(
geometry=td.Box(
center=[0, 0, 0],
size=[td.inf, td.inf, td.inf],
),
medium=background,
name='background')
# waveguide
waveguide = td.Structure(
geometry=td.Box(
center=[0, wg_center_y, 0],
size=[td.inf, wg_width, td.inf],
),
medium=solid,
name='waveguide')
# outside ring
outer_ring = td.Structure(
geometry=td.Cylinder(
center=[0,0,0],
axis=2,
radius=ring_radius+ring_wg_width/2.0,
length=td.inf,
),
medium=solid,
name='outer_ring')
# inside ring fill
inner_ring = td.Structure(
geometry=td.Cylinder(
center=[0,0,0],
axis=2,
radius=ring_radius-ring_wg_width/2.0,
length=td.inf,
),
medium=background,
name='inner_ring')
Compute and visualize the waveguide modes.
[6]:
from tidy3d.plugins import ModeSolver
mode_plane = td.Box(
center=[-wg_insert_x, wg_center_y, 0],
size=[0, 2, td.inf],
)
sim_modesolver = td.Simulation(
size=[x_span, y_span, 0],
grid_spec=td.GridSpec.auto(min_steps_per_wvl=min_steps_per_wvl, wavelength=td.C_0 / freq0),
structures=[background_box, waveguide],
run_time=1e-12,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.Periodic())
)
mode_spec = td.ModeSpec(num_modes=2)
mode_solver = ModeSolver(simulation=sim_modesolver, plane=mode_plane, mode_spec=mode_spec, freqs=[freq0])
mode_data = mode_solver.solve()
WARNING No sources in simulation. simulation.py:494
/Users/twhughes/.pyenv/versions/3.10.9/lib/python3.10/site-packages/numpy/linalg/linalg.py:2139: RuntimeWarning: divide by zero encountered in det
r = _umath_linalg.det(a, signature=signature)
/Users/twhughes/.pyenv/versions/3.10.9/lib/python3.10/site-packages/numpy/linalg/linalg.py:2139: RuntimeWarning: invalid value encountered in det
r = _umath_linalg.det(a, signature=signature)
[7]:
f, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2, 3, tight_layout=True, figsize=(10, 6))
mode_data.Ex.sel(mode_index=0, f=freq0).abs.plot(ax=ax1)
mode_data.Ey.sel(mode_index=0, f=freq0).abs.plot(ax=ax2)
mode_data.Ez.sel(mode_index=0, f=freq0).abs.plot(ax=ax3)
mode_data.Ex.sel(mode_index=1, f=freq0).abs.plot(ax=ax4)
mode_data.Ey.sel(mode_index=1, f=freq0).abs.plot(ax=ax5)
mode_data.Ez.sel(mode_index=1, f=freq0).abs.plot(ax=ax6)
ax1.set_title('|Ex|: mode_index=0')
ax2.set_title('|Ey|: mode_index=0')
ax3.set_title('|Ez|: mode_index=0')
ax4.set_title('|Ex|: mode_index=1')
ax5.set_title('|Ey|: mode_index=1')
ax6.set_title('|Ez|: mode_index=1')
plt.show()
From the above plots, we see that
mode_index=0 corresponds to exciting 0-th order TM mode (E=Ez) and
mode_index=1 corresponds to exciting 0-th order TE mode (E=Ey).
We can therefore switch the mode index accordingly based on our polarization.
Let’s select Ey and create the source for it. As mentioned before, we will use a broadband mode source that will approximate the chosen mode’s frequency dependence in the range (freq0 - 1.5 * fwidth, freq0 + 1.5 * fwidth) using 11 Chebyshev polynomials.
[8]:
mode_source = td.ModeSource(
size=mode_plane.size,
center=mode_plane.center,
source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
mode_spec=td.ModeSpec(num_modes=2),
mode_index=1,
direction='+',
num_freqs=11, # using 11 (Chebyshev) points to approximate frequency dependence
)
In addition, let’s monitor both the fields in plane as well as the output mode amplitudes into the fundamental TE mode.
[9]:
# monitor steady state fields at central frequency over whole domain
field_monitor = td.FieldMonitor(
center=[0, 0, 0],
size=[td.inf, td.inf, 0],
freqs=[freq0],
name='field')
# monitor the mode amps on the output waveguide
lambdas_measure = np.linspace(lambda_beg, lambda_end, 1001)
freqs_measure = td.C_0 / lambdas_measure[::-1]
mode_monitor = td.ModeMonitor(
size=mode_plane.size,
center=mode_plane.center,
freqs=freqs_measure,
mode_spec=td.ModeSpec(num_modes=2),
name='mode'
)
# lets reset the center to the on the right hand side of the simulation though
mode_monitor = mode_monitor.copy(update=dict(center = [+wg_insert_x, wg_center_y, 0]))
Define simulation. (docs) To make the simulation 2D, we can just set the simulation size in one of the dimensions to be 0. However, note that we still have to define a grid size in that direction. Additionally, we lower the shutoff factor to make sure the high-frequency response of this resonant system is accurately resolved.
[10]:
# create simulation
sim = td.Simulation(
size=[x_span, y_span, 0],
grid_spec=td.GridSpec.auto(min_steps_per_wvl=min_steps_per_wvl, wavelength=td.C_0 / freq0),
structures=[background_box, waveguide, outer_ring, inner_ring],
sources=[mode_source],
monitors=[field_monitor, mode_monitor],
run_time = run_time,
boundary_spec=td.BoundarySpec(x=td.Boundary.pml(), y=td.Boundary.pml(), z=td.Boundary.periodic()),
shutoff=1e-9,
)
Visualize structure, source, and modes. (docs)
[11]:
# plot the two simulations
fig, ax = plt.subplots(1, 1, figsize=(6, 6))
sim.plot_eps(z=0.01, ax=ax)
plt.show()
Run Simulation#
Run simulations on our server. (docs)
[12]:
# use function above to run simulation
sim_data = web.run(sim, task_name='ring_resonator', path='data/simulation_data.hdf5')
[13]:
print(sim_data.log)
Simulation domain Nx, Ny, Nz: [1433, 1001, 1]
Applied symmetries: (0, 0, 0)
Number of computational grid points: 1.4393e+06.
Using subpixel averaging: True
Number of time steps: 4.6047e+05
Automatic shutoff factor: 1.00e-09
Time step (s): 2.1717e-17
Compute source modes time (s): 1.4599
Compute monitor modes time (s): 101.3099
Rest of setup time (s): 10.5289
Running solver for 460471 time steps...
- Time step 440 / time 9.56e-15s ( 0 % done), field decay: 1.00e+00
- Time step 18418 / time 4.00e-13s ( 4 % done), field decay: 2.03e-01
- Time step 36837 / time 8.00e-13s ( 8 % done), field decay: 1.89e-02
- Time step 55256 / time 1.20e-12s ( 12 % done), field decay: 4.68e-03
- Time step 73675 / time 1.60e-12s ( 16 % done), field decay: 1.01e-03
- Time step 92094 / time 2.00e-12s ( 20 % done), field decay: 3.39e-04
- Time step 110513 / time 2.40e-12s ( 24 % done), field decay: 1.35e-04
- Time step 128931 / time 2.80e-12s ( 28 % done), field decay: 4.18e-05
- Time step 147350 / time 3.20e-12s ( 32 % done), field decay: 1.88e-05
- Time step 165769 / time 3.60e-12s ( 36 % done), field decay: 6.14e-06
- Time step 184188 / time 4.00e-12s ( 40 % done), field decay: 3.04e-06
- Time step 202607 / time 4.40e-12s ( 44 % done), field decay: 1.54e-06
- Time step 221026 / time 4.80e-12s ( 48 % done), field decay: 7.29e-07
- Time step 239444 / time 5.20e-12s ( 52 % done), field decay: 4.86e-07
- Time step 257863 / time 5.60e-12s ( 56 % done), field decay: 1.62e-07
- Time step 276282 / time 6.00e-12s ( 60 % done), field decay: 1.10e-07
- Time step 294701 / time 6.40e-12s ( 64 % done), field decay: 4.66e-08
- Time step 313120 / time 6.80e-12s ( 68 % done), field decay: 2.62e-08
- Time step 331539 / time 7.20e-12s ( 72 % done), field decay: 1.60e-08
- Time step 349957 / time 7.60e-12s ( 76 % done), field decay: 9.43e-09
- Time step 368376 / time 8.00e-12s ( 80 % done), field decay: 6.07e-09
- Time step 386795 / time 8.40e-12s ( 84 % done), field decay: 2.62e-09
- Time step 405214 / time 8.80e-12s ( 88 % done), field decay: 2.28e-09
- Time step 423633 / time 9.20e-12s ( 92 % done), field decay: 1.23e-09
- Time step 442052 / time 9.60e-12s ( 96 % done), field decay: 7.84e-10
Field decay smaller than shutoff factor, exiting solver.
Solver time (s): 117.2070
[14]:
# visualize normalization run
fig, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(15, 6))
ax1 = sim_data.plot_field('field', 'Ex', val='real', z=0, ax=ax1)
ax2 = sim_data.plot_field('field', 'Ey', val='real', z=0, ax=ax2)
ax1.set_title('Ex')
ax2.set_title('Ey')
plt.show()
Analyze Spectrum#
Now let’s analyze the mode amplitudes in the output waveguide.
First, let’s grab the data to inspect it.
[15]:
sim_data['mode'].amps
[15]:
<xarray.ModeAmpsDataArray (direction: 2, f: 1001, mode_index: 2)>
array([[[-1.27228285e-07+1.93114125e-07j,
-9.90342063e-01-3.09109412e-02j],
[ 6.40936866e-08+1.90399395e-07j,
-9.80814112e-01-1.40976772e-01j],
[ 2.23946689e-07+2.19753956e-07j,
-9.58727946e-01-2.50409607e-01j],
...,
[-1.16812559e-07+6.89954062e-08j,
3.12045632e-01-9.49995321e-01j],
[-3.35518463e-07+2.04260623e-07j,
4.59846595e-01-8.87904341e-01j],
[-4.31400983e-07+1.49033071e-07j,
5.98078371e-01-8.01217452e-01j]],
[[ 1.85592892e-09+7.80862030e-10j,
-1.09921095e-05-8.99535549e-05j],
[ 9.89061559e-09-4.09757933e-10j,
-1.62073811e-05-6.42821364e-05j],
[ 1.62798405e-08+3.22998360e-09j,
-2.37364261e-05-6.40256253e-05j],
...,
[ 9.62290546e-10-4.52698563e-09j,
2.69287517e-04+6.53467658e-05j],
[ 1.70808352e-10-1.14896296e-09j,
3.00599786e-04+1.34636139e-04j],
[-1.30603267e-10+1.03233031e-09j,
2.23401136e-04+1.85317607e-04j]]])
Coordinates:
* direction (direction) <U1 '+' '-'
* f (f) float64 4.997e+14 4.998e+14 5e+14 ... 7.491e+14 7.495e+14
* mode_index (mode_index) int64 0 1
Attributes:
units: sqrt(W)
long_name: mode amplitudesAs we see, the mode amplitude data is complex-valued with three 3 dimensions:
index into the mode order returned by solver (remember, we wanted mode_index=1 for fundamental TE).
direction of the propagation (for decomposition).
frequency.
Let’s select into the first two dimensions to get mode amplitudes as a function of frequency.
[16]:
transmission_amps = sim_data['mode'].amps.sel(mode_index=1, direction='+')
Now let’s plot the data.
[17]:
f, ax = plt.subplots(figsize=(10, 5))
transmission_amps.abs.plot.line(x='f', ax=ax, label='absolute value')
# transmission_amps.real.abs.plot.line(x='f', ax=ax, label='|real|')
ax.legend()
ax.set_title('transmission amplitudes into fundamental TE (forward)')
# ax.set_ylim(0, 1)
ax.set_xlim(freqs_measure[0], freqs_measure[-1])
plt.show()
[ ]: