Radial graded-index fibres

A graded-index fibre varies continuously across its core. The radial provider solves selected bound sectors and their exterior continuation. It is appropriate when radial symmetry and the declared sector approximation describe the guide; it is not a general two-dimensional vector eigensolver. See the solid-fibre workflow for propagation.

Radial graded-index providers

PowerLawPermittivity describes a core law epsilon(r)=epsilon_axis*(1-2*Delta*(r/a)**p), with 0<Delta<1/2 and p>=2. The axis material and the exterior material are explicit. Neither material nor radial law takes a z argument. GradedIndexProvider.radius may be a ZProfile; its parameter_paths declares that geometry path separately, under the name core_radius that every other parameterised mode uses, and at_parameters accepts any subset of it.

import numpy as np

from solprop.constants import C
from solprop.modes import GradedIndexProvider, IndexMaterial, PowerLawPermittivity

grin = GradedIndexProvider(
    permittivity=PowerLawPermittivity(IndexMaterial(1.45), 0.005),
    exterior=IndexMaterial(1.45 * np.sqrt(0.99)),
    radius=25e-6,
    family="TE",
    count=2,
    elements=64,
    extent=2.0,
)
rule = grin.node_rule(2 * np.pi * C / 1.03e-6)
fields = grin.solve(2 * np.pi * C / 1.03e-6, 0.0, nodes=rule)
tail = grin.tail_bound(2 * np.pi * C / 1.03e-6, outer_radius=50e-6)

The solver uses continuous linear finite elements for the radial TE₀ or TM₀ Maxwell sector and matches to a decaying Bessel exterior. It solves the resulting frequency-dependent boundary eigenproblem; it does not multiply transfer matrices through thick evanescent layers. elements and extent are explicit numerical choices. Refine element count, matching radius and quadrature independently. The analytic exterior integral certifies the omitted tail of the discretised mode; it does not certify the finite-element error. boundary_residual includes the one-sided Maxwell trace error at the material interface. Linear elements have first-order derivative traces even when their eigenvalues converge at second order.

The same thing holds pointwise inside an element. The strong-form curl residual of a P1 field does not converge at all – it stays at about \(|\beta^2-k_0^2\epsilon|/(k_0^2\epsilon)\), measured 2.06e-3 relative from 32 to 256 elements – because a linear element has no second derivative. Integrated functionals built from the same field do converge at second order (a full-vector Kerr integral refines by factors 3.95 and 3.99 over the same sequence), so pointwise \(E_z\) near the interface carries about 1 % error at 96 elements while the integrals over it are second-order accurate. Refine for anything that reads the longitudinal field at a point rather than under an integral.

A channel above cutoff is not necessarily bound on this domain. If more than MAX_EXTERIOR_TAIL_FRACTION (1e-4) of a channel’s power lies beyond the matching radius, the exact Bessel continuation rather than the discretised core is carrying the mode, and the provider refuses that channel by name. The floor 1e-4 is a judgement that has not been verified: the reference cases reach 6.6e-6, a factor 15 below it, and a channel that should be refused reads 3.2e-4, a factor 3 above, but no reference yet ties the fraction to the error of an observable. A refusal near the floor is a prompt to widen the matching radius. The fraction is computed exactly, from the analytic modified-Bessel tail integral over the same full-domain norm the fields are normalised by; it is written as a power fraction rather than as the dimensionless exterior decay because the decay also depends on the numerical choice of extent. For reference, a five-channel TE case at 48 elements keeps 5.97e-8 of its last channel outside a 25 um core’s matching radius and 3.24e-4 outside a 22.98 um one.

family="LP" selects a weak-guidance scalar approximation, with one transverse linear polarisation and no longitudinal field. It requires relative permittivity contrast at most 0.02. The supplied tests compare its propagation constants with vector TE/TM sectors and the ideal parabolic ladder; this is not a general bound on nonlinear overlaps or longitudinal-field errors. Hybrid HE/EH GRIN is refused by name. TE₀/TM₀ evidence does not qualify it.

A guide from a constructed provider

A channel’s propagation constant is the generator’s diagonal, so the bridge refuses a generator that couples its channels (MAX_GENERATOR_COUPLING): a coupled basis has no per-channel constant, and its flow is a coupled flow (general modal propagation), not a table of indices. The off-diagonal is measured against the smallest splitting between the channels, \(\min|\beta_i-\beta_j|\), and not against the size of \(\beta\): two channels mix over a length set by the ratio of their coupling to their splitting, so the relevant small parameter is the coupling divided by the splitting. For example, a lossless two-channel model with coupling \(\kappa=10^{-3}\) m\(^{-1}\) and splitting \(\Delta\beta=1\) m\(^{-1}\) has maximum power transfer \(4\kappa^2/(\Delta\beta^2+4\kappa^2)\simeq4\times10^{-6}\), obtained by diagonalising its \(2\times2\) propagation matrix. Comparing the same coupling with a carrier \(\beta=10^7\) m\(^{-1}\) would conceal that finite transfer. The splitting is floored at \(64\,\epsilon\,|\beta|_{\max}\) so that a single channel and an exactly degenerate pair still have a scale.

A provider’s generator is in the internal convention (conventions §2), \(G_{cc} = -i(\omega/c)\,\overline{n_{\rm eff}}\), so its imaginary part is \(-\beta\) and its real part is \(-\alpha/2\): a passive channel has \(\mathrm{Re}\,G \le 0\), and the bridge refuses a positive real part as gain by name. An IndexTable is a record in the user-facing optics convention, where \(\mathrm{Im}\,n_{\rm eff} > 0\) is loss, so the bridge writes \(n_{\rm eff} = \overline{i G c/\omega}\). That conjugation is undone when a run builds its linear operator from the table, the same round trip any optical-convention table makes.

ModalGuide.from_provider(provider, omega) pairs a provider that solves its own fields with the ordered modes the rest of the package needs. The provider owns its sector, channel count, elements, extent and radius path, exactly as from_step_index takes constructed modes: nothing in the constructor invents a radial resolution, and there is no second root solver behind a Mode interface.

from solprop.modes import ModalGuide

omega0 = 2 * np.pi * C / 1.03e-6  # the centre of the band the guide will be used over
omega = np.linspace(0.98 * omega0, 1.02 * omega0, 9)
guide = ModalGuide.from_provider(grin, omega)

With no rule supplied the constructor asks the provider for one at the lowest frequency of the band, where a bound field reaches furthest and a truncating rule discards most power, so that the rule’s certified tail_power_bound is a bound over the whole band. Both providers name the frequency the certificate was taken at in the rule’s identity, as tail_power_bound_omega.

degenerate_blocks names the channels whose degeneracy the provider declares — the two parities of an HE_1m channel, say, as ((0, 1),) in the provider’s own channel order. Tracking then rotates freely inside a block and matches one-to-one outside it, and the generator is rotated with the fields, so a degenerate pair stays physically equivalent instead of being matched by an ordering its equal eigenvalues cannot resolve. It is empty by default, because an arbitrary provider declares no degeneracy, and the blocks are kept on the table so that the propagation-constant checks below align a fresh solve the same way.

# strand: a step-index provider, ModalGuide.from_step_index(...).field_provider
pair = ModalGuide.from_provider(strand, omega, degenerate_blocks=((0, 1),))

Each declared channel becomes a ProviderChannelMode: a dispersion carrier whose propagation constant is the provider’s own generator at the knots, read off the tracked diagonal so the channel is followed by physical overlap rather than by eigenvalue order, and whose beta_1 and beta_2 are the interpolant’s, the dense table’s cubic in frequency differentiated, because a provider declares no analytic frequency derivative. The constructor solves the whole band once, into a single sample_modal_table, and reads every channel out of that one table; solprop.translate.modal_tables.sample_provider_table(provider, channel, omega) is the same bridge for one channel on its own, the sibling of sample_index_table for a provider.

The constructor then asks the provider again, at the entrance frequency and midway across the first, the middle and the last gap between knots (the ends of a band are where the cubic interpolant is worst), and refuses a table whose propagation constant disagrees by more than beta_rtol (default PROVIDER_BETA_RTOL, 1e-6). A provider that names no channels cannot be paired at all. A guide whose geometry does not move is a snapshot at one position, and every ProviderChannelMode reports z_constant. A provider whose geometry moves — a radius path that is not constant, or any longitudinal breakpoint — is sampled along its path as well: its z domain is cut into sections at every declared kink, each divided sections= times, and every section is solved at the Chebyshev–Lobatto points of the degree its own held-out check demands, on the provider’s own rule at each position; each channel then carries its index and area along the path. The checks above are taken at the entrance, and a moving guide’s departure along z is held to path_rtol (by default the frequency interpolation’s own measured error, with a floor). Propagating such a guide is described under a guide whose geometry moves along z.

That check certifies beta and nothing else. validate_modal_table measures a profile, a generator and a flux defect; it measures no frequency derivative, and the guide checks none. beta_1 is the interpolant’s: measured held out against five-point differences of fresh solves on the reference parabolic sector, over a ±2 % band, it is 4.4e-13 relative at nine knots. beta_2 is diagnostic and uncertified: on the same sector and the same held-out frequencies it is 2.8e-2 relative at nine knots, 3.2e-3 at seventeen and 2.1e-3 at thirty-three, and worse than unity over a ±20 % band at every knot count measured. A cubic through the knots reproduces a nearly constant group delay to twelve digits and still carries almost no curvature. Whatever reads a second derivative from such a channel — modes.base.beta2 itself, modes.base.zdw, tools.dispersion_length, tools.soliton_order — gets an order of magnitude, not a certified number.

modes.base.beta2 differences beta1 at 0.3 % of the frequency and its Richardson step doubles that, so it reaches 0.6 % either side of the point asked about. A channel table never extrapolates past its knots, so beta_2 nearer than that to the edge of the sampled band is refused by name rather than guessed.

A ProviderChannelMode has no transverse profile: field refuses by name, because the physical fields of the channel are the provider’s three-component E and H on a node rule. Its aeff is the scalar surrogate (int |E|^2 dA)^2 / int |E|^4 dA summed over the three components, which is what the mode-averaged interfaces mean by an area and not the nonlinear coefficient of a vector channel; that is the full Kerr contraction (modal_kerr_coefficient) over the regions that respond.

That surrogate is one number for the whole band, and it is taken at the median knot, the middle of the band the guide was sampled over. The channel table records which frequency that was, as source_identity["channel_area_omega"]. It matters: on the reference parabolic sector the lowest knot’s area is 2.0 % above the median knot’s over a ±2 % band, 25 % over ±20 % and 66 % over ±40 %. A frequency-dependent area does not reach the scalar route at all; the table carries the area at every knot, and a consumer that wants it as a function of frequency reads effective_area(omega) from the table rather than aeff() from the mode.

To have the area at one chosen frequency – a scalar model’s pump – ask the channel itself: solved_area(omega, z) solves the provider again there on the guide’s own rule (for a moving guide, the rule the provider builds at z, as its path was sampled), and area_at(omega) does so at every sampled position of a moving guide’s path and returns the area as a path along z. Every channel a ModalGuide.from_provider builds carries its provider for this (a ChannelSource, left out of the channel’s serialised record). SpecifiedMode.from_mode(channel, wavelength=...) uses it, and prop_gnlse passes its lambda0 (GNLSE page). Without it from_mode falls back to the median-knot surrogate, which on a graded-index guide can be tens of per cent off the area at the pump.