Getting Waveguide Dispersion Right in MATLAB
I spent about three weeks last year debugging a dispersion simulation where the group velocity dispersion came out negative when it shouldn't have been. Turns out the mesh was too coarse near the core-cladding interface, and the effective index solver was happily converging to garbage. I switched to an adaptive mesh that refined the interface automatically and the numbers fell into place in two runs. That's the kind of thing that eats your time. The basic idea behind waveguide dispersion is that light in a confined structure picks up a wavelength-dependent propagation constant because the mode field distribution changes with wavelength. It's separate from material dispersion, which comes from the refractive index itself varying with wavelength. When you're designing anything from a silicon photonics waveguide to a fiber Bragg grating, you need both pieces, and you need them accurate because they add together.
Waveguide Dispersion Matlab Code
Here's a practical implementation. This code solves for the effective index of the fundamental TE mode in a symmetric planar waveguide using the transcendental eigenvalue equation, then computes the waveguide contribution to group velocity dispersion. clear; clc;
% --- Parameters for a Si-on-insulator strip waveguide (simplified as planar for clarity) ---
lambda = linspace(1.45e-6, 1.60e-6, 150); % wavelengths in meters
n_Si = 3.475; % core refractive index (Si at 1550nm, approximate)
n_SiO2 = 1.444; % cladding (SiO2)
h = 220e-9; % core thickness
% Material dispersion using Sellmeier approximation for Si
% A simple Cauchy form is sufficient here: n(lambda) = A + B/lambda^2
A_Si = 11.75; B_Si = 0.45e-12; % rough fit for Si in IR
n_core = sqrt(A_Si + B_Si ./ lambda.^2);
n_clad = n_SiO2; % silica dispersion is mild in this range
% --- Solve eigenvalue equation for each wavelength ---
neff = zeros(size(lambda));
for i = 1:length(lambda)
V = (2*pi/lambda(i)) * n_core(i) * h;
kappa = sqrt(n_core(i)^2 - n_clad^2);
beta_guess = n_clad * 2*pi/lambda(i);
% Bisection search for beta
beta_min = n_clad * 2*pi/lambda(i);
beta_max = n_core(i) * 2*pi/lambda(i);
for iter = 1:200
beta_mid = (beta_min + beta_max) / 2;
U = h * sqrt(n_core(i)^2 * (2*pi/lambda(i))^2 - beta_mid^2);
W = h * sqrt(beta_mid^2 - n_clad^2 * (2*pi/lambda(i))^2);
% TE mode eigenvalue equation for symmetric slab
lhs = tan(U/2) * (U/W);
rhs = 1;
if lhs > rhs
beta_max = beta_mid;
else
beta_min = beta_mid;
end
end
neff(i) = beta_max / (2*pi/lambda(i));
end
% --- Compute waveguide dispersion ---
% Group index: ng = neff - lambda * d(neff)/d(lambda)
dneff_dlambda = gradient(neff, lambda);
ng = neff - lambda .* dneff_dlambda;
% Group velocity dispersion parameter Dwg (ps/(nm*km))
% D = -(lambda/c) * d^2(neff)/d(lambda)^2
d2neff_dlambda2 = gradient(dneff_dlambda, lambda);
c = 3e8;
Dwg = -(lambda ./ c) .* d2neff_dlambda2 * 1e12; % convert to ps/(nm*km)
% Plot results
figure;
subplot(2,1,1);
plot(lambda*1e9, neff);
xlabel('Wavelength (nm)'); ylabel('Effective Index');
title('Effective Index vs Wavelength');
subplot(2,1,2);
plot(lambda*1e9, Dwg);
xlabel('Wavelength (nm)'); ylabel('D_{wg} [ps/(nm\cdot km)]');
title('Waveguide Dispersion Parameter'); The code above handles the core physics. But here's what most tutorials skip: the bisection method will silently return wrong answers if you pick the wrong initial bracket. For higher-order modes, you need to bracket around the correct root, not just assume the fundamental mode lives near the cladding index. I once had a simulation where the code kept converging to the second TE mode because the eigenvalue crossing wasn't monotonic across my wavelength sweep. The fix was to seed each wavelength step with the previous wavelength's solution, which keeps the solver on the same branch. That single change prevented half a day of confusion.
Another thing to watch out for is the differentiation step. Taking numerical derivatives of effective index versus wavelength amplifies any noise in your solver output. The difference between adjacent points can be tiny, maybe on the order of 1e-5 or less across a 150 nm sweep, so round-off error matters. I usually smooth the neff curve with a Savitzky-Golay filter before differentiating. The sgolayfilt function in MATLAB does this in one line, and using a polynomial order of 3 with a window of about 11 points keeps the dispersion shape intact while killing the jitter. If you need accuracy beyond what a slab model gives you, the same approach works for rectangular strip waveguides, but you'll want to use the effective index method or switch to a full-vectorial solver. The EIM approach is straightforward: solve the vertical slab problem first to get an effective index, then treat that as the core index for a horizontal slab problem. It's fast and usually within a few percent of a full eigenmode solver for aspect ratios up to about 2:1. Beyond that, the error grows because the corner fields aren't captured well. For the full-vectorial case, MATLAB's PDE toolbox can handle it, or you can interface with open-source tools like MPB or Lumerical's FDE solver. If you're doing this repeatedly for optimization, I'd recommend wrapping the slab solver in a function that returns both neff and its derivatives analytically where possible, because analytical derivatives are exact and eliminate the smoothing question entirely.
Get the Full Details
Pitfalls That Will Waste Your Afternoon
Material dispersion model choice matters more than you'd think. The Sellmeier equation for silicon has multiple published parameter sets, and they don't always agree within the 1.4 to 1.6 um range. If your target is a dispersion-flattened waveguide near the zero-dispersion wavelength, even a 0.001 error in n_core can shift that wavelength by several nanometers. I've seen papers where the ZDW position was off by 15 nm between two groups because they used different Si refractive index fits. Check your source. Beware of the thin-oxide substrate trap. On SOI, the bottom cladding is a thick SiO2 layer on silicon, but if your waveguide is etched through the entire silicon film, the mode sees air below. The asymmetry between top (air, n=1) and bottom (SiO2, n=1.44) cladding breaks the symmetry assumption in the eigenvalue equation. You need to use the asymmetric slab equation, which replaces the simple tan(U/2) term with a more complex expression involving the phase shifts at both interfaces. I learned this the hard way when my simulated ZDW didn't match measurements from a foundry PDK, and it turned out I'd been using the symmetric formulation for a fundamentally asymmetric structure. The asymmetric TE eigenvalue equation uses tan(U) = (U*(W1+W2))/(U^2-W1*W2) where W1 and W2 are the decay constants in the top and bottom claddings respectively. It's not much harder to code, and it makes a real difference when the two cladding indices are that different.
One final note on performance: if you're sweeping parameters—core width, height, wavelength—nested loops in MATLAB will be slow. Vectorizing the eigenvalue solve across wavelengths using fzero with vectorized function handles doesn't work directly because fzero expects scalar brackets. The workaround is to precompute an initial guess table using a coarse wavelength grid, then interpolate and refine. This cuts a sweep over 500 wavelengths from about 45 seconds to roughly 3 seconds on a standard laptop. Worth the extra setup if you're doing parameter scans.
