Activate magnetoelectric coupling (ρ1/ρ2) for bianisotropic media

[garrekstemo]

Summary

Add a per-Layer API for the two magnetoelectric coupling tensors (ρ₁, ρ₂ — the off-diagonal 3×3 blocks of the 6×6 constitutive matrix), wire the currently-dead construct_M(ε, μ, ρ1, ρ2) path into the per-layer pipeline, and add bianisotropic / non-reciprocal test cases. This enables modeling of chiral, Tellegen, and other bianisotropic media.

Context

The 6×6 material builder already accepts magnetoelectric coupling: construct_M(ε, μ=…, ρ1=…, ρ2=…) places ρ₁ in the upper-right block and ρ₂ in the lower-left block, returning [ε ρ1; ρ2 μ] (src/matrix_constructors.jl:128-135). This method is unit-tested (test/layer.jl:69-97, test/functions.jl:219-231 build M with non-zero ρ1/ρ2 and assert the blocks land correctly) — but it is never called from the per-layer pipeline. layer_matrices always calls the no-coupling specialized method construct_M(ε, μ) (src/matrix_constructors.jl:117), so ρ blocks are always zero in practice. There is no user-facing way to set ρ₁/ρ₂, so bianisotropic media cannot currently be computed end-to-end.

The interface/Δ stages (construct_a :176-197, construct_Δ :215-243) read generic M[i,j] entries and would pick up non-zero ρ blocks. However, the eigenvector stage calculate_γ (src/matrix_constructors.jl:331-400) is derived for ε/μ only with scalar μ and no magnetoelectric coupling — it builds γ from ε and μ alone and does not see the ρ blocks. So simply un-suppressing the ρ blocks in construct_M is not sufficient: the eigenvector/propagation construction must handle the coupled 6×6 material, which the present γ-based path cannot do.

Why

Bianisotropic media — chiral materials, Tellegen media, moving media, certain metamaterials — are defined precisely by non-zero magnetoelectric coupling. The constitutive infrastructure (construct_M with ρ blocks, tested) is already in place but inert; activating it unlocks a whole material class.

Proposed approach

1. Add a per-Layer API to specify ρ₁ and ρ₂ (3×3 complex tensors, possibly λ-dependent), alongside ε (and μ from the companion μ-tensor issue). Extend Layer (src/layer.jl:48-97) and its predicates accordingly.
2. Route layers carrying ρ blocks through the generic construct_M(ε, μ, ρ1, ρ2) path (src/matrix_constructors.jl:133-135) instead of the no-coupling specialized method called at src/matrix_constructors.jl:117.
3. Replace the γ-eigenvector + dynamical_matrix propagation (src/matrix_constructors.jl:331-426) for coupled layers with a propagator that handles the full coupled 6×6 material — the matrix-exponential propagator from the companion spike issue is the intended mechanism, since the singularity-free γ forms (Xu et al. 2000) assume no magnetoelectric coupling. Re-derive the coupled propagation from Mackay & Lakhtakia (2020); do not transcribe coefficient expressions from any third-party implementation without independent verification.
4. Add validation tests: (a) a reciprocal medium (ρ₂ = -ρ₁ᵀ for a chiral/Pasteur medium) showing the expected optical rotation / circular dichroism; (b) a non-reciprocal (Tellegen-type) case demonstrating distinct forward/backward response; (c) the ρ = 0 limit reproduces the existing no-coupling results exactly (regression guard).
5. Document the bianisotropic constructor and add it to the API summary in CLAUDE.md.

Acceptance criteria

• A Layer can carry non-zero ρ₁/ρ₂; transfer/sweep_*/efield honor them.
• The ρ = 0 limit reproduces current results to machine precision (regression test).
• A chiral (reciprocal) test case shows the expected polarization rotation / circular dichroism via the cross-polarization terms of TransferResult (src/general_TMM.jl:99-108).
• A non-reciprocal test case demonstrates forward ≠ backward response.
• Existing test suite stays green.

Caveats / risks

• Depends on the matrix-exponential propagator issue (or an equivalent generalized-eigenvector propagator). The existing γ-based path cannot represent magnetoelectric coupling, so this issue should not start until that propagator exists — see the companion "Spike + implement a matrix-exponential layer propagator backend" issue. Captured from the original scope assessment's own caveat that ρ1/ρ2 needs the exponential propagator.
• Magnetoelectric coefficient conventions and units vary across the literature (Gaussian vs SI, the placement and sign of ρ₁ vs ρ₂); fix and document the convention up front and validate against an analytic chiral-slab result.
• Re-derive all bianisotropic coefficients from Mackay & Lakhtakia (2020) and validate against analytic / known cases — do not transcribe coefficient expressions from any third-party implementation without independent verification.
• Energy/reciprocity checks differ for non-reciprocal media (R+T=1 still holds for lossless media, but forward/backward symmetry does not); use appropriate validation rather than assuming the existing reciprocal checks.

References

• Mackay & Lakhtakia, The Transfer-Matrix Method in Electromagnetics and Optics, Synthesis Lectures on Electromagnetics, Springer (2020). DOI 10.1007/978-3-031-02022-3 (bianisotropic constitutive relations and propagation)
• Berreman, J. Opt. Soc. Am. 62, 502 (1972). https://doi.org/10.1364/JOSA.62.000502 (4×4 formalism this code builds on)
• Xu et al., Phys. Rev. B 61, 1740 (2000). https://doi.org/10.1103/PhysRevB.61.1740 (current γ eigenvectors — note they assume no magnetoelectric coupling)