Broadband polarization beam splitter using anisotropic metamaterial#
Polarization beam splitters (PBS) play a vital role in polarization diversity circuits to mitigate the strong polarization dependence of silicon nanophotonic devices. Among the various PBS structures, asymmetrical directional couplers are predominantly employed due to their superior overall performance. Notwithstanding, engineering an on-chip silicon PBS that offers a compact footprint, low loss, exceptional extinction ratio, and particularly, a broad bandwidth, all at once, remains a formidable challenge. The working bandwidth of directional couplers is inherently limited.
Metamaterials offer an innovative method of manipulating optical responses and regulating light flow. Subwavelength grating (SWG), an anisotropic metamaterial (AM) formed by a one-dimensional array of deeply subwavelength nano-strips, exemplifies these characteristics. SWGs demonstrate strong anisotropy and flattened dispersion, making them extremely useful for polarization manipulation.
In this notebook, we present an ultra-broadband PBS design that makes effective use of AM. The design is based on the research work Xu, H., Dai, D., & Shi, Y. (2019). Ultra-broadband and ultra-compact on-chip silicon polarization beam splitter by using hetero-anisotropic metamaterials. Laser & Photonics Reviews, 13(4), 1800349. DOI: 10.1002/lpor.201800349.
The design comprises three SWG sections symmetrically arranged to form a hetero-anisotropic slab, with two AM[100] sections sandwiching a central AM[010] section. The design also utilizes an adiabatic taper structure at the input/output ports, which serves as a converter between the standard channel waveguides and the SWGs. The design goal is to have the TE mode transmitted to the through port while the TM mode will be coupled to the cross port.

For more examples utilizing SWGs for on-chip applications, please check out the waveguide bragg gratings, broadband polarizer, and exceptional coupling for waveguide crosstalk reduction.
[1]:
from __future__ import annotations
import matplotlib.pyplot as plt
import numpy as np
import tidy3d as td
import tidy3d.web as web
td.config.simulation.use_local_subpixel = True
Simulation Setup#
We will focus on the broadband capabilities of the designed PBS over a 200 nm spectrum, ranging from 1450 nm to 1650 nm. This range encapsulates the S-band, C-band, and L-band frequencies.
[2]:
lda0 = 1.55 # central wavelength
freq0 = td.C_0 / lda0 # central frequency
ldas = np.linspace(1.45, 1.65, 101) # wavelength range
freqs = td.C_0 / ldas # frequency range
fwidth = 0.5 * (np.max(freqs) - np.min(freqs)) # width of the source frequency range
The three materials utilized are silicon, silicon oxide, and SU8. Both silicon and oxide are readily available in the material library. In contrast, SU8 is modelled as a non-dispersive medium.
[3]:
si = td.material_library["cSi"]["Palik_LowLoss"]
sio2 = td.material_library["SiO2"]["Palik_LowLoss"]
n_su8 = 1.58
su8 = td.Medium(permittivity=n_su8**2)
Determine the geometric parameters. The thickness of the silicon layer is recorded to be 250 nm. The input and output single mode waveguides each have a width \(w_{in}\) of 450 nm. The width of the [100] anisotropic metamaterial (AM[100]) \(w_{100}\) is 650 nm. The junction between the single mode strip waveguide and the AM region takes on a taper structure, incorporating $n_{tp}=$10 grating periods. The periods \(p\) for both the AM[100] and AM[010] are set at 250 nm. The duty cycles \(f_{100}\) and \(f_{010}\) are set to 0.7 for both metamaterials, with an option for fine-tuning later.

[4]:
t = 0.25 # silicon layer thickness
p = 0.25 # pitch of the swg
w_in = 0.45 # input waveguide width
w_100 = 0.65 # width of the AM[100]
n_tp = 10 # number of periods in the taper region
f_0 = 0.7 # initial duty cycle for the anisotropic metamaterial
inf_eff = 1e3 # effective infinity
buffer = 0.6 * lda0 # buffer spacing
Mode Analysis#
To gain insight into the operational mechanism of the PBS design and to assist future design optimization, we start with a mode analysis. The effective refractive indices of the SWG AM are procured via the effective medium theory. The subsequent formulas elaborate on this:
The ordinary refractive index, denoted as \(n_o\), is governed by the formula:
\[n_o^2 = f \cdot n_{Si}^2 + (1 - f) \cdot n_{SU8}^2\]Here, \(f\) represents the SWG duty cycle, \(n_{Si}\) stands for the refractive index of silicon, and \(n_{SU8}\) indicates the refractive index of SU8.
The extraordinary refractive index, denoted as \(n_e\), is demonstrated as:
\[\frac{1}{n_e^2} = \frac{f}{n_{Si}^2} + \frac{1 - f}{n_{SU8}^2}\]The refractive index tensors for various configurations, specifically the anisotropic metamaterial [100] and [010], are depicted as such:
\[n_{[100]} = \text{diag}(n_{[100]}^{xx}, n_{[100]}^{yy}, n_{[100]}^{zz}) = \text{diag}(n_e, n_o, n_o)\]Similarly,
\[n_{[010]} = \text{diag}(n_{[010]}^{xx}, n_{[010]}^{yy}, n_{[010]}^{zz}) = \text{diag}(n_o, n_e, n_o)\]
In this context, \(n_o\) and \(n_e\) are the ordinary and extraordinary refractive indices respectively.
By applying the aforementioned refractive index tensors, a mode analysis is conducted on the effective mediums.
[5]:
n_si = si.nk_model(frequency=freq0)[0] # refractive index of silicon
# calculate effective refractive index tensors
n_o = np.sqrt(f_0 * n_si**2 + (1 - f_0) * n_su8**2)
n_e = np.sqrt(1 / (f_0 / n_si**2 + (1 - f_0) / n_su8**2))
# create effective mediums
n_o_medium = td.Medium(permittivity=n_o**2)
n_e_medium = td.Medium(permittivity=n_e**2)
am_100_medium = td.AnisotropicMedium(xx=n_e_medium, yy=n_o_medium, zz=n_o_medium)
am_010_medium = td.AnisotropicMedium(xx=n_o_medium, yy=n_e_medium, zz=n_o_medium)
# create waveguide structures
n_010_example = 5 # use 5 periods for the AM[010] region width as an example
am_010_waveguide = td.Structure(
geometry=td.Box.from_bounds(
rmin=(-inf_eff, -n_010_example * p / 2, 0), rmax=(inf_eff, n_010_example * p / 2, t)
),
medium=am_010_medium,
)
am_100_waveguide_left = td.Structure(
geometry=td.Box.from_bounds(
rmin=(-inf_eff, -n_010_example * p / 2 - w_100, 0),
rmax=(inf_eff, -n_010_example * p / 2, t),
),
medium=am_100_medium,
)
am_100_waveguide_right = td.Structure(
geometry=td.Box.from_bounds(
rmin=(-inf_eff, n_010_example * p / 2, 0), rmax=(inf_eff, n_010_example * p / 2 + w_100, t)
),
medium=am_100_medium,
)
# create box structure
box = td.Structure(
geometry=td.Box.from_bounds(rmin=(-inf_eff, -inf_eff, -inf_eff), rmax=(inf_eff, inf_eff, 0)),
medium=sio2,
)
# create a simulation for the mode analysis
mode_analysis_sim = td.Simulation(
center=(0, 0, t / 2),
size=(buffer, n_010_example * p + 2 * w_100 + 2 * buffer, t + 2 * buffer),
grid_spec=td.GridSpec.auto(min_steps_per_wvl=25, wavelength=lda0),
structures=[box, am_010_waveguide, am_100_waveguide_left, am_100_waveguide_right],
run_time=1e-12,
medium=su8,
)
# create a mode simulation
mode_simulation = td.ModeSimulation.from_simulation(
simulation=mode_analysis_sim,
plane=td.Box(
center=(0, 0, t / 2), size=(0, n_010_example * p + 2 * w_100 + 2 * buffer, t + 2 * buffer)
),
mode_spec=td.ModeSpec(num_modes=6, target_neff=n_si),
freqs=[freq0],
)
# plot the mode solving cross section
mode_simulation.plot()
plt.show()
To ensure accuracy in our process, we execute the mode analysis locally with local subpixel averaging enabled and subsequently illustrate the mode properties. Utilizing the TE (Ey) fraction as a reference, we can discern that the initial two modes align with the first TE modes. Similarly, the third and fifth modes correspond to the first TM modes.
[6]:
sim_data = mode_simulation.run_local()
mode_data = sim_data.modes
mode_data.to_dataframe()