fix(mode): ensure degenerate modes are orthogonal from mode solver

CHANGELOG.md

@@ -16,6 +16,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added optional `max_results` handling to the kLayout `DRCLoader` and `DRCRunner` and speeded up parsing of big result sets.
 
 ### Fixed
+- Improved degenerate mode handling in the mode solver to ensure modes respect the bi-orthogonality condition.
 
 ## [2.10.0rc3] - 2025-11-26
 

tests/test_plugins/test_mode_solver.py

@@ -13,7 +13,7 @@
 from tidy3d import ScalarFieldDataArray
 from tidy3d.components.data.monitor_data import ModeSolverData
 from tidy3d.components.mode.derivatives import create_sfactor_b, create_sfactor_f
-from tidy3d.components.mode.solver import compute_modes
+from tidy3d.components.mode.solver import TOL_DEGENERATE, EigSolver, compute_modes
 from tidy3d.components.mode_spec import MODE_DATA_KEYS
 from tidy3d.exceptions import DataError, SetupError
 from tidy3d.plugins.mode import ModeSolver
@@ -1503,3 +1503,52 @@ def test_sort_spec_track_freq():
     assert np.allclose(modes_lowest.Ex.abs, modes_lowest_retracked.Ex.abs)
     assert np.all(modes_lowest.n_eff == modes_lowest_retracked.n_eff)
     assert np.all(modes_lowest.n_group == modes_lowest_retracked.n_group)
+
+
+def test_degenerate_mode_processing():
+    """Ensure degenerate modes returned by mode solver are bi-orthogonal."""
+    freq0 = td.C_0
+    sim_size = (0, 2, 2)
+    inf = 10
+    W1 = 0.3
+    n = 1.5
+    num_modes = 4
+    mode_spec = td.ModeSpec(num_modes=num_modes)
+    medium = td.Medium(permittivity=n**2)  # , conductivity=0.1)
+    geom1 = td.Box.from_bounds((-inf, -W1 / 2, -W1 / 2), (inf, W1 / 2, W1 / 2))
+    wg1 = td.Structure(geometry=geom1, medium=medium)
+
+    grid_spec = td.GridSpec.uniform(dl=0.02)
+
+    sim = td.Simulation(
+        size=sim_size,
+        structures=[wg1],
+        grid_spec=grid_spec,
+        run_time=freq0 / 10,
+    )
+
+    ms = ModeSolver(
+        simulation=sim,
+        plane=sim.geometry,
+        mode_spec=mode_spec,
+        freqs=[freq0],
+        direction="+",
+    )
+
+    mode_data = ms.data_raw
+
+    degen_sets = EigSolver._identify_degenerate_modes(
+        mode_data.n_complex.values[0, :], TOL_DEGENERATE
+    )
+    assert len(degen_sets) == 1
+
+    S = mode_data.outer_dot(mode_data, conjugate=True).isel(f=0).values
+    threshold = 1e-7
+    off_diag_mask = ~np.eye(S.shape[0], dtype=bool)
+    large_vals = np.abs(S) > threshold
+    problem_mask = off_diag_mask & large_vals
+
+    indices = np.argwhere(problem_mask)
+    msg = f"Found {len(indices)} off-diagonal values > {threshold}:\n"
+    msg += "\n".join(f"  |S[{i},{j}]| = {np.abs(S[i, j]):.4e}" for i, j in indices)
+    assert not np.any(problem_mask), msg

tidy3d/components/mode/solver.py

@@ -18,6 +18,8 @@
 TOL_COMPLEX = 1e-10
 # Tolerance for eigs
 TOL_EIGS = fp_eps / 10
+# Tolerance to consider modes as degenerate
+TOL_DEGENERATE = TOL_EIGS * 10
 # Tolerance for deciding on the matrix to be diagonal or tensorial
 TOL_TENSORIAL = 1e-6
 # shift target neff by this value, both rel and abs, whichever results in larger shift
@@ -285,6 +287,15 @@ def compute_modes(
         E = E.reshape((3, Nx, Ny, 1, num_modes))
         H = np.sum(jac_h[..., None] * H[:, None, ...], axis=0)
         H = H.reshape((3, Nx, Ny, 1, num_modes))
+
+        # Normalize all modes so self-dot equals 1.0
+        all_modes = list(range(num_modes))
+        E, H = cls._normalize_modes(E, H, dl_f, dl_b, all_modes)
+
+        # Identify and orthogonalize degenerate modes
+        degenerate_modes = cls._identify_degenerate_modes(neff + keff * 1j, TOL_DEGENERATE)
+        E, H = cls._make_orthogonal_basis_for_degenerate_modes(degenerate_modes, E, H, dl_f, dl_b)
+
         fields = np.stack((E, H), axis=0)
 
         neff = neff * np.linalg.norm(kp_to_k)
@@ -1089,6 +1100,172 @@ def mode_plane_contain_good_conductor(material_response) -> bool:
             return False
         return np.any(np.abs(material_response) > GOOD_CONDUCTOR_THRESHOLD * np.abs(pec_val))
 
+    @staticmethod
+    def _identify_degenerate_modes(
+        n_complex: np.ndarray,
+        tol: float,
+    ) -> list[set[int]]:
+        """Inspects the n_complex of modes to find sets of degenerate modes."""
+        num_modes = len(n_complex)
+        ungrouped = set(range(num_modes))
+        degenerate_sets = []
+
+        while ungrouped:
+            # Start a new group with an ungrouped mode
+            seed = ungrouped.pop()
+            current_set = {seed}
+
+            # Find all ungrouped modes similar to the seed
+            to_remove = set()
+            for mode_idx in ungrouped:
+                if np.isclose(n_complex[mode_idx], n_complex[seed], rtol=tol, atol=tol):
+                    current_set.add(mode_idx)
+                    to_remove.add(mode_idx)
+
+            # Remove grouped modes from ungrouped set
+            ungrouped -= to_remove
+
+            # Only keep groups with more than one mode
+            if len(current_set) >= 2:
+                degenerate_sets.append(current_set)
+
+        return degenerate_sets
+
+    @staticmethod
+    def _dot(
+        E: np.ndarray,
+        H: np.ndarray,
+        dl_primal: np.ndarray,
+        dl_dual: np.ndarray,
+        mode_1: int,
+        mode_2: int,
+    ) -> complex:
+        """Dot product based on the bi-orthogonality relationship between E and H."""
+        # Extract field components
+        Ex = E[0, ...]
+        Ey = E[1, ...]
+        Hx = H[0, ...]
+        Hy = H[1, ...]
+
+        # Make the differential area elements
+        Ex_Hy_dS = np.outer(dl_primal[0], dl_dual[1])
+        Ey_Hx_dS = np.outer(dl_dual[0], dl_primal[1])
+
+        term1 = Ex[..., mode_1] * Hy[..., mode_2] + Ex[..., mode_2] * Hy[..., mode_1]
+        term1 *= Ex_Hy_dS[..., np.newaxis]
+        term2 = Ey[..., mode_1] * Hx[..., mode_2] + Ey[..., mode_2] * Hx[..., mode_1]
+        term2 *= Ey_Hx_dS[..., np.newaxis]
+        return np.sum(term1 - term2)
+
+    @staticmethod
+    def _outer_dot(
+        E: np.ndarray,
+        H: np.ndarray,
+        dl_primal: np.ndarray,
+        dl_dual: np.ndarray,
+        mode_indices: set[int],
+    ) -> np.ndarray:
+        """Vectorized modal overlap matrix calculation for a set of modes.
+
+        This overlap is based on the bi-orthogonality relationship between E and H.
+        Returns a matrix S where S[i, j] is the overlap between mode i and mode j.
+        """
+        # Convert set to sorted list for consistent indexing
+        mode_list = sorted(mode_indices)
+
+        # Extract field components
+        Ex_sel = E[0][..., mode_list]  # (Nx, Ny, 1, n)
+        Ey_sel = E[1][..., mode_list]
+        Hx_sel = H[0][..., mode_list]
+        Hy_sel = H[1][..., mode_list]
+
+        # Make the differential area elements
+        Ex_Hy_dS = np.outer(dl_primal[0], dl_dual[1])
+        Ey_Hx_dS = np.outer(dl_dual[0], dl_primal[1])
+
+        # Add broadcasting dimensions: [..., i, j] for outer product
+        Ex_i = Ex_sel[..., :, np.newaxis]  # (Nx, Ny, 1, n, 1)
+        Ex_j = Ex_sel[..., np.newaxis, :]  # (Nx, Ny, 1, 1, n)
+        Hy_i = Hy_sel[..., :, np.newaxis]
+        Hy_j = Hy_sel[..., np.newaxis, :]
+
+        # Compute term1: (Ex[i] * Hy[j] + Ex[j] * Hy[i]) * dS
+        term1 = (Ex_i * Hy_j + Ex_j * Hy_i) * Ex_Hy_dS[..., np.newaxis, np.newaxis, np.newaxis]
+
+        # Compute term2: (Ey[i] * Hx[j] + Ey[j] * Hx[i]) * dS
+        Ey_i = Ey_sel[..., :, np.newaxis]
+        Ey_j = Ey_sel[..., np.newaxis, :]
+        Hx_i = Hx_sel[..., :, np.newaxis]
+        Hx_j = Hx_sel[..., np.newaxis, :]
+        term2 = (Ey_i * Hx_j + Ey_j * Hx_i) * Ey_Hx_dS[..., np.newaxis, np.newaxis, np.newaxis]
+
+        # Sum over spatial dimensions (x, y) to get overlap matrix
+        S = np.sum(term1 - term2, axis=(0, 1))[0]  # (1, n, n) -> (n, n)
+        return S
+
+    @staticmethod
+    def _normalize_modes(
+        E: np.ndarray,
+        H: np.ndarray,
+        dl_primal: np.ndarray,
+        dl_dual: np.ndarray,
+        mode_indices: list[int],
+    ) -> tuple[np.ndarray, np.ndarray]:
+        """Normalize modes so that their self-dot equals 1.0."""
+        for mode_idx in mode_indices:
+            self_dot = EigSolver._dot(E, H, dl_primal, dl_dual, mode_idx, mode_idx)
+            norm_factor = 1.0 / np.sqrt(self_dot)
+            E[..., mode_idx] *= norm_factor
+            H[..., mode_idx] *= norm_factor
+        return E, H
+
+    @staticmethod
+    def _make_orthogonal_basis_for_degenerate_modes(
+        degenerate_mode_sets: list[set[int]],
+        E_vec: np.ndarray,
+        H_vec: np.ndarray,
+        dl_primal: np.ndarray,
+        dl_dual: np.ndarray,
+    ) -> tuple[np.ndarray, np.ndarray]:
+        """Ensures that groups of degenerate modes are orthogonal, which is not guaranteed for eigenvectors with degenerate eigenvalues."""
+        dtype = E_vec.dtype
+        for degenerate_group in degenerate_mode_sets:
+            # Convert set to list for indexing
+            mode_indices = sorted(degenerate_group)
+
+            # S is a symmetric overlap matrix of the degenerate modes.
+            S = EigSolver._outer_dot(E_vec, H_vec, dl_primal, dl_dual, mode_indices)
+
+            # 1. Diagonalize
+            # eigenvalues (w) and right eigenvectors (v)
+            eigvals, eigvecs = np.linalg.eig(S)
+            # 2. Normalize Eigenvectors using UNCONJUGATED dot product
+            # Standard np.linalg.eig normalizes so v.conj().T @ v = 1 (Hermitian)
+            # We need v.T @ v = 1 (Complex Symmetric)
+            for i in range(eigvecs.shape[1]):
+                vec = eigvecs[:, i]
+                # Calculate complex self-dot (no conjugation)
+                self_dot = np.dot(vec, vec)
+
+                # Avoid division by zero
+                if np.isclose(self_dot, 0, atol=np.finfo(dtype).tiny):
+                    raise ValueError(f"Eigenvector {i} is self-orthogonal. Cannot orthogonalize.")
+
+                scale_factor = np.sqrt(self_dot)
+                eigvecs[:, i] = vec / scale_factor
+            # 3. Calculate Inverse Square Root of Eigenvalues
+            lambda_inv_sqrt = np.diag(1.0 / np.sqrt(eigvals))
+            # 4. Construct W
+            # Since V is complex orthogonal (V.T @ V = I), V_inv = V.T
+            # W = V * Lambda^(-1/2) * V.T
+            W = eigvecs @ lambda_inv_sqrt @ eigvecs.T
+            # S_new will be identity
+            # S_new = W.T @ S @ W
+            E_vec[..., mode_indices] = E_vec[..., mode_indices] @ W
+            H_vec[..., mode_indices] = H_vec[..., mode_indices] @ W
+
+        return E_vec, H_vec
+
 
 def compute_modes(*args: Any, **kwargs: Any) -> tuple[Numpy, Numpy, str]:
     """A wrapper around ``EigSolver.compute_modes``, which is used in :class:`.ModeSolver`."""