90 degree optical hybrid#
Note: The cost of running the entire notebook is higher than 1 FlexCredit.
90 degree optical hybrids (also known as quadrature optical hybrids) are essential components in coherent transmission systems. A 90 degree optical hybrid is a six-port device consisting of two inputs and four outputs. In the ideal case, the outputs are the mixed signals of the inputs with 0, \(\pi\)/2, \(\pi\), and 3\(\pi\)/2 relative phase shifts.
This notebook demonstrates the simulation of a compact and low-loss 90 degree optical hybrid based on a silicon-on-insulator platform. The device consists of a Y-branch, three 2x2 MMIs, and four 90 degree waveguide bends. Building those structures natively in Tidy3D is certainly doable but time consuming. Here, we build the device structures in a separate CAD editor and import the stl file into Tidy3D for simulation.
To use this functionality, remember to install Tidy3D as pip install "tidy3d[trimesh]", which will install optional dependencies needed for processing surface meshes.
The device design is adapted from Hang Guan et al., “Compact and low loss 90° optical hybrid on a silicon-on-insulator platform,” Opt. Express 25, 28957-28968 (2017).
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.
[1]:
import numpy as np
import matplotlib.pyplot as plt
import tidy3d as td
import tidy3d.web as web
Simulation of the 2x2 MMI#
To design the optical hybrid, we first design each component individually. For the Y-branch, we will use the same low-loss design demonstrated in the reference and our Y-branch notebook. The waveguide bends can be Euler bends or circular bands, which is also demonstrated in another notebook.
For the 2x2 MMIs, an optimized design is demonstrated here. The design aims to have equal power splitting and a 90 degree phase difference in the through port and cross port. To simulate the MMI, we also built the structure as a stl file and only import it here.
First, define some basic simulation parameters. We simulate at a relatively narrow wavelength band of 1530 nm to 1560 nm.
[2]:
lda0 = 1.55 # central wavelength
freq0 = td.C_0 / lda0 # central frequency
ldas = np.linspace(1.53, 1.56, 101) # wavelength range
freqs = td.C_0 / ldas # frequency range
fwidth = 0.5 * (np.max(freqs) - np.min(freqs))
Define material properties. For simplicity, we only use non-dispersive materials in this simulation.
[3]:
n_si = 3.47 # silicon refractive index
si = td.Medium(permittivity=n_si**2)
n_sio2 = 1.45 # silicon oxide refractive index
sio2 = td.Medium(permittivity=n_sio2**2)
The structures will be imported from stl files, but it is good to define some geometric parameters to help set up the source, monitors, and simulation domains.
[4]:
thickness = 0.22 # si layer thickness
width = 0.5 # waveguide width
Import the MMI geometry from the stl file and define the MMI structure. All stl files used in this notebook can be downloaded from our documentation repo.
[5]:
# import mmi geometry from a stl file
mmi_geometry = td.TriangleMesh.from_stl(
filename="misc/mmi_stl.stl",
)
# define mmi structure
mmi = td.Structure(geometry=mmi_geometry, medium=si)
Define a ModeSource launching the TE0 mode at the top left waveguide. Two ModeMonitors are added to the waveguides on the right to measure the power and phase of the outputs. To visualize the field distribution within the MMI region, we also add a FieldMonitor in the xy plane. Finally, define a Simulation.
[6]:
mode_spec = td.ModeSpec(
num_modes=1, target_neff=n_si
) # define a ModeSpec used in mode source and monitors
# define a mode source
mode_source = td.ModeSource(
center=(-5, 0.45, 0),
size=(0, 2 * width, 6 * thickness),
source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
direction="+",
mode_spec=td.ModeSpec(num_modes=1, target_neff=n_si),
mode_index=0,
)
# define a mode monitor at the through port
mode_monitor_through = td.ModeMonitor(
center=(6, 0.45, 0),
size=(0, 2 * width, 6 * thickness),
freqs=freqs,
name="through",
mode_spec=mode_spec,
)
# define a mode monitor at the cross port
mode_monitor_cross = td.ModeMonitor(
center=(6, -0.45, 0),
size=(0, 2 * width, 6 * thickness),
freqs=freqs,
name="cross",
mode_spec=mode_spec,
)
# define a field monitor at the xp plane
field_monitor = td.FieldMonitor(
center=(0, 0, 0), size=(td.inf, td.inf, 0), freqs=[freq0], name="field"
)
run_time = 1e-12 # simulation run time
sim = td.Simulation(
size=(13, 5, 10 * thickness),
grid_spec=td.GridSpec.auto(min_steps_per_wvl=30, wavelength=lda0),
structures=[mmi],
sources=[mode_source],
monitors=[field_monitor, mode_monitor_through, mode_monitor_cross],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
medium=sio2,
symmetry=(0, 0, 1),
)
Visualize the simulation.
[7]:
ax = sim.plot(z=0)
Submit the simulation job to the server.
[8]:
job = web.Job(simulation=sim, task_name="mmi", verbose=True)
sim_data = job.run(path="data/simulation_data.hdf5")
19:53:45 CEST Created task 'mmi' with task_id 'fdve-a26f1de4-1c45-4278-885a-d72c19338a62' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-a26f1de4-1c 45-4278-885a-d72c19338a62'.
19:53:51 CEST status = queued
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.
19:53:54 CEST status = preprocess
19:53:56 CEST Maximum FlexCredit cost: 0.150. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
starting up solver
running solver
19:55:09 CEST early shutoff detected at 88%, exiting.
status = postprocess
19:55:14 CEST status = success
19:55:15 CEST View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-a26f1de4-1c 45-4278-885a-d72c19338a62'.
19:55:19 CEST loading simulation from data/simulation_data.hdf5
The simulation is complete and the monitor data have been downloaded. First, we plot the field intensity distribution to visualize the power flow. We can observe that the two outputs have a similar intensity.
[9]:
sim_data.plot_field(field_monitor_name="field", field_name="E", val="abs^2", vmax=2000)
plt.show()
To quantitatively investigate the power in the two output waveguides, we compute the power from the mode amplitudes. Compared to the simulation results reported in the reference, the small discrepancy here is likely due to the different material properties used.
Ideally, both power levels are 3 dB. Here, we see the through port has a slightly lower power level than the cross port. However, the small difference is acceptable.
[10]:
# compute power at the through port
P_through = np.abs(sim_data["through"].amps.sel(direction="+").values) ** 2
# compute power at the cross port
P_cross = np.abs(sim_data["cross"].amps.sel(direction="+").values) ** 2
# plot loss
plt.plot(ldas, -10 * np.log10(P_through), label="through port")
plt.plot(ldas, -10 * np.log10(P_cross), label="cross port")
plt.xlabel("Wavelength ($\mu m$)")
plt.ylabel("Insertion loss (dB)")
plt.ylim(2.5, 4)
plt.legend()
plt.show()
Lastly, we check the phase difference in the two outputs. Within the wavelength range, the phase difference is only about 1 degree from the ideal 90 degree.
[11]:
# compute phase at the through port
phase_through = np.angle(sim_data["through"].amps.sel(direction="+").values)
# compute phase at the cross port
phase_cross = np.angle(sim_data["cross"].amps.sel(direction="+").values)
# compute phase difference
delta_phase = np.mod(phase_through - phase_cross, 2 * np.pi) * 180 / np.pi
# plot phase difference
plt.plot(ldas, delta_phase)
plt.xlabel("Wavelength ($\mu m$)")
plt.ylabel("Phase difference")
plt.ylim(85, 95)
plt.show()
Simulation of the Optical Hybrid – LO Input#
Now that we have tested individual components, we can put them together and simulate the whole device. Again, we directly import the geometry from a stl file.
[12]:
# import optical hybrid geometry from a stl file
optical_hybrid_geometry = td.TriangleMesh.from_stl(
filename="misc/optical_hybrid_stl.stl",
)
# define optical hybrid structure
optical_hybrid = td.Structure(geometry=optical_hybrid_geometry, medium=si)
In the first case, define a ModeSource at the LO input and four ModeMonitors at the outputs.
[13]:
# define a mode source at the LO input waveguide on the top y branch
mode_source_LO = td.ModeSource(
center=(0, 1, 0),
size=(2 * width, 0, 4 * thickness),
source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
direction="-",
mode_spec=mode_spec,
mode_index=0,
)
# define mode monitors at the four output waveguides
mode_monitor_Qn = td.ModeMonitor(
center=(-16, -6, 0),
size=(0, 2 * width, 4 * thickness),
freqs=freqs,
name="port_Qn",
mode_spec=mode_spec,
)
mode_monitor_Qp = td.ModeMonitor(
center=(-16, -6.9, 0),
size=(0, 2 * width, 4 * thickness),
freqs=freqs,
name="port_Qp",
mode_spec=mode_spec,
)
mode_monitor_In = td.ModeMonitor(
center=(16, -6, 0),
size=(0, 2 * width, 4 * thickness),
freqs=freqs,
name="port_In",
mode_spec=mode_spec,
)
mode_monitor_Ip = td.ModeMonitor(
center=(16, -6.9, 0),
size=(0, 2 * width, 4 * thickness),
freqs=freqs,
name="port_Ip",
mode_spec=mode_spec,
)
Define the simulation.
[14]:
run_time = 2e-12 # simulation run time
# define simulation
sim = td.Simulation(
center=(0, -10, 0),
size=(33, 23, 10 * thickness),
grid_spec=td.GridSpec.auto(min_steps_per_wvl=20, wavelength=lda0),
structures=[optical_hybrid],
sources=[mode_source_LO],
monitors=[
field_monitor,
mode_monitor_Qn,
mode_monitor_Qp,
mode_monitor_In,
mode_monitor_Ip,
],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
medium=sio2,
symmetry=(0, 0, 1),
)
Visualize the simulation.
[15]:
sim.plot(z=0)
plt.show()
Submit the simulation job to the server.
[16]:
job = web.Job(simulation=sim, task_name="optical_hybrid_lo_input", verbose=True)
sim_data = job.run(path="data/simulation_data.hdf5")
19:55:23 CEST Created task 'optical_hybrid_lo_input' with task_id 'fdve-4a0b5965-8eea-470d-bd06-36a796a1fb69' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-4a0b5965-8e ea-470d-bd06-36a796a1fb69'.
19:55:27 CEST status = queued
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.
19:55:30 CEST status = preprocess
19:55:33 CEST Maximum FlexCredit cost: 0.927. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
starting up solver
19:55:34 CEST running solver
19:58:24 CEST early shutoff detected at 56%, exiting.
19:58:25 CEST status = postprocess
19:58:30 CEST status = success
View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-4a0b5965-8e ea-470d-bd06-36a796a1fb69'.
19:58:47 CEST loading simulation from data/simulation_data.hdf5
Visualize the field intensity distribution.
[17]:
sim_data.plot_field(field_monitor_name="field", field_name="E", val="abs^2", vmax=1500)
plt.show()
Besides the error in the relative phase, there are three additional figures of merit (FOM) the characterize the performance of the 90 degree optical hybrid: insertion loss, common mode rejection ratio (CMRR), and imbalance. We define a function to calculate those FOMs and plot the results.
The insertion loss is calculated as \(-10log_{10}(P_{tot})\), where \(P_{tot}=P_{Ip}+P_{In}+P_{Qp}+P_{Qn}\) is the total power at four output ports.
The CMRR is calculated as CMRR\(_I=-20log_{10}(\left|\frac{P_{Ip}-P_{In}}{P_{Ip}+P_{In}} \right|)\) and CMRR\(_Q=-20log_{10}(\left|\frac{P_{Qp}-P_{Qn}}{P_{Qp}+P_{Qn}} \right|)\).
The imbalance is defined as Imbalance\(_I = -10log_{10}(\frac{P_{Ip}}{P_{In}})\) and Imbalance\(_Q = -10log_{10}(\frac{P_{Qp}}{P_{Qn}})\).
[18]:
def calculate_FOM(sim_data):
# compute power at the Ip port
P_Ip = np.abs(sim_data["port_Ip"].amps.sel(direction="+").values) ** 2
# compute power at the In port
P_In = np.abs(sim_data["port_In"].amps.sel(direction="+").values) ** 2
# compute power at the Qp port
P_Qp = np.abs(sim_data["port_Qp"].amps.sel(direction="-").values) ** 2
# compute power at the Qn port
P_Qn = np.abs(sim_data["port_Qn"].amps.sel(direction="-").values) ** 2
fig, ax = plt.subplots(3, 1, tight_layout=True, figsize=(6, 10))
# compute and plot insertion loss
loss = -10 * np.log10(P_Ip + P_In + P_Qp + P_Qn)
ax[0].plot(ldas, loss)
ax[0].set_xlabel("Wavelength ($\mu m$)")
ax[0].set_ylabel("Insertion loss (dB)")
ax[0].set_ylim(0, 1)
# compute and plot CMRR
CMRR_I = -20 * np.log10(np.abs((P_Ip - P_In) / (P_Ip + P_In)))
CMRR_Q = -20 * np.log10(np.abs((P_Qp - P_Qn) / (P_Qp + P_Qn)))
ax[1].plot(ldas, CMRR_I, label="CMRR$_I$")
ax[1].plot(ldas, CMRR_Q, label="CMRR$_Q$")
ax[1].set_xlabel("Wavelength ($\mu m$)")
ax[1].set_ylabel("CMRR (dB)")
ax[1].set_ylim(25, 45)
ax[1].legend()
# compute and plot imbalance
Imbalance_I = -10 * np.log10(P_Ip / P_In)
Imbalance_Q = -10 * np.log10(P_Qp / P_Qn)
ax[2].plot(ldas, Imbalance_I, label="Imbalance$_I$")
ax[2].plot(ldas, Imbalance_Q, label="Imbalance$_Q$")
ax[2].set_xlabel("Wavelength ($\mu m$)")
ax[2].set_ylabel("Imbalance (dB)")
ax[2].set_ylim(-0.5, 0.5)
ax[2].legend()
Calculate those FOMs and plot the results. Due to symmetry, CMRR\(_I\) and CMRR\(_Q\) are identical. So is Imbalance\(_I\) and Imbalance\(_Q\).
Overall, we see that the design has a good performance. The insertion loss is about 0.3 dB while the CMRRs are above 30 dB and the Imbalances are around 0.1 dB.
[19]:
calculate_FOM(sim_data)
Simulation of the Optical Hybrid – Signal Input#
Finally, we simulate the case where the input signal is from the bottom signal port. To do so, we only need to define a ModeSource at the signal port and update the previous simulation since other settings of the simulation is identical.
[20]:
# define a mode source at the signal input waveguide at the bottom MMI
mode_source_signal = td.ModeSource(
center=(-0.5, -21, 0),
size=(2 * width, 0, 4 * thickness),
source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
direction="+",
mode_spec=mode_spec,
mode_index=0,
)
# copy simulation and change the source
sim = sim.copy(update={"sources": [mode_source_signal]})
sim.plot(z=0)
plt.show()
[21]:
job = web.Job(simulation=sim, task_name="optical_hybrid_signal_input", verbose=True)
sim_data = job.run(path="data/simulation_data.hdf5")
19:58:58 CEST Created task 'optical_hybrid_signal_input' with task_id 'fdve-2e992dcf-3da7-43d4-a23e-5dbf6da5911d' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-2e992dcf-3d a7-43d4-a23e-5dbf6da5911d'.
19:59:02 CEST status = queued
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.
19:59:06 CEST status = preprocess
19:59:08 CEST Maximum FlexCredit cost: 0.927. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
starting up solver
19:59:09 CEST running solver
20:01:53 CEST early shutoff detected at 60%, exiting.
20:01:54 CEST status = postprocess
20:02:00 CEST status = success
20:02:01 CEST View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-2e992dcf-3d a7-43d4-a23e-5dbf6da5911d'.
20:02:18 CEST loading simulation from data/simulation_data.hdf5
Visualize the field intensity distribution.
[22]:
sim_data.plot_field(field_monitor_name="field", field_name="E", val="abs^2", vmax=1500)
plt.show()
Calculate those FOMs and plot the results. Since the signal port is the left waveguide of the bottom MMI, CMRR\(_I\) and CMRR\(_Q\) are expected to be slghtly different. The device performance according to the insertion loss, CMRR, and Imbalance is comparable to the case with input from the LO port. In addition, from the phase simulation, we know the phase error is going to be small too. This means the overall design of the 90 degree optical hybrid works very well.
[23]:
calculate_FOM(sim_data)
[ ]: