Waveguide to ring coupling#
Optical ring resonators are key components in the field of integrated photonics. The unique capability of ring resonators to selectively interact with specific wavelengths of light makes them extremely versatile. They can be utilized in a broad array of applications, including optical filtering, modulating, switching, and sensing. However, simulating a ring resonator can be computationally expensive due to the high-Q resonances. Alternatively, we can investigate the coupling between a straight waveguide to a ring only by simulating the coupling region. The coupling coefficients can be extracted from this much simpler simulation.
Noticeably, in a waveguide-to-ring simulation, one of the simulation domain boundaries intersects the ring structure. Non-translational invariant structures inside PML are known to cause instabilities in a FDTD simulation so such a simulation is likely to diverge. In this notebook, we will try to use the PML boundary first to see if the simulation diverges. When it does, we will apply an effective remedy under this situation: replacing the PML with the adiabatic absorber boundary. The absorber functions similarly to PML such that it absorbs the outgoing radiation to mimic the infinite space. However, the absorber has a slightly higher reflection and requires a bit more computation than PML but it is numerically much more stable.
FDTD simulations can diverge due to various reasons. If you run into any simulation divergence issues, please follow the steps outlined in our troubleshooting guide to resolve it.
For more integrated photonic examples such as the 8-Channel mode and polarization de-multiplexer, the broadband bi-level taper polarization rotator-splitter, and the broadband directional coupler, please visit our examples page.
If you are new to the finite-difference time-domain (FDTD) method, we highly recommend going through our FDTD101 tutorials.
[1]:
import numpy as np
import matplotlib.pyplot as plt
import gdstk
import tidy3d as td
import tidy3d.web as web
Simulation Setup#
Define simulation wavelength range to be 1.5 \(\mu m\) to 1.6 \(\mu m\). .
[2]:
lda0 = 1.55 # central wavelength
ldas = np.linspace(1.5, 1.6, 101) # wavelength range of interest
freq0 = td.C_0 / lda0 # central frequency
freqs = td.C_0 / ldas # frequency range of interest
fwidth = 0.5 * (np.max(freqs) - np.min(freqs)) # frequency width of the source
The ring and waveguide are made of silicon. The top and bottom claddings are made of silicon oxide. Here we use the materials from Tidy3D’s material library directly.
[3]:
# define silicon and silicon dioxide media from the material library
si = td.material_library["cSi"]["Li1993_293K"]
sio2 = td.material_library["SiO2"]["Horiba"]
Define the geometric parameters. The waveguide is 500 nm wide and 220 nm thick. The ring has a radius of 5 \(\mu m\) and the gap size is 50 nm.
[4]:
w = 0.5 # width of the waveguide
h_si = 0.22 # thickness of the silicon layer
gap = 0.05 # gap size between the waveguides and the ring
r = 5 # radius of the ring
inf_eff = 1e2 # effective infinity
# simulation domain size
Lx = 2 * r + 2 * lda0
Ly = r / 2 + gap + 2 * w + lda0
Lz = 9 * h_si
We only need to define two structures: a straight waveguide and a ring. Both are commonly used PIC components introduced in the demonstration notebook. We can simply copy the associated functions here and them use to define the structures quickly. Namely, we copy the straight_waveguide function and the ring_resonator function.
[5]:
def straight_waveguide(x0, y0, z0, x1, y1, wg_width, wg_thickness, medium, sidewall_angle=0):
"""
This function defines a straight strip waveguide and returns the tidy3d structure of it.
Parameters
----------
x0: x coordinate of the waveguide starting position (um)
y0: y coordinate of the waveguide starting position (um)
z0: z coordinate of the waveguide starting position (um)
x1: x coordinate of the waveguide end position (um)
y1: y coordinate of the waveguide end position (um)
wg_width: width of the waveguide (um)
wg_thickness: thickness of the waveguide (um)
medium: medium of the waveguide
sidewall_angle: side wall angle of the waveguide (rad)
"""
cell = gdstk.Cell("waveguide") # define a gds cell
path = gdstk.RobustPath((x0, y0), wg_width, layer=1, datatype=0) # define a path
path.segment((x1, y1))
cell.add(path) # add path to the cell
# define geometry from the gds cell
wg_geo = td.PolySlab.from_gds(
cell,
gds_layer=1,
axis=2,
slab_bounds=(z0 - wg_thickness / 2, z0 + wg_thickness / 2),
sidewall_angle=sidewall_angle,
)
# define tidy3d structure of the bend
wg = td.Structure(geometry=wg_geo[0], medium=medium)
return wg
[6]:
def ring_resonator(
x0,
y0,
z0,
R,
wg_width,
wg_thickness,
medium,
sidewall_angle=0,
):
"""
This function defines a ring and returns the tidy3d structure of it.
Parameters
----------
x0: x coordinate of center of the ring (um)
y0: y coordinate of center of the ring (um)
z0: z coordinate of center of the ring (um)
R: radius of the ring (um)
wg_width: width of the waveguide (um)
wg_thickness: thickness of the waveguide (um)
medium: medium of the waveguide
sidewall_angle: side wall angle of the waveguide (rad)
"""
cell = gdstk.Cell("top") # define a gds cell
# define a path
path_top = gdstk.RobustPath(
(x0 + R, y0), wg_width - wg_thickness * np.tan(np.abs(sidewall_angle)), layer=1, datatype=0
)
path_top.arc(R, 0, np.pi) # make the top half of the ring
cell.add(path_top) # add path to the cell
# the reference plane depends on the sign of the sidewall_angle
if sidewall_angle >= 0:
reference_plane = "top"
else:
reference_plane = "bottom"
# define top half ring geometry
ring_top_geo = td.PolySlab.from_gds(
cell,
gds_layer=1,
axis=2,
slab_bounds=(z0 - wg_thickness / 2, z0 + wg_thickness / 2),
sidewall_angle=sidewall_angle,
reference_plane=reference_plane,
)
# similarly for the bottom half of the ring
cell = gdstk.Cell("bottom")
path_bottom = gdstk.RobustPath(
(x0 + R, y0), wg_width - wg_thickness * np.tan(np.abs(sidewall_angle)), layer=1, datatype=0
)
path_bottom.arc(R, 0, -np.pi)
cell.add(path_bottom)
ring_bottom_geo = td.PolySlab.from_gds(
cell,
gds_layer=1,
axis=2,
slab_bounds=(z0 - wg_thickness / 2, z0 + wg_thickness / 2),
sidewall_angle=sidewall_angle,
reference_plane=reference_plane,
)
# define ring structure
ring = td.Structure(
geometry=td.GeometryGroup(geometries=ring_bottom_geo + ring_top_geo), medium=medium
)
return ring
Use the above functions to define the Structures.
[7]:
# define straight waveguide
waveguide = straight_waveguide(
x0=-inf_eff,
y0=0,
z0=0,
x1=inf_eff,
y1=0,
wg_width=w,
wg_thickness=h_si,
medium=si,
sidewall_angle=0,
)
# define ring
ring = ring_resonator(
x0=0,
y0=w + gap + r,
z0=0,
R=r,
wg_width=w,
wg_thickness=h_si,
medium=si,
sidewall_angle=0,
)
We will use a ModeSource to excite the straight waveguide using the fundamental TE mode. A ModeMonitor is placed at the through to measure the transmission. Another ModeMonitor is placed at the
ring to measure the coupling to the ring. For the monitor at the ring, we need to properly set up angle_theta, angle_theta, and bend_axis in the ModeSpec as demonstrated in the tutorial. Finally, we add a
FieldMonitor to help visualize the field distribution.
[8]:
n_si = 3.47
# mode spec for the source
mode_spec_source = td.ModeSpec(num_modes=1, target_neff=n_si)
# mode spec for the through port
mode_spec_through = mode_spec_source
# angle of the mode at the ring
theta = np.pi / 4
# mode spec for the drop port at the ring
mode_spec_drop = td.ModeSpec(
num_modes=1, target_neff=n_si, angle_theta=theta, bend_radius=-r, bend_axis=1
)
# add a mode source as excitation
mode_source = td.ModeSource(
center=(-r - lda0 / 4, 0, 0),
size=(0, 6 * w, 6 * h_si),
source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
direction="+",
mode_spec=mode_spec_source,
mode_index=0,
)
# add a mode monitor to measure transmission at the through port
mode_monitor_through = td.ModeMonitor(
center=(r + lda0 / 4, 0, 0),
size=mode_source.size,
freqs=freqs,
mode_spec=mode_spec_through,
name="through",
)
# add a mode monitor to measure transmission at the drop port
mode_monitor_drop = td.ModeMonitor(
center=(np.sin(theta) * r, w + gap + r - np.cos(theta) * r, 0),
size=(6 * w, 0, 6 * h_si),
freqs=freqs,
mode_spec=mode_spec_drop,
name="drop",
)
# add a field monitor to visualize the field distribution
field_monitor = td.FieldMonitor(
center=(0, 0, 0), size=(td.inf, td.inf, 0), freqs=[freq0], name="field"
)
Using PML Boundary#
Finally, we define the Simulation. Intuitively, we would use the PML boundary on all sides similar to what we did in most other PIC component simulations. Let’s try it here.
[9]:
run_time = 2e-12 # simulation run time
# construct simulation
sim_pml = td.Simulation(
center=(0, Ly / 4, 0),
size=(Lx, Ly, Lz),
grid_spec=td.GridSpec.auto(min_steps_per_wvl=25, wavelength=lda0),
structures=[waveguide, ring],
sources=[mode_source],
monitors=[mode_monitor_through, mode_monitor_drop, field_monitor],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
medium=sio2,
symmetry=(0, 0, 1),
)
# plot simulation
sim_pml.plot(z=0)
plt.show()
Submit the simulation to the server to run.
[10]:
sim_data = web.run(
simulation=sim_pml, task_name="waveguide_to_ring", path="data/simulation_data.hdf5"
)
16:05:18 PST Created task 'waveguide_to_ring' with task_id 'fdve-f6a84898-dd7d-4874-98c4-87b4d3705495' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-f6a84898-dd7 d-4874-98c4-87b4d3705495'.
16:05:20 PST status = queued
16:05:25 PST status = preprocess
16:05:27 PST Maximum FlexCredit cost: 0.469. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
starting up solver
16:05:28 PST running solver
To cancel the simulation, use 'web.abort(task_id)' or 'web.delete(task_id)' or abort/delete the task in the web UI. Terminating the Python script will not stop the job running on the cloud.
16:07:11 PST early shutoff detected at 52%, exiting.
status = postprocess
16:07:19 PST status = diverged
View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-f6a84898-dd7 d-4874-98c4-87b4d3705495'.
16:07:20 PST loading simulation from data/simulation_data.hdf5
WARNING: The simulation has diverged! For more information, check 'SimulationData.log' or use 'web.download_log(task_id)'.
Now we see that the status of the job is shown to be diverged. As discussed in the introduction, it should not come as a surprise since the ring section in the PML layer is not perpendicular to the boundary and translational invariant as shown in the zoomed-in plot below. We see that we are also given a warning about the divergence and prompted to check the log for more information.
[12]:
print(sim_data.log)
[00:05:25] USER: Simulation domain Nx, Ny, Nz: [759, 280, 82]
USER: Applied symmetries: (0, 0, 1)
USER: Number of computational grid points: 9.2279e+06.
USER: Using subpixel averaging: True
USER: Number of time steps: 6.2950e+04
USER: Automatic shutoff factor: 1.00e-05
USER: Time step (s): 3.1772e-17
USER:
USER: Compute source modes time (s): 0.5671
[00:05:26] USER: Rest of setup time (s): 0.2845
[00:06:15] USER: Compute monitor modes time (s): 46.0400
[00:07:07] USER: Solver time (s): 98.8520
USER: Time-stepping speed (cells/s): 5.98e+09
WARNING: The simulation has diverged!
WARNING: Some structures were found to be spatially varying inside
the PML. Ensure that structures are translationally invariant in the
PML regions in the direction normal to the PML interface.
Alternatively, switching the PML to Absorber boundary should fix the
divergence.
WARNING: A dispersive medium inside the PML regions was detected,
which may be causing the divergence. Consider making the medium
non-dispersive, or fitting a different dispersive model for that
medium. Alternatively, switching the PML to Absorber boundary should
fix the divergence.
USER: Post-processing time (s): 0.5129
====== SOLVER LOG ======
Processing grid and structures...
Building FDTD update coefficients...
Potential divergence: dispersive medium into PML.
Potential divergence: structures are not translationally invariant inside PML.
Solver setup time (s): 1.5782
Running solver for 62950 time steps...
- Time step 2517 / time 8.00e-14s ( 4 % done), field decay: 1.00e+00
- Time step 4010 / time 1.27e-13s ( 6 % done), field decay: 1.00e+00
- Time step 5035 / time 1.60e-13s ( 8 % done), field decay: 1.00e+00
- Time step 7553 / time 2.40e-13s ( 12 % done), field decay: 1.00e+00
- Time step 10071 / time 3.20e-13s ( 16 % done), field decay: 6.82e-02
- Time step 12589 / time 4.00e-13s ( 20 % done), field decay: 2.32e-04
- Time step 15107 / time 4.80e-13s ( 24 % done), field decay: 1.00e+00
- Time step 17625 / time 5.60e-13s ( 28 % done), field decay: 1.00e+00
- Time step 20143 / time 6.40e-13s ( 32 % done), field decay: 1.00e+00
- Time step 22661 / time 7.20e-13s ( 36 % done), field decay: 1.00e+00
- Time step 25179 / time 8.00e-13s ( 40 % done), field decay: 1.00e+00
- Time step 27697 / time 8.80e-13s ( 44 % done), field decay: 1.00e+00
- Time step 30215 / time 9.60e-13s ( 48 % done), field decay: 1.00e+00
- Time step 32733 / time 1.04e-12s ( 52 % done), field decay: 1.00e+00
Field diverged at time step 33860 ( 53 % done), exiting solver.
Time-stepping time (s): 97.0829
Data write time (s): 0.1915
The log indeed tells us that there are two potential causes of divergence that were found in this simulation. One is the presence of dispersive materials in the PML, while the other is the fact that the ring is not translationally invariant. While both of those things can cause divergence, our in-built silicon material model has been fit to perform well for straight waveguides going into the PML, at telecommunication wavelengths. Thus, the more likely cause of the divergence is the fact that the ring is not translationally invariant.
[13]:
# zoom-in plot around the pml
ax = sim_pml.plot(z=0)
ax.set_xlim(4, 5.5)
ax.set_ylim(3.5, 4.5)
plt.show()
This can be further confirmed by plotting the electric field distribution. From the plot, we can see where the field built up and caused the simulation to diverge. As expected, the field built up around the ring to PML area.
[12]:
sim_data.plot_field(field_monitor_name="field", field_name="E", val="abs")
plt.show()
Using Absorber Boundary#
To resolve the diverge, we can switch to Absorber. This can be done by copying the previous simulation and only updating the boundary_spec. Since the ring is only intersecting the boundary in the y direction, we only need to apply Absorber to the positive \(y\) boundary.
In addition, we increased the number of layers in the Absorber to 60 to ensure sufficient absorption and minimal reflection.
[13]:
# copy simulation and update boundary condition
sim_absorber = sim_pml.copy(
update={
"boundary_spec": td.BoundarySpec(
x=td.Boundary.pml(),
y=td.Boundary(plus=td.Absorber(num_layers=60), minus=td.PML()),
z=td.Boundary.pml(),
)
}
)
# run simulation
sim_data = web.run(
simulation=sim_absorber, task_name="waveguide_to_ring", path="data/simulation_data.hdf5"
)
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- f78250d4-0834-42a9-92e0-d47037b0bab1v1'.
After switching to Absorber, the simulation doesn’t run into the divergence issue anymore and thus we can plot the transmission spectra to the through and drop ports.
[14]:
# extract mode amplitude in the through port
t = sim_data["through"].amps.sel(mode_index=0, direction="+")
# extract mode amplitude in the drop port
k = sim_data["drop"].amps.sel(mode_index=0, direction="+")
# plot transmission
plt.plot(ldas, np.abs(t) ** 2, label="Through port")
plt.plot(ldas, np.abs(k) ** 2, label="Drop port")
plt.legend()
plt.xlabel("Wavelength ($\mu m$)")
plt.ylabel("Transmission")
plt.xlim(1.5, 1.6)
plt.ylim(0, 1)
plt.show()
Again we can visualize the field distribution. This time, we see the correct result where the energy is partially coupled to the ring and partially stays in the straight waveguide.
[15]:
sim_data.plot_field(field_monitor_name="field", field_name="E", val="abs")
plt.show()
Additional Notes on Absorber#
The adiabatic absorber is a multilayer system with gradually increasing conductivity. As briefly discussed above, the absorber usually has a larger undesired reflection compared to PML. In practice, this small difference rarely matters, but is important to understand for simulations that require high accuracy. There are two possible sources for the reflection from absorbers. The first, and more common one, is that the ramping up of the conductivity is not sufficiently slow, which can be remedied
by increasing the number of absorber layers (40 by default). The second one is that the absorption is not high enough, such that the light reaches the PEC boundary at the end of the Absorber, travels back through it, and is still not fully attenuated before re-entering the simulation region. If this is the case, increasing the maximum conductivity (see the API reference) can help. In both cases,
changing the order of the scaling of the conductivity (sigma_order) can also have an effect, but this is a more advanced setting that we typically do not recommend modifying.
[ ]: