Eigenmode ports: theory and numerical methods

Start with Eigenmode ports and S-parameter analysis for Python and hash setup, matched feed terminations, modal channel selection, and result interpretation. This page defines the operators and conventions used by the solvers, sources, and monitors. It distinguishes physical-frequency quantities from the discrete symbols needed to reproduce the FDTD update.

The published HTML guide renders the equations and cross-references in full.

Mathematical formulation

The remaining sections describe the component-sampled FDFD eigenproblems, power normalization and modal reconstruction, followed by spectrum synthesis, I/Q injection, Yee staggering, direct DFT reception, and multimode decomposition. The companion Eigenmode tracking theory develops branch assignment, SVD subspace comparison, symmetry decomposition, and basis transport.

Overview

gprMax uses two finite-difference frequency-domain eigenmode solvers:

Time-domain model	Modal cross-section	Solver
2D TM or TE	One physical transverse coordinate	FDFD_1D_mode_solver
3D	Two physical transverse coordinates	FDFD_2D_mode_solver

Both solvers operate directly on component-sampled Yee grids, accept complex permittivity and permeability, enforce PEC and PMC constraints at the corresponding component locations, select the passive effective-index branch, and return either real-power-normalised or E/H-balanced modal fields. Only the propagating, real-power-normalised subset can be used by EigenmodeSource.

Frequency and wavenumber convention

Both FDFD solvers distinguish the physical solve frequency from the symbols of the FDTD differences. At each modal solve frequency,

\[\omega = 2\pi f, \qquad k_0 = \frac{\omega}{c}.\]

These physical quantities are exposed as solver.omega and solver.k0.

Given the owning grid’s time step \(\Delta t\), the eigenproblem and field reconstruction use the leapfrog temporal symbol

\[\Omega = \frac{2}{\Delta t}\sin\left(\frac{\omega\Delta t}{2}\right), \qquad k_{0,\mathrm{operator}} = \frac{\Omega}{c}.\]

These are exposed as solver.operator_omega and solver.operator_k0. The transverse operators retain their component-sampled Yee differences. The eigenvalue determines the longitudinal difference symbol \(K_w\), rather than the phase propagation constant \(\beta\) directly:

\[\begin{split}\lambda&=-n_{\mathrm{operator}}^2, \qquad n_{\mathrm{operator}}=\frac{K_w}{k_{0,\mathrm{operator}}},\\ \beta&=\frac{2}{\Delta w}\sin^{-1}\left(\frac{K_w\Delta w}{2}\right), \qquad n_{\mathrm{eff}}=\frac{\beta}{k_0}.\end{split}\]

Here \(\Delta w\) is the normal cell spacing. The inverse-sine branch is chosen for passive forward propagation, including decay for evanescent modes. solver.beta stores the phase propagation constant. Field reconstruction uses operator_neff; the public complex_neff and modal_real_neff describe its effective index. This distinction keeps the E/H amplitude relationship consistent with the discrete curls while supplying the correct phase to source and monitor staggering.

Lossless modes at or beyond the longitudinal spatial band edge \(|K_w\Delta w/2|\geq1\) are retained for inspection but have power_valid=False and cannot be used for source injection.

Eigenmode sources pass the owning grid’s time step and normal spacing automatically, including a subgrid’s own values. Direct low-level callers enable the same convention with fdtd_dt and propagation_spacing. Omitting fdtd_dt uses \(\Omega=\omega\) and makes operator_k0 equal to k0; omitting propagation_spacing uses \(\beta=K_w\). Omitting both preserves the continuum time and longitudinal conventions of earlier low-level calls. Waveforms, Fourier transforms, and half-time-step phase factors always use the physical frequency \(\omega\).

The source material extraction also includes the midpoint factor in static electric and magnetic conductivity, including static electric conductivity in dispersive media. Bulk material poles and their associated Drude or inclusive conductivity terms are still evaluated from their analytic physical-frequency response, rather than the exact FDTD ADE transfer. Dispersive materials therefore retain a time-discretization mismatch in their pole response that should be checked by convergence.

1D Scalar Solver for 2D Models

fdfd_1d_mode_solver.py supplies the scalar, Yee-staggered mode solve used by eigenmode sources in gprMax 2D TM and TE domains.

1D Coordinates and Yee Shapes

The local basis is (t, a, w):

t

The one physical transverse coordinate.

a

The invariant 2D axis.

w

The in-plane propagation direction and source normal.

For N cells along t, the native staggered shapes are:

Material array	Field component	Shape
eps_r_t	E_t	(N,)
eps_r_a	E_a	(N + 1,)
eps_r_w	E_w	(N + 1,)
mu_r_t	H_t	(N + 1,)
mu_r_a	H_a	(N,)
mu_r_w	H_w	(N,)

The TM reduction solves the node-sampled scalar field E_a and reconstructs H_t and H_w. The TE reduction solves the cell-sampled scalar field H_a and reconstructs E_t and E_w. No derivative is taken through gprMax’s artificial one-cell TM or two-cell TE invariant-axis thickness.

1D Inputs

The constructor signature is:

FDFD_1D_mode_solver(
    frequency,
    dt,
    mode_index,
    polarization,
    eps_r_t,
    eps_r_a,
    eps_r_w,
    mu_r_t,
    mu_r_a,
    mu_r_w,
    pec_t_mask=None,
    pec_a_mask=None,
    pec_w_mask=None,
    pmc_t_mask=None,
    pmc_a_mask=None,
    pmc_w_mask=None,
    guess=None,
    *,
    fdtd_dt=None,
    propagation_spacing=None,
)

frequency

Modal solve frequency in Hz.

dt

Yee-cell spacing along the physical transverse coordinate t, in metres. Despite its name, this is a spatial step rather than a time step.

mode_index

Zero-based requested mode. The solver computes mode_index + 1 modes and exposes the selected one through modal_Et, modal_Ea, modal_Ew, modal_Ht, modal_Ha, modal_Hw and modal_real_neff.

polarization

TM selects the E_a scalar problem. TE selects the H_a scalar problem.

eps_r_* and mu_r_*

Complex relative material arrays sampled at the component locations in 1D Coordinates and Yee Shapes.

pec_*_mask and pmc_*_mask

Optional boolean masks at the matching electric or magnetic component locations. Non-finite values in the corresponding material arrays are also interpreted as constraints.

guess

Optional ARPACK shift. If omitted, the solver derives a shift from the largest finite bulk material magnitude, before SIBC adds surface admittance to the electric coefficients. SIBC remains in the eigenproblem.

fdtd_dt

Optional keyword-only FDTD time step in seconds. A positive finite value enables the leapfrog frequency symbol and requires frequency below temporal Nyquist. This is distinct from the transverse spatial dt.

propagation_spacing

Optional keyword-only positive finite cell spacing along w, in metres. Enables conversion from the longitudinal difference symbol to the phase propagation constant. See Frequency and wavenumber convention.

1D Eigenproblem and Constraints

The solver constructs a node-to-cell difference operator D_nc and its negative adjoint:

D_nc : node fields -> cell fields
D_cn = -D_nc.H : cell fields -> node fields

Both operators are normalised by k0 * dt. For TM, the scalar eigenproblem is assembled on the E_a nodes:

A_TM = -mu_t [D_cn inv(mu_w) D_nc + eps_a]
A_TM E_a = lambda E_a

For TE, it is assembled on the H_a cells:

A_TE = -eps_t [D_nc inv(eps_w) D_cn + mu_a]
A_TE H_a = lambda H_a

In both cases:

operator_neff = sqrt(-lambda)

gprMax uses exp(+j*omega*t - j*beta*w). The square-root branch is therefore chosen with Re(operator_neff) >= 0. For a passive mode Im(operator_neff) <= 0; a purely evanescent mode uses the negative-imaginary branch so that it decays in positive w. The public complex_neff is then recovered as described in Frequency and wavenumber convention.

PEC constraints remove electric scalar degrees of freedom from the TM problem, while PMC constraints remove magnetic scalar degrees of freedom from the TE problem. Longitudinal inverse-material operators are evaluated only on unconstrained degrees of freedom: constrained entries receive a zero inverse. After the reduced sparse eigenproblem is solved, the eigenvectors are expanded back to their full node or cell arrays and every constrained field component is explicitly zeroed.

1D Field Reconstruction

For TM, the selected E_a eigenvector gives:

H_t = -operator_neff E_a / (eta0 mu_t)
H_w = i inv(mu_w) D_nc E_a / eta0

For TE, the selected H_a eigenvector gives:

E_t = eta0 operator_neff H_a / eps_t
E_w = -i eta0 inv(eps_w) D_cn H_a

The other three field components are identically zero for the selected 2D polarization. The reconstructed electric fields are in V/m and magnetic fields are in A/m.

1D Normalisation and Phase Alignment

Each propagating mode is normalised to one watt per metre of invariant-axis length. A non-propagating mode instead uses the positive E/H-balanced scale defined in Mode Selection, Fields, and Power; it is retained as a monitor reference but is not a one-watt power wave. For TM, node-sampled E_a and H_t are first averaged onto transverse cells:

P_TM = 0.5 Re sum(-E_a H_t*) dt

For TE, E_t and H_a already share the cell locations:

P_TE = 0.5 Re sum(E_t H_a*) dt

If a passive branch has negative real power along the solver’s forward axis, the propagation-constant branch is reversed and its dependent fields are reconstructed before normalisation. Each complex mode is then phase-rotated to a deterministic real-profile convention used for tracking and, for a propagating source anchor, FDTD injection.

1D gprMax Integration

EigenmodeSource samples component materials from the mode’s live invariant layer, supplies the corresponding PEC/PMC masks, and maps the returned line profiles back into the thin 3D Yee arrays used by the FDTD source kernels. The TM source uses the single live invariant layer; the TE source uses the shared interior layer of its two-cell invariant thickness. Inactive components and TE outer boundary planes are explicitly zero.

Modes are selected using the same shift-invert convention as the full-vector 2D solver. plot_fields writes one row per computed mode with line plots of all three active fields, including their node- or cell-sampled locations.

2D Full-Vector Solver for 3D Models

fdfd_2d_mode_solver.py contains FDFD_2D_mode_solver, the full-vector solver used when a gprMax eigenmode source has two physical transverse coordinates.

2D Coordinates and Yee Shapes

The solver works in a local coordinate system rather than the global gprMax x, y and z axes:

u

First transverse source-plane axis.

v

Second transverse source-plane axis.

w

Propagation-normal axis.

For example, a source normal to global z uses local u=x, v=y and w=z. A source normal to global x uses local u=y, v=z and w=x.

The solver is built on a true staggered Yee grid. Material arrays supplied to the constructor must already be sampled at the corresponding local field component locations. The solver does not average cell-centred material data and it does not collapse all fields onto a common rectangular array.

For a transverse source region containing Nu by Nv Yee cells, the expected local component shapes are:

Array	Field component	Shape
eps_r_uu	E_u	(Nu, Nv + 1)
eps_r_vv	E_v	(Nu + 1, Nv)
eps_r_ww	E_w	(Nu + 1, Nv + 1)
mu_r_uu	H_u	(Nu + 1, Nv)
mu_r_vv	H_v	(Nu, Nv + 1)
mu_r_ww	H_w	(Nu, Nv)

The returned modal fields use the same native Yee shapes:

Modal field	Shape
modal_Eu	(Nu, Nv + 1)
modal_Ev	(Nu + 1, Nv)
modal_Ew	(Nu + 1, Nv + 1)
modal_Hu	(Nu + 1, Nv)
modal_Hv	(Nu, Nv + 1)
modal_Hw	(Nu, Nv)

Only transverse fields participate in gprMax eigenmode source injection. The longitudinal fields E_w and H_w are still reconstructed because they are part of the full-vector mode solution, but the TF/SF correction kernels use only the tangential/transverse modal components.

2D Inputs

The constructor signature is:

FDFD_2D_mode_solver(
    frequency,
    du,
    dv,
    mode_index,
    eps_r_uu,
    eps_r_vv,
    eps_r_ww,
    mu_r_uu,
    mu_r_vv,
    mu_r_ww,
    pec_u_mask=None,
    pec_v_mask=None,
    pec_w_mask=None,
    pmc_u_mask=None,
    pmc_v_mask=None,
    pmc_w_mask=None,
    guess=None,
    surface_boundary=None,
    *,
    fdtd_dt=None,
    propagation_spacing=None,
)

frequency

Source frequency in Hz.

du, dv

Local transverse cell sizes in metres. The solver normalises finite- difference operators by k0 * du and k0 * dv.

mode_index

Zero-based modal index requested by the caller. The solver computes mode_index + 1 modes, then exposes the requested mode through modal_Eu, modal_Ev, modal_Ew, modal_Hu, modal_Hv, modal_Hw and modal_real_neff.

eps_r_* and mu_r_*

Complex relative permittivity and permeability arrays sampled at the local Yee component locations listed in 2D Coordinates and Yee Shapes.

pec_u_mask, pec_v_mask, pec_w_mask

Optional explicit boolean masks for constrained electric degrees of freedom. They must match the E_u, E_v and E_w shapes.

pmc_u_mask, pmc_v_mask, pmc_w_mask

Optional explicit boolean masks for constrained magnetic degrees of freedom. Non-finite entries in the matching permeability arrays are also interpreted as PMC.

guess

Optional ARPACK shift. If omitted, the solver chooses a conservative shift from the largest finite bulk material magnitude, before SIBC adds surface admittance to the electric coefficients. SIBC remains in the eigenproblem.

surface_boundary

Optional compiled impedance-volume boundary. See Surface impedance and impedance volumes for retained-component masks and clipped curl rows.

fdtd_dt

Optional keyword-only FDTD time step in seconds. A positive finite value enables the leapfrog frequency symbol and requires frequency below temporal Nyquist.

propagation_spacing

Optional keyword-only positive finite cell spacing along w, in metres. Enables conversion from the longitudinal difference symbol to the phase propagation constant. See Frequency and wavenumber convention.

2D Array Ordering

All component arrays are flattened with Fortran order:

flat = array.ravel(order='F')

and modal vectors are reshaped back with:

array = vector.reshape((*shape, num_modes), order='F')

There is no axis-order switch. gprMax must pass local u/v slices in the same native transverse ordering used by the extracted source plane.

2D PEC Handling

PEC is represented as constrained electric degrees of freedom, not as a large finite permittivity approximation.

The solver detects electric PEC in two ways:

1. Explicit electric masks passed through pec_u_mask, pec_v_mask and pec_w_mask.

2. Non-finite values, normally np.inf + 0j, in the electric material arrays.

For example:

eps_r_uu[pec_u_mask] = np.inf + 0j
eps_r_vv[pec_v_mask] = np.inf + 0j
eps_r_ww[pec_w_mask] = np.inf + 0j

Each electric component is treated independently:

• pec_u_mask constrains E_u to zero.

• pec_v_mask constrains E_v to zero.

• pec_w_mask constrains E_w to zero.

Port extraction uses the final native Yee electric material IDs, including samples shared by neighbouring geometry objects. It does not expand cell-centred PEC voxels into extra electric constraints. Consequently an air bore carved with averaging disabled is solved with the same electric clamps as its FDTD update. Averaging and object order can still change the represented geometry and its cutoff.

After masks are built, PEC material entries are replaced by finite placeholders before matrix assembly:

eps_r_uu[self.pec_u_mask] = 1.0 + 0j

The physics is carried by removed/constrained degrees of freedom, not by the placeholder value. Large finite values such as 1e8 or 1e10 are ordinary finite material values and are intentionally not treated as PEC.

2D Eigenproblem

The solver constructs rectangular sparse derivative matrices between true Yee component grids. The core local operators are:

DEU_EW_TO_EU : E_w -> E_u
DEV_EW_TO_EV : E_w -> E_v
DEU_EV_TO_HW : E_v -> H_w
DEV_EU_TO_HW : E_u -> H_w

and their adjoint magnetic-grid counterparts:

DHU_HV_TO_EW = -DEU_EW_TO_EU.H
DHV_HU_TO_EW = -DEV_EW_TO_EV.H
DHU_HW_TO_HU = -DEU_EV_TO_HW.H
DHV_HW_TO_HV = -DEV_EU_TO_HW.H

The transverse electric field vector is:

Euv = [E_u, E_v]^T

The solver forms the standard full-vector FDFD P and Q matrices and solves:

Omega * Euv = eigenvalue * Euv
Omega = P * Q

where the operator index is recovered from:

operator_neff = sqrt(-eigenvalue)

Here the matrix Omega = P * Q is distinct from the scalar temporal symbol \(\Omega\). The branch follows exp(+j*omega*t - j*beta*w): Re(operator_neff) >= 0 and Im(operator_neff) <= 0 for passive forward propagation. When the real part is zero, the negative-imaginary branch gives evanescent decay in positive w. The public complex_neff is then recovered as described in Frequency and wavenumber convention.

Because the operators connect the correct staggered Yee component grids, there is no separate PEC-neighbour spurious-mode rejection heuristic in this solver. The old candidate scoring/filtering path has been removed.

2D Degree-of-Freedom Reduction

PEC constraints are applied by removing constrained transverse electric degrees of freedom from the eigenproblem:

Omega = Omega[self.free_euv_mask, :][:, self.free_euv_mask]

Removing electric samples does not, in general, justify setting their paired transverse magnetic samples to zero. Write \(s=i n_{\mathrm{op}}\) and use the solver’s normalized transverse magnetic vector \(h\), so that the unconstrained equations are \(s e=P h\) and \(s h=Q e\). Split magnetic samples into \(a\), whose paired transverse electric update is live, and \(c\), whose paired electric sample is PEC. The latter have no live Ampere equation, but Faraday’s law at the clamped electric sample still requires

\[P_{ca}h_a+P_{cc}h_c=0,\qquad s h_a=Q_a e,\qquad s h_c=-P_{cc}^{-1}P_{ca}Q_a e.\]

Here the first index of \(P\) selects electric rows (free \(f\) or constrained \(c\)); the second selects magnetic samples. The exact reduced electric operator is therefore

\[\left(P_{fa}Q_a-P_{fc}P_{cc}^{-1}P_{ca}Q_a\right)e =s^2 e=-n_{\mathrm{op}}^2 e.\]

For conventional voxel PEC walls the forcing \(P_{ca}Q_a\) vanishes on the free electric columns, and the usual sparse product suffices. Otherwise the solver applies this Schur complement through a sparse factorization of \(P_{cc}\). Shift-invert iteration factors the sparse block system

\[\begin{split}\begin{pmatrix} P_{fa}Q_a-\sigma I & P_{fc}\\ P_{ca}Q_a & P_{cc} \end{pmatrix},\end{split}\]

whose upper inverse block is the shifted reduced inverse. No dense Schur matrix is formed for the sparse eigensolve. The same static reconstruction supplies the modal H fields and the operator used for tracking residuals. A singular reconstruction is a numerical failure. Actual PMC constraints and impedance-volume retained-sample masks continue to exclude their own magnetic degrees of freedom.

When the reduced matrix has only one more degree of freedom than the requested mode count, the solver uses a dense eigensolve because ARPACK requires k < N - 1. Larger systems use shift-invert ARPACK and retry with a roundoff-scale shift perturbation if the original shift produces a singular factorisation. The reduced eigenvectors are then expanded back to the full transverse field-vector size and constrained fields are explicitly zeroed.

The inverse eps_r_ww and mu_r_ww operators are built only on free longitudinal degrees of freedom. Constrained entries receive zero inverse values so that no division by np.inf or placeholder data affects the reconstructed fields.

2D Field Reconstruction

After solving the eigenproblem, the solver reconstructs:

• E_u and E_v directly from the transverse eigenvector.

• H_u and H_v from Q * Euv / sqrt(eigenvalue), using the branch sqrt(eigenvalue) = +j*operator_neff selected by the propagation convention.

• E_w from transverse magnetic curl terms.

• H_w from transverse electric curl terms.

Magnetic fields are converted to physical A/m using eta0:

H = 1j * H_normalized / eta0

The solver then zeroes all constrained fields to ensure returned modal fields satisfy the enforced constraints exactly.

2D Normalisation and Phase Alignment

Propagating modes are normalised to carry one watt of time-average power. Non-propagating modes instead use the finite E/H-balanced scale defined in Mode Selection, Fields, and Power; they can enter the tracked monitor-reference bank but not the source/power bank. Real power is computed from cell-centred transverse Poynting flux by averaging the staggered transverse fields onto local cells:

P = 0.5 * Re integral((E_u * H_v* - E_v * H_u*) dA)

If a passive branch carries negative real power along the solver’s forward axis, the propagation-constant branch is reversed and all dependent fields are reconstructed before normalisation. Each complex mode is then phase-rotated to a deterministic real-profile convention. This makes plotted, tracked, and (for propagating source anchors) injected fields easier to compare.

2D gprMax Integration

sources.py extracts complex material tensors from G.ID after the Yee grid has been built. This is the correct integration point because G.ID is already sampled at Yee component locations.

For a source plane, sources.py maps global components into local u/v/w components:

local_to_global = (
    self.transverse_axes[0],
    self.transverse_axes[1],
    self.normal_axis,
)

It then extracts six native Yee slices:

• electric local u component: (Nu, Nv + 1)

• electric local v component: (Nu + 1, Nv)

• electric local w component: (Nu + 1, Nv + 1)

• magnetic local u component: (Nu + 1, Nv)

• magnetic local v component: (Nu, Nv + 1)

• magnetic local w component: (Nu, Nv)

For electric materials:

• finite conductivity is folded into complex permittivity,

• se == inf is converted to np.inf + 0j, which the solver treats as PEC.

For magnetic materials:

• finite magnetic conductivity is folded into complex permeability,

• sm == inf is converted to np.inf + 0j, which the solver treats as PMC.

PMC constraints, like PEC constraints, follow the final native Yee component material IDs. Cell-centred PMC voxels are not expanded into additional magnetic constraints, including when reducing the port cross-section to a 1D solve.

After solving, sources.py maps local modal fields back to global component slots. The Cython injection kernels consume the transverse components with their native staggered shapes; longitudinal modal fields are stored but are not used for TF/SF source corrections.

Source synthesis and FDTD injection

This section connects the FDFD solvers described above to the complete eigenmode band/port/excitation workflow. It fixes the phasor and propagation signs, shows how a solved mode becomes real FDTD update terms, and describes the single-frequency, in-phase/quadrature (I/Q), and broadband synthesis paths.

Phasor, Fourier, and Propagation Conventions

gprMax advances real fields in time. Whenever a complex frequency-domain field is used, it follows the electrical-engineering convention

\[\mathbf{E}(u,v,w,t) = \operatorname{Re}\left\{ \widetilde{\mathbf{E}}(u,v,\omega) \exp(+j\omega t-j\beta w) \right\}.\]

The forward Fourier transform uses the complementary negative-sign kernel,

\[X(\omega)=\int_{-\infty}^{\infty}x(t)\exp(-j\omega t)\,dt,\]

so NumPy’s rfft returns the positive-frequency coefficients used by this phasor convention, and irfft reconstructs the real signal with positive frequency time dependence exp(+j*omega*t). The FFT sign is not a second physical convention: it is the analysis kernel paired with the phasor synthesis sign.

With exp(+j*omega*t), the continuum Maxwell curl equations are

\[\begin{split}\nabla\times\widetilde{\mathbf{E}} &= -j\omega\widetilde{\boldsymbol{\mu}} \widetilde{\mathbf{H}},\\ \nabla\times\widetilde{\mathbf{H}} &= +j\omega\widetilde{\boldsymbol{\epsilon}}_c \widetilde{\mathbf{E}}.\end{split}\]

Their continuum electric and magnetic conductivity terms are

\[\begin{split}\epsilon_{r,c}(\omega) &= \epsilon_r(\omega) -j\frac{\sigma}{\omega\epsilon_0},\\ \mu_{r,c}(\omega) &= \mu_r(\omega) -j\frac{\sigma_m}{\omega\mu_0}.\end{split}\]

For nondispersive material extraction on an FDTD grid, the solver instead uses the exact midpoint conductivity terms

\[\begin{split}\epsilon_{r,c}^{\mathrm{Yee}} &= \epsilon_r-j\frac{\sigma\cos(\omega\Delta t/2)}{\Omega\epsilon_0},\\ \mu_{r,c}^{\mathrm{Yee}} &= \mu_r-j\frac{\sigma_m\cos(\omega\Delta t/2)}{\Omega\mu_0},\end{split}\]

with the temporal symbol from Frequency and wavenumber convention. Low-level callers supply their own complex relative material arrays; passing fdtd_dt does not reinterpret those arrays as conductivity parameters.

For a forward passive mode,

\[\beta=k_0n_{\mathrm{eff}}=\beta_r-j\alpha, \qquad \alpha\geq 0,\]

and therefore

\[\exp(-j\beta w) =\exp(-j\beta_r w)\exp(-\alpha w).\]

The selected propagation branch consequently satisfies

\[\begin{split}\operatorname{Re}(n_{\mathrm{eff}})&\geq 0,\\ \operatorname{Im}(n_{\mathrm{eff}})&\leq 0 \quad\text{for passive propagation}.\end{split}\]

If the real part is numerically zero, a purely evanescent mode uses Im(n_eff) < 0 so it decays in positive local w. The imaginary part of the square root is never replaced by its absolute value; its sign contains the loss or gain information.

As a continuum example, at 5 GHz a homogeneous material with epsilon_r=9 and sigma=2 S/m has

\[\epsilon_{r,c}\simeq 9-j7.190, \qquad n\simeq 3.203-j1.122.\]

This gives \(\alpha=-k_0\operatorname{Im}(n)\simeq 117.6\ \mathrm{m}^{-1}\). The field magnitude after 0.5 mm is approximately exp(-alpha * 0.5e-3) = 0.943.

Build-Time Workflow

The source is prepared after geometry construction, when material IDs already refer to their final Yee-component locations. The workflow is:

1. Validate the source plane, requested direction, mode index, waveform, and one or more solve frequencies.

2. Choose local coordinates (u, v, w). The u and v axes lie in the source plane and w is normal to it. In a 2D TM or TE model, one of the transverse axes is the invariant axis and only one is physical.

3. Extract complex relative permittivity and permeability at every native Yee component position on the source plane. PEC and PMC cells become explicit component constraints.

4. At each requested frequency, solve either the 1D scalar TM/TE problem for a 2D FDTD model or the 2D full-vector problem for a 3D FDTD model, using the owning grid’s time step and normal cell spacing.

5. Reconstruct all modal E and H components, zero constrained degrees of freedom, apply either real-power or balanced E/H normalization, and choose a consistent global phase.

6. Map local modal arrays back to global x/y/z component slots. If the local coordinate mapping is left-handed, reverse H so that the Poynting direction remains correct.

7. Select real-only, single-frequency I/Q, or multi-anchor broadband synthesis.

8. During time stepping, apply tangential incident E and H as TF/SF corrections on the appropriate side of the source plane.

The modal solution describes the cross-section at the source reference plane. Propagation away from that plane is performed by the FDTD grid itself. The explicit propagation constant is used in the broadband E/H staggering and in frequency interpolation, not to overwrite fields throughout the guide.

Mode Selection, Fields, and Power

The user-facing mode_index is one-based: mode 1 selects the first eigenpair after the solver orders the computed eigenvalues. The low-level Python FDFD solver classes retain zero-based array indices internally. Shift-invert sparse eigensolving targets modes near a material-derived guess. Near cut-off, degeneracy, or an eigenvalue crossing, the same numerical index can identify a different physical mode at another frequency; broadband sources and receivers therefore check adjacent-anchor overlap.

For a 3D model, transverse modal power is evaluated on local cells as

\[P=\frac{1}{2}\operatorname{Re} \sum_{u,v} \left(E_uH_v^*-E_vH_u^*\right)\Delta u\Delta v.\]

The native Yee components are averaged only as needed to place each product on the same transverse cell. For a 2D model, the equivalent line integral gives power per metre along the invariant axis. If a passive branch initially has negative real power along the solver’s forward axis, gprMax reverses its propagation-constant branch and reconstructs the dependent fields. A propagating mode is then scaled to one watt in 3D or one watt per metre in 2D. A non-propagating mode has no independent forward real power and instead receives a finite balanced E/H field scale for tracking and diagnostics; that scale is not a one-watt normalization.

This normalization defines the scale of the modal profile at the source plane. Multiplying the source waveform amplitude by a factor multiplies both incident E and H by that factor; for a monochromatic mode, time-average power therefore scales with the square of the waveform amplitude only where the mode has valid real-power normalization.

The one-watt propagating profiles form the source-synthesis and real-power bank. Modal monitoring also retains a second, tracked reference bank. Before a reference profile is used in a generalized-only bin, both its E and H fields are multiplied by the same factor so that

\[P_\mathrm{bal}=\frac{1}{4\eta_0}\int \left(|\mathbf E_t|^2+\eta_0^2|\mathbf H_t|^2\right)\,dA=1\]

(or the corresponding 2D line integral). This positive balanced quantity is a field-scale convention, not transported real power. Applying it to both propagating and evanescent monitor references prevents their original solver scales from introducing an artificial coefficient-scale jump. It does not license interpolating through the cutoff singularity: propagating and evanescent references remain in separate interpolation branches.

The solver returns the positive-local-w mode. A source requested in the negative global direction retains the electric profile and reverses the magnetic orientation in the TF/SF updates, as required for the opposite Poynting vector.

Global Phase and the Real-Only Test

An eigenvector has arbitrary complex phase. Before deciding how to inject a single-frequency profile, gprMax finds the phase rotation that makes all tangential E and impedance-scaled H components as real as possible. If \(F_q\) denotes every such component sample, with H samples first scaled by \(\eta_0\), the rotation is

\[\phi=-\frac{1}{2}\arg\left(\sum_q F_q^2\right), \qquad F_q' = F_q\exp(j\phi).\]

Only tangential components enter this test because longitudinal mode components can be intrinsically in quadrature while not participating in the TF/SF correction. After rotation, the normalized imaginary residual is

\[r=\sqrt{ \frac{\sum_q\left[\operatorname{Im}(F_q')\right]^2} {\sum_q\left|F_q'\right|^2} }.\]

If r <= 1e-8, the tangential spatial profiles are effectively real. gprMax stores their real parts and multiplies them directly by the requested real waveform. The H waveform is evaluated with its Yee half-time-step and half-normal-cell phase delay. This path avoids unnecessary FFT preparation for ordinary lossless fixed-profile modes.

Why Complex Modes Need I/Q Injection

A single real scalar waveform cannot reproduce a general complex spatial profile. Let

\[F=F_R+jF_I, \qquad Z(t)=\operatorname{irfft}\{C(\omega)\}.\]

The desired real field is

\[\operatorname{Re}\{F C(\omega)\exp(j\omega t)\}.\]

gprMax realizes this with two real bases:

\[F_R\operatorname{irfft}\{C\} +F_I\operatorname{irfft}\{jC\}.\]

Because irfft(j*C) is the negative quadrature of irfft(C), this sum is exactly the required real part. No complex values are passed to the real FDTD update arrays.

When the single-frequency residual exceeds the tolerance, gprMax uses this I/Q construction with one modal anchor. The same solved profile and n_eff are used for every significant FFT bin of the waveform. This is a fixed-profile approximation: use multiple solve frequencies when the modal shape or propagation constant varies appreciably across the waveform bandwidth.

Spectrum and Piecewise-Linear Modal Interpolation

For the automatic waveform, gprMax first forms the requested smooth spectral magnitude and finds the circular time support of its zero-phase analytic envelope. It then applies the earliest discrete delay that places all envelope samples above the significance threshold inside the causal record. This keeps the requested spectrum while reserving as much of the simulation window as possible for the packet to traverse the model and decay.

The real waveform is sampled at the FDTD interval for the requested number of iterations. It is zero-padded to the next power of two at least twice as long as the simulation record, and transformed with numpy.fft.rfft. For FFT bin \(m\), let its frequency and coefficient be \(f_m\) and \(S_m\). Non-finite samples and a zero or non-finite spectrum are errors; gprMax does not replace the excitation with a zero-valued source.

For source synthesis, piecewise-linear weights \(w_{k,m}\) interpolate between surrounding propagating source anchors and satisfy

\[\sum_k w_{k,m}=1.\]

Below or above the anchor range, the nearest endpoint receives weight one. This avoids a hard spectral truncation. Significant waveform energy outside the anchor range is governed by the spectrum-coverage policy: the default is an error, while an explicit warn policy permits endpoint extrapolation with a warning. The source-interpolated fields and propagation constant are

\[\begin{split}\mathbf{E}_m &= \sum_k w_{k,m}\mathbf{E}_k,\\ \mathbf{H}_m &= \sum_k w_{k,m}\mathbf{H}_k,\\ n_{\mathrm{operator},m} &= \sum_k w_{k,m}n_{\mathrm{operator},k},\\ K_m &= \frac{\Omega(f_m)}{c}n_{\mathrm{operator},m},\\ \beta_m &= \frac{2}{\Delta w}\sin^{-1}\left(\frac{K_m\Delta w}{2}\right).\end{split}\]

The inverse spatial difference is evaluated at each FFT bin after operator index interpolation, including endpoint extrapolation. Single-frequency I/Q sources, banks with only one retained anchor, and downstream solvers without operator-index metadata retain the constant/legacy physical-index convention \(\beta_m=2\pi f_m n_m/c\). In particular, a single-frequency source holds its physical phase index constant across the pulse.

Significant source energy in the longitudinal grid stop band is an error: refine the normal cell spacing or narrow the excitation bandwidth. Sub-threshold tails in that stop band are discarded; the scalar waveform reconstruction error includes their removal. DC and Nyquist bins are excluded before this propagation calculation.

Linear interpolation of individually normalized modes does not in general preserve unit power. gprMax constructs the cross-power matrix

\[P_{kl}=\frac{1}{2}\int \left(\mathbf{E}_k\times\mathbf{H}_l^*\right) \cdot\hat{\mathbf{w}}\,dA\]

(or the corresponding 2D line integral), and calculates

\[p_m=\operatorname{Re} \left\{\sum_{k,l}w_{k,m}P_{kl}w_{l,m}\right\}, \qquad a_m=\frac{1}{\sqrt{p_m}}.\]

Thus the interpolated frequency-bin source mode is renormalized rather than assuming that linear field weights retain one-watt power. Invalid or nearly zero interpolated power is an error because a finite fallback would no longer represent the requested one-watt incident wave.

Yee Time and Space Staggering

Electric and magnetic fields are offset by half a time step and their tangential samples at the TF/SF surface are offset by half a normal cell. For each bin, gprMax applies the magnetic staggering factor

\[M_m=\exp\left[j\left( \frac{\omega_m\Delta t}{2} +\frac{\beta_m\Delta w}{2} \right)\right].\]

Here \(\omega_m=2\pi f_m\) is the physical angular frequency, and \(\beta_m\) is recovered from the longitudinal difference symbol. The temporal difference symbol \(\Omega\) is not a phase frequency.

The local coordinate is defined in the requested propagation direction. The tangential H correction lies half a cell on the incident side of the electric reference plane, giving the positive spatial phase in this relative factor. For a lossy mode this local factor can have magnitude greater than one because the incident-side sample precedes the reference plane; this does not represent growth in the forward direction. Forward propagation over a positive distance d is always tested by exp(-j*beta*d), whose magnitude is below one for a passive mode.

The spectral coefficients assigned to anchor k are therefore

\[\begin{split}C^E_{k,m} &= w_{k,m}S_m a_m,\\ C^H_{k,m} &= w_{k,m}S_m a_m M_m.\end{split}\]

Each anchor field is split into real and imaginary arrays, and each coefficient set is inverse-transformed both normally and after multiplication by j. The FDTD source update then sums the two I/Q bases over all anchors.

DC and, for an even transform length, the Nyquist bin are self-conjugate. They cannot carry a general complex modal coefficient while preserving a real time record, so gprMax discards those two bins. Significant DC or Nyquist energy produces a warning that the requested excitation has been changed. Use a band-limited waveform; for a finite frequency band, EigenmodeExcitation(..., waveform='auto') synthesizes one automatically.

TF/SF Injection into the FDTD Updates

The source plane is a total-field/scattered-field interface. In continuous form, the incident mode corresponds to equivalent surface currents

\[\begin{split}\mathbf{J}_s &= \hat{\mathbf{n}}\times\mathbf{H}_{\mathrm{inc}},\\ \mathbf{M}_s &= -\hat{\mathbf{n}}\times\mathbf{E}_{\mathrm{inc}}.\end{split}\]

The discrete implementation applies the equivalent correction directly to the curl updates adjacent to the plane:

• the magnetic-field update uses the incident tangential E profile;

• the electric-field update uses the incident tangential H profile;

• the side of the plane and correction sign depend on + or - source direction;

• only transverse/tangential components are injected. Longitudinal fields are reconstructed and retained for diagnostics but do not enter these TF/SF corrections.

Because modal fields are stored on their native Yee component grids, the source does not resample every component onto a common rectangle. The update kernels consume the component-specific array shapes documented earlier in this page. This preserves PEC/PMC constraints, tangential staggering, and the mode solver’s discrete curl relationships.

For real-only sources, one real modal array per component is multiplied by the waveform value. For I/Q and broadband sources, each update sums all anchor and quadrature contributions prepared by the inverse FFT. In both cases, source activation is clipped to the configured waveform start and stop times.

Single-frequency sources record single_frequency_iq_reasons and log the individual selection reasons alongside the measured modal-profile residual: complex modal profile, drive phase/delay, and/or complex longitudinal staggering. A small modal residual can therefore coexist with I/Q injection. Negative real propagation is handled by the signed real-only time shift and does not itself require I/Q. Eigenmode solvers, source staggering, and monitor propagation use the active simulation’s em_consts["c"] for the speed of light.

Reduced-mode surface impedance

The 1D solver accepts an optional surface_boundary in its local (t, a, w) basis. gprMax projects this from the compiled SIBC geometry, retaining only the live invariant layer. The TM electric scalar and TE longitudinal electric rows use the clipped Ampere derivative and the exact discrete surface-ADE admittance. Faraday retains its ordinary Yee derivative; the two derivatives need not be negative adjoints at a clipped wall. resistance=float('inf') gives exactly zero surface admittance.

CPU VirtualWaveguide supports all 2D TE/TM orientations. Its modal window spans the full invariant storage dimension. Its artificial PEC rim constrains tangential E only at the ends of the physical transverse coordinate, including where a cropped SIBC row needs a magnetic sample outside the window. All ports constrain complete SIBC rows on that rim to PEC too, whether or not a virtual guide is attached. This matches the auxiliary aperture and PML updates. Opaque padding moves a physical wall inside the window and retains its SIBC equation. The guide and its retained host must remain uniform along propagation through the aperture and PML. Both ordinary and virtual sources apply the surface-row modal forcing and ADE correction. See Waveguides and absorbers for setup and PML and virtual-guide coupling for validation.

Virtual-guide aperture coupling

For a virtual guide the TF/SF incident-mode correction is applied on an internal plane of the auxiliary grid rather than the main grid. Let \(w\) be the port-normal coordinate and \(u,v\) the transverse coordinates. At the main/auxiliary split the normal magnetic sample is shared,

\[H_w^{\mathrm{aux}}\big|_{\Gamma} = H_w^{\mathrm{main}}\big|_{\Gamma},\]

and the aperture updates for \(E_u\) and \(E_v\) use the ordinary Yee curl with the two normal-neighbour magnetic samples taken from different grids. In schematic form,

\[\begin{split}E_u^{n+1}\big|_{\Gamma} &= C_A E_u^n\big|_{\Gamma} + C_v\,\Delta_v H_w - C_w\left(H_v^{\mathrm{aux}}-H_v^{\mathrm{main}}\right),\\ E_v^{n+1}\big|_{\Gamma} &= C_A E_v^n\big|_{\Gamma} + C_w\left(H_u^{\mathrm{aux}}-H_u^{\mathrm{main}}\right) - C_u\,\Delta_u H_w.\end{split}\]

The signs reverse consistently for the opposite port direction. The updated tangential E samples are shared with the main-grid aperture, while the duplicate main-grid continuation behind the aperture is disconnected. The 3D axis/direction variants use compiled Cython kernels. Reduced 2D coupling uses vectorized operations on the live field layer, preserving its native staggering without introducing invariant-axis side walls.

On the non-distributed CPU solver, bulk electric dispersion is supported in both 3D and reduced TE/TM virtual guides. Conductivity and instantaneous polarization storage are included in the electric update coefficients. For each dispersive aperture sample, the curl result above is completed by

\[\begin{split}\Phi^n &= \sum_m \operatorname{Re}(a_m T_m^n),\\ E^{n+1} &= E_{\mathrm{curl}}^{n+1} - C_\Phi\Phi^n,\\ T_m^{n+1} &= f_m T_m^n + b_m(E^n-E^{n+1}).\end{split}\]

These are the same real or complex bulk ADE coefficients used inside the grid. The auxiliary boundary owns this history because the ordinary bulk kernels omit tangential E on its outer plane. The final corrected field is shared with the main grid. For SIBC contact rows, the auxiliary guide instead copies the complete sparse constitutive equation, including every retained bulk pole and independent Foster history, and uses the coupled sparse solve. This also applies inside its longitudinal PML. Dispersive aperture coupling on accelerator and MPI backends remains unsupported.

Modal Receivers, Direct DFT, and S-parameters

Exactly one global band and zero or more distinct port/mode excitations may exist whenever modal ports are used. Excitation may be omitted only when every port is a passive virtual guide; that form writes raw modal spectra but no S matrix. Every EigenmodePort owns one monitor at its reference plane; each selected excitation additionally applies a TF/SF source. Port indices are one-based and unique, and each port carries an explicit tuple of monitored mode indices. Every excitation selects one of the modes listed by its port. All ports accumulate the common DFT bins from EigenmodeBand. Exactly one active channel permits S-parameter normalization; a simultaneous driven state retains its decomposed waves without constructing S. Automatic ports share one candidate-frequency list, but tracking, retained source/reference masks, and fallback policies are resolved independently for each port and mode. Ports with explicit anchors retain their individually requested candidate frequencies. For each requested frequency \(f_q\), a Cython kernel applies the recursive DFT

\[X_q^{(n+1)}=X_q^{(n)} +\Delta t\,x_n\exp(-j2\pi f_q n\Delta t).\]

The phase factor is advanced once per bin and time step. To bound complex64 recurrence drift in long simulations, both oscillators are reconstructed from their physical times every 1024 iterations using float64 argument reduction; the HDF5 PhaseReanchorInterval attribute records this interval. Electric fields use integer FDTD times, while magnetic fields use their half-time-step phase. Transverse Yee components are averaged to common cells only inside the projection. A passive receiver samples the magnetic reference plane half a normal cell upstream of the electric plane. At a TF/SF source, H is sampled half a cell downstream so that both fields are on the total-field side. The final decomposition corrects either offset using each mode’s propagation constant.

The DFT-bin monitor basis is selected independently from source synthesis. A power-wave-valid bin uses the interpolated one-watt propagating bank. A generalized-only bin uses anchor_mode_reference_valid instead. It first selects the applicable contiguous evanescent reference run inside the solved candidate range; interpolation never crosses cutoff or spans disconnected evanescent runs. Outside that range, it uses the nearest tracked reference endpoint. Every selected anchor E/H pair is divided by the square root of its anchor_balanced_power, then E, H, and the operator index are interpolated with identical branch-local weights. The propagation constant is mapped at each DFT frequency using the same discrete symbols as source synthesis; branches retaining only one anchor keep their physical phase index constant. The interpolated cell-centred E/H pair is balanced once more before its Gram matrices are formed. Keeping all three quantities on the same tracked branch is essential below cutoff, where modal admittance becomes reactive; interpolating only the propagation constant while retaining a propagating endpoint E/H pair would not correctly separate forward and backward amplitudes at one plane.

A bin that enters the numerical spatial stop band after interpolation loses power-wave validity. It retains its selected propagating-bank weights, fields, and propagation constant, with final cell-centred balanced normalization for generalized coefficients. It does not select a different physical-cutoff branch. Existing conditioning and separation checks still determine whether those coefficients can be reported.

The HDF5 port group preserves anchor_complex_neff as the physical phase index. When available, anchor_operator_neff stores the corresponding dimensionless operator indices with the same anchor/mode axes. beta records the actual propagation constant used at every DFT bin, in radians per metre, with frequency/mode axes. These datasets also describe runs that reuse cached modal anchors.

For several requested modes, independent overlaps are insufficient when the discrete profiles are not exactly orthogonal. At each frequency gprMax forms electric and magnetic Gram matrices,

\[\begin{split}G^E_{mn} &= \frac{1}{2}\int (\mathbf E_n\times\mathbf H_m^*)\cdot\hat{\mathbf w}\,dA,\\ G^H_{mn} &= \frac{1}{2}\int (\mathbf E_m^*\times\mathbf H_n)\cdot\hat{\mathbf w}\,dA,\end{split}\]

and solves both systems for the total electric and magnetic modal coefficients. If \(x=a+b\) is the electric coefficient and the staggered magnetic coefficient is \(y=p_+a-p_-b\), where \(p_+\) and \(p_-\) are the forward and backward half-cell phase factors, then

\[a=\frac{y+p_-x}{p_++p_-},\qquad b=x-a.\]

Here \(a\) travels in the receiver’s declared direction and \(b\) travels in the opposite direction. With port directions defined into the device, the single-source scattering result is

\[S_{j m,\,1 n}(f)=\frac{b_{j,m}(f)}{a_{1,n}(f)}.\]

Consequently the explicitly numbered source port gives S11 and a downstream multimode port gives one S21 result for every requested destination mode. These are generalized modal-amplitude ratios. reference_basis_valid records only pre-solve reference eligibility. coefficient_valid then marks each coefficient that survives both the electric and magnetic conditioned solves, finite half-cell phase reconstruction, and finite-value checks. coefficient_valid_S additionally requires a usable source coefficient and a -60 dB incident-spectrum floor evaluated separately for power-wave and generalized-only source bins. Bins remain present in the arrays, but unusable S entries are NaN. power_wave_valid further requires destination-mode power-wave support and a valid destination power matrix. power_wave_valid_S includes those destination gates and additionally requires the launched source mode and its power matrix to be physically valid.

For each electric or magnetic Gram matrix, a singular value is retained only when

\[\sigma_i > \max\left( \frac{\epsilon}{10^{-3}}, \frac{\sigma_{\max}}{\kappa_{\mathrm{lim}}} \right), \qquad \kappa_{\mathrm{lim}}=\min\left(10^{10},\frac{10^{-3}}{\epsilon}\right),\]

where \(\epsilon\) is the precision of the stored field components. If all singular values pass, gprMax solves the complete system directly. If not, a truncated full-system SVD fallback is considered only when the active set contains both power-wave and generalized-only coordinates. The discarded right-singular subspace must have Frobenius projection no larger than \(10^{-3}\) onto the power-wave coordinates; otherwise the complete solve is rejected. Among accepted fallbacks, only coordinates whose ambiguity in the discarded subspace is no larger than \(10^{-3}\) in both the electric and magnetic solves become coefficient_valid. This can preserve a power-wave coordinate when only a generalized-only coordinate is singular, but it rejects a nullspace that mixes physical modes.

condition_number reports the larger electric/magnetic full-system condition number for a direct solve. For a successful truncated fallback it instead reports the larger condition number of the two retained singular subspaces, not the singular original system; it is infinite when no coefficient survives. The small Gram systems and SVD are evaluated in complex128 even when the stored FDTD fields and Gram entries use complex64.

A below-cutoff mode can therefore produce finite incident/outgoing coefficients and a continuous S21 amplitude. In a uniform guide its forward amplitude varies as \(\exp(-\alpha L)\) for \(\beta=-j\alpha\); this is the attenuation of a field coefficient, not the transport of real power by an isolated evanescent wave. At exact cutoff, \(\beta=0\) and the true forward/backward eigenmode basis coalesces. An E/H-balanced reference may approach a finite limiting coefficient there, but the directional decomposition is not unique and must be treated as conditioning-sensitive.

The same Gram matrices define the Hermitian forward-wave power form

\[W=\frac{1}{2}\left(G^E+G^H\right),\qquad P(c)=\operatorname{Re}\{c^\mathrm{H}Wc\}.\]

The implementation symmetrizes \(W\) against round-off and checks that it is finite and positive semidefinite. Keeping the off-diagonal terms is essential for degenerate, nearly degenerate, or merely non-orthogonal finite-grid profiles. Individual coefficient magnitudes therefore are not additive modal powers. A pure evanescent mode has zero independent real power, so coefficient_valid_S may be true while power_basis_valid and power_wave_valid_S are false. In particular, \(|S|^2\) must not be used as an evanescent power fraction.

For net accepted power, let \(x=a+b\) be the total electric coefficient and \(y=a-b\) the co-located total magnetic coefficient after the half-cell correction. The direct time-average flux is

\[P_{\mathrm{acc},p} =\operatorname{Re}\{y_p^\mathrm{H}G^E_p x_p\}.\]

This reduces to \(P(a_p)-P(b_p)\) when \(G^E_p\) is Hermitian. For a lossy port, its anti-Hermitian part supplies a forward/backward interference term which must be retained.

At each frequency this quadratic form is evaluated on the valid propagating mode subspace. All off-diagonal terms within that subspace are retained, but generalized-only rows and columns are excluded. An invalid propagating mode invalidates the accepted-power result rather than being silently omitted.

For active port \(p\), let \(D_p\) select its explicitly driven modes. The externally driven incident power is

\[P_{\mathrm{inc}} = \sum_{p\,\mathrm{active}} \operatorname{Re}\!\left\{ a_{p,D_p}^{\mathrm{H}}W_{p,D_pD_p}a_{p,D_p}\right\}.\]

Thus a one-channel run retains the previous single-mode definition, while a simultaneous multimode run keeps the cross terms between non-orthogonal driven modes on the same physical port. Passive modal receivers have zero generator incident power, but their signed accepted power remains in the multiport balance used for gain. This distinction makes realized gain use launched source power while gain uses the net power accepted by the radiating structure. The power adapter applies the power-normalization and power-matrix masks, so generalized below-cutoff coefficients do not enter gain, accepted-power, or energy-balance normalization.

Understanding Lossy-Mode Results

For a passive lossy mode, users should expect:

• positive Re(n_eff) for the selected forward phase branch;

• negative Im(n_eff);

• positive modal power in the requested source direction;

• a forward magnitude factor abs(exp(-j*k0*n_eff*d)) < 1 for d > 0;

• a nonzero complex-profile residual in general, causing I/Q injection;

• downstream attenuation to be produced by the FDTD material updates after the mode has been launched.

The port-mode field plots show tangential E and H vectors for every retained anchor and report the complex effective index. E and H magnitudes are normalised independently in their panels. A nonzero receiver field alone does not validate a lossy mode: sign-sensitive validation must also check Im(n_eff), modal power direction, and forward attenuation.

Accuracy and convergence

The practical accuracy checklist is in Eigenmode ports and S-parameter analysis. Material dispersion is sampled at each anchor; interpolation between anchors remains piecewise linear. A frequency-dependent bulk pole response is not generally identical to the FDTD ADE transfer. Check convergence in time and space, as well as anchor spacing, rather than treating modal overlap as a sufficient accuracy certificate. Significant DC and Nyquist waveform bins are warned about and discarded.

Complete matrices, active reflection, and antenna normalization

For independent cases, assemble the measured incident and outgoing columns into \(A(f)\) and \(B(f)\). The complete matrix satisfies

\[B = S A.\]

Solving this system accounts for measured incident waves at nominally passive channels. A missing or ill-conditioned incident basis invalidates the affected matrix bins. For reciprocal linear, time-invariant materials and consistent modal normalization, phase conventions, and reference planes, \(S_{ij}=S_{ji}\). Reciprocity alone does not require \(S_{11}=S_{22}\); that equality in example 4 additionally follows from the device’s symmetry.

A simultaneous drive q has spectral multiplier

\[w_q(f) = A_q\exp(j\phi_q)\exp(-j2\pi f\tau_q).\]

The API takes \(\phi_q\) in degrees and \(\tau_q\) in seconds. For a coherent incident vector \(a\), active reflection is

\[\Gamma_{\mathrm{active},q} = \frac{b_q}{a_q} = \frac{\sum_j S_{qj}a_j}{a_q}.\]

Its coefficient and power-wave masks also require an adequately excited denominator. It depends on the incident vector, so it cannot replace the independent cases used to determine \(S\).

For the array along y in example 5, spacing d and constant progressive phase \(\Delta\phi\) give the ideal forward xy-plane maximum

\[\sin\phi_{\mathrm{peak}}(f) = -\frac{\Delta\phi}{k(f)d},\qquad k(f)=\frac{2\pi f}{c}.\]

The phase increment is in radians in this equation. Positive observation angle is measured from +x toward +y. For d = 18 mm and an increment of -108 degrees, the predicted angles at 8, 9, 10, 11, and 12 GHz are approximately 38.7, 33.7, 30.0, 27.0, and 24.6 degrees. These ideal array-factor angles do not include the aperture element pattern or mutual coupling.

For radiation intensity U, total radiated power \(P_{\mathrm{rad}}\), accepted power \(P_{\mathrm{acc}}\), and externally driven incident power \(P_{\mathrm{inc}}\), the antenna quantities are

\[D = \frac{4\pi U}{P_{\mathrm{rad}}},\qquad G = \frac{4\pi U}{P_{\mathrm{acc}}},\qquad G_{\mathrm{realized}} = \frac{4\pi U}{P_{\mathrm{inc}}}.\]

The power adapter uses the physical power-wave subspace and its full power matrix. Finite below-cutoff coefficients do not enter gain or energy-balance normalization. The port_power/modal_ports output repeats the modal waves, mode indices, power matrix, and physical masks at the NTFF frequencies so that the normalization can be audited.

Limitations

Automatic mode tracking is experimental and requires further testing. Enable it with tracking="auto" at your discretion; tracking="legacy" remains the default. The additional matching and diagnostic evidence does not remove the following numerical and physical limitations.

• Material tensors are diagonal in the local u/v/w basis.

• The finite-difference operators use first-order sparse Yee-grid differences.

• Bulk dispersive material poles use their analytic physical-frequency response. Matching their exact FDTD ADE transfer remains a separate step beyond the temporal and longitudinal difference compensation.

• A single-frequency source reuses one solved profile across the waveform spectrum. Broadband interpolation remains sensitive to anchor spacing even when tracking confidently identifies the same branch. It does not certify modal-impedance accuracy or a particular S11 floor.

• Automatic tracking can follow changes in eigenvalue order and transport detected degenerate subspaces. Extra candidates and adaptive frequency solves cannot guarantee a unique match in every guide. Unresolved identity gaps invoke the automatic-anchor fallback or an explicit-anchor error; interpolation never bridges them. Resolved eigenvalue splitting retains separate eigenmodes and propagation constants, rather than mixing the fields.

• Confinement and artificial-boundary diagnostics are evidence, not proofs that a mode is physical or spurious. Numerically valid, confidently tracked forward-power profiles remain eligible with warnings when confinement is suspect or unresolved. Failed residuals, nonfinite fields, singular reconstruction, and rank loss remain disqualifying. Evanescent references do not become injectable merely because their transverse fields are bound.

• Full verification concerns the represented voxel geometry. It requests twice the transverse resolution and two larger, PML-free windows for open cross-sections. Missing exterior geometry, unsupported faithful refinement of thin sheets or impedance boundaries, failed solves, and exhausted solve budgets leave diagnostics unresolved. Larger-window verification is currently unavailable with domain-decomposed MPI. verification="fast" omits these extra solves and provides only residual and edge evidence.

See Eigenmode tracking theory for the assignment and diagnostic checks, and Examples 8 and 9 in Eigenmode tracking guide for circular degeneracy and a crossing in a single anisotropic guide.

Tracking and basis selection are developed in Eigenmode tracking theory. See Eigenmode tracking guide for their input arguments and illustrated examples.