Broadband polarizer assisted by anisotropic metamaterial#

Silicon photonic polarizers are specialized optical devices designed to filter, control, and manipulate the polarization state of light in integrated photonic circuits. Comprising carefully engineered silicon waveguide structures, these polarizers play a crucial role in ensuring the efficient and reliable operation of various optical communication and sensing systems. Over the years, various silicon photonic polarizer designs have been proposed. However, most designs suffer from a narrow working bandwidth or a large device footprint.

A polarizer design demonstrated by Hongnan Xu, Daoxin Dai, and Yaocheng Shi, "Anisotropic metamaterial-assisted all-silicon polarizer with 415-nm bandwidth," *Photon. Res.* 7, 1432-1439 (2019), DOI: 10.1364/PRJ.7.001432 utilizes anisotropic subwavelength grating (SWG) metamaterial cladding to a 180 degree waveguide bend. This polarizer design preserves the TE mode and filters out TM mode in a broad wavelength range of 400 nm. The footprint of the device is also more compact than many of the previous designs.

Schematic of the polarizer

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.

[1]:
import numpy as np
import matplotlib.pyplot as plt

import gdstk
import tidy3d as td
import tidy3d.web as web

Effective Anisotropic Refractive Index Tensor of the SWG#

For this device, first we are going to make a regular 180 degree waveguide bend. In the bending region, we are going to add an SWG metamaterial parallel to the bend. The SWG metamaterial acts as an effective medium with an anisotropic refractive index tensor. The ordinary and extraordinary refractive indices are given by the well known Rytov’s approximation: \(1/n_e^2=f/n_{si}^2+(1-f)/n_{air}^2\) and \(n_o^2=fn_{si}^2+(1-f)n_{air}^2\), where \(f\) is the duty cycle of the SWG. We can investigate how the effective indices vary with \(f\).

From the result, we notice that \(n_0\) is much larger than \(n_e\) when \(f\) is not close to 0 or 1. Since the TE waveguide mode has the electric field primarily in the extraordinary direction while the TM mode has the electric field primarily in the ordinary direction, they will feel very different refractive index in the SWG metamaterial. For TE mode, the refractive index is much lower than \(n_{si}\) so it should stay guided inside the waveguide. For TM mode, on the other hand, the effective refractive index is close to \(n_{si}\) so it will be guided into the SWG metamaterial region and radiates out. This way, the waveguide bend together with the SWG metamaterial functions as a polarizer that passes TE mode and suppresses TM mode.

For a more detailed discussion, please refer to the reference.

[2]:
n_si = 3.47

f = np.linspace(0, 1, 30)
n_e = np.sqrt(1 / (f / n_si**2 + (1 - f)))
n_o = np.sqrt(f * n_si**2 + (1 - f))
plt.plot(f, n_e, label="$n_e$")
plt.plot(f, n_o, label="$n_o$")
plt.xlabel("$f$")
plt.ylabel("Effective index")
plt.legend()
plt.show()
../_images/notebooks_SWGBroadbandPolarizer_5_0.png

Structure Setup#

The simulation wavelength range is from 1260 nm to 1675 nm.

[3]:
lda_min = 1.26  # minimal wavelength
lda_max = 1.675  # maximal wavelength
lda0 = (lda_min + lda_max) / 2  # central wavelength
freqs = np.linspace(td.C_0 / lda_max, td.C_0 / lda_min, 301)  # frequency range
freq0 = td.C_0 / lda0  # central frequency

fwidth = 0.4 * (td.C_0 / lda_min - td.C_0 / lda_max)  # width of the frequency range of the source

Next, we define the media used in the simulation. Since the simulation is very broadband, it is better to use dispersive media. We can use the crystalline silicon and silicon dioxide from Tidy3D’s built-in material library.

[4]:
si = td.material_library["cSi"]["Li1993_293K"]
sio2 = td.material_library["SiO2"]["Horiba"]

The silicon layer has a thickness of 250 nm. To stay at single mode operation, the waveguide width is chosen to be 360 nm. The design parameters such as waveguide bend radius, period of the SWG, maximal and minimal duty cycles, and the number of SWG layers need to be optimized for the design. This is explored in the reference. Here we report the optimized parameters.

[5]:
h_si = 0.25  # silicon layer thickness
w_wg = 0.36  # waveguide width
R_0 = 3.5  # radius of the waveguide bend
p = 0.3  # period of the SWG
f_M = 0.75  # maximal duty cycle
f_m = 0.2  # minimal duty cycle
N = 9  # number of SWG layers
inf_eff = 1e3  # effective infinity

