TIDY3D
LEARNING CENTER

Lectures
Tutorials
Applications

Example Library / EME 🌈 / EME convergence for a bent anisotropic waveguide (LiNbO₃)

EME convergence for a bent anisotropic waveguide (LiNbO₃)

Sign up now to simulate .ipynb

Lithium niobate (LiNbO₃) is a key material in integrated photonics due to its large electro-optic coefficient and wide transparency window. However, it is a uniaxial crystal: its permittivity is described by a tensor, not a scalar. This has important consequences when the waveguide is bent.

In this notebook we:

  1. Run an FDTD simulation of a 180° LiNbO₃ waveguide bend to obtain a reference transmittance.
  2. Run EME simulations of the same bend with increasing numbers of cells.
  3. Show that a single EME cell overestimates transmission and explain the physical origin.
  4. Demonstrate convergence of EME to the FDTD result as the cell count increases.

For an introduction to the EME solver in Tidy3D see the EME solver tutorial and the bent waveguides with EME notebook.

from __future__ import annotations

import matplotlib.pyplot as plt
import numpy as np
import tidy3d as td
import tidy3d.web as web
from tidy3d.plugins.mode import ModeSolver

td.config.logging.level = "ERROR"

Physics Background: Why Anisotropy Complicates Bent-Waveguide EME

The Rotating Permittivity Tensor

LiNbO₃ is uniaxial: it has one extraordinary index \(n_e \approx 2.14\) along the optic (c) axis and an ordinary index \(n_o \approx 2.21\) along the other two directions (at 1.55 µm). In the material API, Zelmon1997(1) places the extraordinary axis along the global \(y\) direction (index 2 would place it along \(z\)). The lab-frame tensor used here is therefore

\[\varepsilon_\text{lab} = \begin{pmatrix} n_o^2 & 0 & 0 \\ 0 & n_e^2 & 0 \\ 0 & 0 & n_o^2 \end{pmatrix} \qquad \text{(optic axis along } y\text{)} .\]

The key point is that the optic axis lies in the plane of the bend (the \(x\)–\(y\) plane). In a straight waveguide propagating along \(x\) the tensor is constant — modes are quasi-TE / quasi-TM and propagate without cross-coupling. However, in a bent waveguide (ring in the \(x\)–\(y\) plane) the local propagation direction rotates by an angle \(\phi\) as the wave travels around the bend. From the co-moving frame the effective tensor becomes

\[\varepsilon_\text{local}(\phi) = R_z(\phi)^\top \, \varepsilon_\text{lab} \, R_z(\phi) .\]

Because the optic axis lies in the rotation (\(x\)–\(y\)) plane and \(n_o \neq n_e\), this rotation produces off-diagonal terms that grow with \(\phi\). The local eigenmodes therefore evolve continuously as the wave travels — the quasi-TE and quasi-TM modes of the bent waveguide mix along the propagation path.

Consequence for EME

EME divides the waveguide into \(N\) cells and matches eigenmodes at each cell boundary. At every interface the two adjacent mode bases are slightly misaligned due to the incremental tensor rotation — generating partial back-reflections and mode conversion.

EME cells Physical picture
1 Entire 180° bend is one uniform section. No inter-cell interfaces → no reflections from tensor rotation → transmission is overestimated.
few Tensor rotation is approximated as a staircase; dominant back-reflection is captured.
many Staircase finely resolves the continuous rotation → converges to FDTD.

The loss mechanism is not geometric bending loss (negligible at \(r = 5\) µm for this waveguide) but polarization-state evolution driven by the rotating anisotropy.

# ── Simulation parameters ──────────────────────────────────────────────────
wavelength = 1.55  # free-space wavelength (µm)
freq0 = td.C_0 / wavelength

radius = 5.0  # bend radius (µm)
width = 0.6  # waveguide width (µm)
thickness = 0.4  # waveguide thickness (µm)

num_modes = 5  # number of modes in mode solvers
target_neff = 1.65  # target effective index for mode search
plane_size = (9.0, 9.0)  # cross-section window size (µm)
min_steps_per_wvl = 30  # grid resolution (steps per wavelength)

# ── Materials ───────────────────────────────────────────────────────────────
SiO2 = td.material_library["SiO2"]["Palik_Lossless"]
# (1): extraordinary (optic) axis along global y — in the bend plane, the
# orientation-sensitive case. Use (2) to place the optic axis along z instead.
LiNbO3 = td.material_library["LiNbO3"]["Zelmon1997"](1)

Step 1 — FDTD Reference Simulation

We model the 180° bend as a full ring in the FDTD simulation:

  • Geometry: outer Cylinder minus inner cylinder (LiNbO₃ annulus on SiO₂).
  • Source: ModeSource at the bottom of the ring (0, −r, 0), launching in the +x direction with a negative bend radius (the center of curvature is on the +y side of the source).
  • Monitor: ModeMonitor at the top of the ring (0, +r, 0), measuring the amplitudes of all guided modes travelling in the −x direction (i.e., those that have propagated 180° around the bend).
  • Boundaries: Absorber on ±x to suppress re-entry of any radiation that escapes the open side of the ring; PML on ±y, ±z.

The bend curvature is encoded in the ModeSpec via bend_radius and bend_axis.

The source and monitor sort their modes with ModeSortSpec(sort_key="fill_fraction_box", bounding_box=core_box): in a bent, anisotropic waveguide the default ordering (by effective index) does not reliably put the guided fundamental first, so we instead rank modes by the fraction of their intensity that falls inside a box wrapping the core. The most core-confined (guided) mode then lands at mode_index = 0 regardless of polarization or tensor-rotation mixing. (fill_fraction_box ranks by confinement, not polarization — the polarization is verified separately.)

def build_fdtd_bend(
    radius=5.0,
    width=0.6,
    thickness=0.4,
    num_modes=10,
    target_neff=1.65,
    freq0=td.C_0 / 1.55,
    core_medium=LiNbO3,
    background_medium=SiO2,
    plane_size=(9.0, 9.0),
    min_steps_per_wvl=30,
    mode_index=0,
):
    """Build an FDTD simulation of a 180-degree LiNbO3 waveguide bend.

    The geometry is a ring (annular cylinder) with a ModeSource at the
    bottom and a ModeMonitor at the top, measuring 180-degree transmission.
    """
    wavelength = td.C_0 / freq0
    cross_sec = (0, plane_size[0], plane_size[1])  # cross-section for mode planes

    core_box_src = td.Box(center=(0, -radius, 0), size=(0, width, thickness))
    core_box_mon = td.Box(center=(0, radius, 0), size=(0, width, thickness))

    sort_spec_src = td.ModeSortSpec(sort_key="fill_fraction_box", bounding_box=core_box_src)
    sort_spec_mon = td.ModeSortSpec(sort_key="fill_fraction_box", bounding_box=core_box_mon)

    # Bend mode specs
    # Source at bottom (0, -r, 0) going in +x: center of curvature is at +y → negative radius
    # Monitor at top  (0, +r, 0) receiving in -x: center of curvature is at -y → positive radius
    mode_spec_src = td.ModeSpec(
        num_modes=num_modes,
        target_neff=target_neff,
        bend_radius=-radius,
        bend_axis=1,
        num_pml=(12, 12),
        sort_spec=sort_spec_src,
        precision="double",
    )
    mode_spec_mon = td.ModeSpec(
        num_modes=num_modes,
        target_neff=target_neff,
        bend_radius=radius,
        bend_axis=1,
        num_pml=(12, 12),
        sort_spec=sort_spec_mon,
        precision="double",
    )

    # Ring geometry: outer cylinder minus inner cylinder
    outer = td.Cylinder(center=(0, 0, 0), radius=radius + width / 2, length=thickness, axis=2)
    inner = td.Cylinder(center=(0, 0, 0), radius=radius - width / 2, length=thickness, axis=2)
    ring = td.Structure(geometry=outer - inner, medium=core_medium)

    # Simulation domain — asymmetric: source side open, far side has extra room
    ring_outer = radius + width / 2
    x_min, x_max = -wavelength, ring_outer + 2 * wavelength
    sim_size = (x_max - x_min, 2 * (radius + width + 2 * wavelength), plane_size[1])
    sim_center = ((x_min + x_max) / 2, 0, 0)

    source = td.ModeSource(
        center=(0, -radius, 0),
        size=cross_sec,
        source_time=td.GaussianPulse(freq0=freq0, fwidth=0.1 * freq0),
        mode_spec=mode_spec_src,
        direction="+",
        mode_index=mode_index,
    )

    # Field monitor in the z=0 plane for visualization
    field_xy = td.FieldMonitor(
        center=(0, 0, 0),
        size=(td.inf, td.inf, 0),
        freqs=[freq0],
        name="field_xy",
    )

    # Mode amplitude monitor at the output (top of the ring)
    mode_out = td.ModeMonitor(
        center=(0, radius, 0),
        size=cross_sec,
        freqs=[freq0],
        name="mode_out",
        mode_spec=mode_spec_mon,
    )

    return td.Simulation(
        size=sim_size,
        center=sim_center,
        structures=[ring],
        sources=[source],
        monitors=[field_xy, mode_out],
        medium=background_medium,
        grid_spec=td.GridSpec.auto(min_steps_per_wvl=min_steps_per_wvl),
        run_time=0.6e-12,
        boundary_spec=td.BoundarySpec(
            x=td.Boundary.absorber(num_layers=160),
            y=td.Boundary.pml(),
            z=td.Boundary.pml(),
        ),
    )
fdtd_sim_test = build_fdtd_bend(
    radius=radius,
    width=width,
    thickness=thickness,
    num_modes=num_modes,
    target_neff=target_neff,
    freq0=freq0,
    plane_size=plane_size,
    min_steps_per_wvl=min_steps_per_wvl,
    mode_index=0,  # fundamental quasi-TE bend mode
)

fig, ax = plt.subplots(figsize=(6, 6))
fdtd_sim_test.plot(z=0, ax=ax)
ax.set_title("FDTD simulation — z = 0 cross-section")
plt.tight_layout()
plt.show()

Mode Solver Verification

Before running the FDTD, we use the ModeSolver plugin on the source plane to confirm that the fundamental bend mode is the one being excited.

ms = ModeSolver(
    simulation=fdtd_sim_test,
    mode_spec=fdtd_sim_test.sources[0].mode_spec,
    plane=fdtd_sim_test.sources[0].geometry,
    freqs=fdtd_sim_test.sources[0].source_time.freq0,
)

ms_data = web.run(ms)
23:32:43 UTC Created task 'mode_solver_2026-09-29_23-32-43' with resource_id    
             'mo-4347047f-05b3-4e3f-8988-9067673882cd' and task_type            
             'MODE_SOLVER'.                                                     
             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=mo-4347047f-05b3-
             4e3f-8988-9067673882cd'.                                           
             Task folder: 'default'.                                            
↑ mode_solver.hdf5.gz ━━━━━━━━━━━━━━━━━━━━━━━━ 100.0% • 2.3/2.3 kB • ? • 0:00:00

23:32:45 UTC Estimated FlexCredit cost: 0.013. For this solver type, the        
             estimate is the final billed cost.                                 
23:32:46 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.                                                             
23:32:51 UTC starting up solver                                                 
             running solver                                                     
solver progress ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:00:00

23:33:18 UTC status = success                                                   
23:33:20 UTC View simulation result at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=mo-4347047f-05b3-
             4e3f-8988-9067673882cd'.                                           
↓ mode_solver_data.hdf5.gz ━━━━━━━━━━ 100.0% • 37.1/37.1   • 13.6 MB/s • 0:00:00
                                               MB                               

23:33:24 UTC Loading results from simulation_data.hdf5

Next, we can check the respective TE and TM fundamental modes.

df = ms_data.to_dataframe()
df["Fill factor (%)"] = ms_data.fill_fraction_box.values.squeeze() * 100
df.head(10)
Example library
Previous Tutorial Next Tutorial