Numerical computation of Bragg fiber leaky modes

2.2. Numerical computation of Bragg fiber leaky modes#

The Bragg class computes leaky modes of a Bragg fiber numerically using finite elements and an eigensolver, inheriting from ModeSolver. This is the same eigenproblem introduced in Notebook 2.1, namely either the scalar (Helmholtz) or the vector (Maxwell) leaky-mode eigenproblems with an outgoing radiation condition, but now solved approximately on a truncated computational domain rather than analytically.

We encountered the numerical computation of leaky modes briefly in Notebook 1.3. The same tool we used there can now be used to compute leaky modes of Bragg fibers (and indeed any fiber geometry). The numerical facility is particularly useful when no closed-form solution exists. It gives access to field data on a mesh, enabling direct visualization, error estimation, and confinement loss computations [1].

A key difference from the guided mode case of Notebook 2.1 is that leaky modes do not generally have exponential decay. Hence to truncate the unbounded domain without changing the solution in the region of interest, we now append an outer perfectly matched layer (PML) region, labelled 'Outer' in the material list. The PML transforms the mode in that layer so that it decays to (near) zero at the outer boundary, while leaving the solution untouched in the PML-free interior. The decay in the PML allows us to truncate the eigenproblem to a bounded mesh.

import numpy as np
import ngsolve as ng
from ngsolve.webgui import Draw
import warnings
warnings.filterwarnings('ignore', category=RuntimeWarning)
warnings.filterwarnings('ignore', category=UserWarning)

Constructing Bragg fiber objects#

The Bragg constructor takes the same layer-wise lists as BraggExactScalar and BraggExactVector (see Notebook 2.1), plus mesh-specific parameters, of which two must be highlighted:

  • maxhs — maximum element size per layer, as a fraction of that layer’s outer radius.

  • bcs — boundary condition names. The last entry must be 'OuterCircle' (enforced; this is where the domain is truncated).

The final layer must have mats[-1] = 'Outer' (also enforced) — this tells the mesh builder to treat that region as PML. Its thickness should be large enough for the mode to decay adequately before the outer boundary (and some trial and error may be needed to set this right for each problem).

Scalar leaky mode#

We compute a leaky mode of this fiber at λ = 1.7 µm. The mode is confined to the air core: the field decays rapidly through the glass ring and is negligible in the outer cladding, as we will see.

from fibermode.bragg import Bragg

p = 5  # polynomial degree

B = Bragg(
    ts=[4e-5, 1e-5, 1e-5, 1e-4],   # core | glass ring | air cladding | PML
    mats=['air', 'glass', 'air', 'Outer'],
    ns=[1, 1.44, 1, 1],
    bcs=['r0', 'r1', 'R', 'OuterCircle'],
    maxhs=[0.10, 0.05, 0.2, 0.2],
    wl=1.7e-6,
    curveorder=p,
    scale=1e-6,
)

Draw(B.mesh);

ModeSolver: Checking if mesh has required regions
Mesh has  2554  elements,  1292  points,  and  3845  edges.

We will use our FEAST polynomial eigensolver for computing the Helmholtz leaky modes [2] of this waveguide. Let us recall the non-dimensional complex \(Z\)-plane within which the algorithm searches for eigenvalues. Here \(Z\) is connected to the physical propagation constant \(\beta\) through the relation

\[ Z^2 = L^2\!\left(k_0^2 n_0^2 - \beta^2\right), \]

where \(L\) is the characteristic length (scale), \(k_0 = 2\pi/\lambda\), and \(n_0\) is the outer cladding index. Two conversion methods are available as class methods of Bragg:

  • sqrZfrom(beta) — given a physical \(\beta\) (in rad/m), returns \(Z^2\).

  • betafrom(Z2) — the inverse: given \(Z^2\), returns \(\beta\).

The structure of the spectrum in the \(Z\)-plane is described at the end of Notebook 1.3: guided modes lie on the positive imaginary axis, while leaky modes appear just below the real axis. It is near these leaky-mode locations that we center the FEAST search contour.

We use a frequency-dependent PML [2]. It often gives better convergence in our experience, but requires the additional overhead of solving a polynomial eigenproblem. The leakymode method implements this strategy to solve the frequency-dependent PML eigenproblem for the scalar Helmholtz equation; its key parameters are as follows.

  • p — polynomial degree of the finite elements.

  • ctr, rad — center and radius of the circular search contour in the \(Z\)-plane.

  • alpha — PML absorption strength. Larger values give stronger attenuation but can create unnecessary nonsmoothness in the solution.

  • nspan — dimension of the random initial eigenspace (must be ≥ number of modes sought).

  • npts — number of quadrature points on the FEAST contour.

  • niterations — maximum FEAST iterations.

Z, y, yl, beta, P, info = B.leakymode(p,
    ctr=0.0607,
    rad=0.0001,
    alpha=5,
    seed=1, nspan=4, npts=6, niterations=200, nrestarts=0,
)

ModeSolver.leakymode called on object with these settings:
 <fibermode.bragg.bragg.Bragg object at 0x10b81ea50>
Set freq-dependent PML with p= 5  alpha= 5 and thickness=100.000

SpectralProj: Setting shifted trapezoidal rule quadrature on circular contour
SpectralProj: Radius=0.0001, Center=0.0607+0j
SpectralProjNGPoly: Checking if A[d] is nonzero & Hermitian.
SpectralProjNG: Checking if inverse.H works for umfpack
SpectralProjNGPoly: Computing resolvents along the contour using
                    inverse type umfpack (installed default umfpack)
SpectralProjNGPoly: Factorizing 31996x31996 system at z =  +0.061 +0.000j

SpectralProjNGPoly: Factorizing 31996x31996 system at z =  +0.061 +0.000j

SpectralProjNGPoly: Factorizing 31996x31996 system at z =  +0.061 +0.000j

SpectralProjNGPoly: Factorizing 31996x31996 system at z =  +0.061 -0.000j

SpectralProjNGPoly: Factorizing 31996x31996 system at z =  +0.061 -0.000j

SpectralProjNGPoly: Factorizing 31996x31996 system at z =  +0.061 -0.000j

=========== Starting FEAST iterations ===========
Trying with 4 vectors:

 ITERATION 1 with 4 vectors and 4 left vectors

   During kernel cleaning singular values of B-Gram:
   [2.910e-01 1.646e-16 6.044e-17 3.369e-17]
   Removing ew#[1, 2, 3]
   Real part of computed eigenvalues:
   [0.06042327]
   Imag part of computed eigenvalues:
   [-5.58937774e-05]
   Relative Hausdorff distance from prior iterate: 1.647e+101

 ITERATION 2 with 1 vectors and 1 left vectors

   During kernel cleaning singular values of B-Gram:
   [3.895e-06]
   Real part of computed eigenvalues:
   [0.06042327]
   Imag part of computed eigenvalues:
   [-5.58937774e-05]
   Relative Hausdorff distance from prior iterate: 6.467e-16


***FEAST converged with:
	Relative Hausdorff distance = 6.467e-16

============= FEAST iterations done =============
   Normalizing by M-norms: max=2.8e+02, min=2.8e+02
   Normalizing by M-norms: max=2.7e+00, min=2.7e+00
Results:
 Z: [0.06042327-5.58937774e-05j]
 beta: [3695497.41495736+0.91389179j]
 CL dB/m: [7.93796319]
Mode boundary L² norm = 8.5e-10

For this simple geometry, we can compare against the semi-analytical result from BraggExactScalar (see Notebook 2.1). The exact \(Z\) value for this fiber and wavelength was computed using the transfer-matrix method:

scalar_case_z = Z[0]
exact_scalar_case_z = 0.06069691615738331 - 1.9987492290870344e-05j
exact_beta = B.betafrom(exact_scalar_case_z**2)

print(f'Numerical beta  = {beta[0]}')
print(f'Analytical beta = {exact_beta}')
print(f'Numerical Z     = {scalar_case_z}')
print(f'Analytical Z    = {exact_scalar_case_z}')
print(f'|Z-error|       = {abs(scalar_case_z - exact_scalar_case_z):.3e}')

Numerical beta  = (3695497.4149573636+0.9138917858047603j)
Analytical beta = (3695492.930179003+0.3282861492895721j)
Numerical Z     = (0.060423268719887215-5.5893777406330365e-05j)
Analytical Z    = (0.06069691615738331-1.9987492290870344e-05j)
|Z-error|       = 2.760e-04

The output y is an NGvecs object holding the converged eigenmode. The core-localized mode intensity field is clearly visible below — it is negligible outside the core.

assert y.m == 1, 'Expected only one mode to be found!'

Draw(ng.Norm(y.gridfun())**2, B.mesh);

Vector leaky modes (leakyvecmodes)#

The method leakyvecmodes solves the full Maxwell curl–curl eigenproblem for leaky vector modes. We use the same fiber as in the scalar section above, so the B object and its mesh are already constructed.

The search contour for leakyvecmodes is in the \(Z^2\)-plane (not the \(Z\)-plane). This vector solver uses a standard frequency-independent PML, in contrast to the frequency-dependent PML used for the scalar case. In the scalar solver, \(Z\) enters the problem nonlinearly (as a polynomial eigenvalue), and the FEAST polynomial algorithm [2] finds it directly. In the vector solver, \(Z^2\) appears as an ordinary linear eigenvalue, so no polynomial eigensolver is needed. The practical consequence is that the ctr argument here is a point in the \(Z^2\)-plane, not the \(Z\)-plane.

The search region lies very close to the real axis, where eigenvalues are densely clustered. Such cases are generally difficult for the eigensolver and may require some trial and error to achieve convergence.

betas_num, Zsqrs, Es, phis, R = B.leakyvecmodes(p=p, 
    ctr = scalar_case_z**2, 
    rad=0.0001,
    alpha=5,
    nspan=4, npts=6, seed=1, expected_dim=2, check_contour=10, niterations=100, nrestarts=0)

Set NGSolve automatic PML with p= 5  alpha= 5 and thickness=100.000

Using FEAST to search for vector leaky modes in
circle of radius 0.0001 centered at  (0.003650968278681348-6.75456946398452e-06j)
assuming not more than 4 modes in this interval
System size: 99690  x  99690   Inverse type: umfpack

SpectralProj: Setting shifted trapezoidal rule quadrature on circular contour
SpectralProj: Radius=0.0001, Center=0.00365097-6.75457e-06j
SpectralProjNGR:   Making resolvent at z =  +0.004 +0.000j

SpectralProjNGR:   Making resolvent at z =  +0.004 +0.000j

SpectralProjNGR:   Making resolvent at z =  +0.004 +0.000j

SpectralProjNGR:   Making resolvent at z =  +0.004 -0.000j

SpectralProjNGR:   Making resolvent at z =  +0.004 -0.000j

SpectralProjNGR:   Making resolvent at z =  +0.004 -0.000j

=========== Starting FEAST iterations ===========
Trying with 4 vectors:

 ITERATION 1 with 4 vectors and 4 left vectors

   During kernel cleaning singular values of B-Gram:
   [2.971e-02 3.883e-03 5.334e-09 1.278e-11]
   Removing ew#[3]
   Real part of computed eigenvalues:
   [0.00376649 0.00378879 0.00403383]
   Imag part of computed eigenvalues:
   [-1.05568510e-05 -1.23811578e-05 -1.07711726e-05]
   Relative Hausdorff distance from prior iterate: 2.666e+102

 ITERATION 2 with 3 vectors and 3 left vectors

   During kernel cleaning singular values of B-Gram:
   [1.458e-02 2.394e-04 1.076e-16]
   Removing ew#[2]
   Real part of computed eigenvalues:
   [0.00376649 0.0037888 ]
   Imag part of computed eigenvalues:
   [-1.05568379e-05 -1.23827673e-05]
   Relative Hausdorff distance from prior iterate: 6.533e-02

 ITERATION 3 with 2 vectors and 2 left vectors

   During kernel cleaning singular values of B-Gram:
   [1.152e-01 1.376e-05]
   Real part of computed eigenvalues:
   [0.00376649 0.0037888 ]
   Imag part of computed eigenvalues:
   [-1.05568379e-05 -1.23827688e-05]
   Relative Hausdorff distance from prior iterate: 6.213e-10


***FEAST converged with:
	Relative Hausdorff distance = 6.213e-10

============= FEAST iterations done =============
Results:
 Z²: [0.00376649-1.05568379e-05j 0.0037888 -1.23827688e-05j]
 beta: [3695481.78550748+1.42834392j 3695478.76667676+1.67539439j]
 CL dB/m: [12.40643765 14.55229081]

leakyvecmodes returns:

  • betas_num — physical propagation constants \(\beta\) for each found mode.

  • Zsqrs — the corresponding \(Z^2\) eigenvalues.

  • Es — NGvecs in the HCurl space, holding the transverse electric field \(E_\tau\) of each mode.

  • phis — NGvecs in the H1 space, holding the scaled longitudinal component \(\varphi = i\beta E_z / L\).

The returned numerical approximations to the propagation constants are printed below: they are very close to each other, reflecting the fact that the two found modes approximate the degenerate HE\(_{11}\) pair, which has a single exact propagation constant.

print(f'Numerical betas  = {betas_num}')

Numerical betas  = [3695481.78550748+1.42834392j 3695478.76667676+1.67539439j]

Visualizing the vector leaky mode intensities#

Transverse electric field#

# Intensity of the transverse electric field, with field direction overlaid
for E in Es:
    Draw(ng.Norm(E)**2, B.mesh, vectors={'grid_size': 100});

Longitudinal component#

Although the longitudinal component is much smaller, the next visualization shows tiny ripples in the cladding. Thanks to the studies in [1], we now know that these ripples are not a numerical artifact, but rather a physical feature of the mode that is important to capture to obtain good confinement loss predictions. They will arise again in other notebooks.

# Longitudinal component phi = i*beta*Ez / L
for phi in phis:
    Draw(ng.Norm(phi)**2, B.mesh);

Comparison with exact solution#

We now compute the semi-analytical reference eigenvalue from BraggExactVector (HE\(_{11}\) mode, nu=1) and compare it with the above-found numerical result.

from fibermode.bragg import BraggExactVector
from scipy.optimize import newton

Bexact = BraggExactVector(
    ts=[4e-5, 1e-5, 1e-5],
    ns=[1, 1.44, 1],
    mats=['air', 'glass', 'air'],
    wl=1.7e-6,
    scale=1e-6,
)

nu = 1
outer = 'h1'
guess = B.scale * (3695492.930179003+0.3282861492895721j)
beta2 = newton(Bexact.determinant, guess, args=(nu, outer), tol=1e-15)

print(f'Exact beta    = {beta2 / Bexact.scale}')
print(f'Numerical betas = {betas_num}')
print(f'Relative error  = {abs(betas_num[0] - beta2/Bexact.scale)/abs(beta2/Bexact.scale):.2e}')

Exact beta    = (3695487.7507003495+0.8756591625223714j)
Numerical betas = [3695481.78550748+1.42834392j 3695478.76667676+1.67539439j]
Relative error  = 1.62e-06

References#

[1] P. Vandenberge, J. Gopalakrishnan, and J. Grosek, “Sensitivity of confinement losses in optical fibers to modeling approach,” Optics Express, vol. 31, no. 16, pp. 26735–26760, Jul. 2023. Open Access, DOI: 10.1364/OE.495467

[2] J. Gopalakrishnan, C. Parker, and P. Vandenberge, “Computing leaky modes of optical fibers using a FEAST algorithm for polynomial eigenproblems,” Wave Motion, vol. 101, p. 102826, 2021. DOI: 10.1016/j.wavemoti.2021.102826

Contents