Compact polarization splitter-rotator#

Run this notebook in your browser using Binder.

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).

123c892fd8d041499e4595a5c89edd73

Simulation Setup#

[1]:
import numpy as np
import matplotlib.pyplot as plt
import tidy3d as td
import tidy3d.web as web
from tidy3d.plugins import ModeSolver

[08:03:10] WARNING  This version of Tidy3D was pip installed from the 'tidy3d-beta' repository on   __init__.py:102
                    PyPI. Future releases will be uploaded to the 'tidy3d' repository. From now on,                
                    please use 'pip install tidy3d' instead.                                                       
           INFO     Using client version: 1.8.2                                                     __init__.py:120

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,
)

vertices = np.array(
    [
        (-w3 / 2 - w_gap - w4, L_tp1 + L_tp2 + L_tp3),
        (-w3 / 2 - w_gap - w4, L_tp1 + L_tp2 + L_tp3 + L_dc),
        (-w3 / 2 - w_gap - w4 - shift, L_tp1 + L_tp2 + L_tp3 + L_dc + L_tpout),
        (-w3 / 2 - w_gap - w4 - shift, inf_eff),
        (-w3 / 2 - w_gap - shift, inf_eff),
        (-w3 / 2 - w_gap - shift, L_tp1 + L_tp2 + L_tp3 + L_dc + L_tpout),
        (-w3 / 2 - w_gap, L_tp1 + L_tp2 + L_tp3 + L_dc),
        (-w3 / 2 - w_gap, L_tp1 + L_tp2 + L_tp3),
    ]
)
narrow_wg = td.Structure(
    geometry=td.PolySlab(vertices=vertices, axis=2, slab_bounds=(-h_co / 2, h_co / 2)),
    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
flux_monitor1 = td.FluxMonitor(
    center=(0, Ly - buffer_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, Ly - buffer_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, Ly - buffer_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, Ly - buffer_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")

../_images/notebooks_PolarizationSplitterRotator_11_0.png

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()

/Users/twhughes/.pyenv/versions/3.10.9/lib/python3.10/site-packages/numpy/linalg/linalg.py:2139: RuntimeWarning: divide by zero encountered in det
  r = _umath_linalg.det(a, signature=signature)
/Users/twhughes/.pyenv/versions/3.10.9/lib/python3.10/site-packages/numpy/linalg/linalg.py:2139: RuntimeWarning: invalid value encountered in det
  r = _umath_linalg.det(a, signature=signature)

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()

../_images/notebooks_PolarizationSplitterRotator_15_0.png

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()

../_images/notebooks_PolarizationSplitterRotator_17_0.png

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")
sim_data = job.run(path="data/simulation_data.hdf5")

[08:03:12] INFO     Using Tidy3D credentials from stored file.                                           auth.py:70
[08:03:13] INFO     Authentication successful.                                                           auth.py:30
[08:03:14] INFO     Created task 'polarization_splitter_rotator' with task_id                         webapi.py:120
                    'a42507b6-5605-4c8b-8efe-053ec80631ff'.                                                        
[08:03:16] INFO     Maximum FlexUnit cost: 0.692                                                      webapi.py:252
           INFO     status = queued                                                                   webapi.py:261
[08:03:20] INFO     status = preprocess                                                               webapi.py:273
[08:03:26] INFO     starting up solver                                                                webapi.py:277
[08:03:37] INFO     running solver                                                                    webapi.py:283
[08:04:39] INFO     early shutoff detected, exiting.                                                  webapi.py:294
           INFO     status = postprocess                                                              webapi.py:300
[08:04:53] INFO     status = success                                                                  webapi.py:306
           INFO     Billed FlexUnit cost: 0.361                                                       webapi.py:310
           INFO     downloading file "output/monitor_data.hdf5" to "data/simulation_data.hdf5"        webapi.py:592
[08:05:02] INFO     loading SimulationData from data/simulation_data.hdf5                             webapi.py:414

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", "int", ax=ax, f=freq0)
ax.set_aspect("auto")

../_images/notebooks_PolarizationSplitterRotator_22_0.png

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("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 = np.abs(mode_amp) ** 2 / T1
plt.plot(ldas, mode_power)
plt.xlim(1.5, 1.6)
plt.xlabel("Wavelength ($\mu m$)")
plt.ylabel("Power share at Port1 (%)")
plt.legend(["TE0", "TM0"])

[11]:
<matplotlib.legend.Legend at 0x2839a2410>
../_images/notebooks_PolarizationSplitterRotator_24_1.png

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 = td.ModeSource(
    center=(0, -buffer_y / 2, 0),
    size=(2 * w0, 0, 5 * h_co),
    source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10),
    direction="+",
    mode_spec=mode_spec,
    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")
sim_data = job.run(path="data/simulation_data.hdf5")

[08:05:06] INFO     Created task 'polarization_splitter_rotator' with task_id                         webapi.py:120
                    '8a476982-43ea-430b-9a7c-9998475338aa'.                                                        
[08:05:08] INFO     status = queued                                                                   webapi.py:261
[08:05:14] INFO     Maximum FlexUnit cost: 0.692                                                      webapi.py:252
[08:05:16] INFO     status = preprocess                                                               webapi.py:273
[08:05:22] INFO     starting up solver                                                                webapi.py:277
[08:05:33] INFO     running solver                                                                    webapi.py:283
[08:08:00] INFO     early shutoff detected, exiting.                                                  webapi.py:294
           INFO     status = postprocess                                                              webapi.py:300
[08:08:15] INFO     status = success                                                                  webapi.py:306
[08:08:16] INFO     Billed FlexUnit cost: 0.388                                                       webapi.py:310
           INFO     downloading file "output/monitor_data.hdf5" to "data/simulation_data.hdf5"        webapi.py:592
[08:08:21] INFO     loading SimulationData from data/simulation_data.hdf5                             webapi.py:414

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", "int", ax=ax, f=freq0)
ax.set_aspect("auto")

../_images/notebooks_PolarizationSplitterRotator_29_0.png

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("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 = np.abs(mode_amp) ** 2 / T2
plt.plot(ldas, mode_power)
plt.xlim(1.45, 1.6)
plt.xlabel("Wavelength ($\mu m$)")
plt.ylabel("Power share at Port1 (%)")
plt.legend(["TE0", "TM0"])

[14]:
<matplotlib.legend.Legend at 0x16a7e1cc0>
../_images/notebooks_PolarizationSplitterRotator_31_1.png

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.

[ ]: