4.2. All-Solid Rod-Lattice PBG Fiber

The bandgap fiber of the previous notebook confines light with a lattice of high-index strands. A bandgap can also be opened by a lattice of high-index solid rods embedded in a lower-index background, with a defect core of the same index as the background. There is no index contrast to guide the core mode at all, unlike ordinary step-index fibers. Confinement is instead an antiresonant effect of the rod lattice, closely related to the ARROW (antiresonant reflecting optical waveguide) mechanism seen for the ARF notebook. This design and its resonance structure are studied in [1], and the same fiber is used as a benchmark example in [2].

The rod dictionary in fibermode.pbg.fiber_dicts describes a single ring of 6 dielectric rods (\(n_\text{tube}=1.8\)) around a solid core of index equal to the cladding (\(n_\text{core}=n_\text{clad}=1.44\)), at a wavelength of 825 nm.

import ngsolve as ng
import numpy as np
from ngsolve.webgui import Draw
from fibermode import PBG
from fibermode.pbg.fiber_dicts.rod import params

Constructing the fiber

A = PBG(params)
Draw(A.mesh);

no buffer no polymer
Generating mesh.

Mesh created.
ModeSolver: Checking if mesh has required regions
Mesh has  14932  elements,  7485  points,  and  22416  edges.

Draw(A.N, A.mesh, 'index', settings={"Objects": {"Wireframe": False, "Edges":False}});
A.n_tube, A.n_clad, A.n_core

(1.8, 1.44, 1.44)

Fundamental mode: scalar search

As with earlier notebooks, we can search for the fundamental mode using the scalar (frequency-dependent PML) formulation. Because \(\Delta n = n_\text{tube}-n_\text{clad}=0.36\) is fairly large here, the scalar Helmholtz problem is only a weakly-guiding approximation – useful for locating roughly where the true (vector) mode sits, but not itself the physically correct eigenvalue. The search center below was located by prior numerical experience with this fiber; see also the tip on obtaining a principled guess in the previous notebook.

Z, y, yl, beta, P, extras = A.leakymode(
    2,
    rad=.01,
    ctr=1.73806037 - 0.01388821j,
    alpha=5,
    npts=2,
    nspan=2,
    niterations=100,
    nrestarts=0,
)

ModeSolver.leakymode called on object with these settings:
 <fibermode.pbg.PBG.PBG object at 0x10d826a50>
Set freq-dependent PML with p= 2  alpha= 5 and thickness=2.000

SpectralProj: Setting shifted trapezoidal rule quadrature on circular contour
SpectralProj: Radius=0.01, Center=1.73806-0.0138882j
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 29901x29901 system at z =  +1.738 -0.004j

SpectralProjNGPoly: Factorizing 29901x29901 system at z =  +1.738 -0.024j

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

 ITERATION 1 with 2 vectors and 2 left vectors
   During kernel cleaning singular values of B-Gram:
   [3.782e+02 5.659e-06]
   Real part of computed eigenvalues:
   [0.0723568 1.7398859]
   Imag part of computed eigenvalues:
   [-0.00705228 -0.01438002]
   Relative Hausdorff distance from prior iterate: 5.720e+99

 ITERATION 2 with 2 vectors and 2 left vectors

   During kernel cleaning singular values of B-Gram:
   [9.406e-01 9.389e-09]
   Removing ew#[1]
   Real part of computed eigenvalues:
   [1.73988605]
   Imag part of computed eigenvalues:
   [-0.01437999]
   Relative Hausdorff distance from prior iterate: 9.539e-01

 ITERATION 3 with 1 vectors and 1 left vectors

   During kernel cleaning singular values of B-Gram:
   [0.941]
   Real part of computed eigenvalues:
   [1.73988605]
   Imag part of computed eigenvalues:
   [-0.01437999]
   Relative Hausdorff distance from prior iterate: 3.577e-16


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

============= FEAST iterations done =============
   Normalizing by M-norms: max=1.6e-02, min=1.6e-02
   Normalizing by M-norms: max=1.5e+02, min=1.5e+02
Results:
 Z: [1.73988605-0.01437999j]
 beta: [10950449.88369618+274.0333015j]
 CL dB/m: [2380.22301397]
Mode boundary L² norm = 7.1e-06
*** Mode boundary L2 norm > 1e-6!
PML decay estimates boundary norm ~ 1.1e-03

We should, in principle increase the PML radius to avoid the last warning, but since this scalar search is only a rough guide to the true vector mode, we ignore it for now.

print('Z    =', Z)
print('beta =', beta)
print('CL [dB/m] =', 20 * beta.imag / np.log(10))

Z    = [1.73988605-0.01437999j]
beta = [10950449.88369618+274.0333015j]
CL [dB/m] = [2380.22301397]

Fundamental mode: vector search

Because the scalar approximation isn’t physically trustworthy for this fiber, we solve the full vector (Maxwell) problem next, the physically correct target, as in [2]. We center the search at \((Z_\text{scalar})^2\) anyway, since that’s close enough for FEAST to converge onto the true vector eigenvalues from there; the resulting \(Z^2\) does not match \((Z_\text{scalar})^2\), which is expected given the large index contrast (unlike the ARF notebook, where scalar and vector formulations of the same near-lossless mode agreed closely).

The fundamental mode here is a near-degenerate pair. The two polarization states of the fundamental hybrid mode are split only slightly by the hexagonal lattice’s imperfect rotational symmetry. Resolving both reliably needs nspan comfortably larger than the true multiplicity of 2 (since with nspan set to exactly 2, FEAST’s subspace-cleaning step can, depending on run-to-run roundoff, collapse the subspace to a single vector partway through, after which it can’t converge cleanly to either eigenvalue). So we set nspan=4 to give it enough slack to separate the two reliably and adjust eta_tol.

center2 = Z**2   # for the Z computed above, X^2 = (1.73806037 - 0.01388821j) ** 2

betas, Zsqrs, Es, phis, R = A.leakyvecmodes(
    p=4,
    rad=0.1,
    ctr=center2,
    alpha=A.alpha,
    npts=2,
    nspan=4,
    niterations=100,
    nrestarts=0,
    eta_tol=1e-15  # need this to prevent accidental truncation of the near-degenerate pair
)

Set NGSolve automatic PML with p= 4  alpha= 5 and thickness=2.000

Using FEAST to search for vector leaky modes in
circle of radius 0.1 centered at  [3.02699668-0.05003908j]
assuming not more than 4 modes in this interval
System size: 410720  x  410720   Inverse type: umfpack

SpectralProj: Setting shifted trapezoidal rule quadrature on circular contour
SpectralProj: Radius=0.1, Center=3.027-0.0500391j
SpectralProjNGR:   Making resolvent at z =  +3.027 +0.050j

SpectralProjNGR:   Making resolvent at z =  +3.027 -0.150j

=========== 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.380e-01 8.300e-02 1.787e-04 9.668e-06]

   Real part of computed eigenvalues:
   [0.67417277 3.13156772 3.1326654  3.16511842]
   Imag part of computed eigenvalues:
   [-4.76207074 -0.08499375 -0.08518905 -2.15726393]
   Relative Hausdorff distance from prior iterate: 3.303e+99

 ITERATION 2 with 4 vectors and 4 left vectors

   During kernel cleaning singular values of B-Gram:
   [1.239e+00 3.813e-01 3.612e-09 1.174e-10]

   Real part of computed eigenvalues:
   [2.66353301 3.13187147 3.13187479 4.36315871]
   Imag part of computed eigenvalues:
   [-0.65854416 -0.08505738 -0.08505704 -3.17561973]
   Relative Hausdorff distance from prior iterate: 1.326e+00

 ITERATION 3 with 4 vectors and 4 left vectors

   During kernel cleaning singular values of B-Gram:
   [1.099e+00 9.877e-02 6.817e-16 3.154e-16]
   Removing ew#[2, 3]

   Real part of computed eigenvalues:
   [3.13187354 3.13187471]
   Imag part of computed eigenvalues:
   [-0.08506113 -0.08505557]
   Relative Hausdorff distance from prior iterate: 1.099e+00

 ITERATION 4 with 2 vectors and 2 left vectors

   During kernel cleaning singular values of B-Gram:
   [0.128 0.012]

   Real part of computed eigenvalues:
   [3.13187354 3.13187472]
   Imag part of computed eigenvalues:
   [-0.08506113 -0.08505558]
   Relative Hausdorff distance from prior iterate: 1.197e-09

 ITERATION 5 with 2 vectors and 2 left vectors

   During kernel cleaning singular values of B-Gram:
   [0.09  0.032]

   Real part of computed eigenvalues:
   [3.13187354 3.13187472]
   Imag part of computed eigenvalues:
   [-0.08506113 -0.08505558]
   Relative Hausdorff distance from prior iterate: 2.416e-12


***FEAST converged with:
	Relative Hausdorff distance = 2.416e-12

============= FEAST iterations done =============
Results:
 Z²: [3.13187354-0.08506113j 3.13187472-0.08505558j]
 beta: [10949875.52890159+465.85203106j 10949875.52247826+465.82160895j]
 CL dB/m: [4046.33932946 4046.07508639]

print('Z^2 =', Zsqrs)
print('CL [dB/m]:', 20 * np.array(betas).imag / np.log(10))

Z^2 = [3.13187354-0.08506113j 3.13187472-0.08505558j]
CL [dB/m]: [4046.33932946 4046.07508639]

leakyvecmodes returns the transverse electric fields Es (in the H(curl) space) and the scaled longitudinal components phis (in H1), one pair per polarization state:

for e in Es:
    Draw(e.real, A.mesh, vectors={'grid_size': 40}, 
         settings={"Objects": {"Wireframe": False}})

for phi in phis:
    Draw(phi, A.mesh, settings={"Objects": {"Wireframe": False}})

Both polarization states show the same antiresonant confinement mechanism: the mode is trapped in the solid core purely by the rod lattice’s bandgap effect, with no index contrast between core and cladding to rely on.

References

[1] N. M. Litchinitser, S. C. Dunn, B. Usner, B. J. Eggleton, T. P. White, R. C. McPhedran, and C. M. de Sterke, “Resonances in microstructured optical waveguides,” Optics Express 11(11), 1243-1251 (2003). DOI: 10.1364/OE.11.001243

[2] J. Gopalakrishnan, J. Grosek, G. Pinochet-Soto, and P. Vandenberge, “Adaptive resolution of fine scales in modes of microstructured optical fibers,” SIAM Journal on Scientific Computing, 47(1):B108-B130, 2025. Open Access, DOI: 10.1137/24M1651605