4.1. Modes of photonic bandgap fibers

Photonic bandgap (PBG) fibers guide light by a mechanism distinct from both total internal reflection of step-index fibers and the antiresonance of ARF. A periodic arrangement of high-index strands in the cladding creates photonic bandgaps, which are ranges of effective index in which transverse propagation through the cladding lattice is forbidden. A defect core (here formed by omitting a few strands) then supports modes whose effective indices fall inside a bandgap. Such modes are prevented from fully radiating through the cladding even though the core index is not higher than its surroundings [Litchinitser et al 2002].

The PBG class builds all-solid bandgap fibers from a parameter dictionary specifying the strand lattice (pitch \(\Lambda\), strand radius, number of layers, indices) together with jacket and PML data. Ready-made dictionaries for several published fibers are provided in fibermode.pbg.fiber_dicts.

To contrast photonic bandgap fibers with the broader term photonic crystal fibers (PCF), note that PCF covers any fiber with a periodic lattice microstructure in the cladding, including ordinary index-guided (“holey”) fibers that merely lower the cladding’s average index and guide by a conventional total-internal-reflection-like mechanism (no bandgap involved). Fibers like those we now consider guide via an actual bandgap effect. However, the PBG class’s geometry machinery can build either kind; see demos/pbg/holey_demo.py for a worked index-guided (non-bandgap) example built with the same class.

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

Constructing PBG objects

The lyr6cr2 dictionary describes a fiber with six rings of high-index strands and a two-ring core defect:

params

{'p': 6,
 'layers': 6,
 'skip': 2,
 'pattern': [],
 'Λ': 7e-06,
 'r_tube': 1.4000000000000001e-06,
 'r_core': 8.819999999999999e-06,
 'r_fiber': 6.65e-05,
 'r_poly': 6.65e-05,
 'r_pml': 6.65e-05,
 'r_out': 9.296e-05,
 'scale': 8.819999999999999e-06,
 'n_buffer': 1.45,
 'n_tube': 1.48,
 'n_clad': 1.45,
 'n_poly': 1.45,
 'n_core': 1.45,
 'n_outer': 1.45,
 'n0': 1.45,
 'wavelength': 1.55e-06,
 't_buffer': 0.0,
 't_poly': 0.0,
 't_outer': 2.6459999999999994e-05,
 'alpha': 5,
 'pml_maxh': 3.7698412698412707,
 'buffer_maxh': 2.2619047619047623,
 'tube_maxh': 0.5277777777777779,
 'clad_maxh': 1.5079365079365084,
 'poly_maxh': 6.031746031746033,
 'core_maxh': 0.07539682539682542}

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

no buffer no polymer
Generating mesh.

Mesh created.
ModeSolver: Checking if mesh has required regions
Mesh has  18772  elements,  9403  points,  and  28174  edges.

The refractive index profile is exposed as the coefficient function N (and the derived potential-like coefficient V used in the eigenproblem):

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

(1.48, 1.45)

Two conveniences worth knowing are A.rotate(angle), which rotates the strand lattice (useful for studying orientation sensitivity), and A.reset_mesh() which restores the original. Here A.N may be modified (e.g., thermally perturbed) and restored by A.reset_N():

# example: a smooth thermal perturbation of the index profile
T = ng.exp(-(ng.x**2 / 30 + ng.y**2 / 30))
A.N = A.N + 0.6 * T
Draw(A.N, A.mesh, 'perturbed index', 
     settings={"Objects": {"Wireframe": False, "Edges":True}})
A.reset_N()

Computing leaky modes

Modes of the defect core are leaky: they can tunnel through the finite number of cladding rings. We search the \(Z\)-plane using the polynomial (frequency-dependent) PML formulation, with a thin elliptical contour around the fundamental mode group, following demos/pbg/pbg_mode.py:

Tip: finding a search center systematically. Rather than guessing ctr/rad by trial and error, A.sqrZfrom(beta) converts a target physical propagation constant into the corresponding nondimensional \(Z^2\), the coordinate FEAST searches in. A principled first guess for a core-guided mode is a \(\beta\) just below the cladding’s own index-guided cutoff, e.g. A.sqrZfrom(.9998 * A.k * A.n_core); a wide-radius search centered there (as below) then locates the actual bandgap-guided modes nearby. (Note this assumes the outer/PML region shares the cladding index, n0 == n_clad. This is not the case for every fiber_dict, e.g. lyr6cr2_w_air uses an air outer region instead.)

z, y, yl, beta, P, extras = A.leakymode(
    3,                     # finite element degree
    ctr=2, rad=1,          # elliptic contour in the Z-plane
    rhoinv=0.8,
    quadrule='ellipse_trapez_shift',
    alpha=A.alpha,
    npts=6, nspan=3,
    niterations=12, nrestarts=0,
    stop_tol=1e-8,
    seed=1,
)

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

SpectralProj: Setting shifted trapezoidal rule on elliptical contour
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 84523x84523 system at z =  +2.866 +0.110j

SpectralProjNGPoly: Factorizing 84523x84523 system at z =  +2.000 +0.220j

SpectralProjNGPoly: Factorizing 84523x84523 system at z =  +1.134 +0.110j

SpectralProjNGPoly: Factorizing 84523x84523 system at z =  +1.134 -0.110j

SpectralProjNGPoly: Factorizing 84523x84523 system at z =  +2.000 -0.220j

SpectralProjNGPoly: Factorizing 84523x84523 system at z =  +2.866 -0.110j

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

 ITERATION 1 with 3 vectors and 3 left vectors

   During kernel cleaning singular values of B-Gram:
   [659.506   2.144   0.914]
   Real part of computed eigenvalues:
   [1.73579773 1.87309005 2.31584699]
   Imag part of computed eigenvalues:
   [ 0.00711599 -0.11393957 -0.63633464]
   Relative Hausdorff distance from prior iterate: 5.000e+99

 ITERATION 2 with 3 vectors and 3 left vectors

   During kernel cleaning singular values of B-Gram:
   [0.722 0.212 0.04 ]
   Real part of computed eigenvalues:
   [1.74020886 2.66486935 2.77688019]
   Imag part of computed eigenvalues:
   [-1.04124833e-04  1.45020415e-01 -7.78472366e-02]
   Relative Hausdorff distance from prior iterate: 4.165e-01

 ITERATION 3 with 3 vectors and 3 left vectors

   During kernel cleaning singular values of B-Gram:
   [0.785 0.723 0.705]
   Real part of computed eigenvalues:
   [1.74011405 2.69998728 2.70927472]
   Imag part of computed eigenvalues:
   [-8.03457602e-08  2.47724630e-04  1.72053409e-04]
   Relative Hausdorff distance from prior iterate: 7.449e-02

 ITERATION 4 with 3 vectors and 3 left vectors

   During kernel cleaning singular values of B-Gram:
   [0.821 0.815 0.725]
   Real part of computed eigenvalues:
   [1.74011163 2.70880034 2.70890195]
   Imag part of computed eigenvalues:
   [-2.15792715e-08  2.92410108e-05 -3.79039011e-06]
   Relative Hausdorff distance from prior iterate: 4.408e-03

 ITERATION 5 with 3 vectors and 3 left vectors

   During kernel cleaning singular values of B-Gram:
   [0.842 0.842 0.725]
   Real part of computed eigenvalues:
   [1.74011176 2.70882881 2.7088367 ]
   Imag part of computed eigenvalues:
   [-6.39957286e-09  1.24972600e-06  7.52832758e-06]
   Relative Hausdorff distance from prior iterate: 3.311e-05

 ITERATION 6 with 3 vectors and 3 left vectors

   During kernel cleaning singular values of B-Gram:
   [0.758 0.758 0.725]
   Real part of computed eigenvalues:
   [1.74011177 2.70882813 2.70882842]
   Imag part of computed eigenvalues:
   [-3.38900897e-09  1.62064634e-07 -1.20363476e-07]
   Relative Hausdorff distance from prior iterate: 5.637e-06

 ITERATION 7 with 3 vectors and 3 left vectors

   During kernel cleaning singular values of B-Gram:
   [0.857 0.857 0.725]
   Real part of computed eigenvalues:
   [1.74011177 2.70882805 2.70882843]
   Imag part of computed eigenvalues:
   [-3.37641531e-09 -1.22646702e-07 -1.23696298e-07]
   Relative Hausdorff distance from prior iterate: 1.481e-07

 ITERATION 8 with 3 vectors and 3 left vectors

   During kernel cleaning singular values of B-Gram:
   [0.889 0.889 0.725]
   Real part of computed eigenvalues:
   [1.74011177 2.70882805 2.70882843]
   Imag part of computed eigenvalues:
   [-3.38240379e-09 -1.23445705e-07 -1.23289539e-07]
   Relative Hausdorff distance from prior iterate: 2.330e-09


***FEAST converged with:
	Relative Hausdorff distance = 2.330e-09

============= FEAST iterations done =============
   Normalizing by M-norms: max=3.8e-01, min=1.8e-09
   Normalizing by M-norms: max=7.2e+01, min=4.5e-05
Results:
 Z: [2.70882843-1.23289539e-07j 2.70882805-1.23445705e-07j
 1.74011177-3.38240379e-09j]
 beta: [5869789.23707818+7.31388370e-04j 5869789.2393264 +7.32314686e-04j
 5874506.49031534+1.28793528e-05j]
 CL dB/m: [0.00635276 0.0063608  0.00011187]
Mode boundary L² norm = 5.3e-09

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

Z    = [2.70882843-1.23289539e-07j 2.70882805-1.23445705e-07j
 1.74011177-3.38240379e-09j]
beta = [5869789.23707818+7.31388370e-04j 5869789.2393264 +7.32314686e-04j
 5874506.49031534+1.28793528e-05j]
CL [dB/m] = [0.00635276 0.0063608  0.00011187]

Draw(ng.Norm(y.gridfun())**2, A.mesh, settings={"Objects": {"Wireframe": False}});

The mode intensity is trapped in the low-index core defect, with the strand lattice acting as the bandgap mirror — visibly different from the index-guided modes of earlier notebooks.

Custom strand patterns

The lattice need not be complete. The pattern entry of the parameter dictionary lists, ring by ring, which strand sites are occupied (1) or vacant (0). Here is an example keeping the two innermost rings complete and opening channels in the third ring:

pattern = [
    [1] * 12,          # ring 1: complete
    [1] * 18,          # ring 2: complete
    [1, 1, 0, 1] * 6,  # ring 3: one vacancy per sector
]
params2 = dict(params)
params2['pattern'] = pattern
params2['layers'] = len(pattern)
B = PBG(params2)
Draw(B.mesh);

no buffer no polymer
Generating mesh.

Mesh created.
ModeSolver: Checking if mesh has required regions
Mesh has  7084  elements,  3556  points,  and  10639  edges.

Such patterned claddings model coupling channels, polarization- breaking designs, or fabrication defects; the mode-solving calls are unchanged.

Other ready-made parameter dictionaries in fibermode.pbg.fiber_dicts include lyr7cr1 (seven rings, one-ring core), holey (air-hole fiber), and rod.

References

• N. M. Litchinitser, A. K. Abeeluck, C. Headley, and B. J. Eggleton. Antiresonant reflecting photonic crystal optical waveguides. Optics Letters, 27:1592-1594, 2002.

• F. Luan et al. All-solid photonic bandgap fiber. Optics Letters, 29:2369-2371, 2004.

• A. Argyros, T. Birks, S. Leon-Saval, C. M. B. Cordeiro, and P. St. J. Russell. Guidance properties of low-contrast photonic bandgap fibres. Optics Express, 13:2503-2511, 2005.