Silicon-on-insulator (SOI) devices are known to be polarization sensitive. Devices that can manipulate the polarization of light are important components of an integrated photonic system.

This model demonstrates the design of a compact polarization splitter-rotator that consists of an adiabatic waveguide tapper and an asymmetric directional coupler. When the TM0 mode is launched at the input end, it is efficiently converted into the TE1 mode at the tapper and then coupled to the TE0 mode at the narrow waveguide through the directional coupler. When the TE0 mode is launched at the input end, it propagates through the tapper without polarization change and coupling to the narrow waveguide. That is, the TE and TM polarizations are separated by this device and the output is always TE polarization, as schematically shown below. This model is based on Daoxin Dai and John E. Bowers, "Novel concept for ultracompact polarization splitter-rotator based on silicon nanowires," Opt. Express 19, 10940-10949 (2011).

For more simulation examples, 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.

Simulation Setup¶

[1]:

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

Define geometric parameters. The device consists of a wide tapered waveguide and a narrow waveguide. They are coupled through a directional coupler.

[2]:

w0 = 0.54  # width of the input/output single mode waveguides
w1 = 0.69  # width of the first tapper
w2 = 0.83  # width of the second tapper
w3 = 0.9  # width of the third tapper
w4 = 0.405  # width of the narrow waveguide
w_gap = 0.15  # gap of the directional coupler
L_tp1 = 4  # length of the first tapper
L_tp2 = 44  # length of the second tapper
L_tp3 = L_tp1 * (w3 - w2) / (w1 - w0)  # length of the third tapper
L_dc = 7  # length of the directional coupler
L_tpout = 14  # length of the output tapper
shift = 0.4  # shift of the narrow waveguide output
h_co = 0.22  # thickness of the waveguides
inf_eff = 1000  # effective infinity used to make sure the waveguides extend into pml

Define materials. The waveguides are made of silicon. The upper cladding is silicon nitride and the lower cladding is silicon oxide.

[3]:

si = td.Medium(permittivity=3.455**2)
sio2 = td.Medium(permittivity=1.445**2)
si3n4 = td.Medium(permittivity=2**2)

The silicon structures are defined using PolySlabs. The coordinates of the vertices can be determined by the taper widths and lengths.

[4]:

cladding = td.Structure(
    geometry=td.Box.from_bounds(
        rmin=(-inf_eff, -inf_eff, -h_co / 2), rmax=(inf_eff, inf_eff, inf_eff)
    ),
    medium=si3n4,
)

vertices = np.array(
    [
        (-w0 / 2, -inf_eff),
        (-w0 / 2, 0),
        (-w1 / 2, L_tp1),
        (-w2 / 2, L_tp1 + L_tp2),
        (-w3 / 2, L_tp1 + L_tp2 + L_tp3),
        (-w3 / 2, L_tp1 + L_tp2 + L_tp3 + L_dc),
        (-w0 / 2, L_tp1 + L_tp2 + L_tp3 + L_dc + L_tpout),
        (-w0 / 2, inf_eff),
        (w0 / 2, inf_eff),
        (w0 / 2, L_tp1 + L_tp2 + L_tp3 + L_dc + L_tpout),
        (w3 / 2, L_tp1 + L_tp2 + L_tp3 + L_dc),
        (w3 / 2, L_tp1 + L_tp2 + L_tp3),
        (w2 / 2, L_tp1 + L_tp2),
        (w1 / 2, L_tp1),
        (w0 / 2, 0),
        (w0 / 2, -inf_eff),
    ]
)

wide_wg = td.Structure(
    geometry=td.PolySlab(vertices=vertices, axis=2, slab_bounds=(-h_co / 2, h_co / 2)),
    medium=si,
)

R = 100
cell = gdstk.Cell("bend")  # define a gds cell
bend = gdstk.FlexPath((-w3 / 2 - w_gap - w4 / 2, L_tp1 + L_tp2 + L_tp3), w4, layer=1, datatype=0)
bend.vertical(L_dc, relative=True)
bend.arc(R, 0, np.pi / 50)
bend.arc(R, -np.pi + np.pi / 50, -np.pi)
bend.vertical(inf_eff)
cell.add(bend)

# define the waveguide bend tidy3d geometries
bend_geo = td.PolySlab.from_gds(
    cell,
    gds_layer=1,
    axis=2,
    slab_bounds=(-h_co / 2, h_co / 2),
)[0]

narrow_wg = td.Structure(geometry=bend_geo, medium=si)

Set up a mode source for excitation. First, we launch the TE0 mode with the mode source. Later, we will change the mode source to launch the TM0 mode.

Two flux monitors and two mode monitors are set up at the outputs of the wide and narrow waveguides. A field monitor is added to monitor the field at z=0 plane.

[5]:

lda0 = 1.525  # central wavelength
freq0 = td.C_0 / lda0  # central frequency
ldas = np.linspace(1.45, 1.6, 101)  # wavelength range
freqs = td.C_0 / ldas  # frequency range

# buffer lengths in x and y directions
buffer_x = 1
buffer_y = 2

# simulation domain size
Lx = w3 + w_gap + w4 + shift + 2 * buffer_x
Ly = L_tp1 + L_tp2 + L_tp3 + L_dc + L_tpout + 2 * buffer_y
Lz = 10 * h_co
sim_size = (Lx, Ly, Lz)

# define mode source that launches the TE0 mode (mode_index=0). Later, we will modify it to investigate the TM0 mode case
mode_spec = td.ModeSpec(num_modes=2, target_neff=3)
mode_source = td.ModeSource(
    center=(0, -buffer_y / 2, 0),
    size=(3 * w0, 0, 5 * h_co),
    source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10),
    direction="+",
    mode_spec=mode_spec,
    mode_index=0,
)

# define a field monitor
field_monitor = td.FieldMonitor(
    center=(0, -buffer_y / 2, 0), size=(td.inf, td.inf, 0), freqs=[freq0], name="field"
)

# define two flux monitors at the two outputs to measure transmission
output_y = Ly - 3 * buffer_y / 2
flux_monitor1 = td.FluxMonitor(
    center=(0, output_y, 0), size=(3 * w0, 0, 5 * h_co), freqs=freqs, name="flux1"
)

flux_monitor2 = td.FluxMonitor(
    center=(-w4 / 2 - w3 / 2 - w_gap - shift, output_y, 0),
    size=(3 * w0, 0, 5 * h_co),
    freqs=freqs,
    name="flux2",
)

# define two mode monitors at the two outputs to study output polarization
mode_monitor1 = td.ModeMonitor(
    center=(0, output_y, 0),
    size=(3 * w0, 0, 5 * h_co),
    freqs=freqs,
    mode_spec=mode_spec,
    name="mode1",
)

mode_monitor2 = td.ModeMonitor(
    center=(-w4 / 2 - w3 / 2 - w_gap - shift, output_y, 0),
    size=(3 * w0, 0, 5 * h_co),
    freqs=freqs,
    mode_spec=mode_spec,
    name="mode2",
)

# initialize the Simulation object
sim = td.Simulation(
    center=(-(shift + w_gap) / 2, Ly / 2 - buffer_y, 0),
    size=sim_size,
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=20, wavelength=lda0),
    structures=[cladding, wide_wg, narrow_wg],
    sources=[mode_source],
    monitors=[
        field_monitor,
        flux_monitor1,
        flux_monitor2,
        mode_monitor1,
        mode_monitor2,
    ],
    run_time=2e-12,
    boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
    medium=sio2,
)

# plot the simulation at z=0 to inspect the structure, source, and monitors
fig = plt.figure()
ax = fig.add_subplot(111)
sim.plot(z=0, ax=ax)
ax.set_aspect("auto")

No description has been provided for this image

Before the simulation, we can visualize the mode fields to ensure we are launching the correct mode at the source.

[6]:

mode_solver = ModeSolver(
    simulation=sim,
    plane=td.Box(center=(0, -buffer_y / 2, 0), size=(2 * w0, 0, 5 * h_co)),
    mode_spec=mode_spec,
    freqs=[freq0],
)
mode_data = mode_solver.solve()
mode_data.to_dataframe()

21:07:56 CEST WARNING: Use the remote mode solver with subpixel averaging for   
              better accuracy through 'tidy3d.plugins.mode.web.run(...)'.

[6]:

		wavelength	n eff	k eff	TE (Ex) fraction	wg TE fraction	wg TM fraction	mode area
f	mode_index
1.965852e+14	0	1.525	2.608284	0.0	0.990227	0.845346	0.831149	0.201091
1.965852e+14	1	1.525	2.159161	0.0	0.041372	0.754310	0.908583	0.352458

As expected, the first mode (mode_index=0) is the TE0 mode.

[7]:

f, (ax1, ax2, ax3) = plt.subplots(1, 3, tight_layout=True, figsize=(10, 3))
abs(mode_data.Ex.isel(mode_index=0)).plot(x="x", y="z", ax=ax1, cmap="magma")
abs(mode_data.Ey.isel(mode_index=0)).plot(x="x", y="z", ax=ax2, cmap="magma")
abs(mode_data.Ez.isel(mode_index=0)).plot(x="x", y="z", ax=ax3, cmap="magma")

ax1.set_title("|Ex(x, y)|")
ax2.set_title("|Ey(x, y)|")
ax3.set_title("|Ez(x, y)|")
plt.show()

The second mode (mode_index=1) is the TM0 mode.

[8]:

f, (ax1, ax2, ax3) = plt.subplots(1, 3, tight_layout=True, figsize=(10, 3))
abs(mode_data.Hx.isel(mode_index=1)).plot(x="x", y="z", ax=ax1, cmap="magma")
abs(mode_data.Hy.isel(mode_index=1)).plot(x="x", y="z", ax=ax2, cmap="magma")
abs(mode_data.Hz.isel(mode_index=1)).plot(x="x", y="z", ax=ax3, cmap="magma")

ax1.set_title("|Hx(x, y)|")
ax2.set_title("|Hy(x, y)|")
ax3.set_title("|Hz(x, y)|")
plt.show()

After making sure the simulation setups and mode profiles are correct, submit the simulation to the server.

[9]:

job = web.Job(simulation=sim, task_name="polarization_splitter_rotator", verbose=True)
sim_data = job.run(path="data/simulation_data.hdf5")

21:07:58 CEST Created task 'polarization_splitter_rotator' with task_id         
              'fdve-2ce31ffa-0e8f-4431-a941-9bb41b7c2531' and task_type 'FDTD'.

              View task using web UI at                                         
              'https://tidy3d.simulation.cloud/workbench?taskId=fdve-2ce31ffa-0e
              8f-4431-a941-9bb41b7c2531'.

21:08:01 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.

21:08:11 CEST status = preprocess

21:08:13 CEST Maximum FlexCredit cost: 0.729. Use 'web.real_cost(task_id)' to   
              get the billed FlexCredit cost after a simulation run.

              starting up solver

21:08:14 CEST running solver

21:10:08 CEST early shutoff detected at 52%, exiting.

              status = postprocess

21:10:11 CEST status = success

              View simulation result at                                         
              'https://tidy3d.simulation.cloud/workbench?taskId=fdve-2ce31ffa-0e
              8f-4431-a941-9bb41b7c2531'.

21:10:23 CEST loading simulation from data/simulation_data.hdf5

Case 1: Launching the TE0 Mode at the Input¶

Visualize the field intensity distribution to see the propagation of the input TE0 mode. We can see that it propagates through the wide waveguide with no coupling to the narrow waveguide.

[10]:

fig = plt.figure()
ax = fig.add_subplot(111)
sim_data.plot_field(field_monitor_name="field", field_name="E", val="abs^2", ax=ax, f=freq0)
ax.set_aspect("auto")

Plot the transmission at the two outputs. The transmission through the wide waveguide is nearly 100% with very little coupling to the narrow waveguide.

Then plot the mode composition at the wide waveguide output. We can see that the TE polarization is preserved with no conversion to TM.

[11]:

f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 5))
T1 = sim_data["flux1"].flux
T2 = sim_data["flux2"].flux

plt.sca(ax1)
plt.plot(ldas, T1, ldas, T2)
plt.xlim(1.5, 1.6)
plt.xlabel(r"Wavelength ($\mu m$)")
plt.ylabel("Transmission")
plt.legend(("Wide waveguide", "Narrow waveguide"))

plt.sca(ax2)
mode_amp = sim_data["mode1"].amps.sel(direction="+")
mode_power_share = 100 * np.abs(mode_amp) ** 2 / T1
plt.plot(ldas, mode_power_share)
plt.xlim(1.5, 1.6)
plt.xlabel(r"Wavelength ($\mu m$)")
plt.ylabel("Power share at Port1 (%)")
plt.legend(["TE0", "TM0"])
plt.show()

Case 2: Launching the TM0 Mode at the Input¶

Next, we investigate the case where the TM0 mode is launched at the input. To set up the simulation, we simply need to copy the previous simulation and update the mode source.

[12]:

mode_source = mode_source.copy(update={"mode_index": 1})  # mode_index=1 corresponds to the TM0 mode

sim = sim.copy(update={"sources": [mode_source]})
job = web.Job(simulation=sim, task_name="polarization_splitter_rotator", verbose=True)
sim_data = job.run(path="data/simulation_data.hdf5")

21:10:30 CEST Created task 'polarization_splitter_rotator' with task_id         
              'fdve-a44a3907-eeba-4b48-b151-4da448bd4a55' and task_type 'FDTD'.

              View task using web UI at                                         
              'https://tidy3d.simulation.cloud/workbench?taskId=fdve-a44a3907-ee
              ba-4b48-b151-4da448bd4a55'.

21:10:33 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.

21:10:36 CEST status = preprocess

21:10:39 CEST Maximum FlexCredit cost: 0.729. Use 'web.real_cost(task_id)' to   
              get the billed FlexCredit cost after a simulation run.

              starting up solver

              running solver

21:11:57 CEST early shutoff detected at 84%, exiting.

              status = postprocess

21:11:59 CEST status = success

21:12:00 CEST View simulation result at                                         
              'https://tidy3d.simulation.cloud/workbench?taskId=fdve-a44a3907-ee
              ba-4b48-b151-4da448bd4a55'.

21:12:06 CEST loading simulation from data/simulation_data.hdf5

Visualize the field intensity distribution to see the propagation of the input TM0 mode. We can see that the TM0 mode is efficiently converted to the TE1 mode and then coupled to the narrow waveguide through the directional coupler region.

[13]:

fig = plt.figure()
ax = fig.add_subplot(111)
sim_data.plot_field(field_monitor_name="field", field_name="E", val="abs^2", ax=ax, f=freq0)
ax.set_aspect("auto")

Plot the transmission at the two outputs. At the central wavelength, above 90% of the power is coupled to the narrow waveguide.

Then plot the mode composition at the narrow waveguide output. We can see that the TE polarization is dominant, indicating a good TM to TE mode conversion efficiency.

[14]:

f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 5))
T1 = sim_data["flux1"].flux
T2 = sim_data["flux2"].flux

plt.sca(ax1)
plt.plot(ldas, T1, ldas, T2)
plt.xlim(1.45, 1.6)
plt.xlabel(r"Wavelength ($\mu m$)")
plt.ylabel("Transmission")
plt.legend(("Wide waveguide", "Narrow waveguide"))

plt.sca(ax2)
mode_amp = sim_data["mode2"].amps.sel(direction="+")
mode_power_share = 100 * np.abs(mode_amp) ** 2 / T2
plt.plot(ldas, mode_power_share)
plt.xlim(1.45, 1.6)
plt.xlabel(r"Wavelength ($\mu m$)")
plt.ylabel("Power share at Port2 (%)")
plt.legend(["TE0", "TM0"])
plt.show()

Simulation results from the two cases confirm that we can use this device to separate the TE and TM polarizations as well as convert the TM polarization to TE.

[ ]: