TIDY3D
LEARNING CENTER

Lectures
Tutorials
Applications

Example Library / Active Photonic Integrated Circuit Components / Thin film lithium niobate electro-optic modulator

Thin film lithium niobate electro-optic modulator

Sign up now to simulate .ipynb

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.

Schematic of the grating coupler

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.

# 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 

# 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

# # 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()

mode_data_opt = web.run(mode_solver_opt, task_name="TFLN-opt_mode", folder_name="TFLN-VPI")
17:45:51 UTC Created task 'TFLN-opt_mode' with resource_id                      
             'mos-df4e866d-03c2-4432-a530-796a7311e069' and task_type 'MODE'.   







17:50:15 UTC Estimated FlexCredit cost: 0.007. For this solver type, the        
             estimate is the final billed cost.                                 







17:50:16 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.                                                             


















17:51:01 UTC starting up solver                                                 







             running solver                                                     







17:51:06 UTC status = success                                                   
























17:51:08 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.






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 the Flex RF CPW in RF Photonics 1: Transmission Line Basics example.






# 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.






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 size is chosen from the CPW cross-section. For AutoImpedanceSpec, 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.






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()











# 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()

























# 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(...).






mzm_mode_data = web.run(mzm_solver, task_name="MZM mode solver", path="data/MZM_mode_data.hdf5")







17:51:10 UTC Created task 'MZM mode solver' with resource_id                    
             'mos-1e2a39a9-1358-48d2-a3b8-e090b5a09357' and task_type 'MODE'.   


















17:51:12 UTC Estimated FlexCredit cost: 0.015. For this solver type, the        
             estimate is the final billed cost.                                 







17:51:13 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.                                                             


















17:51:30 UTC starting up solver                                                 







             running solver                                                     







17:51:40 UTC status = success                                                   
























17:52:07 UTC Loading results from data/MZM_mode_data.hdf5                       







Visualizing the results:






# 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.






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






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











batch_data = web.run(
    simulation=sims,
    folder_name="mzm_solver",
    path="batch_mzm_solver",
)


















17:53:36 UTC Started working on Batch containing 100 tasks.                     







17:55:03 UTC Maximum FlexCredit cost: 1.993 for the whole batch.                







             Use 'Batch.real_cost()' to get the billed FlexCredit cost after    
             completion.                                                        












17:56:26 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.






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











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.






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"
)







18:25:49 UTC Created task 'MZM mode solver' with resource_id                    
             'mos-91b66dcc-3910-4da0-937d-778cc7461993' and task_type 'MODE'.   


















18:26:19 UTC Estimated FlexCredit cost: 0.019. For this solver type, the        
             estimate is the final billed cost.                                 







18:26:21 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.                                                             


















18:26:56 UTC starting up solver                                                 







             running solver                                                     







18:27:10 UTC status = success                                                   
























18:27:27 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.






# 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:


\[ \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)






# 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.






# 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.






# 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.






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()

























# Run all simulations
batch_data = web.run(
    simulation=sims,
    folder_name="perturbed_mode",
    path="batch_perturbed_mode",
)


















18:28:13 UTC Started working on Batch containing 5 tasks.                       







18:28:18 UTC Maximum FlexCredit cost: 0.036 for the whole batch.                







             Use 'Batch.real_cost()' to get the billed FlexCredit cost after    
             completion.                                                        












18:28:36 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.






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()

























# 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












          

          

            Example library
          

          
          

          

            

              


                

                  More Examples
                



              Read More
            

            

              

              

                
                
                

                  
                    Photo-thermal optical control in silicon waveguides
                    

                    Photo-thermal optical control in silicon waveguides
                  
                

                
                

                  
                    Forward-bias PIN phase shifter with thermal analysis
                    

                    Forward-bias PIN phase shifter with thermal analysis
                  
                

                
                

                  
                    Avalanche photodiode
                    

                    Avalanche photodiode
                  
                

                
                
              

              

              

                
No matching results found.

              

            

          

          

          
          

          

            
            
              
              Previous Tutorial
            
            
            
            
              Next Tutorial