The radius of the \(i\)-th SWG layer is given by \(R_i = R_0 + w_wg/2 + p(i-1/2)\) and the duty cycle of the layer is also given by a linear ramping \(f_i=f_M+(f_m-f_M)(i-1)/(N-1)\). The width of the waveguide in each SWG layer is thus given by \(pf_i\).

[6]:
R_is = [R_0 + w_wg / 2 + p * (i - 0.5) for i in range(1, N + 1)]  # radius of each SWG layer
f_is = [
    f_M + (f_m - f_M) * (i - 1) / (N - 1) for i in range(1, N + 1)
]  # duty cycle of each SWG layer
w_is = [p * f for f in f_is]  # waveguide width of each SWG layer

# insert the circular waveguide band into the list
R_is.insert(0, R_0)
w_is.insert(0, w_wg)

With all the parameters determined, we can set up the geometries. One easy way is to make use of the gdstk library as demonstrated below.

[7]:
cell = gdstk.Cell("SWG")  # define a gds cell

# define a path of each swg layer
for R_i, w_i in zip(R_is, w_is):
    points_x = np.linspace(-R_i, R_i, 1001)  # x coordinates of the waveguide in the swg layer
    points_y = np.sqrt(R_i**2 - points_x**2)  # y coordinates

    # if it's the innermost waveguide bend, add two points to the coordinates to include the straight waveguides
    if R_i == R_0:
        points_x = np.insert(points_x, 0, -R_0)
        points_y = np.insert(points_y, 0, -inf_eff)
        points_x = np.append(points_x, R_0)
        points_y = np.append(points_y, -inf_eff)

    points = points_x + 1j * points_y

    # add the path to the cell
    cell.add(gdstk.FlexPath(points, w_i, layer=1, datatype=0))
[8]:
# define geometry from the paths in the cell
device_geo = td.PolySlab.from_gds(
    cell,
    gds_layer=1,
    axis=2,
    slab_bounds=(0, h_si),
)

To make sure the geometries are defined correctly, we can plot them first. From the plot, it indeed looks like the correct device.

[9]:
l_straight = 4
pad = lda0

fig, ax = plt.subplots()
for geo in device_geo:
    geo.plot(z=h_si / 2, ax=ax)
ax.set_ylim(-l_straight, R_is[-1] + pad)
ax.set_xlim(-R_is[-1] - pad, R_is[-1] + pad)
plt.show()
../_images/notebooks_SWGBroadbandPolarizer_19_0.png

Lastly, we define Tidy3D Structures from the geometries. A silicon dioxide substrate structure is also defined.

[10]:
device_structure = [td.Structure(geometry=geo, medium=si) for geo in device_geo]
substrate = td.Structure(
    geometry=td.Box.from_bounds(rmin=(-inf_eff, -inf_eff, -inf_eff), rmax=(inf_eff, inf_eff, 0)),
    medium=sio2,
)

TE Input Simulation#

To fully characterize the device performance, we need to simulate two cases: TE input and TM input. First, let’s set up the simulation for the TE input. A ModeSource is defined to inject the TE mode at the input waveguide by setting mode_index=0. Note that since the simulation is quite broadband, we use a broadband source by setting num_freqs = 10, which accounts for the frequency-dependent mode profile.

[11]:
# define a mode source that injects the fundamental te mode
mode_spec = td.ModeSpec(num_modes=2, target_neff=n_si)
mode_source = td.ModeSource(
    center=(-R_0, -l_straight / 2, h_si / 2),
    size=(4 * w_wg, 0, 4 * h_si),
    source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
    direction="+",
    mode_spec=mode_spec,
    mode_index=0,
    num_freqs=10,
)

# define a mode monitor to measure the transmission of the mode
mode_monitor = td.ModeMonitor(
    center=(R_0, -l_straight * 0.9, h_si / 2),
    size=mode_source.size,
    freqs=freqs,
    mode_spec=mode_spec,
    name="mode",
)

lda_field = 1.55  # frequency at which field is monitored

# define a field monitor to visualize mode propagation
field_monitor = td.FieldMonitor(
    center=(0, 0, h_si / 2), size=(td.inf, td.inf, 0), freqs=[td.C_0 / lda_field], name="field"
)
[16:32:58] 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.

[12]:
# simulation domain size
Lx = 2 * (R_is[-1] + pad)
Ly = l_straight + R_is[-1] + pad
Lz = 8 * h_si

run_time = 1e-12  # simulation run time

# define a simulation
sim_te = td.Simulation(
    center=(0, 0.5 * (R_is[-1] + pad - l_straight), 0),
    size=(Lx, Ly, Lz),
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=20, wavelength=lda0),
    structures=device_structure + [substrate],
    sources=[mode_source],
    monitors=[mode_monitor, field_monitor],
    run_time=run_time,
)

Visualize the simulation setup before submitting it to the server.

[13]:
sim_te.plot(z=h_si / 2)
plt.show()
../_images/notebooks_SWGBroadbandPolarizer_28_0.png

Submit the simulation to the server.

[14]:
job = web.Job(simulation=sim_te, task_name="TE_mode")
sim_data_te = job.run(path="data/simulation_data.hdf5")
[16:32:59] Created task 'TE_mode' with task_id                     webapi.py:188
           'fdve-9f8e016c-918d-4504-a6c7-b4fa9229596dv1'.                       
[16:33:01] status = queued                                         webapi.py:361
[16:33:15] status = preprocess                                     webapi.py:355
[16:33:21] Maximum FlexCredit cost: 0.361. Use                     webapi.py:341
           'web.real_cost(task_id)' to get the billed FlexCredit                
           cost after a simulation run.                                         
           starting up solver                                      webapi.py:377
           running solver                                          webapi.py:386
           To cancel the simulation, use 'web.abort(task_id)' or   webapi.py:387
           '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.                                            
[16:36:30] early shutoff detected, exiting.                        webapi.py:404
[16:36:31] status = postprocess                                    webapi.py:419
[16:36:53] status = success                                        webapi.py:426
[16:37:04] loading SimulationData from data/simulation_data.hdf5   webapi.py:590

TM Input Simulation#

Before analyzing the simulation result for the TE input case, let’s run the TM input case first. To have a TM input simulation, we only need to update the mode_index to 1 in the ModeSource. This can be done easily by copying the previous simulation.

[15]:
mode_source = mode_source.copy(
    update={"mode_index": 1}
)  # copy the mode source and update mode_index
sim_tm = sim_te.copy(update={"sources": [mode_source]})  # copy the simulation and update the source

job = web.Job(simulation=sim_tm, task_name="TM_mode")
sim_data_tm = job.run(path="data/simulation_data.hdf5")
[16:37:08] Created task 'TM_mode' with task_id                     webapi.py:188
           'fdve-d44f5e5e-8cfa-4ffb-801f-4754e74e0fcav1'.                       
[16:37:10] status = queued                                         webapi.py:361
[16:37:19] status = preprocess                                     webapi.py:355
[16:37:25] Maximum FlexCredit cost: 0.361. Use                     webapi.py:341
           'web.real_cost(task_id)' to get the billed FlexCredit                
           cost after a simulation run.                                         
           starting up solver                                      webapi.py:377
           running solver                                          webapi.py:386
           To cancel the simulation, use 'web.abort(task_id)' or   webapi.py:387
           '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.                                            
[16:41:09] early shutoff detected, exiting.                        webapi.py:404
           status = postprocess                                    webapi.py:419
[16:41:32] status = success                                        webapi.py:426
[16:41:37] loading SimulationData from data/simulation_data.hdf5   webapi.py:590

Result Analysis#

Now that both simulations are complete, let’s visualize the results. First, we plot the transmission spectra. We see that in the TE input, transmission is close to 0 dB while in the TM input case, transmission is below -20 dB. This is the desirable outcome for a polarizer.

[16]:
amp_te = sim_data_te["mode"].amps.sel(mode_index=0, direction="-")
T_te = np.abs(amp_te) ** 2

amp_tm = sim_data_tm["mode"].amps.sel(mode_index=1, direction="-")
T_tm = np.abs(amp_tm) ** 2

plt.plot(td.C_0 / freqs, 10 * np.log10(T_te), label="TE input")
plt.plot(td.C_0 / freqs, 10 * np.log10(T_tm), label="TM input")
plt.xlabel("Wavelength ($\mu m$)")
plt.ylabel("Transmission (dB)")
plt.ylim(-35, 0)
plt.xlim(lda_min, lda_max)
plt.legend()
plt.show()
../_images/notebooks_SWGBroadbandPolarizer_36_0.png

