EME convergence for a bent anisotropic waveguide (LiNbOβ)#
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:
Run an FDTD simulation of a 180Β° LiNbOβ waveguide bend to obtain a reference transmittance.
Run EME simulations of the same bend with increasing numbers of cells.
Show that a single EME cell overestimates transmission and explain the physical origin.
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.
[1]:
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 :math:`y` direction (index 2 would place it along \(z\)). The lab-frame tensor used here is therefore
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
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.
[2]:
# ββ 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.)
[3]:
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(),
),
)
[4]:
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.
[5]:
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)
01:03:18 UTC Created task 'mode_solver_2026-07-11_01-03-18' with resource_id 'mo-fa081187-bfa6-47ed-a9b3-c2bac0b3fbcd' and task_type 'MODE_SOLVER'.
View task using web UI at 'https://tidy3d.simulation.cloud/workbench?taskId=mo-fa081187-bfa6- 47ed-a9b3-c2bac0b3fbcd'.
Task folder: 'default'.
01:03:20 UTC Estimated FlexCredit cost: 0.013. For this solver type, the estimate is the final billed cost.
01:03:21 UTC status = success
01:03:24 UTC Loading results from simulation_data.hdf5
Next, we can check the respective TE and TM fundamental modes.
[6]:
df = ms_data.to_dataframe()
df["Fill factor (%)"] = ms_data.fill_fraction_box.values.squeeze() * 100
df.head(10)
