🔬 optics

Practical Vectorial Mode Solver for Dielectric Waveguides Based on Finite Differences

This paper introduces a finite-difference-based numerical solver that formulates a generalized eigenvalue problem using all three electric field components to accurately compute propagation constants and mode profiles in dielectric waveguides while minimizing numerical artifacts at permittivity discontinuities.

Original authors: Ergun Simsek

Published 2026-04-15

📖 4 min read☕ Coffee break read

CC BY 4.0

Original authors: Ergun Simsek

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to design a tiny, invisible highway for light. This highway is a dielectric waveguide—basically a microscopic glass or plastic pipe that traps light and guides it from one place to another, just like a fiber optic cable.

To build these highways, engineers need to know exactly how the light will behave inside them. They need to calculate two main things:

The Speed: How fast does the light travel down the pipe? (In physics terms, this is the "propagation constant").
The Shape: What does the light wave look like as it moves? Does it hug the center, or does it spill over the edges? (This is the "mode profile").

The Problem: The "Spooky" Math

For a long time, scientists used complex math to predict this light behavior. However, when light hits the sharp edges where two different materials meet (like the glass core meeting the air outside), the math gets messy.

Think of it like trying to draw a smooth curve on a pixelated screen. If you aren't careful, the curve looks jagged, or you get "ghost pixels" (called spurious solutions) that don't actually exist in the real world. These errors make it hard to design reliable devices for things like internet cables or medical lasers.

The Solution: A New "Grid" Solver

The author of this paper, Ergun Simsek, has built a new, practical tool (a numerical solver) to fix this. Here is how it works, using simple analogies:

1. The Three-Dimensional Puzzle

Light is an electromagnetic wave. It has three dimensions of movement: up/down, left/right, and forward/backward.

Old way: Some methods tried to solve for just two dimensions and guess the third. It's like trying to bake a cake by only looking at the top and hoping the inside is right. This often led to errors at the "crust" (the boundaries).
This paper's way: The author decided to solve for all three dimensions at the same time. Imagine a 3D jigsaw puzzle where you fit all the pieces together perfectly before you even start. By treating the light's electric field as a complete 3D team, the math stays consistent, and the "ghost pixels" disappear.

2. The Digital Grid (The "Pixel" Map)

To solve this, the computer breaks the waveguide into a tiny grid of squares (like a spreadsheet or a Minecraft world).

The Mesh: The computer looks at every single square in this grid.
The Rules: It uses a clever set of rules (called Finite Differences) to estimate how the light changes from one square to the next.
The Secret Sauce: The author uses a "high-definition" rule. Instead of just looking at the immediate neighbor (a blurry view), the computer looks two squares away in every direction to get a super-sharp, accurate picture of the light's shape.

3. The "Magic" Equation

The author turned the physics of light into a giant math problem called a Generalized Eigenvalue Problem.

Analogy: Imagine a giant drum. If you hit it, it vibrates at specific notes (frequencies). The shape of the drum determines the notes.
In this case, the "drum" is the waveguide, and the "notes" are the specific ways light can travel through it. The computer solves a massive equation to find these "notes" (the modes) and how fast they travel.

Did It Work?

The author tested this new tool on two different types of light highways:

A Rectangular Pipe: A standard glass waveguide.
A Trapezoid Pipe: A waveguide with slanted sides (like a trapezoid shape).

They compared their results with a very expensive, industry-standard software called Tidy3D.

The Result: Their new, simpler tool was almost identical to the expensive one (less than 1% difference).
The Benefit: It's faster, easier to understand, and doesn't require a supercomputer to run. It's like using a reliable, well-calibrated ruler instead of a complex laser scanner for a job that doesn't need nanometer precision.

Why Does This Matter?

This paper provides a practical, reliable, and easy-to-use tool for engineers and researchers.

For Students: It's a clear way to learn how light behaves without getting lost in overly complex math.
For Engineers: It offers a quick way to test designs before building them.
For the Future: Because the code is simple, it can be easily combined with AI and machine learning to design even better optical devices in the future.

In a nutshell: The author built a better, simpler way to simulate how light travels through tiny glass pipes, ensuring that the "map" of the light is accurate and free of digital glitches, making it easier to build the next generation of high-speed internet and laser technology.

1. Problem Statement

The paper addresses the challenge of accurately computing the guided modes (propagation constants and field profiles) of optically linear, non-magnetic dielectric waveguides.

Core Difficulty: Solving the full vector wave equation for the electric field ( $\mathbf{E}$ ) in structures with permittivity discontinuities (material interfaces).
Specific Issues:
- Spurious Solutions: Many numerical methods generate non-physical modes.
- Boundary Condition Enforcement: Standard finite-difference methods often struggle to enforce continuity conditions correctly at material interfaces, leading to numerical artifacts, particularly when calculating the longitudinal field component ( $E_z$ ) from transverse components ( $E_x, E_y$ ).
- Complexity vs. Accuracy: High-order, full-vectorial solvers exist but can be computationally complex and difficult to interpret or integrate with data-driven approaches.

2. Methodology

The author proposes a finite-difference-based numerical solver that formulates the problem as a generalized eigenvalue problem using the electric field formulation.

A. Theoretical Formulation

Governing Equation: The solver starts with the vector wave equation for the electric field in a source-free region:
$\nabla^2 \mathbf{E} + \nabla \left( \frac{1}{\varepsilon} \nabla \varepsilon \cdot \mathbf{E} \right) + k_0^2 \varepsilon_r \mathbf{E} = 0$
Field Decomposition: Assuming a time-harmonic variation $e^{j(\omega t - \beta z)}$ , the field is decomposed into transverse ( $E_x, E_y$ ) and longitudinal ( $E_z$ ) components.
Self-Consistent Approach: Unlike methods that solve for $E_x$ and $E_y$ first and derive $E_z$ via a derivative relation (which introduces errors at interfaces), this method simultaneously solves for all three components ( $E_x, E_y, E_z$ ).
Matrix Construction: The coupled scalar equations are cast into a matrix eigenvalue problem of the form:
$\mathbf{M} \mathbf{E} + \beta \mathbf{L} \mathbf{E} - \beta^2 \mathbf{I} \mathbf{E} = 0$
To linearize this quadratic eigenvalue problem, an auxiliary variable $\mathbf{F} = \beta \mathbf{E}$ is introduced, converting it into a linear generalized eigenvalue problem:
$\begin{pmatrix} \mathbf{M} & \mathbf{L} \\ \mathbf{0} & \mathbf{I} \end{pmatrix} \begin{pmatrix} \mathbf{E} \\ \mathbf{F} \end{pmatrix} = \beta \begin{pmatrix} \mathbf{0} & \mathbf{I} \\ \mathbf{I} & \mathbf{0} \end{pmatrix} \begin{pmatrix} \mathbf{E} \\ \mathbf{F} \end{pmatrix}$
This ensures that boundary conditions and derivatives are handled consistently across the entire field vector.

B. Numerical Implementation

Discretization: A rectangular grid is used where electric field components are evaluated at cell edges and relative permittivity ( $\varepsilon_r$ ) at cell centers.
Derivative Approximation: While first-order differences are mentioned for conceptual clarity, the implementation uses fourth-order central difference stencils for spatial derivatives to improve accuracy:
$\frac{\partial A}{\partial \rho} \approx \frac{-A(\rho+2h) + 8A(\rho+h) - 8A(\rho-h) + A(\rho-2h)}{12h}$
Solver: The resulting large sparse matrix system is solved using MATLAB's eigs function.
Post-Processing:
- Magnetic fields ( $H_x, H_y, H_z$ ) are calculated directly from the electric field derivatives using Maxwell's curl equations (which do not involve $\varepsilon$ derivatives, thus avoiding interface errors).
- Modes are filtered to retain only physically meaningful solutions ( $n_{eff} < n_{cladding}$ ).
- Fields are normalized to a power density of $1 \text{ W/m}^2$ .

3. Key Contributions

Unified Vectorial Formulation: The primary contribution is the inclusion of all three electric field components in a single, self-consistent eigenvalue formulation. This mitigates the numerical inaccuracies associated with calculating $E_z$ post-hoc at permittivity discontinuities.
Practicality and Simplicity: The method avoids the extreme complexity of high-order schemes while maintaining sufficient accuracy for exploratory analysis and integration with data-driven workflows. It offers a "sweet spot" between computational efficiency and physical accuracy.
Robustness at Interfaces: By treating the full vector field simultaneously, the solver effectively reduces spurious solutions and numerical artifacts at material boundaries.

4. Numerical Results

The solver was validated against Tidy3D (a commercial FDTD solver) using two distinct waveguide geometries at a wavelength of 1550 nm.

Case 1: Rectangular Si $_3$ N $_4$ Waveguide
- Geometry: 1.6 $\mu$ m $\times$ 0.7 $\mu$ m core, surrounded by SiO $_2$ .
- Mesh: 20 points per wavelength (PPW).
- Result: The effective refractive indices ( $n_{eff}$ ) for the first four modes showed a difference of less than 0.23% compared to Tidy3D (used as ground truth with 80 PPW).
- Field Profiles: The magnitude distributions of $|E_x|$ , $|E_y|$ , $|E_z|$ , and magnetic fields matched expected physical behavior, with energy concentrated in the core.
Case 2: Trapezoidal Waveguide on Thin Film
- Geometry: A trapezoidal waveguide (width varying from 1.2 $\mu$ m to 1.0 $\mu$ m) mimicking thin-film lithium niobate, on a substrate.
- Mesh: 30 PPW to capture the varying width.
- Result: The difference in $n_{eff}$ compared to Tidy3D was less than 0.88% for all four modes.
- Observation: The solver successfully captured the field variations near the waveguide edges and the symmetric profile of the fundamental mode.

5. Significance

Accessible Tool: The paper provides a "practical" tool that is easier to implement and interpret than state-of-the-art high-order solvers, making it suitable for rapid prototyping and educational purposes.
Data-Driven Integration: The simplified, interpretable nature of the code makes it highly suitable for integration with machine learning or optimization loops where computational speed and code transparency are prioritized over machine-precision accuracy.
Validation of Approach: The results demonstrate that a carefully constructed finite-difference scheme, which respects the vector nature of the field at interfaces, can achieve sub-1% accuracy without the overhead of complex integral-equation or high-order finite-element methods.

In conclusion, this work offers a reliable, efficient, and mathematically sound approach for the modal analysis of dielectric waveguides, effectively balancing the trade-off between numerical complexity and physical accuracy.