Skip to content

Eigenmode expansion (EME)

Mode-matching propagation through a sequence of uniform sections. Use the convergence helpers to assess the retained modes and interfaces.

Run

run_eme

run_eme(
    sections: Sequence[Section],
    n_modes: Optional[int] = None,
    *,
    interface_rcond: float = 1e-10,
    enforce_passivity: bool = False,
    interface: str = "galerkin",
) -> EMEResult

Cascade a sequence of Sections into one device S-matrix.

Parameters

sections: Cross-sections in propagation order. The first and last define the input/output ports (typically length-0 leads). All must share the same transverse grid (ny, nx) and spacing. n_modes: Optional cap applied independently to every section. Adjacent sections may retain different counts; their interface blocks are rectangular. A cap may not split a degenerate set of modes (a validated Yee beta group, or the TE/TM pair of a square or round core). enforce_passivity: Clip the gain of each ordinary real guided interface (see interface_smatrix). Off by default: the clip biases the fundamental transmission low by 0.1 to 0.4 %, and more modes or slabs do not remove the bias. Without it, a cascade of ordinary guided ports whose energy_balance(0) exceeds 1 by more than 1e-3 warns: that gain is a basis that has not converged.

interface

"galerkin" (default) or "projection", the interface model of interface_smatrix. Use "projection" for an open device modelled with a truncated (for example guided-only) basis, where power that a section's modes cannot represent radiates away: gentle bends, lateral offsets and other junctions whose scattered field does not return to the guide. Keep "galerkin" with a complete basis; the two agree once the basis is complete, which is a useful check.

Returns

EMEResult The device scattering matrix and convenience power ratios.

run_eme_band

run_eme_band(
    base_sections: Sequence[Tuple[VectorModeSolver, float]],
    wlens_um: Sequence[float],
    *,
    num_modes: int,
    neff_margin: float = 0.0,
    n_clad: Optional[float] = None,
    n_modes: Optional[int] = None,
    solve_kwargs: Optional[dict] = None,
    basis_kwargs: Optional[dict] = None,
    workers: Optional[int] = None,
    backend: str = "auto",
    enforce_passivity: bool = False,
    interface: str = "galerkin",
) -> Dict[float, EMEResult]

Cascade a staircased device across a band of wavelengths, solving the independent per-(wavelength × section) FDE eigenmodes in parallel.

Parameters

base_sections: Cross-sections in propagation order as (base_solver, length_um) pairs, e.g. from rectangular_base_section. Each cross-section is sampled on the mesh once and reused across the whole band. length_um == 0 marks a port lead (as in Section). All sections must share one transverse grid. wlens_um: The free-space wavelengths (microns) to sweep. num_modes: Guided modes solved per cross-section. Passed to VectorModeSolver.solve on the ordinary path, or used as num_guided when basis_kwargs selects a complete EME basis. neff_margin: Drop modes whose n_eff is within this margin of n_clad (near-cutoff, see waveguide_section). Requires n_clad. n_clad: Cladding index, needed only when neff_margin > 0. n_modes: Cap on the cascade basis size, forwarded to run_eme. solve_kwargs: Extra keyword arguments forwarded to VectorModeSolver.solve (e.g. n_guess); must be picklable for the process backend. basis_kwargs: When provided, solve a complete basis with VectorModeSolver.solve_eme_basis instead. num_modes supplies num_guided and this dictionary supplies options such as num_radiation, num_evanescent and pml_cells_xy. It is mutually exclusive with solve_kwargs and must be picklable for the process backend. workers: Worker count for the eigensolves. None → min(#solves, os.cpu_count). 1 (or a single solve) runs serially with no pool. backend: "auto" (default), "process" on fork platforms (Linux/CI/cloud, ~3.5x), "thread" on spawn platforms (macOS/Windows, ~1.7x without the process-startup tax or __main__-guard requirement). Force "process" for a large spawn-platform sweep, or "serial" to disable parallelism. enforce_passivity: Forwarded to run_eme (off by default). interface: Forwarded to run_eme ("galerkin" by default).

Returns

dict[float, EMEResult] {wavelength_um: EMEResult}, one cascaded device S-matrix per wavelength, identical to a per-wavelength run_eme cascade.

EMEResult dataclass

EMEResult(
    s11: ndarray,
    s12: ndarray,
    s21: ndarray,
    s22: ndarray,
    n_modes: int,
    n_left_modes: int = 0,
    n_right_modes: int = 0,
    section_mode_counts: Tuple[int, ...] = (),
    left_modes: Tuple[VectorMode, ...] = (),
    right_modes: Tuple[VectorMode, ...] = (),
    left_power_metric: Optional[ndarray] = None,
    right_power_metric: Optional[ndarray] = None,
    interface_diagnostics: Tuple[
        InterfaceDiagnostics, ...
    ] = (),
    interface: str = "galerkin",
)

Device scattering matrix from run_eme.

Blocks may be rectangular: S21 is N_right x N_left. Physical power is evaluated with Hermitian Poynting metrics, so evanescent/PML amplitudes are never mistaken for |a|^2 power channels. n_modes retains its legacy meaning as the minimum section-basis count; use n_left_modes, n_right_modes, and section_mode_counts for the rectangular model.

s11 instance-attribute

s11: ndarray

s12 instance-attribute

s12: ndarray

s21 instance-attribute

s21: ndarray

s22 instance-attribute

s22: ndarray

n_modes instance-attribute

n_modes: int

n_left_modes class-attribute instance-attribute

n_left_modes: int = 0

n_right_modes class-attribute instance-attribute

n_right_modes: int = 0

section_mode_counts class-attribute instance-attribute

section_mode_counts: Tuple[int, ...] = ()

left_modes class-attribute instance-attribute

left_modes: Tuple[VectorMode, ...] = ()

right_modes class-attribute instance-attribute

right_modes: Tuple[VectorMode, ...] = ()

left_power_metric class-attribute instance-attribute

left_power_metric: Optional[ndarray] = None

right_power_metric class-attribute instance-attribute

right_power_metric: Optional[ndarray] = None

interface_diagnostics class-attribute instance-attribute

interface_diagnostics: Tuple[InterfaceDiagnostics, ...] = ()

interface class-attribute instance-attribute

interface: str = 'galerkin'

transmission property

transmission: float

Fundamental-to-fundamental physical power transmission.

Index-based, so it follows n_eff ordering and NOT polarization, see _fundamental_index. Warns when the two ports disagree on the polarization at index 0. Prefer transmission_of (result.transmission_of("TE")) whenever the ports can reorder, which is the normal case for a taper.

reflection property

reflection: float

Fundamental-to-fundamental physical power reflection.

transmission_of

transmission_of(
    polarization: Optional[str] = None,
) -> float

Power transmission between the fundamental polarization modes of the two ports, resolved by mode identity rather than by index.

Metric-aware like transmission (physical power ratio); polarization=None reproduces it exactly.

transmitted_power

transmitted_power(input_mode: int = 0) -> float

Total physical output-port power from input_mode.

reflected_power

reflected_power(input_mode: int = 0) -> float

Total physical reflected power from input_mode.

energy_balance

energy_balance(input_mode: int = 0) -> float

T_total + R_total for input_mode.

Unlike the former forced-unitary interface, this is a physical flux diagnostic: radiation/PML attenuation or an incomplete trace basis may produce a deficit. A value above one by more than about 1e-3 is a passivity/convergence failure (run_eme warns); below that it is the residual of guided traces that are not exactly orthogonal.

Sections

Section dataclass

Section(
    modes: Sequence[VectorMode], length_um: float = 0.0
)

One z-invariant cross-section of an EME device.

Attributes

modes: The eigenmodes of this cross-section. A guided-only tuple from VectorModeSolver.solve remains valid; a complete EME basis from VectorModeSolver.solve_eme_basis may additionally contain experimental radiation and evanescent modes. length_um: Physical length of the section along propagation (microns). 0.0 marks a port / semi-infinite lead, it contributes its interface with its neighbour but no propagation phase.

modes instance-attribute

modes: Sequence[VectorMode]

length_um class-attribute instance-attribute

length_um: float = 0.0

waveguide_section

waveguide_section(
    *,
    wlen_um: float,
    dl_um: float,
    core_w_um: float,
    core_h_um: float,
    n_core: float,
    n_clad: float,
    window_w_um: float,
    window_h_um: float,
    num_modes: int,
    num_radiation_modes: int = 0,
    num_evanescent_modes: int = 0,
    pml_cells_xy: Optional[Tuple[int, int]] = None,
    pml_imaginary_thickness_um: Tuple[float, float] = (
        0.1,
        0.1,
    ),
    length_um: float = 0.0,
    neff_margin: float = 0.0,
    subpixel: bool = True,
    subpixel_method: str = "tensor",
) -> Section

Solve a centered rectangular-core cross-section and wrap it as a Section.

Pass the same wlen_um, dl_um, window_w_um and window_h_um for every section of a device so they share one transverse grid (required by run_eme). Only core_w_um / core_h_um should vary between sections.

neff_margin drops modes whose n_eff is within this margin of the cladding index, i.e. near-cutoff modes. Those are poorly resolved (large evanescent tails into the window walls) and break the within-section orthonormality the interface relies on, so a small margin (e.g. 0.05) keeps the modal basis clean. Default 0.0 keeps every guided mode the solver returns.

rectangular_base_section

rectangular_base_section(
    *,
    core_w_um: float,
    length_um: float = 0.0,
    dl_um: float,
    core_h_um: float,
    n_core: float,
    n_clad: float,
    window_w_um: float,
    window_h_um: float,
    ref_wlen_um: float = 1.31,
    subpixel: bool = True,
    subpixel_method: str = "tensor",
) -> Tuple[VectorModeSolver, float]

Sample a centered rectangular-core cross-section once and pair it with its length, ready for a band/length sweep via run_eme_band.

Mirrors the geometry arguments of waveguide_section, but returns the unsolved base solver instead of solved modes: the permittivity sampled grid is wavelength-independent (non-dispersive prototype), so ref_wlen_um is only a placeholder, run_eme_band rebases every solve to the swept wavelength with VectorModeSolver.at_wavelength. Pass the same dl_um / window_* for every section of a device so they share one transverse grid (required by run_eme).

Scattering algebra

interface_smatrix

interface_smatrix(
    left_modes: Sequence[VectorMode],
    right_modes: Sequence[VectorMode],
    dl_x_um: float,
    dl_y_um: float,
    *,
    rcond: float = 1e-10,
    enforce_passivity: bool = False,
    return_diagnostics: bool = False,
    method: str = "galerkin",
) -> Union[SMatrix, Tuple[SMatrix, InterfaceDiagnostics]]

Generalized two-sided Petrov-Galerkin match at a waveguide step.

Independent E- and H-continuity projections are tested from both sides with reaction-adjoint modal traces. The resulting tall system supports unequal mode counts and complex/PML bases. Uniform node-collocated and validated full-box uniform Yee-staggered bases are both supported, provided every mode on both sides uses the same registration. (On the Yee lattice the two reaction-product pairs Ex*Hy and Ey*Hx are natively co-located.) A rank-revealing SVD exposes truncation through continuity residuals; it never fabricates exact unitarity.

method chooses how the scattering matrix is built from those projections:

  • "galerkin" (default) solves the tall continuity system in the least-squares sense. With a complete basis this is the exact interface. With a truncated basis it cannot place power in channels the basis lacks, so it returns that power to the represented modes: a guided-only basis reports almost no loss at a junction whose modes are offset or reshaped.
  • "projection" matches each incidence only on the far side's modes: from the left, the E- and H-continuity equations tested with the right modes, solved for the transmitted and reflected amplitudes together (for one mode per side, t = 2cb / (c + b) and r = (b - c) / (b + c), with c and b the E and H projections); from the right, the same with the left modes. With a complete basis this is the exact interface, like the default. With a truncated basis, power that the far basis cannot represent is lost from the cascade, as radiation leaves an open device: a guided-only basis then keeps the mode-overlap loss of a lateral offset or a straight-to-bend junction. It drops radiation that would return to the guide further on, so it reads low wherever that matters, such as a sharp curvature reversal, a closed box or a very short abrupt step; there a complete basis is needed. Identical bases give the identity in both methods.

enforce_passivity=True clips any gain of an ordinary real guided interface (singular values of the power-normalized S-matrix above 1). It is off by default: the guided traces are not exactly orthogonal, so raw interfaces are slightly active (about 1e-3 on an adiabatic taper, 1e-2 at an abrupt step), and the clip does not remove only that gain. It also pulls the fundamental channel down, by 0.1 % (num_modes=2) to 0.4 % (num_modes=4) on an adiabatic taper at dl 0.04 um, a bias that does not shrink with more modes or slabs (it does shrink with a finer cross-section mesh), while the raw cascade stays lossless to 3e-4. The raw violation is always reported in the diagnostics.

propagation_smatrix

propagation_smatrix(
    modes: Sequence[VectorMode], length_um: float
) -> SMatrix

Diagonal propagation S-matrix for a uniform section of length_um.

Forward and backward both pick up Φ = diag(exp(-i k0 n_eff L - k0 k_eff L)) with no inter-mode coupling; reflection blocks are zero. Equivalently this is exp(-i*beta_eme*L) for beta_eme = k0*mode.n_eff_complex.conjugate().

The generalized interface transforms only its dual test basis, so these raw primal eigenmodes retain their diagonal propagation constants. For an evanescent mode n_eff≈0, k_eff>0 the same expression is the causal decay.

cascade

cascade(segments: Sequence[SMatrix]) -> SMatrix

Cascade a left-to-right sequence of S-matrices with the star product.

star_product

star_product(sa: SMatrix, sb: SMatrix) -> SMatrix

Redheffer star product Sa ⋆ Sb (Sa on the left, Sb on the right), connecting Sa's right port to Sb's left port.

Sa may connect n0 -> n1 and Sb n1 -> n2; only their internal port counts must match. Linear solves are used rather than explicit matrix inverses.

Convergence

eme_convergence_report

eme_convergence_report(
    results: Sequence[EMEResult],
    *,
    labels: Optional[Sequence[str]] = None,
    port_left: int = 1,
    port_right: int = 1,
    tolerance: float = 0.001,
    stability_steps: int = 2,
) -> EMEConvergenceReport

Compare a nested sweep on one fixed propagating external-port subspace.

Radiation eigenvectors reorder when the window or PML changes, so comparing their individual columns is meaningless. Propagating eigenvectors can likewise phase-flip or rotate inside a degenerate subspace. On a common external-port grid this report least-squares aligns each current port basis to the previous one before comparing the four complex S blocks. A requested guided subset may not split a degenerate beta cluster.

energy_balances is the worst-case generalized output/input power ratio over the selected left port subspace, not the coordinate-dependent result of exciting whichever degenerate vector happened to be returned first. Field-continuity (both all-channel and guided-input maxima), reaction orthogonality, and generalized reciprocity are exposed separately; they are not folded into smatrix_stable because useful thresholds depend on the problem and retained channel types.

stability_steps prevents one accidental flat pair from being declared stable. The default requires two consecutive sub-tolerance deltas, hence at least three results. Use stability_steps=1 only for an exact algebraic identity check, not a numerical convergence claim.

EMEConvergenceReport dataclass

EMEConvergenceReport(
    labels: Tuple[str, ...],
    section_mode_counts: Tuple[Tuple[int, ...], ...],
    port_smatrix_deltas: Tuple[float, ...],
    port_subspace_errors: Tuple[float, ...],
    energy_balances: Tuple[float, ...],
    max_interface_residuals: Tuple[float, ...],
    max_reciprocity_errors: Tuple[float, ...],
    tolerance: float,
    stability_steps: int,
    smatrix_stable: bool,
    max_guided_interface_residuals: Tuple[float, ...] = (),
    max_reaction_orthogonality_errors: Tuple[
        float, ...
    ] = (),
)

Successive fixed-port stability evidence for EME basis/PML sweeps.

smatrix_stable means that the requested number of consecutive, basis-aligned propagating-port S-matrix deltas and port-subspace errors are within tolerance. It deliberately does not mean that the underlying field expansion is physically converged. Interface continuity, reciprocity, passivity/energy balance, mesh, window, and PML sweeps remain independent gates.

labels instance-attribute

labels: Tuple[str, ...]

section_mode_counts instance-attribute

section_mode_counts: Tuple[Tuple[int, ...], ...]

port_smatrix_deltas instance-attribute

port_smatrix_deltas: Tuple[float, ...]

port_subspace_errors instance-attribute

port_subspace_errors: Tuple[float, ...]

energy_balances instance-attribute

energy_balances: Tuple[float, ...]

max_interface_residuals instance-attribute

max_interface_residuals: Tuple[float, ...]

max_reciprocity_errors instance-attribute

max_reciprocity_errors: Tuple[float, ...]

tolerance instance-attribute

tolerance: float

stability_steps instance-attribute

stability_steps: int

smatrix_stable instance-attribute

smatrix_stable: bool

max_guided_interface_residuals class-attribute instance-attribute

max_guided_interface_residuals: Tuple[float, ...] = ()

max_reaction_orthogonality_errors class-attribute instance-attribute

max_reaction_orthogonality_errors: Tuple[float, ...] = ()

InterfaceDiagnostics dataclass

InterfaceDiagnostics(
    n_left: int,
    n_right: int,
    rank: int,
    equations: int,
    condition_number: float,
    raw_residual_by_input: ndarray,
    residual_by_input: ndarray,
    raw_passivity_violation: float,
    passivity_violation: float,
    passivity_correction_applied: bool,
    reciprocity_error: float,
    left_gram_condition: float,
    right_gram_condition: float,
    left_gram_asymmetry: float,
    right_gram_asymmetry: float,
    left_reaction_orthogonality_error: float = 0.0,
    right_reaction_orthogonality_error: float = 0.0,
    max_guided_input_residual: float = 0.0,
    method: str = "galerkin",
)

Numerical evidence for one generalized mode-matched interface.

*_reaction_orthogonality_error is the largest self-normalized reaction overlap between modes with distinct propagation constants. It complements Gram condition/asymmetry: a large symmetric off-diagonal overlap can expose an inaccurate high-order E/H reconstruction while leaving both of those older diagnostics deceptively modest. max_guided_input_residual filters the continuity residual to incident guided channels on either side, so a poor continuum test channel cannot hide the accuracy of the physical ports. method records the interface model that built the S-matrix; rank, condition_number and equations describe that model's own system.

n_left instance-attribute

n_left: int

n_right instance-attribute

n_right: int

rank instance-attribute

rank: int

equations instance-attribute

equations: int

condition_number instance-attribute

condition_number: float

raw_residual_by_input instance-attribute

raw_residual_by_input: ndarray

residual_by_input instance-attribute

residual_by_input: ndarray

raw_passivity_violation instance-attribute

raw_passivity_violation: float

passivity_violation instance-attribute

passivity_violation: float

passivity_correction_applied instance-attribute

passivity_correction_applied: bool

reciprocity_error instance-attribute

reciprocity_error: float

left_gram_condition instance-attribute

left_gram_condition: float

right_gram_condition instance-attribute

right_gram_condition: float

left_gram_asymmetry instance-attribute

left_gram_asymmetry: float

right_gram_asymmetry instance-attribute

right_gram_asymmetry: float

left_reaction_orthogonality_error class-attribute instance-attribute

left_reaction_orthogonality_error: float = 0.0

right_reaction_orthogonality_error class-attribute instance-attribute

right_reaction_orthogonality_error: float = 0.0

max_guided_input_residual class-attribute instance-attribute

max_guided_input_residual: float = 0.0

method class-attribute instance-attribute

method: str = 'galerkin'