1.2. Numerical computation of guided modes of a step-index fiber#
In this notebook we see how to use finite elements to solve the following Helmholtz eigenproblem for guided modes of a step-index fiber.
The eigenproblem#
Recall from Notebook 1.1 that guided modes are the eigenfunctions arising as solutions to the following eigenproblem: find a nontrivial \(\varphi(x, y)\) on the fiber cross section together with an eigenvalue \(\beta^2\) such that
The time-independent Schrödinger equation from quantum mechanics, \(-\frac{\hbar^2}{2m}\Delta\varphi + V\varphi = E\varphi, \) has an important analogy with that Helmholtz eigenproblem. This is evident when we rewrite the Helmholtz eigenproblem as
with \(Z=r_{\text{core}}\sqrt{k^2n_{\text{clad}}^2-\beta^2}\), which shows that the index well \(\left(k^2n_{\text{clad}}^2-k^2n^2(x,y)\right)\) is like a Schrödinger potential well, and that \(\left({Z}/{r_{core}}\right)^2\) is like a Schrödinger eigenvalue. Note also that this parameter \(Z\) is connected to \(X\) by the relation
Beyond the immediate similarity of the Laplacian, these operators have quite a bit in common; their other two terms are a potential and an eigenvalue. The potential \(V\) in quantum mechanics is the topography of the potential energy landscape experienced by the particle, either from an atomic nucleus or an electric field. This usually entails some form of “potential well” that traps the particles. The resulting spectrum can be decomposed into the bound states (those few eigenstates which are incapable of leaving the well), and unbound states (the infinitude of states that can leave the well). The eigenpairs (eigenfunctions and eigenvalues) are thus the bound/unbound wavefunctions and the corresponding eigenenergies of the system.
The Helmholtz equivalent works similarly. Instead of a wavefunction, we have an electric field, and instead of a potential well, we have a refractive index well. This well uses mechanisms like total internal reflection to trap certain guided modes (bound states) of the waveguide. The remaining modes of the waveguide are radiation modes (unbound states). The eigenpairs of the usual Helmholtz operator are the modes and their propagation constants, in this case we have modes with a nondimensional parameter \(Z\).
Where to search for \(\beta^2\)#
To find these eigenvalues \(Z\) and their corresponding guided modes, we utilize the FEAST eigensolver. The matrix system of our eigenproblem could be very large, which means the number of eigenvalues is also large. Given that we only care about a select few eigenvalues, we need something to filter through these solutions and find only the eigenvalues we want. Utilizing contour integration techniques from quantum mechanics, FEAST iteratively solves the eigenproblem for a selected cluster of propagation constants within a “search contour” of our choosing.
We can choose this contour to contain only the part of the spectrum that we want, the guided propagation constants. While the search contour can be defined over the complex plane, the propagation constants of guided modes are strictly real. We thus only need an interval on the real line for FEAST to search.
From the theory of optical fibers, we know that the propagation constants are bounded by
With the substitution of \(Z\) for \(\beta\), we instead obtain the inequality
We can thus automatically search the interval \((-V^2,0)\) with FEAST to find all of our guided propagation constants and their corresponding guided modes. This interval is our “refractive index well”, in which the bound states are trapped.
An example#
We use the StepIndex class from fibermode.
from fibermode import StepIndex
import warnings
warnings.filterwarnings('ignore', category=RuntimeWarning)
warnings.filterwarnings('ignore', category=UserWarning)
from ngsolve.webgui import Draw
fiber = StepIndex(fibername='Nufern_Yb')
print(fiber)
ModeSolver: Checking if mesh has required regions
Mesh has 414 elements, 224 points, and 637 edges.
StepIndex (ModeSolver) Object:------------------------------------------
Nondimensional Computational Parameters:
Geometry consists of circular core (radius = 1), an annular
cladding layer 1<r<R=8.5, and an outer layer R<r<Rout=16.
Max mesh sizes: 0.3 (core), 3 (cladding), 3 (outer)
Mesh curved by order 3
Physical Parameters:
Wavelength = 1.064e-06 meters
Refractive indices: 1.44973 (cladding), 1.45097 (core)
No curvature set
No temperature set
Included fiber object:
STEP INDEX FIBER ATTRIBUTES: -------------------------------------------
ks: 5.90525e+06 signal wavenumber
wavelength: 1.064e-06 signal wavelength
V: 4.42894 V-number of fiber
ncore: 1.45097 core refractive index
nclad: 1.44973 cladding refractive index
NA: 0.06 numerical aperture
rcore: 1.25e-05 core radius
rclad: 0.0002 cladding radius
------------------------------------------------------------------------
betas, zsqrs, Y = fiber.guidedmodes(p=3, verbose=False, nspan=10)
Running selfadjoint FEAST to capture guided modes in (-19.6154830530724,0)
assuming not more than nspan=10 modes in this interval
Computed approximation of physical propagation constants
betas
array([8566927.4595491 , 8564823.64168799, 8564823.64054636,
8562235.29694509, 8562235.28199795, 8561636.68876595])
Note that there are propagation constants with multiplicity two. These degenerate eigenvalues are related to the mode profiles that have angular variation, like \(LP_{11}\) and \(LP_{21}\). The same propagation constants correspond to the different symmetries of the mode shape, and we denote these variations with subscript \(a\) or \(b\).
Computed non-dimensional Z-squared Schrödinger eigenvalues
zsqrs
array([-15.81869904, -10.18712343, -10.18712037, -3.2604464 ,
-3.26040641, -1.65880736])
Clearly, these \(Z^2\) values are of much more manageable order of magnitude than the large dimensional \(\beta\) values. Nondimensionalization allows FEAST to avoid numerical roundoff errors and compute the eigenvalues more effectively.
Corresponding mode functions
All mode functions are computed together in the Y object. A query reveals how many modes were found:
len(Y)
6
Displaying numerically found modes
for y in Y:
Draw(y)
Comparing numerical modes with exact modes#
The following facility uses StepIndexExact class to get exact modes and attempts to make a correspondence between each computed mode and an exact mode through the LP naming scheme.
n2i, exactbetas = fiber.name2indices(betas)
n2i
{'LP01': 0, 'LP02': 5, 'LP11_a': 2, 'LP11_b': 1, 'LP21_a': 4, 'LP21_b': 3}
exactbetas
array([8566927.47748206, 8564823.6961377 , 8564823.6961377 ,
8562235.54866228, 8562235.54866228, 8561636.86660481])
betas
array([8566927.4595491 , 8564823.64168799, 8564823.64054636,
8562235.29694509, 8562235.28199795, 8561636.68876595])
We can compare the errors with respect to the semi-analytical (exact) eigenvalues, and see that the numerical eigenvalues are in good agreement with the exact values.
(exactbetas-betas) / betas
array([2.09327822e-09, 6.35736446e-09, 6.49065779e-09, 2.93985372e-08,
3.11442423e-08, 2.07715958e-08])