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 ModeSolver.
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
import numpy as np
import tidy3d as td
# RF utilities for impedance and port definitions
import tidy3d.plugins.microwave as mw
import tidy3d.plugins.smatrix as sm
from matplotlib import pyplot as plt
from tidy3d import web
from tidy3d.plugins.mode import ModeSolver
# Suppress verbose logs to keep the notebook output clear
td.config.logging_level = "ERROR"
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_Lossless"]
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 ModeSolver object and running the mode simulation¶
[3]:
## Creating the ModeSolver
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 = ModeSolver(simulation=sim_opt, freqs=[freq_opt], mode_spec=mode_spec, plane=plane)
_ = mode_solver_opt.plot()
plt.show()
[4]:
mode_data_opt = web.run(mode_solver_opt, task_name="TFLN-opt_mode", folder_name="TFLN-VPI")
15:11:49 -03 Created task 'TFLN-opt_mode' with resource_id 'mo-95e96998-051a-4a00-8363-5ce569d1a4f6' and task_type 'MODE_SOLVER'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=mo-95e96998-051a- 4a00-8363-5ce569d1a4f6'.
Task folder: 'TFLN-VPI'.
Output()
15:11:53 -03 Estimated FlexCredit cost: 0.004. Minimum cost depends on task execution details. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
15:11:55 -03 status = success
Output()
15:12:04 -03 loading simulation 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_solver_opt.plot_field("E", "abs^2", mode_index=0)
mode_data_opt.to_dataframe()
plt.show()
[5]:
| wavelength | n eff | k eff | TE (Ey) fraction | wg TE fraction | wg TM fraction | mode area | group index | dispersion (ps/(nm km)) | ||
|---|---|---|---|---|---|---|---|---|---|---|
| f | mode_index | |||||||||
| 1.934145e+14 | 0 | 1.55 | 1.807221 | 0.0 | 0.997124 | 0.961891 | 0.870859 | 0.799926 | 2.214488 | -545.866993 |
| 1 | 1.55 | 1.721202 | 0.0 | 0.088426 | 0.800950 | 0.951599 | 1.230860 | 2.299366 | -1871.424868 | |
| 2 | 1.55 | 1.689922 | 0.0 | 0.941347 | 0.952167 | 0.875359 | 3.595095 | 2.092426 | -2536.536988 | |
| 3 | 1.55 | 1.682458 | 0.0 | 0.999967 | 0.995621 | 0.857858 | 3.013763 | 2.011008 | -954.413416 | |
| 4 | 1.55 | 1.678206 | 0.0 | 0.985192 | 0.969271 | 0.865669 | 2.993368 | 2.054942 | 135.879602 |
RF CPW Transmission Line ¶
Next, we define the CPW geometry and again use ModeSolver to calculate the RF mode.
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()
[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,
fields=["Ex", "Ey", "Ez", "Hx", "Hy", "Hz"],
): # 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],
)
# Integration paths for impedance calculation
# Define current integration path (loop around CPW signal)
I_int = mw.CurrentIntegralAxisAligned(
center=cpw_center.geometry.center,
size=(
0,
cpw_width + cpw_gap,
5 * cpw_thickness,
),
sign="+",
)
# Define voltage integration path (line from signal to ground)
V_int = mw.VoltageIntegralAxisAligned(
center=(0, cpw_width / 2 + cpw_gap / 2, cpw_center.geometry.center[2]),
size=(0, cpw_gap, 0),
sign="-",
)
# 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,
)
# Create and visualize the mode solver
plane = td.Box(
center=cpw_center.geometry.center,
size=(0, sim_vpi.size[1], sim_vpi.size[2]),
)
mzm_solver = ModeSolver(
simulation=sim_vpi,
plane=plane,
mode_spec=td.ModeSpec(num_modes=1, target_neff=2.2),
freqs=rf_freqs,
fields=fields,
)
return sim_vpi, mzm_solver, I_int, V_int
sim_vpi, mzm_solver, I_int, V_int = get_sim()
[9]:
# Visualize the mesh
ax = mzm_solver.plot()
mzm_solver.plot_grid(ax=ax)
center = mzm_solver.simulation.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 V and I integration paths
fig, ax = plt.subplots(figsize=(10, 3))
mzm_solver.plot(ax=ax)
I_int.plot(x=0, ax=ax)
V_int.plot(x=0, ax=ax)
ax.set_xlim(-20, 20)
ax.set_ylim(-10, 10)
plt.show()
Simulation Definition¶
Now, we define a function to can create the Simulation and ModeSolver objects.
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.
For impedance calculation, we define a current integration path around the CPW signal, and a voltage integration path from the signal to ground.
[11]:
mzm_mode_data = web.run(mzm_solver, task_name="MZM mode solver", path="data/MZM_mode_data.hdf5")
15:12:42 -03 Created task 'MZM mode solver' with resource_id 'mo-f356db97-e9ab-44ad-a652-38f14457d53c' and task_type 'MODE_SOLVER'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=mo-f356db97-e9ab- 44ad-a652-38f14457d53c'.
Task folder: 'default'.
Output()
15:12:45 -03 Estimated FlexCredit cost: 0.008. Minimum cost depends on task execution details. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
15:12:47 -03 status = success
Output()
15:13:09 -03 loading simulation from data/MZM_mode_data.hdf5
Visualizing the results:
[12]:
# Visualize the RF mode
ax = mzm_solver.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()
We can now use the built-in ImpedanceCalculator to calculate and visualize the impedance results.
[13]:
def get_impedance(data, I_int, V_int, plot=True):
n_eff = data.n_eff.isel(mode_index=0)
k_eff = data.k_eff.isel(mode_index=0)
# Calculate Z0
Z0_mzm = np.conjugate(
mw.ImpedanceCalculator(voltage_integral=V_int, current_integral=I_int).compute_impedance(
data
)
).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_solver.data, I_int, V_int)
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)
dict = {}
dic_integrals = {}
for g in g_values:
for w1 in w1_values:
sim_vpi, mzm_solver, I_int, V_int = get_sim(
cpw_gap=g,
cpw_width=w1,
h4=2,
h1=1,
optical_waveguide_gap=w1 + g / 2,
fields=["Ey", "Hy", "Hz"],
) # Fields for calculating impedance
dict[f"({g:.2f},{w1:.2f})"] = mzm_solver
dic_integrals[f"({g:.2f},{w1:.2f})"] = I_int, V_int
batch = web.Batch(simulations=dict, folder_name="mzm_solver")
[15]:
batch_data = batch.run(path_dir="batch_mzm_solver")
Output()
15:13:57 -03 Started working on Batch containing 100 tasks.
15:16:13 -03 Maximum FlexCredit cost: 0.895 for the whole batch.
Use 'Batch.real_cost()' to get the billed FlexCredit cost after the Batch has completed.
Output()
15:18:09 -03 Batch complete.
Output()
[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():
I_int, V_int = dic_integrals[task_name]
Z0_mzm, n_eff, k_eff = get_impedance(sim_data, I_int, V_int, 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()
[18]:
cpw_width = 18
cpw_gap = 5
cpw_thickness = 1
h1 = 0.5
optical_waveguide_gap = 25
sim_final, mzm_solver_final, I_int_final, V_int_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,
optical_waveguide_gap=optical_waveguide_gap,
refinement=20,
resolution=60,
) # g
mzm_solver_final_data = web.run(
mzm_solver_final, task_name="MZM mode solver", path="data/MZM_mode_data.hdf5"
)
16:23:18 -03 Created task 'MZM mode solver' with resource_id 'mo-e411ff34-43aa-489b-815b-5be0324c5b0e' and task_type 'MODE_SOLVER'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=mo-e411ff34-43aa- 489b-815b-5be0324c5b0e'.
Task folder: 'default'.
Output()
16:23:21 -03 Estimated FlexCredit cost: 0.009. Minimum cost depends on task execution details. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost after a simulation run.
16:23:22 -03 status = success
Output()
16:23:30 -03 loading simulation 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]:
# Compute voltage across the gap and normalize electric field
v0 = V_int_final.compute_voltage(mzm_solver_final_data)
ey_norm = (mzm_solver_final_data.Ey / (0.5 * v0.real)).isel(f=0, mode_index=0, drop=True).abs
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:
$$ \Delta_n = -0.5 n_e^3 r_{33} E $$
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 ModeSolver 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 = mode_solver_opt.updated_copy(
simulation=perturbed_sim, plane=td.Box(size=size, center=center)
)
return ms
Next, we create a Batch and run all mode simulations in parallel.
[22]:
# Define voltages to apply and create a Batch of simulations
voltages = np.arange(5)
sims = {str(v): perturbed_solver(v) for v in voltages}
batch = web.Batch(simulations=sims, folder_name="perturbed_mode")
We can now visualize the index change for 4V bias.
[23]:
eps0 = sims["0"].simulation.epsilon(sims["0"].plane)
eps4 = sims["4"].simulation.epsilon(sims["4"].plane)
delta_eps = eps4 - eps0
delta_eps.T.real.plot()
plt.show()
[24]:
# Run all simulations in parallel
batch_data = batch.run(path_dir="batch_perturbed_mode")
Output()
16:38:31 -03 Started working on Batch containing 5 tasks.
16:38:38 -03 Maximum FlexCredit cost: 0.020 for the whole batch.
Use 'Batch.real_cost()' to get the billed FlexCredit cost after the Batch has completed.
Output()
16:38:42 -03 Batch complete.
Output()
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.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 = 0.5 * opt_wavelength / dneff_dV # in V·μm"
print(f"Vπ·L (push-pull) = {vpil * 1e-4 / 2:.2f} V·cm")
Vπ·L (push-pull) = 1.60 V·cm
TERMS & CONDITIONS
How the process works:
Try this real URL to see an example: