fix(modes): bound |n_eff| to reject spurious metal modes

src/gsim/common/modes.py

@@ -292,13 +292,15 @@ def propagating_modes[ModeT](
     min_index: float = 1.0,
     max_loss_ratio: float = 1.0,
     max_gain_ratio: float = MAX_GAIN_RATIO,
+    max_index: float | None = None,
 ) -> list[ModeT]:
     """Keep the Modes that propagate, dropping evanescent and spurious ones.
 
     A Mode is kept when it is guided (``Re(n_eff)`` above *min_index*),
     loses no more than *max_loss_ratio* per unit length than it advances
     in phase, and does not grow: ``Im(n_eff)`` positive beyond
-    *max_gain_ratio* is gain, which a passive line cannot have.
+    *max_gain_ratio* is gain, which a passive line cannot have. With
+    *max_index* set, a Mode whose ``|n_eff|`` exceeds it is dropped too.
 
     Args:
         modes: Every solved Mode.
@@ -311,6 +313,13 @@ def propagating_modes[ModeT](
         max_gain_ratio: Upper bound on ``+Im(n_eff) / Re(n_eff)``, the
             room left for a lossless Mode's numerical noise (default:
             :data:`MAX_GAIN_RATIO`).
+        max_index: Upper bound on ``|n_eff|``, or ``None`` for no bound
+            (default). A solve asked for many Modes around a lossy
+            conductor also returns Modes living inside the metal, with
+            ``|n_eff|`` near ``|sqrt(eps_metal)|``, some of them inside
+            the loss bound; the line Mode sits orders of magnitude below.
+            ``gsim.femwell.conductor_index_bound`` derives a bound from
+            the permittivity the solve was given.
 
     Returns:
         The propagating Modes, in the order they were solved.
@@ -321,6 +330,7 @@ def propagating_modes[ModeT](
         if (
             n_eff.real > min_index
             and -max_loss_ratio * n_eff.real < n_eff.imag <= max_gain_ratio * n_eff.real
+            and (max_index is None or abs(n_eff) <= max_index)
         ):
             kept.append(mode)
     return kept
@@ -340,6 +350,7 @@ def select_line_mode[ModeT](
     rule: LineModeRule | None = None,
     min_index: float = 1.0,
     max_loss_ratio: float = 1.0,
+    max_index: float | None = None,
     degeneracy_rtol: float = 0.03,
 ) -> ModeT:
     """Select the physical line Mode from a set of solved Modes.
@@ -365,6 +376,14 @@ def select_line_mode[ModeT](
         min_index: Lower bound on ``Re(n_eff)`` for the default rule.
         max_loss_ratio: Upper bound on ``|Im(n_eff)| / Re(n_eff)`` for
             the default rule.
+        max_index: Upper bound on ``|n_eff|`` for the default rule, or
+            ``None`` for no bound (default). The default rule picks the
+            highest ``Re(n_eff)``, so a spurious Mode inside a lossy
+            conductor, which carries an index near
+            ``|sqrt(eps_metal)|``, wins over the line Mode when the solve
+            returns one; set this when solving more than a few Modes
+            around such a conductor (``gsim.femwell.conductor_index_bound``
+            gives a bound from the solve's permittivity).
         degeneracy_rtol: Relative spread in ``Re(n_eff)`` within which two
             candidates count as ambiguous and a warning is issued.
 
@@ -381,7 +400,12 @@ def select_line_mode[ModeT](
             choice between them is not physically meaningful.
     """
     candidates: Sequence[ModeT] = (
-        propagating_modes(modes, min_index=min_index, max_loss_ratio=max_loss_ratio)
+        propagating_modes(
+            modes,
+            min_index=min_index,
+            max_loss_ratio=max_loss_ratio,
+            max_index=max_index,
+        )
         if rule is None
         else list(rule(modes))
     )

src/gsim/femwell/__init__.py

@@ -13,6 +13,7 @@
 
 from gsim.femwell.adapter import (
     boundary_field_ratio,
+    conductor_index_bound,
     electrode_current,
     elementwise_epsilon,
     epsilon_by_region,
@@ -25,6 +26,7 @@
 
 __all__ = [
     "boundary_field_ratio",
+    "conductor_index_bound",
     "electrode_current",
     "elementwise_epsilon",
     "epsilon_by_region",

src/gsim/femwell/adapter.py

@@ -24,7 +24,9 @@
 Marks-Williams power-current impedance it implies
 (:func:`z0_power_current`), and — for an RF line — the one reading of a
 selected Mode an RF line solve asks for: index, impedance and whether it is
-the wall Mode (:func:`line_reading`). A conductor is named to every
+the wall Mode (:func:`line_reading`). :func:`conductor_index_bound` reads, from the
+permittivity a solve was given, the index above which a Mode lives inside a
+conductor rather than on the line. A conductor is named to every
 current integral the same way, as a
 :class:`~gsim.common.modes.Conductor`.
 
@@ -34,6 +36,7 @@
 
 from __future__ import annotations
 
+import warnings
 from pathlib import Path
 from typing import TYPE_CHECKING, Any, cast
 
@@ -211,6 +214,93 @@ def elementwise_epsilon(
     return np.asarray(sampled, dtype=np.complex128)
 
 
+def conductor_index_bound(
+    epsilon: dict[str, complex] | ArrayLike,
+    *,
+    fraction: float = 0.1,
+    margin: float = 10.0,
+) -> float | None:
+    """Upper bound on ``|n_eff|`` that rejects Modes living inside a conductor.
+
+    A mode solve asked for more than a few Modes around a lossy
+    conductor also returns Modes confined to the metal, with ``|n_eff|``
+    close to ``|sqrt(eps_metal)|`` (about 2.6e4 for aluminium at 1 GHz).
+    Some lose less per unit length than they advance, so the loss bound
+    keeps them, and the slowest-Mode rule of
+    :func:`~gsim.common.modes.select_line_mode` picks one over the line
+    Mode. The line Mode's index is orders of magnitude lower: over a good
+    conductor ``n_line**2 / eps_metal`` is about
+    ``eps_d R' eps0 / (L' sigma)``, independent of frequency and of
+    order 1e-7 for an aluminium strip in oxide. Pass the bound returned
+    here as ``max_index`` to ``select_line_mode``.
+
+    A region or element is a conductor when conduction dominates its
+    permittivity, ``|Im(eps)| > |Re(eps)|``. That alone does not make
+    it a safe source for the bound: a weak conductor, such as moderately
+    doped silicon, has ``|sqrt(eps)|`` close to the indices a line Mode
+    can reach. A conductor sets the bound only when
+    ``fraction * |sqrt(eps)|`` is at least *margin* times the largest
+    ``|sqrt(eps)|`` of the other regions (and above the light line);
+    the others are skipped with a warning. A slow-wave line Mode can sit
+    several times above the dielectric index (3.5 times on an aluminium
+    strip in oxide at 1 GHz), hence the default margin of 10.
+
+    Args:
+        epsilon: The permittivity given to :func:`solve_modes`: either
+            ``{region_name: eps}`` or per-element values.
+        fraction: Bound as a fraction of a conductor's ``|sqrt(eps)|``
+            (default 0.1; spurious metal Modes have been seen from 0.9
+            times that value up).
+        margin: Smallest ratio between a conductor's bound and the
+            largest non-conductor index for the conductor to count
+            (default 10).
+
+    Returns:
+        The smallest bound over the conductors that count, or ``None``
+        when none does (no bound is then given from this source).
+
+    Raises:
+        ValueError: If *fraction* or *margin* is not positive.
+
+    Warns:
+        UserWarning: When a conductor is skipped because its bound would
+            not clear the dielectric index by *margin*.
+    """
+    if fraction <= 0.0:
+        raise ValueError(f"fraction must be positive, got {fraction}.")
+    if margin <= 0.0:
+        raise ValueError(f"margin must be positive, got {margin}.")
+    values = (
+        np.fromiter(epsilon.values(), dtype=np.complex128)
+        if isinstance(epsilon, dict)
+        else np.asarray(epsilon, dtype=np.complex128).ravel()
+    )
+    is_conductor = np.abs(values.imag) > np.abs(values.real)
+    if not is_conductor.any():
+        return None
+    dielectric = values[~is_conductor]
+    dielectric_index = (
+        float(np.sqrt(np.abs(dielectric)).max()) if dielectric.size else 1.0
+    )
+    floor = max(margin * dielectric_index, 1.0)
+    bounds = fraction * np.sqrt(np.abs(values[is_conductor]))
+    usable = bounds >= floor
+    if not usable.all():
+        skipped = np.unique(values[is_conductor][~usable])
+        shown = ", ".join(f"{eps:.4g}" for eps in skipped[:3])
+        more = f" and {skipped.size - 3} more" if skipped.size > 3 else ""
+        warnings.warn(
+            f"conductor_index_bound skipped conductor permittivity {shown}{more}: "
+            f"{fraction:g} * |sqrt(eps)| does not reach {margin:g} times the "
+            f"largest dielectric index ({dielectric_index:.4g}), so it cannot "
+            "tell a Mode inside it from the line Mode.",
+            stacklevel=2,
+        )
+    if not usable.any():
+        return None
+    return float(bounds[usable].min())
+
+
 def solve_modes(
     msh_path: str | Path,
     *,
@@ -252,6 +342,9 @@ def solve_modes(
             for RF materials with big conductive ``|Im(eps)|`` can put the
             shift far from the physical quasi-TEM mode; pass an explicit
             guess (e.g. the expected slow-wave index) there.
+            Asked for more than a few Modes around a lossy conductor,
+            the solve also returns Modes confined to the metal; select
+            the line Mode with ``max_index=conductor_index_bound(epsilon)``.
 
     Returns:
         The femwell ``Modes`` result (each mode carries ``n_eff``).

tests/common/test_mode_selection.py

@@ -130,6 +130,45 @@ def test_a_bound_nothing_survives_is_reported(self):
             select_line_mode([3.0 - 2.0j], max_loss_ratio=0.5)
 
 
+class TestIndexBound:
+    """A Mode living inside a lossy conductor is not the line Mode (#498).
+
+    Asked for 8 to 16 Modes on an aluminium strip in oxide, the solve also
+    returns Modes with ``|n_eff|`` near ``|sqrt(eps_metal)|`` (about 2.6e4
+    at 1 GHz), some inside the loss bound. The slowest-Mode rule then picks
+    one of them over the line Mode, 5.13-4.68j on that geometry.
+    """
+
+    LINE = 5.13 - 4.68j
+    METAL = 23174.0 - 12017.0j  # -Im/Re = 0.52, inside the loss bound
+
+    def test_without_a_bound_the_metal_mode_wins(self):
+        """Unchanged default: the bound is the caller's to give."""
+        assert select_line_mode([mode(self.LINE), mode(self.METAL)]).n_eff == (
+            self.METAL
+        )
+
+    def test_the_bound_drops_the_metal_mode(self):
+        line = mode(self.LINE)
+        picked = select_line_mode([line, mode(self.METAL)], max_index=2.6e3)
+        assert picked is line
+
+    def test_the_bound_reaches_the_candidates(self):
+        assert propagating_modes([self.METAL, self.LINE], max_index=2.6e3) == [
+            self.LINE
+        ]
+
+    def test_the_bound_is_on_the_magnitude(self):
+        """A Mode at the bound stays; one just above it on |n_eff| goes."""
+        at_bound = 4.0 - 3.0j  # |n_eff| = 5, inside the loss bound
+        assert propagating_modes([at_bound], max_index=5.0) == [at_bound]
+        assert propagating_modes([at_bound], max_index=4.99) == []
+
+    def test_a_bound_nothing_survives_is_reported(self):
+        with pytest.raises(NoLineModeError, match="No propagating line mode"):
+            select_line_mode([self.METAL], max_index=2.6e3)
+
+
 class TestGainBound:
     """A Mode that grows along the line is a numerical artefact, not a candidate.
 

tests/femwell/test_adapter.py

@@ -6,9 +6,11 @@
 import numpy as np
 import pytest
 
+from gsim.common.modes import select_line_mode
 from gsim.common.stack.extractor import Layer, LayerStack
 from gsim.common.stack.materials import MaterialProperties
 from gsim.femwell.adapter import (
+    conductor_index_bound,
     elementwise_epsilon,
     epsilon_by_region,
 )
@@ -144,3 +146,66 @@ def test_user_override_wins_over_database(self, tmp_path):
             overrides={"si": MaterialProperties(permittivity=9.0)},
         )
         assert eps["core"].real == pytest.approx(9.0)
+
+
+class TestConductorIndexBound:
+    """The |n_eff| above which a Mode lives inside a conductor (#498)."""
+
+    def test_a_tenth_of_the_metal_index(self):
+        eps_metal = 1.0 - 6.8e8j  # aluminium at 1 GHz
+        bound = conductor_index_bound({"metal": eps_metal, "oxide": 3.9 + 0j})
+        assert bound == pytest.approx(0.1 * abs(np.sqrt(eps_metal)))
+
+    def test_the_least_conductive_strong_conductor_sets_it(self):
+        bound = conductor_index_bound(
+            {"al": 1.0 - 6.8e8j, "w": 1.0 - 3.3e8j, "oxide": 3.9 + 0j}
+        )
+        assert bound == pytest.approx(0.1 * abs(np.sqrt(1.0 - 3.3e8j)))
+
+    def test_a_weak_conductor_alone_gives_no_bound(self):
+        """11.7-12j conducts, but 0.1|sqrt(eps)| = 0.41 is below the light line."""
+        with pytest.warns(UserWarning, match="skipped conductor"):
+            assert conductor_index_bound({"si": 11.7 - 12j, "oxide": 3.9 + 0j}) is None
+
+    def test_a_weak_conductor_does_not_lower_the_metal_bound(self):
+        """Doped silicon next to aluminium leaves the aluminium bound in place."""
+        eps = {"al": 1.0 - 6.8e8j, "si": 11.7 - 1000j, "oxide": 3.9 + 0j}
+        with pytest.warns(UserWarning, match="skipped conductor"):
+            bound = conductor_index_bound(eps)
+        assert bound == pytest.approx(0.1 * abs(np.sqrt(1.0 - 6.8e8j)))
+        line = 5.13 - 4.68j
+        assert select_line_mode([line, 2.5 - 0.01j], max_index=bound) == line
+
+    def test_the_margin_is_against_the_largest_dielectric_index(self):
+        """A conductor clears oxide by 10x but not silicon is skipped."""
+        eps_cond = 1.0 - 9e4j  # 0.1|sqrt(eps)| = 30
+        assert conductor_index_bound({"c": eps_cond, "ox": 3.9 + 0j}) == (
+            pytest.approx(30.0, rel=1e-4)
+        )
+        with pytest.warns(UserWarning, match="skipped conductor"):
+            assert (
+                conductor_index_bound({"c": eps_cond, "si": 11.7 + 0j, "ox": 3.9 + 0j})
+                is None
+            )
+
+    def test_per_element_values(self):
+        eps = np.array([3.9, 3.9, 1.0 - 6.8e8j, 3.9], dtype=complex)
+        assert conductor_index_bound(eps) == pytest.approx(
+            0.1 * abs(np.sqrt(1.0 - 6.8e8j))
+        )
+
+    def test_a_lossy_dielectric_is_not_a_conductor(self):
+        """Loss below the real permittivity leaves the Mode free to sit there."""
+        assert conductor_index_bound({"si": 11.7 - 5.0j, "oxide": 3.9 + 0j}) is None
+
+    def test_the_fraction_is_a_parameter(self):
+        eps = {"metal": 1.0 - 6.8e8j}
+        default = conductor_index_bound(eps)
+        assert default is not None
+        assert conductor_index_bound(eps, fraction=0.5) == pytest.approx(5 * default)
+
+    def test_the_fraction_and_margin_must_be_positive(self):
+        with pytest.raises(ValueError, match="fraction"):
+            conductor_index_bound({"metal": 1.0 - 6.8e8j}, fraction=0.0)
+        with pytest.raises(ValueError, match="margin"):
+            conductor_index_bound({"metal": 1.0 - 6.8e8j}, margin=0.0)