Subpixel averaging for dispersive (anisotropic) susceptibilities

[oskooi]

Status: proposal / design note
Area: src/anisotropic_averaging.cpp, src/susceptibility.cpp, src/update_pols.cpp
Motivation: permanently fix the anisotropic-Lorentzian boundary instability that #666 worked around with a runtime guard, and remove a known accuracy limitation.

Background: the #666 guard

lorentzian_susceptibility::update_P (src/susceptibility.cpp) guards its anisotropic (2×2 and 3×3) branches with a per-point test that skips the polarization update where the diagonal sigma is zero:

meep/src/susceptibility.cpp

Lines 229 to 238 in ea11e6c

	// s[i] != 0 check is a bit of a hack to work around
	// some instabilities that occur near the boundaries
	// of materials; see PR #666
	if (s[i] != 0) {
	realnum pcur = p[i];
	p[i] = gamma1inv * (pcur * (2 - omega0dtsqr_denom) - gamma1 * pp[i] +
	omega0dtsqr * (s[i] * w[i] + OFFDIAG(s1, w1, is1, is) +
	OFFDIAG(s2, w2, is2, is)));
	pp[i] = pcur;
	}

It exists to suppress an instability "near the boundaries of materials."

Root cause (empirically established)

A controlled 2D TE reproducer (anisotropic Lorentzian medium, Gaussian pulse, no PML, watch total energy) isolates the trigger precisely:

Off-diagonal sigma support vs. diagonal support	No-guard result
Off-diagonal exceeds diagonal	diverges
Off-diagonal strictly inside diagonal	stable (decays)

With the guard, the diverging case becomes stable, and for the well-posed case the guard is a no-op (bit-identical results).

The mechanism:

1. The off-diagonal coupling in update_P uses the 4-point OFFDIAG averaging stencil, which reads sigma/field at i and i+is (one cell along the component's own direction).
2. The dispersive sigma is point-sampled in structure_chunk::add_susceptibility (src/anisotropic_averaging.cpp) — it is not subpixel-averaged. At a material boundary the off-diagonal sigma is therefore nonzero one stencil-step beyond where the diagonal sigma has dropped to zero.
3. At such a point P_c is driven with no diagonal sigma to anchor it. The local discrete susceptibility tensor is effectively [[0, σ], [σ, 0]], whose eigenvalues are ±σ — not positive-semidefinite, i.e. non-passive — so it pumps energy and the simulation diverges.

By contrast, the instantaneous permittivity chi1inv is built with subpixel (Kottke) averaging (structure_chunk::set_chi1inv → material_function::eff_chi1inv_row), which yields a smooth, PSD-consistent tensor whose diagonal and off-diagonal taper together at boundaries. That path is stable. The dispersive sigma lacks this treatment — that is the root cause.

A naive construction-time fix ("zero off-diagonal sigma where the diagonal is zero, same pixel") was prototyped and fails: a diagnostic showed it touches zero points, because the mismatch lives one stencil-step (i+is) away, not at the same pixel.

Proposed fix

Give the dispersive susceptibility tensor the same subpixel averaging the instantaneous tensor already receives, so diagonal and off-diagonal sigma taper together (PSD-preserving) and the off-diagonal can never "stick out" past the diagonal. Then the update_P guard becomes unnecessary and can be removed, and dispersive interfaces gain the same second-order accuracy that non-dispersive interfaces already have.

Approach

1. Reuse the Kottke machinery. eff_chi1inv_row (and the libctl/MPB-backed override) computes the effective inverse-permittivity tensor row for a pixel from the material normal vector and the harmonic/arithmetic mean across the interface. The dispersive sigma tensor should be averaged with the same normal-vector + projection construction so the smoothed tensor is symmetric PSD by construction (the projection from P·a + (I-P)·b preserves PSD when a,b ≥ 0).
2. New material-function hook. Add an eff_sigma_row(component c, double sigrow[3], const volume &v, ...) analogous to eff_chi1inv_row, returning the subpixel-averaged sigma row for the pixel volume v. Default implementation: numerical surface/volume average mirroring eff_chi1inv_row: the libctl interface (python/meep/geom.py / libctlgeom) can override with a semi-analytic average, as it already does for eff_chi1inv_row.
3. Wire it into add_susceptibility. In structure_chunk::add_susceptibility (src/anisotropic_averaging.cpp), replace the point-sampling loop (the sigma_row calls at here and here - shift1) with calls to eff_sigma_row over the pixel volume when subpixel averaging is enabled, falling back to point sampling otherwise (for backward-compatible, averaging disabled runs).
4. Remove the guard. Once sigma is PSD-consistent at boundaries, delete the s[i] != 0 guard from lorentzian_susceptibility::update_P. Keep the non-PSD-tensor warning (added in a separate PR) for genuinely non-passive user input.

Subtleties

• Per-direction storage. sigma[c][d] stores one scalar per (component, direction); the averaged tensor must be expressible in this staggered layout. The diagonal and off-diagonal are sampled at different Yee points, so the averaging must respect the existing half-pixel shift1 convention used for off-diagonals.
• Drude / multi-pole / gyrotropic. The averaging applies to the sigma prefactor only (the per-pole omega_0, gamma are scalars), so Lorentzian, Drude, and noisy-Lorentzian share the treatment. Gyrotropic media use a separate gyro_tensor and a different update; out of scope.
• Backward compatibility. Results for anisotropic dispersive media at boundaries will change (they become more accurate). Gate behind the existing use_anisotropic_averaging flag so users can reproduce old behavior; document the change.

Validation plan.

1. Stability: the reproducer above (off-diagonal exceeding diagonal) must remain stable without the guard once averaging is on.
2. No regression: python/tests/test_faraday_rotation.py, test_material_dispersion.py, test_dispersive_eigenmode.py, test_multilevel_atom.py; C++ tests/pml.cpp (offdiag case).
3. Accuracy / convergence: add a convergence test for an anisotropic dispersive interface (e.g. transmission through a rotated Lorentzian slab) showing second-order error vs. resolution with averaging on, first-order with it off — mirroring the non-dispersive subpixel tests.
4. Energy conservation: a closed cavity with a lossless (gamma=0) anisotropic dispersive inclusion should conserve energy (no secular growth) without the guard.

Effort / risk

Substantial: it introduces a new averaging path for dispersive materials (a long-standing Meep limitation, not just a bug fix), touches the C++ core and the libctl/Python material interface, and changes numerical results for anisotropic dispersive media. Recommend landing behind the averaging flag with the convergence + energy tests above, and keeping the guard until averaging is validated to make it removable.

Near-term status

Until this lands, the repository ships:

• the branchless guard (to be added in a separate PR; efficient, correct passitivity enforcement; no-op for well-posed media), and
• a non-PSD sigma warning in structure::add_susceptibility that flags genuinely non-passive user tensors (which neither the guard nor subpixel averaging can fix).

[stevengj]

It doesn't make sense to me to do the more complicated "Kottke" averaging here, as opposed to simple fill-fraction averaging, since Kottke averaging of σ alone (not taking into account the ω-dependent denominator) will not increase accuracy.

A separate issue is that the timestepping looks at diagonal and off-diagonal σ components from different points in space (different Yee-grid points for different components), so I'm not sure that subpixel averaging alone will solve the issue.

If you want to play around with this empirically, try some test case that first reproduces the problem without the guard — e.g. a cylinder (in 2d) or a sphere (in 3d) of a polarizable material with an anisotropic σ surrounded by vacuum.