We go on to plot the field distributions at 1550 nm in the two simulations. The TE mode (left plot) stays guided inside the waveguide bend and thus the transmission is high. The TM mode (right plot), on the other hand, is leaked into the SWG metamaterial region and radiates to free space eventually, leading to a low transmission.

[17]:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 3))
sim_data_te.plot_field(field_monitor_name="field", field_name="Hz", vmin=-0.7, vmax=0.7, ax=ax1)
sim_data_tm.plot_field(field_monitor_name="field", field_name="Ez", vmin=-40, vmax=40, ax=ax2)
ax1.set_title("Hz at 1550 nm (TE input)")
ax2.set_title("Ez at 1550 nm (TM input)")
plt.show()
../_images/notebooks_SWGBroadbandPolarizer_38_0.png

Reference Simulations without SWG#

As a comparison, we can run a baseline reference and simulate the same waveguide bend without the SWG metamaterial cladding. We only need to copy the previous simulations and update the structure such that they don’t contain the SWG.

To further speed up the simulation, we put the TE and TM input cases into a batch and run them on the server simultaneously.

[18]:
# copy previous simulations and update the structures
sim_te_ref = sim_te.copy(update={"structures": [device_structure[0], substrate]})
sim_tm_ref = sim_tm.copy(update={"structures": [device_structure[0], substrate]})

# make a batch and run it
sims = {
    "TE": sim_te_ref,
    "TM": sim_tm_ref,
}
batch = web.Batch(simulations=sims, verbose=True)
batch_results = batch.run(path_dir="data")
[16:41:44] Created task 'TE' with task_id                          webapi.py:188
           'fdve-cf53b5fa-6795-45b2-87ee-a8304cddea46v1'.                       
[16:41:45] Created task 'TM' with task_id                          webapi.py:188
           'fdve-2db2e385-f319-40ab-8f2b-aaae78dd4a3bv1'.                       
[16:41:47] Started working on Batch.                            container.py:475
[16:41:56] Maximum FlexCredit cost: 0.418 for the whole batch.  container.py:479
           Use 'Batch.real_cost()' to get the billed FlexCredit                 
           cost after the Batch has completed.                                  
[16:43:40] Batch complete.                                      container.py:522

Lastly, we plot the same results. Without SWG, we see that both the TE and TM mode transmissions are high. Adding the SWG indeed couples the TM mode out of the waveguide.

[19]:
amp_te_ref = batch_results["TE"]["mode"].amps.sel(mode_index=0, direction="-")
T_te_ref = np.abs(amp_te_ref) ** 2

amp_tm_ref = batch_results["TM"]["mode"].amps.sel(mode_index=1, direction="-")
T_tm_ref = np.abs(amp_tm_ref) ** 2

plt.plot(td.C_0 / freqs, 10 * np.log10(T_te_ref), label="TE input")
plt.plot(td.C_0 / freqs, 10 * np.log10(T_tm_ref), label="TM input")
plt.xlabel("Wavelength ($\mu m$)")
plt.ylabel("Transmission (dB)")
plt.ylim(-35, 0)
plt.xlim(lda_min, lda_max)
plt.legend()
plt.show()
[16:43:46] loading SimulationData from                             webapi.py:590
           data/fdve-cf53b5fa-6795-45b2-87ee-a8304cddea46v1.hdf5                
[16:43:49] loading SimulationData from                             webapi.py:590
           data/fdve-2db2e385-f319-40ab-8f2b-aaae78dd4a3bv1.hdf5                
../_images/notebooks_SWGBroadbandPolarizer_43_8.png
[20]:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 3))
batch_results["TE"].plot_field(
    field_monitor_name="field", field_name="Hz", vmin=-0.7, vmax=0.7, ax=ax1
)
batch_results["TM"].plot_field(
    field_monitor_name="field", field_name="Ez", vmin=-40, vmax=40, ax=ax2
)
ax1.set_title("Hz at 1550 nm (TE input)")
ax2.set_title("Ez at 1550 nm (TM input)")
plt.show()
[16:43:50] loading SimulationData from                             webapi.py:590
           data/fdve-cf53b5fa-6795-45b2-87ee-a8304cddea46v1.hdf5                
[16:43:51] loading SimulationData from                             webapi.py:590
           data/fdve-2db2e385-f319-40ab-8f2b-aaae78dd4a3bv1.hdf5                
../_images/notebooks_SWGBroadbandPolarizer_44_2.png
[ ]: