Eigenmode Sensitivity in Step-Index Fibers: A Sturm–Liouville Approach to Fabrication Tolerances

Optical fibers and integrated photonic waveguides are fundamental to modern telecommunications and sensing systems. Their performance depends critically on fabrication precision—small deviations in refractive index, core radius, or material purity can shift propagation characteristics and degrade system performance. Understanding these sensitivities is essential for manufacturing optimization, robust device design, and yield analysis.

This post presents a rigorous analysis of the fundamental LP₀₁ mode in step-index optical fibers, combining Sturm–Liouville operator theory with first-order perturbation methods and validated numerical solvers. We derive analytical formulas for eigenvalue sensitivity and validate them against exact numerical recomputation across realistic fabrication parameter ranges.

This article is an entry point to the full implementation; the complete code, data, and figures are available in the GitHub repository: mmaharebi/fiber-perturbation.

1. Physical Model and Eigenmode Structure

Step-Index Fiber Geometry

A step-index fiber consists of a cylindrical core (radius $a$, refractive index $n_1$) surrounded by an infinite cladding (refractive index $n_2 < n_1$). Under the weakly guiding approximation ($n_1 \approx n_2$), the transverse field envelope $\Psi(r)$ satisfies the scalar Helmholtz equation:

$$\frac{1}{r}\frac{d}{dr}\left( r\frac{d\Psi}{dr} \right) + \left( k_0^2 n^2(r) - \beta^2 \right) \Psi(r) = 0$$

where $k_0 = 2\pi/\lambda_0$ is the free-space wavenumber, $\beta$ is the propagation constant, and the index profile is:

$$n^2(r) = \begin{cases} n_1^2, & r < a \\ n_2^2, & r > a \end{cases}$$

LP₀₁ Mode Solutions

For the fundamental azimuthally symmetric mode (LP₀₁), the field profile consists of Bessel functions in each region:

Core ($r < a$):

$$\Psi_{\text{core}}(r) = A \, J_0\left( \frac{ur}{a} \right), \quad u = a\sqrt{k_0^2 n_1^2 - \beta^2}$$

Cladding ($r > a$):

$$\Psi_{\text{cladding}}(r) = C \, K_0\left( \frac{wr}{a} \right), \quad w = a\sqrt{\beta^2 - k_0^2 n_2^2}$$

where $J_0$ is the Bessel function of the first kind (regular at the origin) and $K_0$ is the modified Bessel function of the second kind (exponentially decaying).

Dispersion Relation

Continuity of $\Psi$ and $d\Psi/dr$ at the core-cladding interface ($r = a$) yields the eigenvalue equation:

$$\frac{J_1(u)}{u J_0(u)} + \frac{K_1(w)}{w K_0(w)} = 0$$

The parameters $u$ and $w$ are constrained by the normalized frequency (V-number):

$$V = \sqrt{u^2 + w^2} = k_0 a \sqrt{n_1^2 - n_2^2}$$

For a given wavelength $\lambda_0$ and fiber parameters $(a, n_1, n_2)$, the dispersion relation is solved numerically using Brent's root-finding method to obtain $u$, from which the propagation constant follows:

$$\beta = \sqrt{k_0^2 n_1^2 - \left(\frac{u}{a}\right)^2}$$

The effective index is defined as:

$$n_{\text{eff}} = \frac{\beta}{k_0}$$

which satisfies the physical bound $n_2 < n_{\text{eff}} < n_1$ for guided modes.

Single-mode regime: The fiber supports only the LP₀₁ mode when $V < 2.405$. This criterion is fundamental for telecommunications applications requiring dispersion-free transmission.

2. Sturm–Liouville Operator Formulation

Self-Adjoint Form

The Helmholtz equation can be recast as a Sturm–Liouville eigenvalue problem:

$$\mathcal{L}\Psi = -\frac{1}{r}\frac{d}{dr}\left( r\frac{d\Psi}{dr} \right) + \frac{m^2}{r^2}\Psi + k_0^2 n^2(r)\Psi = \beta^2 \Psi$$

where $\mathcal{L}$ is a self-adjoint operator on the Hilbert space $L^2([0,\infty), r\,dr)$ with inner product:

$$\langle \Psi_1, \Psi_2 \rangle = \int_0^\infty \Psi_1^*(r) \Psi_2(r) \, 2\pi r \, dr$$

Self-adjointness ensures:

• Real eigenvalues: $\beta^2 \in \mathbb{R}$ (physical propagation constants).
• Orthogonality: Modes with different $\beta$ are orthogonal under the weighted inner product.
• Completeness: Guided and radiation modes form a complete basis for field expansions.

Mode Normalization

We adopt the unit-power normalization:

$$\int_0^\infty |\Psi(r)|^2 \, 2\pi r \, dr = 1$$

This convention simplifies perturbation integrals and ensures dimensional consistency.

3. First-Order Perturbation Theory

Perturbation of Refractive Index

Consider a small perturbation to the index profile:

$$n^2(r) \to n^2(r) + \Delta n^2(r), \quad |\Delta n^2(r)| \ll n^2(r)$$

The perturbed operator is $\mathcal{L} + \Delta\mathcal{L}$, where:

$$\Delta\mathcal{L} = k_0^2 \Delta n^2(r) \cdot \text{(multiplication operator)}$$

Using first-order perturbation theory for self-adjoint operators, the eigenvalue shift is:

$$\Delta(\beta^2) = \langle \Psi | \Delta\mathcal{L} | \Psi \rangle = k_0^2 \int_0^\infty \Delta n^2(r) |\Psi(r)|^2 \, 2\pi r \, dr$$

where $\Psi(r)$ is the unperturbed normalized eigenfunction.

Since $\beta^2 = (k_0 n_{\text{eff}})^2$, and $\Delta(\beta^2) = 2\beta \Delta\beta$ for small perturbations:

$$\Delta\beta = \frac{k_0^2}{2\beta} \int_0^\infty \Delta n^2(r) |\Psi(r)|^2 \, 2\pi r \, dr$$

Effective Index Shift

The shift in effective index is:

$$\boxed{\Delta n_{\text{eff}} = \frac{1}{2n_{\text{eff}}} \int_0^\infty \Delta n^2(r) |\Psi(r)|^2 \, 2\pi r \, dr}$$

This formula reveals key physical insights:

• Linear dependence: $\Delta n_{\text{eff}}$ scales linearly with $\Delta n^2(r)$, valid only for small perturbations.
• Spatial weighting: Perturbations are weighted by the mode intensity $|\Psi(r)|^2$. Since most power is confined to the core, core perturbations dominate.
• Inverse scaling: Lower $n_{\text{eff}}$ modes are more sensitive to perturbations.

4. Perturbation Scenarios

Localized Index Perturbation (Doping Diffusion)

A common fabrication imperfection is dopant diffusion near the core-cladding interface, creating a thin shell of perturbed index:

$$\Delta n^2(r) = \begin{cases} \Delta n_0^2, & a - \delta < r < a + \delta \\ 0, & \text{otherwise} \end{cases}$$

For $\delta \ll a$, the mode profile varies slowly across the shell, so $|\Psi(r)|^2 \approx |\Psi(a)|^2$. The integral simplifies:

$$\Delta n_{\text{eff}} \approx \frac{2\pi a \delta \, \Delta n_0^2 |\Psi(a)|^2}{n_{\text{eff}}}$$

This shows that sensitivity is proportional to the shell thickness $\delta$, the index change $\Delta n_0^2$, and the mode intensity at the interface $|\Psi(a)|^2$.

Core Radius Variation

A change in core radius $a \to a + \Delta a$ affects the V-number and indirectly shifts $\beta$. This can be modeled as an equivalent index perturbation or computed via numerical differentiation:

$$\Delta n_{\text{eff}} \approx \frac{d n_{\text{eff}}}{d a} \Delta a$$

The derivative $dn_{\text{eff}}/da$ is computed from the dispersion relation and quantifies geometric sensitivity.

Weak Material Absorption

Absorption introduces a complex refractive index $n \to n + i\kappa$, where $\kappa$ is the extinction coefficient. For weak absorption ($\kappa \ll n$):

$$\Delta n^2 = (n + i\kappa)^2 - n^2 \approx 2in\kappa$$

The resulting attenuation coefficient (dB/km) is:

$$\alpha_{\text{dB}} = \frac{40\pi}{\lambda_0 \ln 10} \cdot \text{Im}(\Delta n_{\text{eff}})$$

This is critical for fiber loss budgets in telecommunications (target: $\alpha < 0.2$ dB/km at 1550 nm).

5. Numerical Implementation and Validation

Solver Pipeline

The numerical framework consists of:

1. Dispersion solver: Brent root-finding with automatic bracket detection to solve the transcendental eigenvalue equation. Achieves $10^{-12}$ convergence residual.

2. Mode normalization: Compute $\Psi(r)$ and normalize to unit power via quadrature integration.

3. Perturbation integrals: Evaluate $\Delta n_{\text{eff}}$ via adaptive Simpson quadrature for each perturbation scenario.

4. Validation: Compare first-order predictions against exact numerical recomputation (resolve dispersion relation for perturbed fiber).

Validation Results

For radius perturbations $\Delta a \in [-0.2, +0.2]\,\mu$m (±5% of core radius):

• Mean relative error: 1.76%
• Max relative error: 4.3%
• Validity: Excellent agreement for $|\Delta a| < 0.15\,\mu$m.

The perturbation theory breaks down when $\Delta a/a > 0.05$, requiring second-order corrections or direct numerical solution.

6. Applications and Relevance

Photonics Manufacturing

The framework directly applies to:

• Fiber fabrication yield: Predict performance degradation from dopant diffusion or geometric tolerances.
• Integrated waveguide design: Silicon nitride and silicon photonic waveguides exhibit similar sensitivity to process variations (sidewall roughness, etch depth).
• Sensor design: Evanescent-field sensors rely on precise mode overlap with the sensing region; perturbation theory quantifies detection sensitivity.

Microwave Dielectric Waveguides

The Sturm–Liouville formulation is mathematically identical to dielectric rod waveguides at microwave frequencies. The analysis extends to:

• Millimeter-wave antennas
• Substrate-integrated waveguides (SIW)
• Metamaterial resonators

Connection to Broader Estimation Theory

The weighted integral form of $\Delta n_{\text{eff}}$ is analogous to linear functionals in signal estimation. Just as the Wiener filter optimally weights signal contributions based on their power spectral density, perturbation theory weights index changes by the mode intensity distribution—a spatial analogue of matched filtering.

7. Computational Reproducibility

All results are generated from a validated Python implementation:

• 9/9 test cases passing: Dispersion monotonicity, V-scaling, normalization, physical bounds, perturbation accuracy.
• Data provenance: Dispersion curves, mode profiles, and sensitivity maps are stored as JSON datasets and regenerated via run_all.sh.
• Figure generation: Publication-quality PDFs produced from numerical data.

The framework is extensible to graded-index fibers, planar waveguides, and higher-order modes with minimal modification.

8. Conclusion

This work demonstrates how rigorous operator-theoretic methods combined with validated numerical solvers provide deep insight into fabrication sensitivities in optical fibers. The Sturm–Liouville formulation justifies first-order perturbation theory, while numerical validation establishes accuracy envelopes for practical parameter ranges.

Key takeaways:

• Mathematical rigor: Self-adjoint operators ensure real eigenvalues, orthogonality, and valid perturbation expansions.
• Physical insight: Mode intensity weighting reveals which spatial regions dominate sensitivity.
• Numerical validation: First-order theory is accurate within ±5% for typical fabrication tolerances.
• Broad applicability: The framework extends to integrated photonics, microwave waveguides, and sensing applications.

Understanding these fundamental principles is essential for designing robust photonic systems and predicting performance under realistic manufacturing conditions. The intersection of electromagnetic theory, operator methods, and numerical analysis provides a powerful toolkit for modern waveguide engineering.

References

Optical Fiber Theory:

• Snyder, A. W., & Love, J. D. (1983). Optical Waveguide Theory. Chapman and Hall.
• Marcuse, D. (1991). Theory of Dielectric Optical Waveguides (2nd ed.). Academic Press.
• Okamoto, K. (2006). Fundamentals of Optical Waveguides (2nd ed.). Academic Press.

Mathematical Methods:

• Zettl, A. (2005). Sturm-Liouville Theory. American Mathematical Society.
• Yariv, A., & Yeh, P. (2007). Photonics: Optical Electronics in Modern Communications (6th ed.). Oxford University Press.

Perturbation Theory:

• Marcuse, D. (1972). "Mode conversion caused by surface imperfections of a dielectric slab waveguide." Bell System Technical Journal, 48(10), 3187–3215.
• Snyder, A. W. (1970). "Asymptotic expressions for eigenfunctions and eigenvalues of a dielectric or optical waveguide." IEEE Transactions on Microwave Theory and Techniques, 17(12), 1130–1138.