Note: the cost of running the entire notebook is larger than 1 FlexCredit.
Efficient manipulation of the polarization state of light is essential for many photonic applications including optical communications, quantum information processing, and sensing. Polarization beam splitters and rotators (PSRs) are key devices that can split an input beam into two TE and TM output beams and rotate the TM beam into TE polarization. This notebook presents the simulation of a compact and highly efficient PSR based on 90 degree bends in a silicon-on-insulator platform. The design is proposed in the work Kang Tan, Ying Huang, Guo-Qiang Lo, Chengkuo Lee, and Changyuan Yu, "Compact highly-efficient polarization splitter and rotator based on 90° bends," Opt. Express 24, 14506-14512 (2016) DOI: 10.1364/OE.24.014506 and experimentally demonstrated in Kang Tan, Ying Huang, Guo-Qiang Lo, Changyuan Yu, and Chengkuo Lee, "Experimental realization of an O-band compact polarization splitter and rotator," Opt. Express 25, 3234-3241 (2017) DOI: 10.1364/OE.25.003234.
The proposed PSR consists of two parallel 90 degree bend waveguides with slightly different cross-sections designed for phase matching between the TM mode of the inner bend and the TE mode of the outer bend. The inner bend is designed to strongly confine the TE mode while allowing the TM mode to couple to the outer bend. By optimizing the geometry, efficient polarization splitting and rotating with low loss and high extinction ratio is achieved in a compact device footprint.
import gdstk
import matplotlib.pyplot as plt
import numpy as np
import tidy3d as td
import tidy3d.web as web
from tidy3d.plugins.mode import ModeSolver
from tidy3d.plugins.mode.web import run as run_mode_solver
Simulation Setup¶
Define simulation wavelength range to be 1.25 μm to 1.37 μm, covering the entire O-band.
lda0 = 1.31 # central wavelength
freq0 = td.C_0 / lda0 # central frequency
ldas = np.linspace(1.25, 1.37, 21) # wavelength range
freqs = td.C_0 / ldas # frequency range
fwidth = 0.5 * (np.max(freqs) - np.min(freqs)) # width of the source frequency
In this device, the PSR is made of silicon. The cladding and substrate are made of silicon oxide. We will simply use the popularly used refractive index data from Palik by calling Tidy3D's built-in material library.
# define silicon medium
si = td.material_library["cSi"]["Palik_Lossless"]
# define silicon oxide medium
sio2 = td.material_library["SiO2"]["Palik_Lossless"]
Define the geometric parameters. The inner bend is fully etched while the outer bend is partially etched. The cross section is schematically shown below.
R1 = 10 # radius of the inner bend
H1 = 0.22 # thickness of the fully etched waveguide
W1 = 0.4 # width of the inner bend waveguide
Wg = 0.2 # width of the gap
W2 = 0.21 # width of the fully etched outer waveguide
W3 = 0.285 # width of the partially etched outer waveguide
H2 = 0.11 # thickness of the partially etched waveguide
buffer = 5 # buffer spacing
Since the PSR contains various bending parts, the easiest way to define the geometry is to use the gdstk library.
# define a gds cell
cell = gdstk.Cell("psi")
# radius of the s bend
s_bend_radius = 6
# angle of the s bend
s_bend_angle = np.pi / 3
# define the inner bend path
inner_path = gdstk.RobustPath(
initial_point=(-buffer, 0), width=W1, tolerance=1e-4, layer=1, datatype=0
)
inner_path.horizontal(x=0)
inner_path.arc(radius=R1, initial_angle=-np.pi / 2, final_angle=0)
inner_path.arc(radius=s_bend_radius, initial_angle=0, final_angle=s_bend_angle)
inner_path.arc(radius=s_bend_radius, initial_angle=-np.pi + s_bend_angle, final_angle=-np.pi)
inner_path.vertical(y=buffer, relative=True)
cell.add(inner_path)
# starting coordinates of the outer bend path
x0 = -2
y0 = -2
# define the fully etched outer bend path
outer_path = gdstk.RobustPath(
initial_point=(x0, -buffer), width=W2, tolerance=1e-4, layer=1, datatype=0
)
outer_path.vertical(y=-W1 / 2 - Wg - W2 / 2 + y0)
outer_path.arc(radius=2, initial_angle=np.pi, final_angle=np.pi / 2)
outer_path.arc(radius=R1 + W1 / 2 + Wg + W2 / 2, initial_angle=-np.pi / 2, final_angle=0)
outer_path.vertical(y=buffer * 3, relative=True)
cell.add(outer_path)
# define the partially etched outer bend path
partially_etched_path = gdstk.RobustPath(
initial_point=(-2 + W3 / 2, -buffer), width=W3 - W2, tolerance=1e-4, layer=2, datatype=0
)
partially_etched_path.vertical(y=-W1 / 2 - Wg - W2 / 2 + y0)
partially_etched_path.arc(radius=2 - W3 / 2, initial_angle=np.pi, final_angle=np.pi / 2)
partially_etched_path.arc(
radius=R1 + W1 / 2 + Wg + W2 / 2 + W3 / 2, initial_angle=-np.pi / 2, final_angle=0
)
partially_etched_path.vertical(y=buffer * 3, relative=True)
cell.add(partially_etched_path)
# define structures from the gds
inner_bend = td.Structure(
geometry=td.PolySlab.from_gds(
cell,
gds_layer=1,
axis=2,
slab_bounds=(0, H1),
)[0],
medium=si,
)
outer_bend = td.Structure(
geometry=td.PolySlab.from_gds(
cell,
gds_layer=1,
axis=2,
slab_bounds=(0, H1),
)[1],
medium=si,
)
partially_etched_bend = td.Structure(
geometry=td.PolySlab.from_gds(
cell,
gds_layer=2,
axis=2,
slab_bounds=(0, H2),
)[0],
medium=si,
)
# plot the structures to visualize
ax = inner_bend.plot(z=H2 / 2)
ax = outer_bend.plot(z=H2 / 2, ax=ax)
partially_etched_bend.plot(z=H2 / 2, ax=ax)
# define simulation domain bounds
x_min = -3
x_max = 12
y_min = -3
y_max = 22
z_min = -1
z_max = 1
ax.set_xlim(x_min, x_max)
ax.set_ylim(y_min, y_max)
plt.show()
Define a ModeSource first to launch the TM mode at the input waveguide. Later, we will also simulate a TE mode source. mode_index=1 specifies that we are selecting the TM0 mode while mode_index=0 corresponds to the TE0 mode since they are ordered by their effective index values in a decreasing order.
To measure mode conversion efficiencies at the through and cross ports, we define two ModeMonitors. We also define a FieldMonitor to help visualize the mode propagation and conversion.
# add a mode source as excitation
mode_spec = td.ModeSpec(num_modes=2, target_neff=3.5)
mode_source = td.ModeSource(
center=(-lda0, 0, H1 / 2),
size=(0, 4 * W1, 6 * H1),
source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
direction="+",
mode_spec=mode_spec,
mode_index=1,
)
# add a mode monitor to measure transmission at the cross port
mode_monitor_cross = td.ModeMonitor(
center=(R1 + W1 / 2 + Wg + W2 / 2 + W3 / 2, R1 + 2 * buffer, H1 / 2),
size=(8 * W1, 0, 6 * H1),
freqs=freqs,
mode_spec=mode_spec,
name="cross",
)
# add a mode monitor to measure transmission at the through port
mode_monitor_through = td.ModeMonitor(
center=(4, R1 + 2 * buffer, H1 / 2),
size=(8 * W1, 0, 6 * H1),
freqs=freqs,
mode_spec=mode_spec,
name="through",
)
# add a field monitor to visualize field distribution at z=t/2
field_monitor = td.FieldMonitor(
center=(0, 0, H2 / 2), size=(td.inf, td.inf, 0), freqs=[freq0], name="field"
)
Define a Tidy3D simulation. Remember to set medium=sio2, which makes the background medium silicon oxide instead of air. To ensure good accuracy, we use an automatic nonuniform grid with min_steps_per_wvl=20.
run_time = 1e-12 # simulation run time
# define a simulation domain box
sim_box = td.Box.from_bounds(rmin=(x_min, y_min, z_min), rmax=(x_max, y_max, z_max))
# construct simulation
sim_tm = td.Simulation(
center=sim_box.center,
size=sim_box.size,
grid_spec=td.GridSpec.auto(min_steps_per_wvl=20, wavelength=lda0),
structures=[inner_bend, outer_bend, partially_etched_bend],
sources=[mode_source],
monitors=[mode_monitor_cross, mode_monitor_through, field_monitor],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
medium=sio2,
)
# plot simulation
sim_tm.plot(z=H2 / 2)
plt.show()
To have a better visualization, we can also plot the simulation in 3D.
sim_tm.plot_3d()
Mode Analysis¶
To analyze the working principle of the design, we need to perform mode solving for the outer bend waveguide. To ensure sufficient accuracy, we will run the mode solving in the server, where subpixel averaging is applied. First, solve for two modes and plot the mode profiles.
mode_solver = ModeSolver(
simulation=sim_tm,
plane=td.Box(center=mode_monitor_cross.center, size=mode_monitor_cross.size),
mode_spec=mode_spec,
freqs=freqs,
)
mode_data_outer = run_mode_solver(mode_solver)
def plot_field_intensity(mode_data):
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 3), tight_layout=True)
cmap = "magma"
mode_data.intensity.sel(mode_index=0, f=freq0).plot(x="x", y="z", cmap=cmap, ax=ax1)
mode_data.intensity.sel(mode_index=1, f=freq0).plot(x="x", y="z", cmap=cmap, ax=ax2)
ax1.set_title("First mode")
ax2.set_title("Second mode")
plot_field_intensity(mode_data_outer)
12:33:48 CEST Mode solver created with task_id='fdve-c57b6d8d-ff2f-46e8-aa53-d40725e5ead9', solver_id='mo-09590491-ffa2-406a-a4fe-1c9e900e083b'.
12:33:52 CEST Mode solver status: queued
12:33:54 CEST Mode solver status: running
12:34:14 CEST Mode solver status: success
We can also plot the field components to better visualize the polarizations of the modes. The first mode is TE-like as the dominant field component is $E_x$ while the second mode is TM-like with the dominant component being $E_z$
def plot_field_components(mode_data):
f, ax = plt.subplots(2, 3, figsize=(10, 5), tight_layout=True)
cmap = "coolwarm"
abs(mode_data.Ex.sel(mode_index=0, f=freq0)).plot(x="x", y="z", ax=ax[0][0], cmap=cmap)
abs(mode_data.Ey.sel(mode_index=0, f=freq0)).plot(x="x", y="z", ax=ax[0][1], cmap=cmap)
abs(mode_data.Ez.sel(mode_index=0, f=freq0)).plot(x="x", y="z", ax=ax[0][2], cmap=cmap)
abs(mode_data.Ex.sel(mode_index=1, f=freq0)).plot(x="x", y="z", ax=ax[1][0], cmap=cmap)
abs(mode_data.Ey.sel(mode_index=1, f=freq0)).plot(x="x", y="z", ax=ax[1][1], cmap=cmap)
abs(mode_data.Ez.sel(mode_index=1, f=freq0)).plot(x="x", y="z", ax=ax[1][2], cmap=cmap)
ax[0][0].set_title("|Ex|")
ax[0][1].set_title("|Ey|")
ax[0][2].set_title("|Ez|")
ax[1][0].set_title("|Ex|")
ax[1][1].set_title("|Ey|")
ax[1][2].set_title("|Ez|")
plot_field_components(mode_data_outer)
Efficient mode conversion from the inner bend to the outer bend is achieved when the phase matching condition is met. For waveguide bends, the phase matching condition is given by the equal optical path length (OPL). That is, $N_1k_0R_1\theta = N_2k_0R_2\theta$, where $N_1$ and $N_2$ are the effective indices of the inner bend mode and outer bend mode, $k_0$ is the free space wave vector, $R_1$ and $R_2$ are the inner bend radius and outer bend radius, and $\theta$ is the angle of the bend. Let's check the OPLs of the second mode (TM) at the inner waveguide to the first mode (TE) at the outer waveguide to see if they match. To do so, first compute the OPLs for the outer waveguide as we have solved for the modes.
# extract effective indices
n_eff_outer = mode_data_outer.n_eff.values
# compute OPLs
OPL_outer = n_eff_outer * (R1 + W1 / 2 + Wg + W3 / 2)
Then we will do the same for the inner bend.
mode_solver = ModeSolver(
simulation=sim_tm,
plane=td.Box(center=mode_monitor_through.center, size=mode_monitor_through.size),
mode_spec=mode_spec,
freqs=freqs,
)
mode_data_inner = run_mode_solver(mode_solver)
plot_field_intensity(mode_data_inner)
12:34:20 CEST Mode solver created with task_id='fdve-9778b3df-a6c3-421a-9231-26a3f1e91328', solver_id='mo-4399ab2f-8719-4f75-9e4e-5c4f0977f010'.
12:34:24 CEST Mode solver status: queued
12:34:38 CEST Mode solver status: running
12:35:02 CEST Mode solver status: success
plot_field_components(mode_data_inner)
Calculate the OPLs for the inner waveguide and plot them as well as the OPLs for the outer waveguide. We see that the OPL of the TM mode on the inner bend (red curve) closely matches the OPL of the TE mode on the outer bend (blue dashed curve).
n_eff_inner = mode_data_inner.n_eff.values
OPL_inner = n_eff_inner * R1
# plot calculated OPLs
plt.plot(ldas, OPL_outer[:, 0], linewidth=3, linestyle="--", label="TE mode on outer bend")
plt.plot(ldas, OPL_outer[:, 1], linewidth=3, linestyle="--", label="TM mode on outer bend")
plt.plot(ldas, OPL_inner[:, 0], linewidth=3, label="TE mode on inner bend")
plt.plot(ldas, OPL_inner[:, 1], linewidth=3, label="TM mode on inner bend")
plt.legend()
plt.ylim(16, 30)
plt.xlabel(r"Wavelength ($\mu m$)")
plt.ylabel("OPL (a.u.)")
plt.show()
Running the TM Simulation¶
Now that we have confirmed the design from a mode analysis perspective, we need to verify the device performance using 3D FDTD. We are ready to submit the simulation job to the server. Before running the simulation, we can get a cost estimation using estimate_cost. This prevents us from accidentally running large jobs that we set up by mistake. The estimated cost is the maximum cost corresponding to running all the time steps.
job_tm = web.Job(simulation=sim_tm, task_name="psr_tm", verbose=True)
estimated_cost = web.estimate_cost(job_tm.task_id)
12:35:08 CEST Created task 'psr_tm' with task_id 'fdve-c6bcbef1-0fcb-45d6-9b6a-cbdd217b11b8' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-c6bcbef1-0f cb-45d6-9b6a-cbdd217b11b8'.
Task folder: 'default'.
12:35:10 CEST Maximum FlexCredit cost: 1.889. Minimum cost depends on task execution details. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
12:35:11 CEST Maximum FlexCredit cost: 1.889. Minimum cost depends on task execution details. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
The cost is reasonable so we can run the simulation.
sim_data_tm = job_tm.run(path="data/simulation_data.hdf5")
12:35:12 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.
12:35:22 CEST status = preprocess
12:35:26 CEST starting up solver
running solver
12:37:28 CEST early shutoff detected at 64%, exiting.
status = postprocess
12:37:33 CEST status = success
12:37:35 CEST View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-c6bcbef1-0f cb-45d6-9b6a-cbdd217b11b8'.
12:37:43 CEST loading simulation from data/simulation_data.hdf5
TM Result Visualization¶
After the simulation is complete, we first plot the field distribution. From the field distribution, we observe an efficient coupling from the inner bend to the outer bend.
sim_data_tm.plot_field(
field_monitor_name="field", field_name="E", val="abs", f=freq0, vmin=0, vmax=70
)
plt.show()
More quantitatively, we plot the mode conversion efficiencies. Indeed, a high TM to TE mode conversion at the cross port is observed as desired.
def plot_conversion_efficiency(sim_data, pol):
# extract the transmission data from the mode monitor
T_cross_te = np.abs(sim_data["cross"].amps.sel(mode_index=0, direction="+")) ** 2
T_cross_tm = np.abs(sim_data["cross"].amps.sel(mode_index=1, direction="+")) ** 2
T_through_te = np.abs(sim_data["through"].amps.sel(mode_index=0, direction="+")) ** 2
T_through_tm = np.abs(sim_data["through"].amps.sel(mode_index=1, direction="+")) ** 2
# plot transmission
plt.plot(ldas, 10 * np.log10(T_cross_te), linewidth=3, label=f"{pol}->TE (cross)")
plt.plot(ldas, 10 * np.log10(T_cross_tm), linewidth=3, label=f"{pol}->TM (cross)")
plt.plot(
ldas, 10 * np.log10(T_through_te), linewidth=3, linestyle="--", label=f"{pol}->TE (through)"
)
plt.plot(
ldas, 10 * np.log10(T_through_tm), linewidth=3, linestyle="--", label=f"{pol}->TM (through)"
)
plt.legend()
plt.xlabel(r"Wavelength ($\mu m$)")
plt.ylabel("Mode conversion efficiency (dB)")
plt.title(f"Mode conversion efficiency for {pol} mode")
plt.show()
plot_conversion_efficiency(sim_data_tm, "TM")
Running the TE Simulation¶
Now that we have confirmed the TM to TE polarization rotation of the device, we need to confirm the TE-TE conversion at the through port. The TE mode at the inner waveguide will not couple to the outer waveguide due to the phase matching condition is not met as shown above.
To run the TE simulation, we only need to update the mode_index from 1 to 0 in the ModeSource.
mode_source = mode_source.copy(update={"mode_index": 0})
sim_te = sim_tm.copy(update={"sources": [mode_source]})
sim_data_te = web.run(simulation=sim_te, task_name="psr_te", path="data/simulation_data.hdf5")
12:37:47 CEST Created task 'psr_te' with task_id 'fdve-b07a8a68-4c2b-4452-a101-6d8fa8db3453' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-b07a8a68-4c 2b-4452-a101-6d8fa8db3453'.
Task folder: 'default'.
12:37:49 CEST Maximum FlexCredit cost: 1.889. Minimum cost depends on task execution details. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
12:37:50 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.
12:37:57 CEST status = preprocess
12:38:02 CEST starting up solver
running solver
12:39:39 CEST early shutoff detected at 96%, exiting.
status = postprocess
12:39:44 CEST status = success
12:39:46 CEST View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-b07a8a68-4c 2b-4452-a101-6d8fa8db3453'.
12:39:50 CEST loading simulation from data/simulation_data.hdf5
TE Result Visualization¶
Field distribution shows minimal coupling to the outer waveguide.
sim_data_te.plot_field(
field_monitor_name="field", field_name="E", val="abs", f=freq0, vmin=0, vmax=70
)
plt.show()
As expected, the input TE mode stays at the inner waveguide as TE polarization. Both simulations confirm the good performance of the PSR design.
plot_conversion_efficiency(sim_data_te, "TE")
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. 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.
TERMS & CONDITIONS
How the process works:
Try this real URL to see an example: