Designing a polarization splitter/rotator on thin-film lithium niobate#
Note: the cost of running the total notebook is almost 27.5 FlexCredits. The cost of running the FDTD validation is just over 26 FlexCredits.
Here we design a polarization splitter/rotator (PSR) on the lithium niobate on insulator (LNOI) platform based on Xuanhao Wang, An Pan, Tingan Li, Cheng Zeng, and Jinsong Xia, "Efficient polarization splitter-rotator on thin-film lithium niobate," Opt. Express 29, 38044-38052 (2021) DOI: 10.1364/OE.443798 and follow the same design process. This will operate in a bandwidth of 1.5 μm to 1.6 μm.
PSRs are useful devices that control optical input. By splitting and rotating its input, a PSR can isolate different modes and transform them according to our needs. This PSR consists of an adiabatic taper to convert the TM0 mode to TE1 mode and an asymmetric Y-junction that separates TE0 and TE1 modes. In this way, our PSR is able to split TE0 and TM0 modes, and rotate the TM0 mode into a TE0 mode. Devices like these allow us to modulate inputs and interface with other optical devices in a controlled way.
In this design process, we will optimize the taper and Y-junction separately using a combination of Tidy3D’s mode solver and eigenmode expansion (EME) solver. Finally we validate our optimized designs with FDTD to further ensure the performance.
[1]:
import numpy as np
from matplotlib import pyplot as plt
import tidy3d as td
import tidy3d.web as web
from tidy3d import material_library
from tidy3d.plugins.mode import ModeSolver
from tidy3d.plugins.mode.web import run_batch
Simulation Setup#
First we define some values in our bandwidth and LNOI setup. Our LNOI setup consists of a lithium niobate (LN) structure with a sidewall angle of 20 degree on top of a thin film that is also LN. All of our LNOI structures will be defined in this way.
[2]:
wvl0 = 1.55 # central wavelength for device
freq0 = td.C_0 / wvl0 # central frequency
wvls = np.linspace(1.5, 1.6, 41) # bandwidth we'll consider for device
freqs = td.C_0 / wvls # corresponding frequency range
sidewall_angle = np.deg2rad(20)
etch_depth = 0.26
film_thickness = 0.5
LN is an anisotropic medium, with different refractive indices in the ordinary and extraordinary directions. We incorporate this information using Tidy3D’s AnisotropicMedium object. Since the propagation direction of the PSR will be defined along the x-direction, this will have the ordinary refractive index.
NOTE: Tidy3D offers lithium niobate in its material library. However, the values there differ slightly from the values used in the paper on which this notebook is based, so we will use the latter.
[3]:
n_e = 2.169 # extraordinary refractive index of LN
n_o = 2.234 # ordinary refractive index of LN
n_SiO2 = 1.46 # refractive index of SiO2
# define LN medium
LN_e = td.Medium(permittivity=n_e**2)
LN_o = td.Medium(permittivity=n_o**2)
LN = td.AnisotropicMedium(xx=LN_o, yy=LN_e, zz=LN_o)
# define SiO2 medium
SiO2 = td.Medium(permittivity=n_SiO2**2)
Next we define a function that takes in a set of points (like a polyslab) and produces the corresponding LNOI structure:
[4]:
def make_ridge_waveguide(pts):
"""
Define the sidewall angled waveguide that sits on top of a Lithium Niobate film:
"""
LN_slab_thickness = film_thickness - etch_depth
waveguide_geo = td.PolySlab(
sidewall_angle=sidewall_angle,
vertices=pts,
reference_plane="bottom",
slab_bounds=[-0.005, etch_depth],
)
waveguide = td.Structure(geometry=waveguide_geo, medium=LN)
"""
Define the film that the waveguide sits on:
"""
LN_slab_geo = td.Box(
center=(0, 0, -LN_slab_thickness / 2), # ensure overlap
size=(td.inf, td.inf, LN_slab_thickness),
)
LN_slab = td.Structure(geometry=LN_slab_geo, medium=LN)
return [LN_slab, waveguide]
Mode Analysis#
To design our adiabatic taper, we want to determine the LNOI cross section at which a TM0 mode hybridizes with the TE1 mode. To do this, we will use Tidy3D’s mode solver for cross sections at increasing LNOI taper widths. We will do this by creating a single taper structure on LNOI and placing mode solver planes along the taper.
[5]:
# We will define the taper by a (somewhat arbitrary) difference in widths over a length.
width_test_w0 = 2 # starting width of taper
width_test_w1 = 3 # ending width of taper
width_test_length = 10 # length of taper
"""
Define points that take the shape of a taper. Note that we extend the inital and final
widths past the boundaries to avoid the sidewall angles from interfering with the ends
of the taper.
"""
width_test_pts = [
(-1, -width_test_w0 / 2),
(0, -width_test_w0 / 2),
(width_test_length, -width_test_w1 / 2),
(width_test_length + 1, -width_test_w1 / 2),
(width_test_length + 1, width_test_w1 / 2),
(width_test_length, width_test_w1 / 2),
(0, width_test_w0 / 2),
(-1, width_test_w0 / 2),
]
width_test_structures = make_ridge_waveguide(width_test_pts)
# create simulation environment of the taper
width_test_sim = td.Simulation(
center=(width_test_length / 2, 0, 0),
size=(width_test_length + 0.2, 2 * width_test_w0, 3 * film_thickness),
structures=width_test_structures,
grid_spec=td.GridSpec.auto(min_steps_per_wvl=30, wavelength=wvl0),
medium=SiO2,
run_time=1e-15,
)
# plot the simulation to check the taper and LNOI geometry
_, ax = plt.subplots(1, 2, figsize=(12, 4))
width_test_sim.plot(z=etch_depth / 2, ax=ax[0])
width_test_sim.plot(x=0, ax=ax[1])
plt.show()
To find the relationship between modes and width, we will place mode solver planes along the taper (along the x direction). These planes will solve for modes at each cross-section, which will be increasing in thickness as we go along the taper. This will tell us at which thickness a mode conversion will occur.
NOTE: Since this cell runs the remote mode solver in series, it may take a couple minutes to run.
[6]:
plane_xs = np.linspace(0, width_test_length, 31) # x values of the mode solver planes
# modespec for the mode solver planes
width_test_mode_spec = td.ModeSpec(num_modes=3, target_neff=n_o)
mode_solvers = []
for i, m in enumerate(plane_xs):
# define plane at x coordinate
plane = td.Box(center=(m, 0, 0), size=(0, td.inf, td.inf))
# create mode solver from plane
mode_solver = ModeSolver(
simulation=width_test_sim, plane=plane, mode_spec=width_test_mode_spec, freqs=[freq0]
)
# add modesolver for this mode solver plane
mode_solvers.append(mode_solver)
[7]:
# run mode solvers in parallel
mode_batch = run_batch(mode_solvers=mode_solvers, task_name="mode_conversion")
08:49:13 PDT Running a batch of 31 mode solvers.
08:51:26 PDT A batch of `ModeSolver` tasks completed successfully!
[8]:
# collect relevant data from batch solve
modes = np.reshape([mode_data.n_eff.data[0] for mode_data in mode_batch], (len(plane_xs), 3))
te_fracs = np.reshape([mode_data.pol_fraction.te.data[0] for mode_data in mode_batch], (len(plane_xs), 3))
[9]:
from matplotlib.patches import Ellipse
import matplotlib.colors as mcolors
ridge_waveguide_widths = (
width_test_w1 - width_test_w0
) / width_test_length * plane_xs + width_test_w0
_, ax = plt.subplots(1, 1, figsize=(7, 4))
conversion = [2.52, 2.7]
ax.plot(ridge_waveguide_widths, modes[:, 0], label="TE0")
ax.scatter(ridge_waveguide_widths, modes[:, 0], c=te_fracs[:, 0], cmap="cool")
ax.plot(ridge_waveguide_widths, modes[:, 1], label="TM0")
ax.scatter(ridge_waveguide_widths, modes[:, 1], c=te_fracs[:, 1], cmap="cool")
ax.plot(ridge_waveguide_widths, modes[:, 2], label="TE1")
ax.scatter(ridge_waveguide_widths, modes[:, 2], c=te_fracs[:, 2], cmap="cool")
circle = Ellipse((2.61, 1.887), width=0.1, height=0.013, ec="blue", fill=False, ls="--")
ax.add_patch(circle)
ax.axvline(x=conversion[0], color="gray", ls="--")
ax.axvline(x=conversion[1], color="gray", ls="--")
ax.axvspan(conversion[0], conversion[1], color="gray", alpha=0.2)
ax.set_xlabel("Ridge waveguide width (μm)")
ax.set_ylabel("Effective refractive indices")
norm = mcolors.Normalize(vmin=0, vmax=1)
plt.colorbar(plt.cm.ScalarMappable(cmap="cool", norm=norm), ax=ax, label="TE fraction")
plt.show()
To show the TM0 to TE1 hybridization, we will plot the effective indices vs taper width, as well as the TE fraction as the color of the curves. We can see that around a taper width of 2.55 μm is when the modes cross.
This mode conversion behavior differs from what is shown in Wang et. al. This is likely due to sub-pixel averaging. When the mode solving is done locally (and thus without subpixel averaging) the result is very similar to what is given in Wang et. al.
Adiabatic Taper Design#
Efficient mode conversion can only happen if the taper changes gradually. Therefore, we will set our conversion taper to be made up of 3 tapers. The first taper will start with a width of 1 μm and end with a width of 2.52 μm. This will be for the single-mode input. The next taper is from 2.52 μm to 2.7 μm for mode conversion. The final taper from 2.7 μm to 3 μm is the connection to the asymmetric Y-junction.
Now that we’ve established the widths of the tapers, the next question is how long the tapers should be. According to the paper this notebook is based on, the choice of lengths of the first and third tapers do not have a significant effect (as the length of the second taper can be tuned for these choices), and thus will be chosen as 100 μm for each.
[10]:
l1 = 100 # length of first taper
l3 = 100 # length of third taper
w0 = 1 # starting width of first taper
w1 = conversion[0] # starting width of second taper
w2 = conversion[1] # starting width of third taper
w3 = 3 # ending width of third taper
To determine the length of the second taper (\(L_2\)) that is optimal for mode conversion efficiency, we will run an EME parameter sweep for different simulations with different \(L_2\) values. To do this efficiently, we will use the EMELengthSweep feature in Tidy3D’s EME simulations, requiring us to only specify the scale factors we’ll assign to the EME cells we’d like to vary.
EME is an excellent choice for investigations of this kind. This type of simulation is mainly about solving for modes at different cross sections, and seeing how these modes change from one cross section to another (in this case, along the taper). Since these changes are adiabatic, it is necessary to check for very long lengths in the propagation direction. For this reason, the device can get very large, making FDTD solving very slow and expensive. On the other hand, mode solving in the device’s cross section is impacted a lot less by extending the length normal to the cross-section. Since this is mainly what EME relies on, EME will be a much more efficient option.
[11]:
l2_sweep = 1 # value to construct EME sweep
l2_test_pts = [
(-10, -w0 / 2), # extend taper into -x boundary
(0, -w0 / 2),
(l1, -w1 / 2),
(l1 + l2_sweep, -w2 / 2),
(l1 + l2_sweep + l3, -w3 / 2),
(l1 + l2_sweep + l3 + 10, -w3 / 2), # extend taper into +x boundary
(l1 + l2_sweep + l3 + 10, w3 / 2),
(l1 + l2_sweep + l3, w3 / 2),
(l1 + l2_sweep, w2 / 2),
(l1, w1 / 2),
(0, w0 / 2),
(-10, w0 / 2),
]
l2_test_structures = make_ridge_waveguide(l2_test_pts)
"""
Here we define the cells we'll use in our EME simulation. We will segment our cells into the three
tapers. Since we care most about our L2 taper, we will divide this into more cells than the other tapers.
"""
num_modes = 4
l1_grid_spec = td.EMEUniformGrid(num_cells=8, mode_spec=td.EMEModeSpec(num_modes=num_modes))
l2_grid_spec = td.EMEUniformGrid(num_cells=15, mode_spec=td.EMEModeSpec(num_modes=num_modes))
l2_total_grid = td.EMECompositeGrid(
subgrids=[l1_grid_spec, l2_grid_spec, l1_grid_spec], subgrid_boundaries=[l1, l1 + l2_sweep]
)
# define regular mesh grid. This is required for mode solving
auto_x = td.AutoGrid(min_steps_per_wvl=10)
auto_yz = td.AutoGrid(min_steps_per_wvl=30)
l2_scales = np.linspace(1 / l2_sweep, 500 / l2_sweep, 51)
l2_scale_factors = np.ones((len(l2_scales), 2 * l1_grid_spec.num_cells + l2_grid_spec.num_cells))
l2_scale_factors[:, l1_grid_spec.num_cells : l1_grid_spec.num_cells + l2_grid_spec.num_cells] = (
l2_scales[:, None]
)
# define EME simulation
l2_sweep_sim = td.EMESimulation(
center=((l1 + l2_sweep + l3) / 2, 0, 0),
size=(l1 + l2_sweep + l3, w3 + 3 * w0, 5 * film_thickness),
structures=l2_test_structures,
grid_spec=td.GridSpec(grid_x=auto_x, grid_y=auto_yz, grid_z=auto_yz, wavelength=wvl0),
medium=SiO2,
eme_grid_spec=l2_total_grid,
freqs=[freq0],
axis=0, # EME will be solved along the x axis
sweep_spec=td.EMELengthSweep(scale_factors=l2_scale_factors),
symmetry=(0, 1, 0), # the adiabatic taper has symmetry in y
)
ERROR: Expected 1 dimensions for ArrayLike, got 2.
We will plot this simulation just to make sure everything looks correct:
[12]:
_, ax = plt.subplots(1, 2, figsize=(12, 4))
l2_sweep_sim.plot(z=0.005, ax=ax[0])
l2_sweep_sim.plot(x=l2_sweep, ax=ax[1])
ax[0].set_aspect(6)
plt.show()
Now we will run the EME simulation (sweep):
[13]:
l2_EME_data = web.run(l2_sweep_sim, task_name="L2 EME sweep")
08:51:27 PDT Created task 'L2 EME sweep' with task_id 'eme-5ce796bc-d2fe-4828-a5ea-46a76d7b261e' and task_type 'EME'.
Tidy3D's EME solver is currently in the beta stage. Cost of EME simulations is subject to change in the future.
08:51:29 PDT 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.
08:51:31 PDT Maximum FlexCredit cost: 0.231. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
starting up solver
running solver
08:54:14 PDT status = success
08:54:28 PDT loading simulation from simulation_data.hdf5
ERROR: Expected 1 dimensions for ArrayLike, got 2.
Plotting the mode conversion efficiency vs \(L_2\), we can determine the optimal \(L_2\) value:
[14]:
l2_TM0_TE1 = l2_EME_data.smatrix.S21.sel(f=freq0, mode_index_in=0, mode_index_out=0).abs.data ** 2
l2_TM0_TM0 = l2_EME_data.smatrix.S21.sel(f=freq0, mode_index_in=0, mode_index_out=1).abs.data ** 2
arg_min = np.argmin(l2_TM0_TM0[:25])
l2 = l2_scales[arg_min] * l2_sweep
[15]:
plt.plot(l2_scales * l2_sweep, l2_TM0_TE1, c="purple", linewidth=2, label="TM0 → TE1")
plt.plot(l2_scales * l2_sweep, l2_TM0_TM0, c="green", linewidth=2, label="TM0 → TM0")
plt.xlim(0, 500)
plt.ylim(0, 1)
plt.xlabel("L2 length (μm)")
plt.ylabel("Mode conversion efficiency")
plt.title("Mode conversion efficiency vs L2")
plt.axvline(x=l2, color="gray", ls="--")
plt.legend()
plt.show()
print(f"Optimal length: {l2}")
Optimal length: 170.66
To verify our results, we will run FDTD validations of \(L_2=170.66\) and \(L_2=250\):
FDTD Validation of EME L2 Sweep#
The FDTD method, while much slower and more expensive than the EME method, is more robust. It is always a good idea to double-check our simulation results with multiple methods. Thus we will define a function that creates an FDTD simulation with a given \(L_2\) value. This way we can easily create multiple FDTD simulations with multiple \(L_2\) values.
[16]:
def make_l2_FDTD_sim(l2_sweep):
# create taper structure
l2_test_pts = [
(-10, -w0 / 2),
(0, -w0 / 2),
(l1, -w1 / 2),
(l1 + l2_sweep, -w2 / 2),
(l1 + l2_sweep + l3, -w3 / 2),
(l1 + l2_sweep + l3 + 10, -w3 / 2),
(l1 + l2_sweep + l3 + 10, w3 / 2),
(l1 + l2_sweep + l3, w3 / 2),
(l1 + l2_sweep, w2 / 2),
(l1, w1 / 2),
(0, w0 / 2),
(-10, w0 / 2),
]
l2_test_structures = make_ridge_waveguide(l2_test_pts)
# add mode monitor
l2_mode_monitor = td.ModeMonitor(
name="l2 mode",
size=(0, td.inf, td.inf),
center=(l1 + l2_sweep + l3 * 0.95, 0, 0),
freqs=[freq0],
mode_spec=td.ModeSpec(num_modes=2, target_neff=n_o),
)
# specify FDTD grid
auto_x = td.AutoGrid(min_steps_per_wvl=10)
auto_yz = td.AutoGrid(min_steps_per_wvl=30)
# create simulation with these specs
l2_sweep_sim = td.Simulation(
center=((l1 + l2_sweep + l3) / 2, 0, 0),
size=(l1 + l2_sweep + l3, w3 + 2 * w0, 5 * film_thickness),
sources=[],
structures=l2_test_structures,
grid_spec=td.GridSpec(grid_x=auto_x, grid_y=auto_yz, grid_z=auto_yz, wavelength=wvl0),
medium=SiO2,
monitors=[l2_mode_monitor],
symmetry=(0, 1, 0),
run_time=5e-12,
)
"""
To inject the right mode into the FDTD simulation, we will need to solve for this correct
mode in FDTD. We will convert this to a source and add it to the simulation.
"""
fdtd_mode_solver = ModeSolver(
simulation=l2_sweep_sim,
plane=td.Box(center=(1, 0, 0), size=(0, td.inf, td.inf)),
mode_spec=td.ModeSpec(num_modes=3, target_neff=n_o),
freqs=[freq0],
)
fdtd_source = fdtd_mode_solver.to_source(
mode_index=0, direction="+", source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10)
)
return l2_sweep_sim.updated_copy(sources=[fdtd_source])
[17]:
# create FDTD simulations where L2 = l2, 250 μm
l2_FDTD_opt_sim = make_l2_FDTD_sim(l2)
l2_FDTD_250_sim = make_l2_FDTD_sim(250)
[18]:
# validate our FDTD-creating function
_, ax = plt.subplots(1, 2, figsize=(12, 4))
l2_FDTD_opt_sim.plot(z=0.005, ax=ax[0])
l2_FDTD_opt_sim.plot(x=l1 + 125 + l3, ax=ax[1])
ax[0].set_aspect(6)
plt.show()
[19]:
l2_FDTD_opt_data = web.run(simulation=l2_FDTD_opt_sim, task_name="L2 FDTD opt")
l2_FDTD_250_data = web.run(simulation=l2_FDTD_250_sim, task_name="L2 FDTD 250")
Created task 'L2 FDTD opt' with task_id 'fdve-d1931ecc-7a98-4561-a467-4748cba61b61' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-d1931ecc-7a9 8-4561-a467-4748cba61b61'.
08:54:30 PDT 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.
08:54:51 PDT status = preprocess
08:55:00 PDT Maximum FlexCredit cost: 2.081. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
starting up solver
running solver
08:56:00 PDT early shutoff detected at 68%, exiting.
status = postprocess
08:56:01 PDT status = success
View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-d1931ecc-7a9 8-4561-a467-4748cba61b61'.
08:56:02 PDT loading simulation from simulation_data.hdf5
Created task 'L2 FDTD 250' with task_id 'fdve-6367aed1-9fe9-4315-b5c9-144862fa1f7c' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-6367aed1-9fe 9-4315-b5c9-144862fa1f7c'.
08:56:03 PDT 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.
08:56:18 PDT status = preprocess
08:56:21 PDT Maximum FlexCredit cost: 2.511. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
starting up solver
08:56:22 PDT running solver
08:58:40 PDT early shutoff detected at 80%, exiting.
status = postprocess
status = success
View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-6367aed1-9fe 9-4315-b5c9-144862fa1f7c'.
08:58:41 PDT loading simulation from simulation_data.hdf5
Now that we have run our equivalent FDTD simulations, we will check these values against our EME values.
[20]:
TM0_TM0_opt_FDTD = (
l2_FDTD_opt_data["l2 mode"].amps.sel(direction="+", mode_index=1).abs.data[0] ** 2
)
TM0_TE1_opt_FDTD = (
l2_FDTD_opt_data["l2 mode"].amps.sel(direction="+", mode_index=0).abs.data[0] ** 2
)
TM0_TM0_opt_EME = l2_TM0_TM0[arg_min]
TM0_TE1_opt_EME = l2_TM0_TE1[arg_min]
TM0_TM0_250_FDTD = (
l2_FDTD_250_data["l2 mode"].amps.sel(direction="+", mode_index=1).abs.data[0] ** 2
)
TM0_TE1_250_FDTD = (
l2_FDTD_250_data["l2 mode"].amps.sel(direction="+", mode_index=0).abs.data[0] ** 2
)
TM0_TM0_250_EME = l2_TM0_TM0[25]
TM0_TE1_250_EME = l2_TM0_TE1[25]
lengths_l2 = np.array([l2, 250])
TM0_FDTD = np.array([TM0_TM0_opt_FDTD, TM0_TM0_250_FDTD])
TM0_EME = np.array([TM0_TM0_opt_FDTD, TM0_TM0_250_EME])
TE1_FDTD = np.array([TM0_TE1_opt_FDTD, TM0_TE1_250_FDTD])
TE1_EME = np.array([TM0_TE1_opt_EME, TM0_TE1_250_EME])
_, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].scatter(lengths_l2, TM0_FDTD, marker="s", c="blue", label="FDTD")
ax[0].scatter(lengths_l2, TM0_EME, marker="^", c="red", label="EME")
ax[0].legend()
ax[0].set_ylim(0, 1)
ax[0].set_xticks(lengths_l2, lengths_l2)
ax[0].margins(0.25)
ax[0].set_xlabel("L2 (μm)")
ax[0].set_ylabel("Conversion Efficiency")
ax[0].set_title("TM0 to TM0")
ax[1].scatter(lengths_l2, TE1_FDTD, marker="s", c="blue", label="FDTD")
ax[1].scatter(lengths_l2, TE1_EME, marker="^", c="red", label="EME")
ax[1].legend()
ax[1].set_ylim(0, 1)
ax[1].set_xticks(lengths_l2, lengths_l2)
ax[1].margins(0.25)
ax[1].set_xlabel("L2 (μm)")
ax[1].set_title("TM0 to TE1")
plt.show()
These results coincide nicely, verifying the accuracy of our EME sweep. Now let’s check the field profile from our EME simulation where \(L_2\) = 170.66 μm, as this is where maximum efficiency occurs. We will do this by adding a field monitor to this simulation.
[21]:
# create taper structure
l2_ideal_pts = [
(-10, -w0 / 2),
(0, -w0 / 2),
(l1, -w1 / 2),
(l1 + l2, -w2 / 2),
(l1 + l2 + l3, -w3 / 2),
(l1 + l2 + l3 + 10, -w3 / 2),
(l1 + l2 + l3 + 10, w3 / 2),
(l1 + l2 + l3, w3 / 2),
(l1 + l2, w2 / 2),
(l1, w1 / 2),
(0, w0 / 2),
(-10, w0 / 2),
]
l2_ideal_structures = make_ridge_waveguide(l2_ideal_pts)
# create field monitor
l2_field_monitor = td.EMEFieldMonitor(
name="l2 field",
size=(td.inf, td.inf, 0),
center=(l1 + l2 + l3 * 0.95, 0, 0),
freqs=[freq0],
fields=["Ex", "Ey", "Ez"],
num_modes=4,
)
# add field monitor to EME simulation
optimal_l2_sim = l2_sweep_sim.updated_copy(
structures=l2_ideal_structures, monitors=[l2_field_monitor], symmetry=(0, 0, 0), sweep_spec=None
)
[22]:
ideal_l2_data = web.run(simulation=optimal_l2_sim, task_name="Ideal L2 field validation")
Created task 'Ideal L2 field validation' with task_id 'eme-2b665774-73df-41a3-8e92-970bce64dc36' and task_type 'EME'.
Tidy3D's EME solver is currently in the beta stage. Cost of EME simulations is subject to change in the future.
08:58:43 PDT 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.
08:59:01 PDT Maximum FlexCredit cost: 0.420. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
starting up solver
running solver
09:01:11 PDT status = success
09:17:10 PDT loading simulation from simulation_data.hdf5
[23]:
_, ax = plt.subplots(1, 2, figsize=(12, 4))
ideal_l2_data.plot_field("l2 field", "E", "abs", f=freq0, eme_port_index=0, mode_index=0, ax=ax[0])
ideal_l2_data.plot_field("l2 field", "E", "abs", f=freq0, eme_port_index=1, mode_index=1, ax=ax[1])
ax[0].set_aspect(10)
ax[1].set_aspect(10)
ax[0].set_title("TE0 to TE0")
ax[1].set_title("TM0 to TE1")
plt.show()
Design of Asymmetric Y-Junction#
We now move forward with the design of the Y-junction. The purpose of this junction is to separate TE0 and TE1 modes - in the wider arm of the junction, a TE0 mode will stay a TE0 mode since the effective refractive indices of the modes will be closest, and, in the thinner arm, a TE1 mode will evolve to a TE0 mode for the same reason. The widths we will choose in order to enable this will be 1.7 μm and 1.3 μm, respectively. We will choose the final distance between the arms to be 1 μm.
This mode conversion will only happen adiabatically, so the angle splitting the two arms must be very small. The central question of our design, therefore, is how long the Y-junction will be. To find the optimal length (which we’ll denote as \(L_4\)), we will use Tidy3D’s EME sweep feature to do screen various \(L_4\) lengths. Tidy3D’s EME sweep feature is a fast and inexpensive way to run sweeps where the entire simulation is scaled, which is ideal for testing this geometry.
We will measure the TE0 and TM0 conversion efficiencies in both arms. The upper, wider arm will be labeled ‘port 1’ and the lower, thinner arm will be labeled ‘port 2.’
[24]:
w4 = 1.7 # width of upper arm
w5 = 1.3 # width of lower arm
split = 1 # ending distance between two arms
# the starting length for the Y-junction. This will change as we change the scale
l4_sweep = 100
# create Y-junction geometry and structure
l4_test_pts = [
(-2 * l3, -w2 / 2), # extend beyond boundary
(-l3, -w2 / 2),
(0, -w3 / 2),
(l4_sweep * 2, 2 * (-split / 2) - w5), # extend beyond boundary
(l4_sweep * 2, 2 * (-split / 2)), # extend beyond boundary
(0, -w3 / 2 + w5),
(l4_sweep * 2, 2 * (split / 2)), # extend beyond boundary
(l4_sweep * 2, 2 * (split / 2) + w4), # extend beyond boundary
(0, w3 / 2),
(-l3, w2 / 2),
(-2 * l3, w2 / 2) # extend beyond boundary
]
l4_test_structures = make_ridge_waveguide(l4_test_pts)
num_modes = 5
l4_port1 = td.ModeSolverMonitor(
name="l4 port 1",
size=(0, w4 + 1.5 * split, td.inf),
center=(l4_sweep - 0.1, (-w3 / 2 + w5) + split / 2 + w4 / 2, 0),
freqs=[freq0],
mode_spec=td.ModeSpec(
num_modes=num_modes,
target_neff=n_o,
angle_theta=np.arctan((split / 2 - (-w3 / 2 + w5)) / l4_sweep),
),
)
l4_port2 = td.ModeSolverMonitor(
name="l4 port 2",
size=(0, w5 + 1.5 * split, td.inf),
center=(l4_sweep - 0.1, -split / 2 - w5 / 2, 0),
freqs=[freq0],
mode_spec=td.ModeSpec(
num_modes=num_modes,
target_neff=n_o,
angle_theta=np.arctan((-split / 2 - w5 / 2 + w3 / 4) / l4_sweep),
),
)
# again owing to the adiabatic nature of this problem, we need only use a unifrom grid throughout the simulation
l3_grid_spec = td.EMEUniformGrid(num_cells=16, mode_spec=td.EMEModeSpec(num_modes=num_modes))
l4_grid_spec = td.EMEUniformGrid(num_cells=41, mode_spec=td.EMEModeSpec(num_modes=num_modes))
l4_total_grid = td.EMECompositeGrid(subgrids=[l3_grid_spec, l4_grid_spec], subgrid_boundaries=[0])
auto_x = td.AutoGrid(min_steps_per_wvl=31)
auto_yz = td.AutoGrid(min_steps_per_wvl=41)
refine_box = td.MeshOverrideStructure(
geometry=td.Box(center=(0, -w3 / 2 + w5, 0), size=(l4_sweep + l3, split / 2, 2.5 * etch_depth)),
dl=[0.1, 0.01, 0.01],
)
l4_scales = np.linspace(1 / l4_sweep, 200 / l4_sweep, 61)
l4_scale_factors = np.ones((len(l4_scales), l3_grid_spec.num_cells + l4_grid_spec.num_cells))
l4_scale_factors[:, l3_grid_spec.num_cells : l3_grid_spec.num_cells + l4_grid_spec.num_cells] = (
l4_scales[:, None]
)
l4_EME = td.EMESimulation(
center=((l4_sweep - l3) / 2, 0, 0),
size=(l4_sweep + l3, 3 * split + w4 + w5, 5 * film_thickness),
structures=l4_test_structures,
grid_spec=td.GridSpec(grid_x=auto_x, grid_y=auto_yz, grid_z=auto_yz, wavelength=wvl0),
medium=SiO2,
monitors=[l4_port1, l4_port2],
eme_grid_spec=l4_total_grid,
freqs=[freq0],
sweep_spec=td.EMELengthSweep(scale_factors=list(l4_scale_factors)),
axis=0,
)
09:17:12 PDT ERROR: Expected 1 dimensions for ArrayLike, got 2.
[25]:
x_plot = 0
_, ax = plt.subplots(1, 2, figsize=(12, 4))
l4_EME.plot(z=0.005, ax=ax[0])
l4_EME.plot(x=x_plot, ax=ax[1])
l4_EME.plot_grid(x=x_plot, ax=ax[1], lw=0.4, colors="r")
ax[0].set_aspect(6)
ax[1].set_xlim(-0.7, 0.3)
ax[1].set_ylim(-0.25, 0.32)
plt.show()
[26]:
l4_EME_data = web.run(l4_EME, task_name="L4 EME sweep")
09:17:14 PDT Created task 'L4 EME sweep' with task_id 'eme-2b82840a-5bf9-4859-a4d4-f296cedb2ca6' and task_type 'EME'.
Tidy3D's EME solver is currently in the beta stage. Cost of EME simulations is subject to change in the future.
09:17:18 PDT 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.
09:17:56 PDT Maximum FlexCredit cost: 0.708. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
starting up solver
09:17:57 PDT running solver
09:24:54 PDT status = success
09:25:44 PDT loading simulation from simulation_data.hdf5
ERROR: Expected 1 dimensions for ArrayLike, got 2.
We want to examine the scattering matrix in the bases of modes in each junction arm. The EMESimulationData.smatrix_in_basis data is thus what we want to examine.
[27]:
TE0_TE0_port1 = (
l4_EME_data.smatrix_in_basis(modes2=l4_EME_data["l4 port 1"])
.S21.sel(f=freq0, mode_index_in=0, mode_index_out=0)
.abs.data
** 2
)
TE0_TE0_port2 = (
l4_EME_data.smatrix_in_basis(modes2=l4_EME_data["l4 port 2"])
.S21.sel(f=freq0, mode_index_in=0, mode_index_out=0)
.abs.data
** 2
)
TE0_TM0_port1 = (
l4_EME_data.smatrix_in_basis(modes2=l4_EME_data["l4 port 1"])
.S21.sel(f=freq0, mode_index_in=0, mode_index_out=1)
.abs.data
** 2
)
TE0_TM0_port2 = (
l4_EME_data.smatrix_in_basis(modes2=l4_EME_data["l4 port 2"])
.S21.sel(f=freq0, mode_index_in=0, mode_index_out=1)
.abs.data
** 2
)
TE1_TE0_port1 = (
l4_EME_data.smatrix_in_basis(modes2=l4_EME_data["l4 port 1"])
.S21.sel(f=freq0, mode_index_in=1, mode_index_out=0)
.abs.data
** 2
)
TE1_TE0_port2 = (
l4_EME_data.smatrix_in_basis(modes2=l4_EME_data["l4 port 2"])
.S21.sel(f=freq0, mode_index_in=1, mode_index_out=0)
.abs.data
** 2
)
TE1_TM0_port1 = (
l4_EME_data.smatrix_in_basis(modes2=l4_EME_data["l4 port 1"])
.S21.sel(f=freq0, mode_index_in=1, mode_index_out=1)
.abs.data
** 2
)
TE1_TM0_port2 = (
l4_EME_data.smatrix_in_basis(modes2=l4_EME_data["l4 port 2"])
.S21.sel(f=freq0, mode_index_in=1, mode_index_out=1)
.abs.data
** 2
)
[28]:
_, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].plot(l4_scales * l4_sweep, TE0_TE0_port1, c="red", label="TE0 in port1")
ax[0].plot(l4_scales * l4_sweep, TE0_TE0_port2, c="blue", label="TE0 in port2")
ax[0].plot(l4_scales * l4_sweep, TE0_TM0_port1, c="purple", label="TM0 in port1")
ax[0].plot(l4_scales * l4_sweep, TE0_TM0_port2, c="green", label="TM0 in port2")
ax[0].set_xlim((0, 200))
ax[0].set_ylim((0, 1))
ax[0].legend()
ax[0].set_xlabel("L4 (μm)")
ax[0].set_ylabel("Mode conversion efficiency")
ax[0].set_title("TE0 Input")
ax[1].plot(l4_scales * l4_sweep, TE1_TE0_port1, c="red", label="TE0 in port1")
ax[1].plot(l4_scales * l4_sweep, TE1_TE0_port2, c="blue", label="TE0 in port2")
ax[1].plot(l4_scales * l4_sweep, TE1_TM0_port1, c="purple", label="TM0 in port1")
ax[1].plot(l4_scales * l4_sweep, TE1_TM0_port2, c="green", label="TM0 in port2")
ax[1].set_xlim((0, 200))
ax[1].set_ylim((0, 1))
ax[1].legend()
ax[1].set_xlabel("L4 (μm)")
ax[1].set_ylabel("Mode conversion efficiency")
ax[1].set_title("TE1 Input")
plt.show()
l4 = l4_sweep * l4_scales[np.argmin(TE0_TE0_port2[:40])]
print(f"Optimal L4: {l4}")
Optimal L4: 117.08333333333331
Full FDTD Device Simulation#
Now we will combine our optimized taper and Y-junction to create a full polarizer/rotator in FDTD and examine the results:
[29]:
# Add buffer to extend structure into the boundary to avoid sidewall angle trouble
buffer = 10
# Create full polarization rotator/splitter structure in LNOI
full_pts = [
(-buffer, -w0 / 2),
(0, -w0 / 2),
(l1, -w1 / 2),
(l1 + l2, -w2 / 2),
(l1 + l2 + l3, -w3 / 2),
(l1 + l2 + l3 + 2 * l4, 2 * (-split / 2) - w5),
(l1 + l2 + l3 + 2 * l4, 2 * (-split / 2)),
(l1 + l2 + l3, -w3 / 2 + w5),
(l1 + l2 + l3 + 2 * l4, 2 * (split / 2)),
(l1 + l2 + l3 + 2 * l4, 2 * (split / 2) + w4),
(l1 + l2 + l3, w3 / 2),
(l1 + l2, w2 / 2),
(l1, w1 / 2),
(0, w0 / 2),
(-buffer, w0 / 2),
]
full_structure = make_ridge_waveguide(full_pts)
# Create mode monitors
full_port1 = td.ModeMonitor(
name="Full Port 1",
size=(0, w4 + 1.5 * split, td.inf),
center=(l1 + l2 + l3 + l4 - 1, (-w3 / 2 + w5) + split / 2 + w4 / 2, 0),
freqs=freqs[::2],
mode_spec=td.ModeSpec(
num_modes=num_modes,
target_neff=n_o,
angle_theta=np.arctan((split / 2 - (-w3 / 2 + w5)) / l4),
),
)
full_port2 = td.ModeMonitor(
name="Full Port 2",
size=(0, w5 + 1.5 * split, td.inf),
center=(l1 + l2 + l3 + l4 - 1, -split / 2 - w5 / 2, 0),
freqs=freqs[::2],
mode_spec=td.ModeSpec(
num_modes=num_modes,
target_neff=n_o,
angle_theta=np.arctan((-split / 2 - w5 / 2 + w3 / 4) / l4),
),
)
full_field = td.FieldMonitor(name="Full Field", size=(td.inf, td.inf, 0), freqs=[freq0])
auto_x = td.AutoGrid(min_steps_per_wvl=11)
# Create simulation
full_sim = td.Simulation(
center=((l1 + l2 + l3 + l4) / 2, 0, 0),
size=(l1 + l2 + l3 + l4, 4 * split + w4 + w5, 5 * film_thickness),
structures=full_structure,
grid_spec=td.GridSpec(grid_x=auto_x, grid_y=auto_yz, grid_z=auto_yz, wavelength=wvl0),
medium=SiO2,
sources=[],
monitors=[full_port1, full_port2, full_field],
run_time=5e-12,
)
"""
Solve for modes to be able to launch them as a source. We will create two simulations, one with the TE0 mode as a source, and one
with a TM0 mode as a source.
"""
full_mode_solver = ModeSolver(
simulation=full_sim,
plane=td.Box(center=(1, 0, 0), size=(0, td.inf, td.inf)),
mode_spec=td.ModeSpec(num_modes=3, target_neff=n_o),
freqs=[freq0],
)
full_source_TE0 = full_mode_solver.to_source(
mode_index=0, direction="+", source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10)
)
full_source_TM0 = full_source_TE0.updated_copy(mode_index=1)
full_sim_TE0 = full_sim.updated_copy(sources=[full_source_TE0])
full_sim_TM0 = full_sim.updated_copy(sources=[full_source_TM0])
Before we run we’ll check that we’ve constructed the simulation correctly.
[30]:
_, ax = plt.subplots(1, 2, figsize=(12, 4))
full_sim_TE0.plot(z=0.005, ax=ax[0])
full_sim_TE0.plot(x=l1 + l2 + l3 + l4, ax=ax[1])
ax[0].set_aspect(6)
plt.show()
[31]:
full_TE0_data = web.run(simulation=full_sim_TE0, task_name="Full FDTD pol/rot TE0")
full_TM0_data = web.run(simulation=full_sim_TM0, task_name="Full FDTD pol/rot TM0")
09:25:55 PDT Created task 'Full FDTD pol/rot TE0' with task_id 'fdve-01f39cf2-8570-4f78-ae70-c440c349bb0b' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-01f39cf2-857 0-4f78-ae70-c440c349bb0b'.
09:25:56 PDT 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.
09:26:01 PDT status = preprocess
09:26:17 PDT Maximum FlexCredit cost: 16.718. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
starting up solver
running solver
09:30:26 PDT early shutoff detected at 84%, exiting.
status = postprocess
09:30:33 PDT status = success
View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-01f39cf2-857 0-4f78-ae70-c440c349bb0b'.
09:35:11 PDT loading simulation from simulation_data.hdf5
Created task 'Full FDTD pol/rot TM0' with task_id 'fdve-75c4177d-9a1e-49c3-8a9e-debe4f7fa3de' and task_type 'FDTD'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-75c4177d-9a1 e-49c3-8a9e-debe4f7fa3de'.
09:35:13 PDT 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.
09:35:18 PDT status = preprocess
09:35:34 PDT Maximum FlexCredit cost: 16.718. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
starting up solver
running solver
09:39:41 PDT early shutoff detected at 84%, exiting.
09:39:42 PDT status = postprocess
09:39:48 PDT status = success
09:39:49 PDT View simulation result at 'https://tidy3d.simulation.cloud/workbench?taskId=fdve-75c4177d-9a1 e-49c3-8a9e-debe4f7fa3de'.
09:55:17 PDT loading simulation from simulation_data.hdf5
Full Device Results#
Plotting our results, we can see that we get high contrast and low insertion loss for both TE0 and TM0 inputs.
[32]:
full_TE0_TE0_port1 = np.abs(full_TE0_data["Full Port 1"].amps.sel(mode_index=0, direction="+")) ** 2
full_TE0_TE0_port2 = np.abs(full_TE0_data["Full Port 2"].amps.sel(mode_index=0, direction="+")) ** 2
full_TE0_TM0_port1 = np.abs(full_TE0_data["Full Port 1"].amps.sel(mode_index=1, direction="+")) ** 2
full_TE0_TM0_port2 = np.abs(full_TE0_data["Full Port 2"].amps.sel(mode_index=1, direction="+")) ** 2
full_TM0_TE0_port1 = np.abs(full_TM0_data["Full Port 1"].amps.sel(mode_index=0, direction="+")) ** 2
full_TM0_TE0_port2 = np.abs(full_TM0_data["Full Port 2"].amps.sel(mode_index=0, direction="+")) ** 2
full_TM0_TM0_port1 = np.abs(full_TM0_data["Full Port 1"].amps.sel(mode_index=1, direction="+")) ** 2
full_TM0_TM0_port2 = np.abs(full_TM0_data["Full Port 2"].amps.sel(mode_index=1, direction="+")) ** 2
_, ax = plt.subplots(2, 2, figsize=(12, 8))
full_TE0_data.plot_field("Full Field", "E", "abs^2", ax=ax[0][0])
full_TM0_data.plot_field("Full Field", "E", "abs^2", ax=ax[0][1])
ax[0][0].set_aspect(10)
ax[0][1].set_aspect(10)
ax[0][0].set_title("TE0 Input")
ax[0][1].set_title("TM0 Input")
ax[1][0].plot(wvls[::2], 10 * np.log10(full_TE0_TE0_port1), linewidth=2, label="TE0 in port 1")
ax[1][0].plot(wvls[::2], 10 * np.log10(full_TE0_TE0_port2), linewidth=2, label="TE0 in port 2")
ax[1][0].plot(wvls[::2], 10 * np.log10(full_TE0_TM0_port1), linewidth=2, label="TM0 in port 1")
ax[1][0].plot(wvls[::2], 10 * np.log10(full_TE0_TM0_port2), linewidth=2, label="TM0 in port 2")
ax[1][0].set_title("TE0 Input")
ax[1][0].set_xlabel("Wavelength (μm)")
ax[1][0].set_ylabel("Transmission (dB)")
ax[1][0].legend()
ax[1][1].plot(wvls[::2], 10 * np.log10(full_TM0_TE0_port1), linewidth=2, label="TE0 in port 1")
ax[1][1].plot(wvls[::2], 10 * np.log10(full_TM0_TE0_port2), linewidth=2, label="TE0 in port 2")
ax[1][1].plot(wvls[::2], 10 * np.log10(full_TM0_TM0_port1), linewidth=2, label="TM0 in port 1")
ax[1][1].plot(wvls[::2], 10 * np.log10(full_TM0_TM0_port2), linewidth=2, label="TM0 in port 2")
ax[1][1].set_title("TM0 Input")
ax[1][1].set_xlabel("Wavelength (μm)")
ax[1][1].set_ylabel("Transmission (dB)")
ax[1][1].legend()
plt.show()
[ ]: