Thin film lithium niobate electro-optic modulator#
In this example, we will simulate a thin-film lithium niobate electro-optic modulator based on the Pockels effect. The device is based on the work of *Ying Li et al.*, “High-Performance Mach–Zehnder Modulator Based on Thin-Film Lithium Niobate with Low Voltage-Length Product,” *ACS Omega* 2023, 8 (10), 9644–9651. DOI: 10.1021/acsomega.3c00310.
We start by calculating the propagating modes of the TFLN waveguide and the coplanar (CPW) transmission line using mode simulation.
Finally, we calculate the electro-optic overlap and predict the Vπ·L figure of merit.
You can check here the same model, along with the full chip layout integrated with a foundry PDK, using Photonforge, our photonic design automation tool.
Workflow Overview#
1. Define the thin-film lithium niobate waveguide geometry and solve for the optical TE mode using the ModeSimulation.
2. Build the CPW transmission line in the same cross-section and compute the RF mode and its overlap with the optical waveguide.
3. Normalize the RF field to 1 V across the electrodes, evaluate the electro-optic overlap, and estimate Vπ·L.
[1]:
# Core simulation and plotting libraries
from __future__ import annotations
import numpy as np
import tidy3d as td
# RF utilities for impedance and port definitions
import tidy3d.rf as rf
from matplotlib import pyplot as plt
from tidy3d import web
# Suppress verbose logs to keep the notebook output clear
td.config.logging.level = "ERROR"
# Use permittivity staircasing for local epsilon queries
td.config.simulation.use_local_subpixel = False
Optical Waveguide Geometry#
[2]:
# Simulation parameters
eff_inf = 1e3 # Large extent to emulate semi-infinite regions
opt_wavelength = 1.55
freq_opt = td.C_0 / opt_wavelength
# Optical materials: SiO2 cladding and TFLN core
SiO2 = td.material_library["SiO2"]["Palik_LowLoss"]
LiNbO3 = td.material_library["LiNbO3"]["Zelmon1997"](1)
# Ridge geometry parameters from the reference design
sidewall_angle = 17
core_width = 1.1 # w0
slab_thickness = 0.22 # h3
h3 = slab_thickness
core_thickness = 0.4 - slab_thickness
plane_size = 8
plane_limits = (-1.5, 1.9)
plane_height = 3.4
# Define waveguide structures
ridge = td.Structure(
name="ridge",
geometry=td.Box(
center=(0.0, 0.0, -slab_thickness / 2), size=(eff_inf, eff_inf, slab_thickness)
),
medium=LiNbO3,
)
core = td.Structure(
geometry=td.PolySlab(
sidewall_angle=sidewall_angle * np.pi / 180,
reference_plane="top",
slab_bounds=[0, core_thickness],
vertices=[
[-eff_inf, -core_width / 2],
[eff_inf, -core_width / 2],
[eff_inf, core_width / 2],
[-eff_inf, core_width / 2],
],
),
name="core",
medium=LiNbO3,
)
# Grid specification
grid_spec = td.GridSpec.auto(min_steps_per_wvl=20, wavelength=opt_wavelength)
# Create optical simulation
sim_opt = td.Simulation(
size=(10.0, plane_size, plane_height),
run_time=1e-12,
medium=SiO2,
structures=[ridge, core],
grid_spec=grid_spec,
)
Creating the ModeSimulation object and running the mode simulation#
[3]:
# # Creating the ModeSimulation
mode_spec = td.ModeSpec(num_modes=5, group_index_step=True)
plane_size = (0, sim_opt.bounding_box.size[1], sim_opt.bounding_box.size[2])
plane = sim_opt.bounding_box.updated_copy(size=plane_size)
mode_solver_opt = td.ModeSimulation.from_simulation(
simulation=sim_opt, plane=plane, mode_spec=mode_spec, freqs=[freq_opt]
)
_ = mode_solver_opt.plot()
plt.show()
[4]:
mode_data_opt = web.run(mode_solver_opt, task_name="TFLN-opt_mode", folder_name="TFLN-VPI")
22:06:30 UTC Created task 'TFLN-opt_mode' with resource_id 'mos-76a1461e-961c-45ff-ae87-bd26bcc29493' and task_type 'MODE'.
22:06:58 UTC Estimated FlexCredit cost: 0.007. For this solver type, the estimate is the final billed cost.
22:07:00 UTC 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.
22:07:07 UTC starting up solver
running solver
22:07:10 UTC status = success
22:07:11 UTC Loading results from simulation_data.hdf5
We can now visualize and inspect the optical modes. We are interested in the TE-like fundamental mode (with mode_index == 0), which has an effective index of 1.85 and a group index of 2.20. This information is very important to ensure velocity matching between the optical and RF modes, and hence optimize the electro-optical effect.
[5]:
mode_data_opt.plot_field("E", "abs^2", mode_index=0)
mode_data_opt.modes.to_dataframe()
plt.show()
RF CPW Transmission Line#
Next, we define the CPW geometry and use a WavePort to solve the RF mode and automatically extract its characteristic impedance.
For more information on CPW mode calculation, one can refer to this example.
[6]:
# RF frequency range
rf_freqs = np.linspace(10e9, 90e9, 21)
# RF materials
si_rf = td.Medium(permittivity=11.7)
sio2_rf = td.Medium(permittivity=3.9)
air_rf = td.Medium()
tfln_rf = td.AnisotropicMedium(
xx=td.Medium(permittivity=43), yy=td.Medium(permittivity=27.9), zz=td.Medium(permittivity=43)
)
metal = td.LossyMetalMedium(
frequency_range=(100000000.0, 100000000000.0),
conductivity=41,
fit_param=td.SurfaceImpedanceFitterParam(
max_num_poles=16,
),
)
Since we will sweep over many variables, it is convenient to create a function to return the CPW and optical waveguide geometries.
[7]:
def get_structures(
h4=4.7,
h2=0.4,
h1=1,
cpw_thickness=1, # CPW geometry parameters
cpw_width=18, # w1
cpw_gap=5, # g
optical_waveguide_gap=25,
): # G
sio2_thickness = h1 + h3 + h4
si_substrate = td.Structure(
geometry=td.Box.from_bounds(
rmin=(-eff_inf, -eff_inf, -eff_inf), rmax=(eff_inf, eff_inf, -sio2_thickness / 2)
),
medium=si_rf,
)
# Define substrate and layer structures
sio2_slab = td.Structure(
geometry=td.Box(center=(0, 0, 0), size=(eff_inf, eff_inf, sio2_thickness)), medium=sio2_rf
)
# Optical waveguide
center_slab = sio2_slab.geometry.bounds[0][2] + h4 + slab_thickness / 2
tfln_slab = td.Structure(
geometry=td.Box(center=(0, 0, center_slab), size=(eff_inf, eff_inf, slab_thickness)),
medium=tfln_rf,
)
# Calculate vertical offset to position waveguides at TFLN slab level
delta_z = tfln_slab.geometry.center[2] - ridge.geometry.center[2]
# Define conformal cladding for waveguides
delta_z_conformal = sio2_slab.geometry.bounds[1][2]
waveguide_conformal_cladding = td.PolySlab(
sidewall_angle=sidewall_angle * np.pi / 180,
reference_plane="top",
slab_bounds=[delta_z_conformal, delta_z_conformal + core_thickness],
vertices=[
[-eff_inf, -core_width / 1.5],
[eff_inf, -core_width / 1.5],
[eff_inf, core_width / 1.5],
[-eff_inf, core_width / 1.5],
],
)
# Define right waveguide positioned under CPW gap
waveguide_r = [
ridge.updated_copy(
geometry=ridge.geometry.translated(
x=0, y=optical_waveguide_gap / 2 + core_width / 2, z=delta_z
),
name="ridge_r",
medium=tfln_rf,
),
core.updated_copy(
geometry=core.geometry.translated(
x=0, y=optical_waveguide_gap / 2 + core_width / 2, z=delta_z
),
name="core_r",
medium=tfln_rf,
),
td.Structure(
geometry=waveguide_conformal_cladding.translated(
x=0, y=optical_waveguide_gap / 2 + core_width / 2, z=0
),
medium=sio2_rf,
),
]
# Define left waveguide positioned under CPW gap
waveguide_l = [
ridge.updated_copy(
geometry=ridge.geometry.translated(
x=0, y=-optical_waveguide_gap / 2 - core_width / 2, z=delta_z
),
name="ridge_l",
medium=tfln_rf,
),
core.updated_copy(
geometry=core.geometry.translated(
x=0, y=-optical_waveguide_gap / 2 - core_width / 2, z=delta_z
),
name="core_l",
medium=tfln_rf,
),
td.Structure(
geometry=waveguide_conformal_cladding.translated(
x=0, y=-optical_waveguide_gap / 2 - core_width / 2, z=0
),
medium=sio2_rf,
),
]
cpw_center_pos = sio2_slab.geometry.bounds[1][2] + cpw_thickness / 2
# Define CPW transmission line structures
cpw_left = td.Structure(
geometry=td.Box.from_bounds(
rmin=(-eff_inf, -eff_inf, cpw_center_pos - cpw_thickness / 2),
rmax=(eff_inf, -cpw_width / 2 - cpw_gap, cpw_center_pos + cpw_thickness / 2),
),
medium=metal,
)
cpw_right = td.Structure(
geometry=td.Box.from_bounds(
rmin=(-eff_inf, cpw_width / 2 + cpw_gap, cpw_center_pos - cpw_thickness / 2),
rmax=(eff_inf, eff_inf, cpw_center_pos + cpw_thickness / 2),
),
medium=metal,
)
cpw_center = td.Structure(
geometry=td.Box.from_bounds(
rmin=(-eff_inf, -cpw_width / 2, cpw_center_pos - cpw_thickness / 2),
rmax=(eff_inf, cpw_width / 2, cpw_center_pos + cpw_thickness / 2),
),
medium=metal,
)
air_gap = td.Structure(
geometry=td.Box.from_bounds(
rmin=(-eff_inf, -eff_inf, sio2_thickness / 2), rmax=(eff_inf, eff_inf, eff_inf)
),
medium=air_rf,
)
return (
[sio2_slab, si_substrate, air_gap, cpw_left, cpw_right, cpw_center, tfln_slab]
+ waveguide_l
+ waveguide_r
)
structures = get_structures()
The `WavePort <https://docs.flexcompute.com/projects/tidy3d/en/latest/api/microwave/_autosummary/tidy3d.rf.WavePort.html>`__ size is chosen from the CPW cross-section. For `AutoImpedanceSpec <https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.rf.AutoImpedanceSpec.html>`__, the lateral port boundaries must intersect conducting regions, so the port should include the signal electrode, the gaps, and part of the ground electrodes. Here, using the full transverse
simulation size works because the ground electrodes extend to the lateral simulation boundaries.
[8]:
def get_sim(
h4=4.7, # top SiO2 cladding thickness (µm)
h2=0.4,
h1=1, # lower layer thickness beneath the electrodes (µm)
cpw_thickness=1, # CPW geometry parameters (t)
cpw_width=18, # w1
cpw_gap=5,
optical_waveguide_gap=23 - core_width / 2, # G
refinement=10,
resolution=30,
): # g)
structures = get_structures(
h4, h2, h1, cpw_thickness, cpw_width, cpw_gap, optical_waveguide_gap
)
# Define layer refinement specification for enhanced mesh at CPW corners
cpw_center = structures[5]
cpw_left = structures[3]
cpw_right = structures[4]
tfln_slab = structures[6]
lr_spec = td.LayerRefinementSpec.from_structures(
structures=[cpw_center, cpw_left, cpw_right],
min_steps_along_axis=5,
corner_refinement=td.GridRefinement(dl=cpw_thickness / refinement, num_cells=5),
refinement_inside_sim_only=False,
)
# Mesh override for waveguide region (fine mesh in z-direction)
override_region = td.MeshOverrideStructure(
geometry=td.Box(
size=(eff_inf, eff_inf, slab_thickness + core_thickness),
center=(0, 0, tfln_slab.geometry.center[2] + slab_thickness / 2),
),
dl=(None, None, 0.05),
)
# Calculate gap center position
gap_center = cpw_center.geometry.center[1] + cpw_width / 2 + cpw_gap / 2
# Mesh override for left CPW gap (fine mesh in y-direction)
override_region_gap_l = td.MeshOverrideStructure(
geometry=td.Box(size=(eff_inf, cpw_gap, eff_inf), center=(0, gap_center, 0)),
dl=(None, 0.05, None),
)
# Mesh override for right CPW gap (fine mesh in y-direction)
override_region_gap_r = td.MeshOverrideStructure(
geometry=td.Box(size=(eff_inf, cpw_gap, eff_inf), center=(0, -gap_center, 0)),
dl=(None, 0.05, None),
)
# Overall grid specification
grid_spec = td.GridSpec.auto(
min_steps_per_sim_size=resolution,
wavelength=td.C_0 / max(rf_freqs),
layer_refinement_specs=[lr_spec],
override_structures=[override_region, override_region_gap_l, override_region_gap_r],
)
# Simulation domain size
sim_rf_size = (25, 138, 150)
# Create RF simulation
sim_vpi = td.Simulation(
size=sim_rf_size,
structures=structures,
grid_spec=grid_spec,
run_time=1e-12,
)
# Port boundaries must intersect the ground electrodes for AutoImpedanceSpec.
wave_port = rf.WavePort(
center=(0, 0, 0),
size=(0, sim_vpi.size[1], sim_vpi.size[2]),
direction="+",
name="rf_port",
num_grid_cells=None, # keep the simulation grid unchanged
mode_spec=rf.MicrowaveModeSpec(
num_modes=1,
target_neff=2.2,
impedance_specs=rf.AutoImpedanceSpec(), # automatic V, I and Z0
),
)
# Build the ModeSimulation directly from the WavePort
mzm_solver = wave_port.to_mode_simulation(simulation=sim_vpi, freqs=rf_freqs)
return sim_vpi, mzm_solver, wave_port
sim_vpi, mzm_solver, wave_port = get_sim()
[9]:
# Visualize the mesh
ax = mzm_solver.plot(x=0)
mzm_solver.plot_grid(ax=ax, x=0)
center = mzm_solver.structures[0].geometry.center
ax.set_xlim(center[0] - 20, center[0] + 20)
ax.set_ylim(center[1] - 10, center[1] + 10)
plt.show()
[10]:
# Plot the WavePort plane (translucent) over the cross-section
fig, ax = plt.subplots(figsize=(10, 3))
mzm_solver.plot(ax=ax)
wave_port.plot(x=0, ax=ax, alpha=0.3)
ax.set_xlim(-20, 20)
ax.set_ylim(-10, 10)
plt.show()
Simulation Definition#
Now, we define a function to create the Simulation and the RF port.
Since the RF wavelength is much larger than the geometric features, we use a LayerRefinementSpec at the CPW to automatically enhance the mesh and better resolve its corners.
For the rest of the simulation, we use the automatic grid with min_steps_per_sim_size = 100, which ensures at least 100 grid cells along the longest dimension of the simulation domain.
We will also add mesh override regions to properly resolve the thickness of the optical waveguide, as well as the gaps, to ensure accurate results.
To excite and analyze the CPW, we use a WavePort. Combined with an AutoImpedanceSpec inside its MicrowaveModeSpec, the port automatically determines the voltage and
current integration paths from the conductors that intersect its plane and computes the characteristic impedance — so there is no need to define integration paths or run an ImpedanceCalculator by hand. We obtain the 2D mode simulation directly from the port via wave_port.to_mode_simulation(...).
[11]:
mzm_mode_data = web.run(mzm_solver, task_name="MZM mode solver", path="data/MZM_mode_data.hdf5")
22:07:13 UTC Created task 'MZM mode solver' with resource_id 'mos-04ea42d9-3fa1-472b-becd-97f1e4729dd4' and task_type 'MODE'.
22:07:14 UTC Estimated FlexCredit cost: 0.015. For this solver type, the estimate is the final billed cost.
22:07:15 UTC 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.
22:07:25 UTC starting up solver
running solver
22:07:35 UTC status = success
22:07:41 UTC Loading results from data/MZM_mode_data.hdf5
Visualizing the results:
[12]:
# Visualize the RF mode
ax = mzm_mode_data.plot_field(
"E", "abs", f=rf_freqs[0], mode_index=0, robust=True, shading="gouraud"
)
ax.set_ylim(1, 5)
ax.set_xlim(5, 18)
plt.show()
The characteristic impedance is computed automatically by the WavePort and stored in the mode solver results. We can read it directly from transmission_line_data.Z0, together with the effective index and attenuation, without any manual integration paths.
[13]:
def get_impedance(data, plot=True):
n_eff = data.n_eff.isel(mode_index=0)
k_eff = data.k_eff.isel(mode_index=0)
# The characteristic impedance is computed automatically by the WavePort.
# Use np.conjugate to convert from the physics phase convention (Tidy3D
# default) to the standard electrical engineering convention.
Z0_mzm = np.conjugate(data.transmission_line_data.Z0.isel(mode_index=0)).squeeze()
if plot:
fig, ax = plt.subplots(1, 2, figsize=(8, 3), tight_layout=True)
ax[0].plot(rf_freqs * 1e-9, n_eff, ".-")
ax[0].set(ylabel="Effective Index", xlabel="Frequency (GHz)")
ax[1].plot(rf_freqs * 1e-9, np.real(Z0_mzm), ".-")
_ = ax[1].set(ylabel="Re{Z₀} (Ω)", xlabel="Frequency (GHz)")
return Z0_mzm, n_eff, k_eff
Z0_mzm, n_eff, k_eff = get_impedance(mzm_mode_data.modes)
plt.show()
Sweep Over Different Parameters#
Now, we can sweep over different g and w1 parameters, and calculate impedance, effective index and attenuation constant
[14]:
g_values = np.linspace(2, 16, 10)
w1_values = np.linspace(3, 25, 10)
sims = {}
for g in g_values:
for w1 in w1_values:
sim_vpi, mzm_solver, wave_port = get_sim(
cpw_gap=g,
cpw_width=w1,
h4=2,
h1=1,
optical_waveguide_gap=w1 + g / 2,
)
sims[f"({g:.2f},{w1:.2f})"] = mzm_solver
[15]:
batch_data = web.run(
simulation=sims,
folder_name="mzm_solver",
path="batch_mzm_solver",
)
22:09:07 UTC Started working on Batch containing 100 tasks.
22:10:31 UTC Maximum FlexCredit cost: 1.993 for the whole batch.
Use 'Batch.real_cost()' to get the billed FlexCredit cost after completion.
22:14:10 UTC Batch complete.
Data Processing#
Now we can calculate impedance, RF effective index, RF losses, and reproduce Figures 4a, 4b, and 4c from the reference paper.
[16]:
import pandas as pd
df_z0 = pd.DataFrame()
df_z0.index.name = "h1"
df_z0.columns.name = "cpw_thickness"
df_neff = pd.DataFrame()
df_neff.index.name = "h1"
df_neff.columns.name = "cpw_thickness"
df_k_eff = pd.DataFrame()
df_k_eff.index.name = "h1"
df_k_eff.columns.name = "cpw_thickness"
for task_name, sim_data in batch_data.items():
Z0_mzm, n_eff, k_eff = get_impedance(sim_data.modes, plot=False)
h1, cpw_thickness = eval(task_name)
df_z0.loc[h1, cpw_thickness] = Z0_mzm.sel(f=70e9, method="nearest").values
df_neff.loc[h1, cpw_thickness] = n_eff.sel(f=70e9, method="nearest").values
df_k_eff.loc[h1, cpw_thickness] = k_eff.sel(f=70e9, method="nearest").values
[17]:
fig, Ax = plt.subplots(ncols=3, figsize=(17, 5))
pm1 = Ax[0].pcolormesh(
df_z0.index, df_z0.columns, df_z0.values.T.real, cmap="jet", shading="gouraud"
)
cb1 = plt.colorbar(pm1)
cb1.set_label("Impedance (Ω)")
pm2 = Ax[1].pcolormesh(
df_neff.index, df_neff.columns, df_neff.values.T, cmap="jet", shading="gouraud"
)
cb2 = plt.colorbar(pm2)
cb2.set_label("Effective Index")
df_alpha = 2 * np.pi * (70e9 / td.C_0) * df_k_eff * 10**6
pm3 = Ax[2].pcolormesh(
df_alpha.index, df_alpha.columns, df_alpha.values.T, cmap="jet", shading="gouraud"
)
cb3 = plt.colorbar(pm3)
cb3.set_label("Attenuation Constant $\\alpha$ (1/m)")
for ax in Ax:
ax.set_ylim(Ax[2].get_ylim()[::-1])
ax.set_xlabel("h1 (µm)")
plt.show()
Computation of Vπ·L#
Now we can use the RF fields to calculate the Pockels effect and Vπ·L.
First, we create a final simulation with the optimized parameters, as reported in the reference paper.
[18]:
cpw_width = 18
cpw_gap = 5
cpw_thickness = 1
h1 = 0.5
optical_waveguide_gap = 25
sim_final, mzm_solver_final, wave_port_final = get_sim(
h4=4.7, # top SiO2 cladding thickness (µm)
h2=0.4,
h1=h1, # lower layer thickness beneath the electrodes (µm)
cpw_thickness=cpw_thickness, # CPW geometry parameters (t)
cpw_width=cpw_width, # w1
cpw_gap=cpw_gap, # g
optical_waveguide_gap=optical_waveguide_gap,
refinement=20,
resolution=60,
)
mzm_solver_final_data = web.run(
mzm_solver_final, task_name="MZM mode solver", path="data/MZM_mode_data.hdf5"
)
22:21:43 UTC Created task 'MZM mode solver' with resource_id 'mos-75b38df9-86a7-431c-ba14-76f65c525dc9' and task_type 'MODE'.
22:22:14 UTC Estimated FlexCredit cost: 0.019. For this solver type, the estimate is the final billed cost.
22:22:15 UTC 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.
22:22:22 UTC starting up solver
running solver
22:22:32 UTC status = success
22:22:43 UTC Loading results from data/MZM_mode_data.hdf5
To calculate the Pockels effect, we apply an electric field along the LiNbO₃ cut direction to perturb the optical medium. First, we normalize the field by the applied voltage on the push-pull configuration.
[19]:
# Modal voltage computed automatically by the WavePort, used to normalize the field
v0 = mzm_solver_final_data.modes.transmission_line_data.voltage_coeffs.isel(mode_index=0)
ey_norm = (mzm_solver_final_data.modes.Ey / v0.real).isel(f=0, mode_index=0, drop=True).real
waveguide_center = sim_final.structures[8].geometry.bounding_box.center
ax = ey_norm.isel(x=0).transpose("z", "y").plot(robust=True).axes
ax.set_xlim(waveguide_center[1] - 2, waveguide_center[1] + 2)
ax.set_ylim(waveguide_center[2] - 1.5, waveguide_center[2] + 1.5)
ax.set_aspect("equal")
plt.show()
The normalized electric field is applied to the LiNbO₃ crystal along its z-axis, following the Pockels effect model:
Where: - \(\Delta_n\): Index variation - \(n_e\): Extraordinary refractive index - \(r_{33}\): Pockels coefficient (30.9 pm/V) - \(E\): Normalized electric field from CPW (1V across gap)
[20]:
# Extract LiNbO3 material properties at optical wavelength
eps_o, eps_e, _ = (e.real for e in LiNbO3.eps_diagonal(td.C_0 / opt_wavelength))
n_e = eps_e**0.5
r33 = 30.9e-6 # Pockels coefficient in μm/V
# Calculate refractive index variation from Pockels effect
Δn = -0.5 * n_e**3 * r33 * ey_norm
# Plot index variation (perturbation will only take effect over the LiNbO₃ regions)
ax = Δn.isel(x=0).transpose("z", "y").plot(robust=False).axes
ax.set_xlim(waveguide_center[1] - 2, waveguide_center[1] + 2)
ax.set_ylim(waveguide_center[2] - 1.5, waveguide_center[2] + 1.5)
ax.set_aspect("equal")
plt.show()
Now, assuming that the field scales with the applied voltage, we can define an auxiliary function to create the perturbed medium using a CustomAnisotropicMedium object as a function of the voltage.
This incorporates the index change due to the RF fields and returns a ModeSimulation object.
[21]:
# Use a single data point for the homogeneous directions
eps_o_array = td.SpatialDataArray(np.full((1, 1, 1), eps_o), coords={"x": [0], "y": [0], "z": [0]})
def perturbed_solver(voltage):
perturbed_ln = td.CustomAnisotropicMedium(
xx=td.CustomMedium(permittivity=eps_o_array, subpixel=True),
yy=td.CustomMedium(permittivity=(n_e + voltage * Δn) ** 2, subpixel=True),
zz=td.CustomMedium(permittivity=eps_o_array, subpixel=True),
)
# Changing the RF mediums to optical mediums
si02_structures_index = [0, 9, 12]
ln_structures_index = [7, 8, 10, 11]
perturbed_structures = list(sim_final.structures)
for i in si02_structures_index:
perturbed_structures[i] = sim_final.structures[i].updated_copy(medium=SiO2)
for i in ln_structures_index:
perturbed_structures[i] = sim_final.structures[i].updated_copy(medium=perturbed_ln)
# Remove metals
del perturbed_structures[2:5]
center = waveguide_center
size = (0, 8, 3)
perturbed_sim = sim_final.updated_copy(
structures=perturbed_structures, size=size, center=center, grid_spec=sim_opt.grid_spec
)
ms = td.ModeSimulation.from_simulation(
simulation=perturbed_sim,
plane=td.Box(size=size, center=center),
mode_spec=mode_spec,
freqs=[freq_opt],
)
return ms
Next, we run the voltage-dependent mode simulations with web.run.
[22]:
# Define voltages to apply and create mode simulations
voltages = np.arange(5)
sims = {str(v): perturbed_solver(v) for v in voltages}
We can now visualize the index change for 4V bias.
[23]:
eps0 = sims["0"].epsilon(sims["0"].plane)
eps4 = sims["4"].epsilon(sims["4"].plane)
delta_eps = eps4 - eps0
delta_eps.T.real.plot()
plt.show()
[24]:
# Run all simulations
batch_data = web.run(
simulation=sims,
folder_name="perturbed_mode",
path="batch_perturbed_mode",
)
22:23:24 UTC Started working on Batch containing 5 tasks.
22:23:29 UTC Maximum FlexCredit cost: 0.036 for the whole batch.
Use 'Batch.real_cost()' to get the billed FlexCredit cost after completion.
22:24:26 UTC Batch complete.
Finally, we can observe the variation of the waveguide mode effective index as a function of the applied voltage and calculate the Vπ·L figure of merit.
[25]:
n_eff = [sim_data.modes.n_eff.isel(mode_index=0, f=0).item() for sim_data in batch_data.values()]
_ = plt.plot(voltages, n_eff, ".-")
_ = plt.ticklabel_format(style="plain", useOffset=False)
plt.gca().set(xlabel="Voltage (V)", ylabel="Optical Effective Index")
plt.show()
[26]:
# The variation is quite linear, so we simply use the largest interval
dneff_dV = abs(n_eff[-1] - n_eff[0]) / (voltages[-1] - voltages[0])
vpil = (1 / 2) * opt_wavelength / dneff_dV # in V·μm"
print(f"Vπ·L (push-pull) = {vpil * 1e-4 / 2:.2f} V·cm")
Vπ·L (push-pull) = 3.25 V·cm