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 numpy as np
import matplotlib.pyplot as plt
import gdstk
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"
)
[13:59:43] WARNING: Default value for the field monitor monitor.py:261 'colocate' setting has changed to 'True' in Tidy3D 2.4.0. All field components will be colocated to the grid boundaries. Set to 'False' to get the raw fields on the Yee grid instead.
Define a Tidy3D simulation. Rememeber 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)
Output()
Output()
Output()
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)