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.

diagram

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

Initial setup#

Our ring resonator will include a ring centered at (0,0) with a waveguide just above the ring spanning the x direction.

schematic

[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.0 + couple_width + wg_width / 2.0

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.mode 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()

[20:00:22] WARNING: No sources in simulation.                                                             log.py:50

Tidy3D is warning us that the simulation does not contain a source. However, since this simulation is used to construct the mode solver and will not be run directly, we can ignore this warning.

[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()

../_images/notebooks_RingResonator_13_0.png

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()

../_images/notebooks_RingResonator_21_0.png

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", verbose=True
)

 monitor_data.hdf5 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100.0%69.5/69.5 MB3.0 MB/s0:00:00
[20:03:52] loading SimulationData from data/simulation_data.hdf5                                      webapi.py:512
[13]:
# 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()

../_images/notebooks_RingResonator_24_0.png

Analyze Spectrum#

Now let’s analyze the mode amplitudes in the output waveguide.

First, let’s grab the data to inspect it.

[14]:
sim_data["mode"].amps

[14]:
<xarray.ModeAmpsDataArray (direction: 2, f: 1001, mode_index: 2)>
array([[[-1.27942002e-07+1.93495888e-07j,
         -9.90677167e-01-1.80035454e-02j],
        [ 6.44060089e-08+1.87396504e-07j,
         -9.82550124e-01-1.28286490e-01j],
        [ 2.24544422e-07+2.14037291e-07j,
         -9.61820596e-01-2.38181328e-01j],
        ...,
        [-1.16989442e-07+6.87071635e-08j,
          3.09820209e-01-9.50716471e-01j],
        [-3.34940083e-07+2.04879566e-07j,
          4.57415412e-01-8.89149400e-01j],
        [-4.30859483e-07+1.49696720e-07j,
          5.95632682e-01-8.03027192e-01j]],

       [[ 1.81772549e-09+6.25308137e-10j,
         -1.15496427e-05-9.02740434e-05j],
        [ 9.80917427e-09-2.88172801e-10j,
         -1.69948681e-05-6.49245322e-05j],
        [ 1.65801478e-08+2.97276175e-09j,
         -2.62162557e-05-6.96966250e-05j],
        ...,
        [ 9.69312680e-10-4.50085732e-09j,
          2.70893699e-04+6.52740963e-05j],
        [ 1.88242534e-10-1.19755894e-09j,
          2.99134015e-04+1.36260126e-04j],
        [-1.73317010e-10+1.04108480e-09j,
          2.22787967e-04+1.82471699e-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 amplitudes

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

[15]:
transmission_amps = sim_data["mode"].amps.sel(mode_index=1, direction="+")

Now let’s plot the data.

[16]:
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()

../_images/notebooks_RingResonator_30_0.png
[ ]: