Mode Solvers

For standard all-glass optical fibers, which usually exhibit circular symmetry and a low refractive index contrast, one can use the simplest kind of scalar paraxial wave equation and compute so-called LP modes. “LP” means “linearly polarized”, but actually in this context one does not consider the polarization of light at all, regarding this aspect as separable from mode profiles and unrelated to mode properties like effective refractive indices. One considers only one electric field component, neglecting any coupling between such components, as occurs in cases with high index contrast.

One exploits the circular symmetry by writing the mode function as a product of a radial function ($F_{lm}(r)$) and an azimuthal part ($\cos l \varphi$) or ($\exp i l \varphi$). While the azimuthal part is simple to handle, the radial part at least means numerically solving only a single second-order differential equation, which can be transformed into two coupled first-order differential equations.

The radial equation contains a parameter ($\beta$), the initially unknown propagation constant (phase constant). Modes require solutions of the radial equation which converge to 0 for ($r \to \infty$). Such solutions exist only for specific values ($\beta = \beta_{lm}$), where ($l$) is the chosen azimuthal index and ($m$) is a counter index; for other ($\beta$) values, the function rapidly goes to infinity after some number of oscillations. (The number of possible ($m$) values is also initially not known.) So the numerical algorithm must reliably identify all ($\beta$) values for which the mentioned convergence to zero is achieved. It can start the propagation at ($r = 0$), going up to the core radius, i.e., the radius beyond which the refractive index stays constant (the cladding index). As the solutions in the area of the cladding are known to have the form ($K_l(w r)$), where ($w$) depends on ($\beta$), one can consider a boundary condition at the core–cladding interface (also including first-order derivatives) to indirectly check convergence to 0 for ($r \to \infty$). The number of zero crossings of the radial function can be used to ensure that no mode is overlooked.

Although the technical details are altogether certainly not trivial, highly robust and efficient LP mode solvers have been developed which can calculate the details of thousands of guided modes of a large-core multimode fiber within a few seconds, even on a standard PC. Fig. 1 shows an example with a fiber having only a few modes.

For more details, see the article on LP modes.

Only for rough estimates can LP modes be used in the context of index-guiding photonic crystal fibers. In such cases, the effective index method replaces the structured cladding with an average refractive index, treating the fiber as a step-index fiber. This gives reasonable qualitative results but is too approximate for precise properties such as chromatic dispersion or confinement loss.

If the refractive index contrast is not small — as is the case, for example, in photonic crystal fibers with air holes — the scalar wave equation cannot be used. In addition, such air holes break the assumption of radial symmetry — a second reason why LP mode solvers cannot be used. Similar conditions are found for waveguides on photonic integrated circuits, which often exhibit high index contrast (especially when the optical field reaches the interface to air, but also at many semiconductor interfaces) and reduced symmetry.

The polarization of light cannot be neglected in that regime. There can be substantial coupling e.g. between ($E_x$), ($E_y$) and ($E_z$), so that one cannot consider just one transverse field component. There is a possible hybridization between TE-like and TM-like modes where symmetry and separability are lost. The refractive index now depends on both ($x$) and ($y$), the boundary conditions mix components of electric and magnetic field, and as a result there are no longer pure TE or TM solutions; one has non-negligible longitudinal components ($E_z$) and/or ($H_z$). This leads to hybrid modes of type HE (hybrid electric, TE-like with dominantly transverse electric field) and EH (hybrid magnetic, TM-like with dominantly transverse ($H$) field). For example, the HE₁₁ mode in an optical fiber corresponds to the LP₀₁ approximation; it is mostly TE-like, but with nonzero ($E_z$) and ($H_z$). Some higher-order modes like EH₁₂ are more TM-like. Obviously, some understanding of these more complicated types of modes is a prerequisite for using a vectorial mode solver.

Note that a scalar model would not, for example, be able to predict birefringence of a fiber, introduced by an asymmetric design. For that, a vectorial treatment is essential.

So a vectorial mode solver must fully account for all three electric field components ($E_x$), ($E_y$) and ($E_z$). It must respect electromagnetic boundary conditions at material interfaces and correctly handle polarization coupling and hybridization. The goal is to find fields ($E(x,y)$) and propagation constants ($\beta$) (or effective refractive indices ($n_\textrm{eff} = \beta / k_0$)) such that the fields satisfy Maxwell's equations without making the scalar (LP) approximation or paraxial approximation. The starting point is Maxwell's curl equations, often in the frequency domain, i.e., for a given optical frequency ($\nu$) (or angular frequency ($\omega$)):

where ($\epsilon_\textrm{r} = n^2$) is generally a function of the position coordinates ($x$) and ($y$) (and can be a tensor for non-isotropic media). One assumes that there is no ($z$) dependence of the structure (as usual when seeking modes). Guided modes have

Plugging this into Maxwell's equations eliminates the ($z$) derivative (replaced by ($-i\beta$)), leaving two-dimensional coupled equations in ($x$) and ($y$). By eliminating the ($H$) field, one obtains a generalized eigenvalue equation for ($\beta^2$):

​ where ($\nabla_{\perp}$) is the transverse differential operator.

For finding the eigenvalues and eigenvectors (field distributions ($\vec{e}$) and ($\vec{h}$), and finally ($\vec{E}$) and ($\vec{H}$)), one applies some kind of discretization. Different numerical schemes exist:

• Finite-difference (FD) method: One uses a regular ($x-y$) grid, approximates derivatives by central differences, and applies Maxwell's equations on this grid (often using the Yee cell scheme) [9]. That leads to a system of linear equations containing a sparse matrix, to which an efficient algorithm (often based on an iterative method) can be applied.
• This approach is relatively simple to implement, but is not very accurate for curved boundaries or nonrectangular shapes, as encountered e.g. with photonic crystal fibers — while for photonic integrated circuits that solution can be fully satisfactory.
• Finite-element (FE) method: One uses triangular or tetrahedral elements, represents fields inside each element using vector basis functions, and assembles the global matrix enforcing field continuity across elements [6, 10].
• Such an algorithm is difficult to implement, but can well handle curved and complex geometries with high accuracy.
• Fourier / plane wave expansion (PWE) methods: One expands fields and dielectric function as Fourier series (plane wave superpositions). This leads to a matrix eigenvalue equation in the Fourier domain [7]. This works well for periodic photonic crystal structures, but requires extra methods for handling the central “defect” (e.g., a missing air hole). For example, one may assume periodically repeated defects as an approximation, which is reasonable provided that the period is long enough. Convergence can be slow for sharp boundaries or high index contrast.

Different boundary conditions can be applied depending on the physical situation:

• Perfectly Matched Layers (PMLs) absorb fields, suitable for simulating open boundaries for leaky modes.
• Metallic walls are used for confined or symmetric waveguides.
• Symmetry planes can reduce the computational domain for structures with mirror symmetry.

For each computed field, one can calculate additional quantities like the effective mode area, the group index (using solutions for slightly different frequencies) and chromatic dispersion parameters, and figures for polarization (e.g. field amplitude ratios). Propagation losses may also be estimated, e.g. for lossy materials (with complex refractive index) or leaky modes. Note that photonic crystal fiber modes are inherently leaky due to a limited area being covered with air holes, although the leakage loss can be very small.

Substantial technical difficulties can be encountered, and cannot always be fully solved:

• The matrices obtained from discretization can have millions of unknowns. Solving for all eigenvalues would require full diagonalization, which is computationally infeasible. Therefore, practical solvers generally use iterative algorithms (e.g. Arnoldi or Lanczos) that find only the lowest few eigenvalues near a target effective index.
• Guided modes often cluster near certain effective indices (especially near the cladding index). Higher-order or weakly guided modes may be very close together, and numerical searches can then miss some nearly degenerate modes.
• Some structures (e.g. circular fibers) have degenerate mode pairs. Small numerical asymmetries can cause the solver to find only one of the degenerate pair, or to mix them into a single mode. Reliably detecting and distinguishing degenerate or near-degenerate modes requires extra care.
• Leaky modes have complex effective indices, and their fields extend beyond the computational domain. Depending on boundary conditions (PMLs, absorbing layers), these may not appear as true eigenmodes or may appear distorted. Thus, some solvers capture only bound modes, missing leaky or radiation-type solutions.
• Using symmetry planes (to reduce the computational load) can restrict the solver to certain symmetry classes; modes of other symmetry are then not computed.
• Discretization errors or poor mesh resolution can lead to spurious eigenvalues (non-physical solutions) or numerical “pollution” that hides true modes.
• Poorly chosen initial guesses or convergence thresholds can cause missing weakly confined modes, whose fields are mostly outside the core.

Of course, missing modes can have major implications. For example, the simulated modal content, coupling, and losses are incomplete, which affects e.g. coupled-mode simulations. To minimize such risks, a user may need to perform multiple searches with different initial guesses or solver parameters, check convergence as mesh resolution is improved, and use mode overlap integrals or field orthogonality to confirm completeness.

Recent research has begun exploring alternative paradigms beyond classical discretization. These include meshless deep-learning eigenmode solvers [11], physics-informed neural-networks that solve full-vector waveguide modes [13], and physics-aware surrogate neural models for rapid prediction of effective indices [15]. Parallel to these, hybrid analytical/mode-matching solvers have advanced high-index-contrast waveguide modeling [14]. These methods are not yet as widely used as FEM/FD solvers, but they point toward future directions such as real-time mode computation, inverse design integration, and mesh-free geometry handling.

Generally, vectorial mode solvers are not only difficult to make but also slow to operate — although the computational load depends substantially on (a) the quality of used algorithms and their implementation, and (b) on the devices to be modeled. Obviously, it also matters whether one needs solutions only for a specific optical frequency, or scanning of all modes (in multimode regimes) over a substantial frequency range. Fortunately, computers (even personal computers) have become far more powerful in recent decades, and substantial progress has also been made on numerical algorithms. Many utilize special SIMD instructions of modern CPUs which result in much enhanced speed.

There are also simpler versions, called semi-vectorial mode solvers, which utilize additional simplifying approximations for lowering the computational load. They can be appropriate for waveguides exhibiting a moderate refractive index contrast (e.g. <0.3 between different semiconductor materials, but not between solids and air) and are only mildly asymmetric.

Essentially, one treats only the dominant transverse field component (for TE-like or TM-like modes, treated separately), assuming that the coupling between components is weak. Polarization hybridization cannot be investigated. One still accounts for the full refractive index profile, using correction factors for treating discontinuities of refractive index with approximations.

The resulting accuracy can be good if the conditions are reasonably fulfilled, and the CPU time can be substantially lower than for a full vectorial solver.

As mentioned above, new classes of mode solvers have emerged that leverage machine-learning or meshless formulations rather than traditional grid or element discretization. Such novel approaches are still maturing but point toward future trends in photonic simulation.

Even analytical (or semi-analytical) techniques are evolving for high-index-contrast and complex geometries [14], so besides numerical solvers, analytical and hybrid methods remain active research areas.