diff --git a/tests/test_data/test_monitor_data.py b/tests/test_data/test_monitor_data.py
index dedaa1d494..271d1b15ca 100644
--- a/tests/test_data/test_monitor_data.py
+++ b/tests/test_data/test_monitor_data.py
@@ -1243,3 +1243,182 @@ def test_symmetry_expansion_no_interpolation_warning():
 
     with AssertLogLevel(None):
         _ = field_data.symmetry_expanded_copy
+
+
+@pytest.mark.parametrize("colocate", [True, False])
+@pytest.mark.parametrize("sim_2d", [False, True])
+def test_diff_area_elements(colocate, sim_2d):
+    """Test differential area elements for different colocate and simulation dimension settings.
+
+    Tests that:
+    1. All 4 area elements are returned with correct shapes
+    2. For colocate=True, all 4 elements are identical
+    3. Total integrated area is consistent between methods
+    4. For 2D simulations (1 grid point in y), the unit handling (size=1.0) works correctly
+    5. Non-uniform mesh produces different primal/dual cell sizes (for non-colocated)
+    """
+    # Set up simulation and monitor sizes
+    if sim_2d:
+        # 2D simulation: y=0 means only 1 grid point in y direction
+        # Monitor is in x=0 plane with tangential dims y and z
+        # The y dimension will have only 1 grid point, triggering the 2D handling
+        sim_size = (4.0, 0, 3.0)
+        monitor_size = (0, td.inf, 3.0)  # inf in y to cover full sim
+    else:
+        # 3D simulation: both y and z dimensions are nonzero
+        sim_size = (4.0, 2.0, 3.0)
+        monitor_size = (0, 2.0 - np.pi / 10, 3.0 - np.pi / 10)
+    monitor_center = (0, 0, 0)
+
+    # Create a structure to induce non-uniform meshing
+    # A small box near the edge will cause mesh refinement in that region
+    structure = td.Structure(
+        geometry=td.Box(center=(0, 0.3, 0.5), size=(0.5, 0.3, 0.3)),
+        medium=td.Medium(permittivity=4.0),
+    )
+
+    # Use GridSpec.auto with override to create non-uniform mesh
+    # The structure will cause finer mesh near it, coarser elsewhere
+    grid_spec = td.GridSpec.auto(
+        wavelength=3.0,  # Large wavelength for coarse base mesh
+        min_steps_per_wvl=6,
+        override_structures=[
+            td.Structure(
+                geometry=td.Box(center=(0, 0.3, 0.5), size=(0.6, 0.4, 0.4)),
+                medium=td.Medium(permittivity=4.0),
+            )
+        ],
+    )
+
+    # Create simulation with non-uniform grid
+    sim = td.Simulation(
+        size=sim_size,
+        run_time=1e-12,
+        grid_spec=grid_spec,
+        structures=[structure],
+        sources=[
+            td.PointDipole(
+                source_time=td.GaussianPulse(freq0=1e14, fwidth=1e13),
+                polarization="Ez",
+                center=(0.1, 0, 0),
+            )
+        ],
+        boundary_spec=td.BoundarySpec.all_sides(boundary=td.Periodic()),
+    )
+
+    # Create field monitor
+    monitor = td.FieldMonitor(
+        size=monitor_size,
+        center=monitor_center,
+        freqs=[1e14],
+        name="field_monitor",
+        colocate=colocate,
+    )
+
+    # Get grid for monitor
+    grid = sim.discretize_monitor(monitor)
+
+    # Create field data arrays with appropriate coordinates
+    def make_field_array(field_name):
+        x, y, z = grid[field_name].to_list
+        shape = (len(x), len(y), len(z), 1)
+        data = np.random.rand(*shape) + 1j * np.random.rand(*shape)
+        return td.ScalarFieldDataArray(data, coords={"x": x, "y": y, "z": z, "f": [1e14]})
+
+    field_data = FieldData(
+        monitor=monitor,
+        Ex=make_field_array("Ex"),
+        Ey=make_field_array("Ey"),
+        Ez=make_field_array("Ez"),
+        Hx=make_field_array("Hx"),
+        Hy=make_field_array("Hy"),
+        Hz=make_field_array("Hz"),
+        symmetry=(0, 0, 0),
+        symmetry_center=(0, 0, 0),
+        grid_expanded=grid,
+    )
+
+    # Get differential areas
+    diff_area = field_data._diff_area
+    assert len(diff_area) == 4, "Should return 4 area elements"
+    dS_EuHv, dS_EvHu, dS_Ez, dS_Hz = diff_area
+
+    # Get tangential dimensions
+    tan_dims = field_data._tangential_dims
+    assert len(tan_dims) == 2
+
+    # Check that all areas have correct dimensions
+    for dS in diff_area:
+        assert set(dS.dims) == set(tan_dims), f"Area dimensions should be {tan_dims}"
+    # Additional check: all area values should be non-negative
+    for dS in diff_area:
+        assert np.all(dS.values >= 0), "All area elements should be non-negative"
+
+    # For colocated monitors, all 4 elements should be identical
+    if colocate:
+        np.testing.assert_allclose(dS_EuHv.values, dS_EvHu.values, rtol=1e-10)
+        np.testing.assert_allclose(dS_EuHv.values, dS_Ez.values, rtol=1e-10)
+        np.testing.assert_allclose(dS_EuHv.values, dS_Hz.values, rtol=1e-10)
+
+    # Test that _diff_area_at_boundaries returns all identical elements
+    dS_boundaries = field_data._diff_area_at_boundaries()
+
+    for i in range(1, 4):
+        np.testing.assert_allclose(dS_boundaries[0].values, dS_boundaries[i].values, rtol=1e-10)
+
+    # Test total integrated area consistency
+    # For a rectangular monitor, the total area should be approximately monitor_size[1] * monitor_size[2]
+    # (allowing for grid discretization effects)
+    # For 2D sims (y=0), the y dimension has size 1.0 for unit handling (W/um instead of W)
+    if sim_2d:
+        # 2D sim: y dim uses 1.0 multiplier for unit handling
+        expected_area = 1.0 * monitor_size[2]  # y=1.0, z=monitor_size[2]
+    else:
+        expected_area = monitor_size[1] * monitor_size[2]
+
+    # Total area from boundaries method
+    total_area_boundaries = float(dS_boundaries[0].sum())
+
+    # Total area from yee positions method
+    dS_yee = field_data._diff_area_at_yee_positions(truncate_to_monitor_bounds=True)
+    # For Yee grid, sum of cell×dual and dual×cell should both give consistent total areas
+    # when summed over the integration domain
+    total_area_EuHv = float(dS_yee[0].sum())
+    total_area_EvHu = float(dS_yee[1].sum())
+    total_area_Ez = float(dS_yee[2].sum())
+    total_area_Hz = float(dS_yee[3].sum())
+
+    # All total areas be the expected monitor area
+    rtol = 1e-12  # Allow 15% tolerance for discretization effects
+    np.testing.assert_allclose(total_area_boundaries, expected_area, rtol=rtol)
+    np.testing.assert_allclose(total_area_EuHv, expected_area, rtol=rtol)
+    np.testing.assert_allclose(total_area_EvHu, expected_area, rtol=rtol)
+    np.testing.assert_allclose(total_area_Ez, expected_area, rtol=rtol)
+    np.testing.assert_allclose(total_area_Hz, expected_area, rtol=rtol)
+
+    np.testing.assert_allclose(total_area_EuHv, total_area_boundaries, rtol=rtol)
+    np.testing.assert_allclose(total_area_EvHu, total_area_boundaries, rtol=rtol)
+    np.testing.assert_allclose(total_area_Ez, total_area_boundaries, rtol=rtol)
+    np.testing.assert_allclose(total_area_Hz, total_area_boundaries, rtol=rtol)
+
+    # For non-colocated monitors, the _diff_area should return the Yee grid areas
+    if not colocate:
+        np.testing.assert_allclose(dS_EuHv.values, dS_yee[0].values, rtol=rtol)
+        np.testing.assert_allclose(dS_EvHu.values, dS_yee[1].values, rtol=rtol)
+        np.testing.assert_allclose(dS_Ez.values, dS_yee[2].values, rtol=rtol)
+        np.testing.assert_allclose(dS_Hz.values, dS_yee[3].values, rtol=rtol)
+
+        # For non-colocated 3D monitors with non-uniform mesh, the 4 differential areas
+        # should be different from each other (cell×dual ≠ dual×cell ≠ cell×cell ≠ dual×dual)
+        if not sim_2d:
+            # Verify the mesh is actually non-uniform by checking areas differ
+            assert not np.allclose(dS_EuHv.values, dS_EvHu.values), (
+                "Non-uniform mesh should produce different areas for EuHv vs EvHu"
+            )
+            assert not np.allclose(dS_Ez.values, dS_Hz.values), (
+                "Non-uniform mesh should produce different areas for Ez vs Hz"
+            )
+
+    # For colocated monitors, _diff_area should return boundary areas
+    if colocate:
+        np.testing.assert_allclose(dS_EuHv.values, dS_boundaries[0].values, rtol=rtol)
diff --git a/tests/test_plugins/test_mode_solver.py b/tests/test_plugins/test_mode_solver.py
index cdc412513d..6d52f05667 100644
--- a/tests/test_plugins/test_mode_solver.py
+++ b/tests/test_plugins/test_mode_solver.py
@@ -196,8 +196,16 @@ def compare_colocation(ms):
     data_at_boundaries = ms_nocol.sim_data.at_boundaries(MODE_MONITOR_NAME)
 
     for key, field in data_col.field_components.items():
-        # Check the colocated data is the same
-        assert np.allclose(data_at_boundaries[key], field, atol=1e-7)
+        # Normalize both fields per-mode and per-frequency to the same peak value
+        # (colocate=True and colocate=False may have different normalizations)
+        field_at_boundaries = data_at_boundaries[key]
+        # Get dims to reduce over (all except mode_index and f)
+        reduce_dims = [d for d in field.dims if d not in ("mode_index", "f")]
+        max_field = np.abs(field).max(dim=reduce_dims)
+        max_at_boundaries = np.abs(field_at_boundaries).max(dim=reduce_dims)
+        field_normalized = field / max_field
+        field_at_boundaries_normalized = field_at_boundaries / max_at_boundaries
+        assert np.allclose(field_at_boundaries_normalized, field_normalized, atol=1e-7)
 
         # Also check coordinates
         for dim, coords1 in field.coords.items():
diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py
index b2bda4e653..92c0fb0939 100644
--- a/tidy3d/components/data/monitor_data.py
+++ b/tidy3d/components/data/monitor_data.py
@@ -16,6 +16,12 @@
 
 from tidy3d.components.base import cached_property, skip_if_fields_missing
 from tidy3d.components.base_sim.data.monitor_data import AbstractMonitorData
+from tidy3d.components.data.utils import (
+    _complex_power_flow_numpy,
+    _dot_numpy,
+    _instantaneous_power_flow_numpy,
+    _outer_dot_numpy,
+)
 from tidy3d.components.grid.grid import Coords, Grid
 from tidy3d.components.medium import Medium, MediumType
 from tidy3d.components.mode_spec import ModeSortSpec, ModeSpec
@@ -439,6 +445,15 @@ def _grid_correction_dict(self):
             "grid_dual_correction": self.grid_dual_correction,
         }
 
+    @property
+    def _normal_dim(self) -> str:
+        """For a 2D monitor data, return the name of the normal dimension. Raise if cannot
+        confirm that the associated monitor is 2D."""
+        if len(self.monitor.zero_dims) != 1:
+            raise DataError("Data must be 2D to get normal dimension.")
+        normal_dim = "xyz"[self.monitor.size.index(0)]
+        return normal_dim
+
     @property
     def _tangential_dims(self) -> list[str]:
         """For a 2D monitor data, return the names of the tangential dimensions. Raise if cannot
@@ -450,6 +465,145 @@ def _tangential_dims(self) -> list[str]:
 
         return tangential_dims
 
+    def _find_enclosing_boundary(
+        self, field_coord: float, boundaries: np.ndarray, side: str
+    ) -> float:
+        """Find the grid boundary just outside a field coordinate.
+
+        Parameters
+        ----------
+        field_coord : float
+            The field coordinate to find enclosing boundary for.
+        boundaries : np.ndarray
+            Array of grid boundaries from grid_expanded.
+        side : str
+            "lower" to find boundary below/at, "upper" to find boundary above/at.
+
+        Returns
+        -------
+        float
+            The grid boundary just outside the field coordinate.
+        """
+        if side == "lower":
+            idx = np.searchsorted(boundaries, field_coord, side="right") - 1
+            return boundaries[max(0, idx)]
+        else:  # "upper"
+            idx = np.searchsorted(boundaries, field_coord, side="left")
+            return boundaries[min(len(boundaries) - 1, idx)]
+
+    def _diff_area_at_yee_positions(
+        self, truncate_to_monitor_bounds: bool = False
+    ) -> tuple[DataArray, DataArray, DataArray, DataArray]:
+        """Return differential area elements at Yee grid stagger positions.
+
+        The four returned DataArrays correspond to differential areas at the four
+        staggered Yee grid locations for field components:
+        - dS_EuHv: area elements at Eu/Hv positions = cell_dim1 x dual_dim2
+        - dS_EvHu: area elements at Ev/Hu positions = dual_dim1 x cell_dim2
+        - dS_Ez: area elements at Ez positions = cell_dim1 x cell_dim2
+        - dS_Hz: area elements at Hz positions = dual_dim1 x dual_dim2
+
+        Parameters
+        ----------
+        truncate_to_monitor_bounds : bool = False
+            If True, clamp integration region to monitor bounds (for flux calculations).
+            If False, use grid_expanded bounds enclosing field data (for dot/outer_dot).
+            When monitor bounds are np.inf, always uses grid_expanded fallback.
+
+        Returns
+        -------
+        tuple[DataArray, DataArray, DataArray, DataArray]
+            A tuple of (dS_EuHv, dS_EvHu, dS_Ez, dS_Hz) DataArrays with differential
+            area elements for Yee grid integration. Each DataArray has dimensions
+            corresponding to the two tangential dimensions of the monitor plane.
+        """
+        if not self.grid_expanded:
+            raise DataError(
+                "Monitor data requires 'grid_expanded' to compute Yee grid integration sizes."
+            )
+
+        _, plane_inds = self.monitor.pop_axis([0, 1, 2], self.monitor.size.index(0.0))
+        dims = ["x", "y", "z"]
+        mnt_bounds = np.array(self.monitor.bounds)
+
+        cell_sizes = {}
+        dual_sizes = {}
+        # Get boundaries for determining integration limits
+        # Use full grid boundaries (not colocation_boundaries which truncates for non-colocated)
+        # to properly find enclosing boundaries when monitor bounds are infinite
+        full_bounds = self.grid_expanded.boundaries.to_dict
+
+        for axis in plane_inds:
+            dim = dims[axis]
+            # Use full grid boundaries (without colocation truncation) for finding integration limits
+            grid_boundaries = full_bounds[dim][:-1]  # Drop last boundary like colocation_boundaries
+
+            # Do not apply the spurious dl along a dimension where the simulation is 2D.
+            # Instead, we just set the boundaries such that the cell size along the zero dimension is 1,
+            # such that quantities like flux will come out in units of W / um.
+            if grid_boundaries.size == 1:
+                dual_sizes[dim] = np.array([1.0])
+                cell_sizes[dim] = np.array([1.0])
+                continue
+
+            mnt_min = mnt_bounds[0, axis]
+            mnt_max = mnt_bounds[1, axis]
+
+            # Get field coordinates (E at centers, H at boundaries along tangential axis)
+            e_field_centers = self._tangential_fields["E" + dim].coords[dim].to_numpy()
+            h_field_boundaries = self._tangential_fields["H" + dim].coords[dim].to_numpy()
+
+            # Determine integration bounds
+            if truncate_to_monitor_bounds:
+                # Use monitor bounds, handling inf by finding grid boundary outside field data
+                if np.isinf(mnt_min):
+                    integration_min = self._find_enclosing_boundary(
+                        e_field_centers[0], grid_boundaries, "lower"
+                    )
+                else:
+                    integration_min = mnt_min
+
+                if np.isinf(mnt_max):
+                    integration_max = self._find_enclosing_boundary(
+                        e_field_centers[-1], grid_boundaries, "upper"
+                    )
+                else:
+                    integration_max = mnt_max
+            else:
+                # Use grid_expanded bounds that enclose all field data
+                integration_min = self._find_enclosing_boundary(
+                    e_field_centers[0], grid_boundaries, "lower"
+                )
+                integration_max = self._find_enclosing_boundary(
+                    e_field_centers[-1], grid_boundaries, "upper"
+                )
+
+            # Build coordinate arrays for size calculation
+            centers = np.concatenate([[integration_min], e_field_centers])
+            boundaries = np.concatenate([h_field_boundaries, [integration_max]])
+
+            # Dual sizes: distances between cell centers, clamped
+            dual_coords = np.clip(centers, integration_min, integration_max)
+            dual_sizes[dim] = dual_coords[1:] - dual_coords[:-1]
+
+            # Cell/primal sizes: distances between boundaries, clamped
+            cell_coords = np.clip(boundaries, integration_min, integration_max)
+            cell_sizes[dim] = cell_coords[1:] - cell_coords[:-1]
+
+        dim1 = self._tangential_dims[0]
+        dim2 = self._tangential_dims[1]
+        dS_EuHv = np.outer(cell_sizes[dim1], dual_sizes[dim2])
+        dS_EvHu = np.outer(dual_sizes[dim1], cell_sizes[dim2])
+        dS_Ez = np.outer(cell_sizes[dim1], cell_sizes[dim2])
+        dS_Hz = np.outer(dual_sizes[dim1], dual_sizes[dim2])
+
+        return (
+            DataArray(dS_EuHv, dims=self._tangential_dims),
+            DataArray(dS_EvHu, dims=self._tangential_dims),
+            DataArray(dS_Ez, dims=self._tangential_dims),
+            DataArray(dS_Hz, dims=self._tangential_dims),
+        )
+
     @property
     def colocation_boundaries(self) -> Coords:
         """Coordinates to be used for colocation of the data to grid boundaries."""
@@ -502,13 +656,18 @@ def _plane_grid_centers(self) -> tuple[Coords1D, Coords1D]:
         """For 2D monitor data, return the centers of the in-plane grid"""
         return [(bs[1:] + bs[:-1]) / 2 for bs in self._plane_grid_boundaries]
 
-    @property
-    def _diff_area(self) -> DataArray:
-        """For a 2D monitor data, return the area of each cell in the plane, for use in numerical
-        integrations. This assumes that data is colocated to grid boundaries, and uses the
-        difference in the surrounding grid centers to compute the area.
-        """
+    def _diff_area_at_boundaries(self) -> tuple[DataArray, DataArray, DataArray, DataArray]:
+        """Return differential area elements for colocated fields at grid boundaries.
+
+        This assumes that data is colocated to grid boundaries, and uses the difference in
+        the surrounding grid centers to compute the area. All four returned elements are
+        identical since fields are at the same positions when colocated.
 
+        Returns
+        -------
+        tuple[DataArray, DataArray, DataArray, DataArray]
+            A tuple of (dS_EuHv, dS_EvHu, dS_Ez, dS_Hz) DataArrays, all identical.
+        """
         # Monitor values are interpolated to bounds
         bounds = self._plane_grid_boundaries
         # Coords to compute cell sizes around the interpolation locations
@@ -519,14 +678,14 @@ def _diff_area(self) -> DataArray:
         coords[0] = np.array([bounds[0][0], *coords[0].tolist(), bounds[0][-1]])
         coords[1] = np.array([bounds[1][0], *coords[1].tolist(), bounds[1][-1]])
 
-        """Truncate coords to monitor boundaries. This implicitly makes extra pixels which may be
-        present have size 0 and so won't be included in the integration. For pixels intersected
-        by the monitor edge, the size is truncated to the part covered by the monitor. When using
-        the differential area sizes defined in this way together with integrand values
-        defined at cell boundaries, the integration is equivalent to trapezoidal rule with the first
-        and last values interpolated to the exact monitor start/end location, if the integrand
-        is zero outside of the monitor geometry. This should usually be the case for flux and dot
-        computations"""
+        # Truncate coords to monitor boundaries. This implicitly makes extra pixels which may be
+        # present have size 0 and so won't be included in the integration. For pixels intersected
+        # by the monitor edge, the size is truncated to the part covered by the monitor. When using
+        # the differential area sizes defined in this way together with integrand values
+        # defined at cell boundaries, the integration is equivalent to trapezoidal rule with the
+        # first and last values interpolated to the exact monitor start/end location, if the
+        # integrand is zero outside of the monitor geometry. This should usually be the case for
+        # flux and dot computations.
         mnt_bounds = np.array(self.monitor.bounds)
         mnt_bounds = mnt_bounds[:, plane_inds].T
         coords[0][np.argwhere(coords[0] < mnt_bounds[0, 0])] = mnt_bounds[0, 0]
@@ -540,7 +699,26 @@ def _diff_area(self) -> DataArray:
         sizes_dim0 = coords[0][1:] - coords[0][:-1] if bounds[0].size > 1 else [1.0]
         sizes_dim1 = coords[1][1:] - coords[1][:-1] if bounds[1].size > 1 else [1.0]
 
-        return DataArray(np.outer(sizes_dim0, sizes_dim1), dims=self._tangential_dims)
+        area = DataArray(np.outer(sizes_dim0, sizes_dim1), dims=self._tangential_dims)
+        return (area, area, area, area)
+
+    @cached_property
+    def _diff_area(self) -> tuple[DataArray, DataArray, DataArray, DataArray]:
+        """Return differential area elements for flux/integration calculations.
+
+        Returns
+        -------
+        tuple[DataArray, DataArray, DataArray, DataArray]
+            (dS_EuHv, dS_EvHu, dS_Ez, dS_Hz) where:
+            - dS_EuHv: areas for Eu × Hv products
+            - dS_EvHu: areas for Ev × Hu products
+            - dS_Ez: areas for Ez integration
+            - dS_Hz: areas for Hz integration
+            For colocated monitors, all four elements are identical.
+        """
+        if self.monitor.colocate:
+            return self._diff_area_at_boundaries()
+        return self._diff_area_at_yee_positions(truncate_to_monitor_bounds=True)
 
     def _tangential_corrected(self, fields: dict[str, DataArray]) -> dict[str, DataArray]:
         """For a 2D monitor data, extract the tangential components from fields and orient them
@@ -559,7 +737,7 @@ def _tangential_corrected(self, fields: dict[str, DataArray]) -> dict[str, DataA
         tan_dims = self._tangential_dims
         components = [fname + dim for fname in "EH" for dim in tan_dims]
 
-        normal_dim = "xyz"[self.monitor.size.index(0)]
+        normal_dim = self._normal_dim
 
         tan_fields = {}
         for component in components:
@@ -637,7 +815,7 @@ def grid_corrected_copy(self) -> ElectromagneticFieldData:
         if len(self.monitor.zero_dims) != 1:
             return field_data
 
-        normal_dim = "xyz"[self.monitor.zero_dims[0]]
+        normal_dim = self._normal_dim
         update = {"grid_primal_correction": 1.0, "grid_dual_correction": 1.0}
         for field_name, field in field_data.field_components.items():
             eig_val = self.symmetry_eigenvalues[field_name](normal_dim)
@@ -652,7 +830,7 @@ def intensity(self) -> ScalarFieldDataArray:
         """Return the sum of the squared absolute electric field components."""
         self._check_fields_stored(["Ex", "Ey", "Ez"])
 
-        normal_dim = "xyz"[self.monitor.size.index(0)]
+        normal_dim = self._normal_dim
         fields = self._colocated_fields
         components = ("Ex", "Ey", "Ez")
         if any(cmp not in fields for cmp in components):
@@ -686,25 +864,69 @@ def poynting(self) -> ScalarFieldDataArray:
         projected to the direction normal to the monitor plane."""
         return self.complex_poynting.real
 
-    def package_flux_results(self, flux_values: DataArray) -> Any:
-        """How to package flux based on the coordinates present in the data."""
-        # Choose appropriate data array type based on coordinates
-        if "mode_index" in flux_values.dims:
-            return FreqModeDataArray(flux_values)
-        return FluxDataArray(flux_values)
+    def _prepare_data_arrays_for_flux(
+        self,
+        fields: dict[str, DataArray],
+    ) -> tuple[dict, dict[str, DataArray]]:
+        """Make DataArrays transposed for use in flux calculations.
+
+        Final arrays have spatial dimensions last: (..., Nu, Nv).
+        """
+        tangential_dims = self._tangential_dims
+        normal_dim = self._normal_dim
+        coords = list(fields.values())[0].coords
+
+        def prepare_one_array(data_array: DataArray) -> DataArray:
+            """Prepare array: remove normal dim and put spatial dims last."""
+            prepared = data_array
+            if normal_dim in data_array.dims:
+                prepared = prepared.squeeze(dim=normal_dim, drop=True)
+            return prepared.transpose(..., *tangential_dims)
+
+        prepped_fields = {key: prepare_one_array(field) for key, field in fields.items()}
+
+        # Build final_coords dict
+        final_coords = {"f": coords["f"].to_numpy()}
+        if "mode_index" in coords:
+            final_coords["mode_index"] = coords["mode_index"].to_numpy()
+
+        return final_coords, prepped_fields
 
     @cached_property
     def complex_flux(self) -> Union[FluxDataArray, FreqModeDataArray]:
-        """Flux for data corresponding to a 2D monitor."""
+        """Complex flux for data corresponding to a 2D monitor."""
+        # Get fields based on colocate setting
+        if self.monitor.colocate:
+            fields = self._colocated_tangential_fields
+        else:
+            fields = self._tangential_fields
 
-        # Compute flux by integrating Poynting vector in-plane
-        d_area = self._diff_area
-        poynting = self.complex_poynting
+        # Get differential area elements (handles both colocated and non-colocated)
+        dS_EuHv, dS_EvHu, _, _ = self._diff_area
+        dS_numpy = (dS_EuHv.to_numpy(), dS_EvHu.to_numpy())
 
-        flux_values = poynting * d_area
-        flux_values = flux_values.sum(dim=d_area.dims)
+        # Prepare fields for numpy calculation
+        final_coords, prepped_fields = self._prepare_data_arrays_for_flux(fields)
 
-        return self.package_flux_results(flux_values)
+        # Extract tangential components as tuples of numpy arrays
+        u, v = self._tangential_dims
+        E = (
+            prepped_fields["E" + u].to_numpy(),
+            prepped_fields["E" + v].to_numpy(),
+        )
+        H = (
+            prepped_fields["H" + u].to_numpy(),
+            prepped_fields["H" + v].to_numpy(),
+        )
+
+        # Compute complex flux using the numpy helper
+        flux_result = _complex_power_flow_numpy(E, H, dS_numpy)
+
+        # Return appropriate data array type
+        if "mode_index" in final_coords:
+            return FreqModeDataArray(flux_result, coords=final_coords)
+        else:
+            return FluxDataArray(flux_result, coords=final_coords)
 
     @cached_property
     def flux(self) -> Union[FluxDataArray, FreqModeDataArray]:
@@ -721,7 +943,7 @@ def mode_area(self) -> FreqModeDataArray:
         """
         intensity = self.intensity
         # integrate over the plane
-        d_area = self._diff_area
+        d_area = self._diff_area_at_boundaries()[0]
         num = (intensity * d_area).sum(dim=d_area.dims) ** 2
         den = (intensity**2 * d_area).sum(dim=d_area.dims)
 
@@ -731,9 +953,118 @@ def mode_area(self) -> FreqModeDataArray:
 
         return FreqModeDataArray(area)
 
+    def _prepare_data_arrays_for_dot(
+        self,
+        fields_self: dict[str, DataArray],
+        fields_other: dict[str, DataArray],
+    ) -> tuple[
+        dict, Union[FieldData, ModeData, ModeSolverData], Union[FieldData, ModeData, ModeSolverData]
+    ]:
+        """Make DataArrays transposed for use in "dot" overlap calculations.
+
+        Final arrays have shape (N_f, N_modes, Nu, Nv) where Nu and Nv are the spatial
+        dimensions in the transverse plane.
+
+        Broadcasting rules (matching develop's xarray behavior):
+        - Normal dimension: always dropped
+        - Frequency/mode_index: if other has length 1, squeeze and broadcast to self's coords;
+          otherwise use intersection
+        """
+        normal_dim = self._normal_dim
+        tangential_dims = self._tangential_dims
+
+        # Grab the first E field as a template
+        test_field = "E" + tangential_dims[0]
+        self_coords = fields_self[test_field].coords
+        other_coords = fields_other[test_field].coords
+
+        def get_broadcast_selection(
+            coord_name: str,
+        ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+            """Get final values and selection indices for a coordinate.
+
+            Broadcasting rules (matching develop's xarray behavior):
+            - If other has size 1: squeeze it out, use self's coords (broadcast)
+            - Otherwise: use intersection
+
+            Returns (final_values, self_sel, other_sel) or (None, None, None) if
+            coordinate is not present in either.
+            """
+            has_self = coord_name in self_coords
+            has_other = coord_name in other_coords
+
+            if not has_self and not has_other:
+                return None, None, None
+
+            if has_self and has_other:
+                idx_self = Index(self_coords[coord_name].values)
+                idx_other = Index(other_coords[coord_name].values)
+                size_self = len(idx_self)
+                size_other = len(idx_other)
+
+                if size_other == 1:
+                    # Other is squeezed out, broadcast to self's coords
+                    return idx_self.to_numpy(), np.arange(size_self), np.zeros(size_self, dtype=int)
+                else:
+                    # Use intersection (includes case where self has size 1)
+                    common = idx_self.intersection(idx_other, sort=False)
+                    return (
+                        common.to_numpy(),
+                        idx_self.get_indexer(common),
+                        idx_other.get_indexer(common),
+                    )
+            elif has_self:
+                vals = self_coords[coord_name].values
+                return vals, np.arange(len(vals)), None
+            else:
+                vals = other_coords[coord_name].values
+                return vals, None, np.arange(len(vals))
+
+        final_freqs, self_f_sel, other_f_sel = get_broadcast_selection("f")
+        final_modes, self_m_sel, other_m_sel = get_broadcast_selection("mode_index")
+
+        def prepare_one_array(
+            input_coords: xr.Coordinates,
+            freq_sel: np.ndarray,
+            mode_sel: np.ndarray,
+            data_array: DataArray,
+        ) -> DataArray:
+            """Prepare array for numpy-based dot calculation."""
+            # Select frequencies
+            prepared = data_array.isel(f=freq_sel)
+            # Remove the normal dimension if it exists as a dimension (not just scalar coord)
+            if normal_dim in prepared.dims:
+                prepared = prepared.squeeze(dim=normal_dim, drop=True)
+            # Handle mode_index
+            if "mode_index" not in input_coords:
+                # Add mode_index dimension if not present
+                prepared = prepared.expand_dims(dim={"mode_index": [0]}, axis=1)
+            elif mode_sel is not None:
+                # Select mode indices
+                prepared = prepared.isel(mode_index=mode_sel)
+            # Transpose to standard order: (f, mode_index, u, v)
+            prepared = prepared.transpose("f", "mode_index", *tangential_dims)
+            return prepared
+
+        prepped_fields_self = {
+            key: prepare_one_array(self_coords, self_f_sel, self_m_sel, field)
+            for key, field in fields_self.items()
+        }
+        prepped_fields_other = {
+            key: prepare_one_array(other_coords, other_f_sel, other_m_sel, field)
+            for key, field in fields_other.items()
+        }
+
+        # Build final_coords dict
+        final_coords = {"f": final_freqs}
+        if final_modes is not None:
+            final_coords["mode_index"] = final_modes
+
+        return final_coords, prepped_fields_self, prepped_fields_other
+
     def dot(
-        self, field_data: Union[FieldData, ModeData, ModeSolverData], conjugate: bool = True
-    ) -> ModeAmpsDataArray:
+        self, field_data: FieldData | ModeData | ModeSolverData, conjugate: bool = True
+    ) -> FreqDataArray | FreqModeDataArray:
         r"""Dot product (modal overlap) with another :class:`.FieldData` object. Both datasets have
         to be frequency-domain data associated with a 2D monitor. Along the tangential directions,
         the datasets have to have the same discretization. Along the normal direction, the monitor
@@ -747,15 +1078,23 @@ def dot(
 
         .. math:
 
-           \frac{1}{4} \int \left( E_0 \times H_1^* + H_0^* \times E_1 \) \, {\rm d}S
+           \frac{1}{4} \int \left( E_0^* \times H_1 + H_0^* \times E_1 \right) \, {\rm d}S
 
         Parameters
         ----------
-        field_data : :class:`ElectromagneticFieldData`
+        field_data : :class:`.FieldData` | :class:`.ModeData` | :class:`.ModeSolverData`
             A data instance to compute the dot product with.
         conjugate : bool, optional
             If ``True`` (default), the dot product is defined as above. If ``False``, the definition
-            is similar, but without the complex conjugation of the $H$ fields.
+            is similar, but without the complex conjugation of the fields.
+
+        Returns
+        -------
+        :class:`.FreqDataArray` | :class:`.FreqModeDataArray`
+            Data array with the complex-valued modal overlaps.
+
+            - If neither dataset has ``mode_index``: returns :class:`.FreqDataArray`.
+            - If either dataset has ``mode_index``: returns :class:`.FreqModeDataArray`.
 
         Note
         ----
@@ -767,48 +1106,50 @@ def dot(
             In the non-conjugated definition, modes are orthogonal, but the interpretation of the
             dot product power carried by a given mode is no longer valid.
         """
-
-        # Tangential fields for current and other field data
-        fields_self = self._colocated_tangential_fields
-
-        if conjugate:
-            fields_self = {key: field.conj() for key, field in fields_self.items()}
-
-        fields_other = field_data._interpolated_tangential_fields(self._plane_grid_boundaries)
         dim1, dim2 = self._tangential_dims
-        d_area = self._diff_area
 
-        # After interpolation, the tangential coordinates should match. However, the two arrays
-        # may either have the same shape along other dimensions, or be broadcastable.
-        if (
-            fields_self[next(iter(fields_self))].shape
-            == fields_other[next(iter(fields_other))].shape
-        ):
-            # Arrays are same shape, so we can use numpy
-            e_self_x_h_other = fields_self["E" + dim1].values * fields_other["H" + dim2].values
-            e_self_x_h_other -= fields_self["E" + dim2].values * fields_other["H" + dim1].values
-            h_self_x_e_other = fields_self["H" + dim1].values * fields_other["E" + dim2].values
-            h_self_x_e_other -= fields_self["H" + dim2].values * fields_other["E" + dim1].values
-            integrand = xr.DataArray(
-                e_self_x_h_other - h_self_x_e_other, coords=fields_self["E" + dim1].coords
-            )
-            integrand *= d_area
+        # Get the two sets of fields on the same grid, either colocated to boundaries or on the original
+        # uncolocated grid positions
+        use_interpolation = self.monitor.colocate or field_data.monitor.colocate
+        if use_interpolation:
+            fields_self = self._colocated_tangential_fields
+            fields_other = field_data._interpolated_tangential_fields(self._plane_grid_boundaries)
+            dS_EuHv, dS_EvHu, _, _ = self._diff_area_at_boundaries()
+            dS_numpy = (dS_EuHv.to_numpy(), dS_EvHu.to_numpy())
         else:
-            # Broadcasting is needed, which may be complicated depending on the dimensions order.
-            # Use xarray to handle robustly.
+            fields_self = self._tangential_fields
+            fields_other = field_data._tangential_fields
+            # Sanity check
+            assert fields_self["E" + dim1].coords[dim1] == fields_other["E" + dim1].coords[dim1]
+            assert fields_self["E" + dim2].coords[dim2] == fields_other["E" + dim2].coords[dim2]
+            assert fields_self["E" + dim1].coords[dim1] == fields_self["H" + dim2].coords[dim1]
+            assert fields_self["E" + dim2].coords[dim2] == fields_self["H" + dim1].coords[dim2]
+            # Get boundary-adjusted differential areas for proper integration weights
+            dS_EuHv, dS_EvHu, _, _ = self._diff_area_at_yee_positions(
+                truncate_to_monitor_bounds=False
+            )
+            dS_numpy = (dS_EuHv.to_numpy(), dS_EvHu.to_numpy())
+
+        # Determine broadcast behavior and final dimensions
+        final_coords, prepped_fields_self, prepped_fields_other = self._prepare_data_arrays_for_dot(
+            fields_self, fields_other
+        )
 
-            # Drop size-1 dimensions in the other data
-            fields_other = {key: field.squeeze(drop=True) for key, field in fields_other.items()}
+        # Extract tangential components as tuples
+        u, v = self._tangential_dims
+        E1 = (prepped_fields_self["E" + u].to_numpy(), prepped_fields_self["E" + v].to_numpy())
+        H1 = (prepped_fields_self["H" + u].to_numpy(), prepped_fields_self["H" + v].to_numpy())
+        E2 = (prepped_fields_other["E" + u].to_numpy(), prepped_fields_other["E" + v].to_numpy())
+        H2 = (prepped_fields_other["H" + u].to_numpy(), prepped_fields_other["H" + v].to_numpy())
 
-            # Cross products of fields
-            e_self_x_h_other = fields_self["E" + dim1] * fields_other["H" + dim2]
-            e_self_x_h_other -= fields_self["E" + dim2] * fields_other["H" + dim1]
-            h_self_x_e_other = fields_self["H" + dim1] * fields_other["E" + dim2]
-            h_self_x_e_other -= fields_self["H" + dim2] * fields_other["E" + dim1]
-            integrand = (e_self_x_h_other - h_self_x_e_other) * d_area
+        dot_result = _dot_numpy(E1, H1, E2, H2, dS_numpy, conjugate=conjugate)
 
-        # Integrate over plane
-        return ModeAmpsDataArray(0.25 * integrand.sum(dim=d_area.dims))
+        # Squeeze out mode_index dimension (axis 1) if not in final_coords
+        if "mode_index" not in final_coords:
+            dot_result = dot_result.squeeze(axis=1)
+            return FreqDataArray(dot_result, coords=final_coords)
+        else:
+            return FreqModeDataArray(dot_result, coords=final_coords)
 
     def _tangential_fields_match_coords(self, coords: ArrayFloat2D) -> bool:
         """Check if the tangential fields already match given coords in the tangential plane."""
@@ -855,35 +1196,126 @@ def _interpolated_tangential_fields(self, coords: ArrayFloat2D) -> dict[str, Dat
 
         return fields
 
+    def _prepare_data_arrays_for_outer_dot(
+        self,
+        fields_self: dict[str, DataArray],
+        fields_other: dict[str, DataArray],
+    ) -> tuple[
+        dict, Union[FieldData, ModeData, ModeSolverData], Union[FieldData, ModeData, ModeSolverData]
+    ]:
+        """Make DataArray that are transposed for use in "outer dot" overlap calculations. Final
+        DataArray to all be shape (N_f, N_modes, Nu, Nv), where Nu and Nv are the spatial dimensions in the transverse plane.
+
+        If "f" or "mode_index" dimensions do not exist in the input fields, they are added as dummy dims so that we always
+        get the same number of dimensions.
+        """
+
+        # Drop normal direction completely
+        normal_dim = self._normal_dim
+        tangential_dims = self._tangential_dims
+
+        # Grab the first E field as a template for the remaining data arrays.
+        test_field = "E" + tangential_dims[0]
+        self_coords = fields_self[test_field].coords
+        other_coords = fields_other[test_field].coords
+
+        # Common frequencies to both data arrays
+        freq_self = Index(self_coords["f"].values)
+        freq_other = Index(other_coords["f"].values)
+        common_freqs = freq_self.intersection(freq_other, sort=False)
+
+        # Keep frequency order consistent with the current data while aligning the other dataset.
+        self_f_sel = freq_self.get_indexer(common_freqs)
+        other_f_sel = freq_other.get_indexer(common_freqs)
+
+        def prepare_one_array(
+            input_coords: xr.Coordinates,
+            normal_dim: str,
+            tangential_dims: tuple[str, str],
+            freq_sel: np.ndarray,
+            data_array: DataArray,
+        ) -> DataArray:
+            """Modifies coordinates and dimensions of arrays for use with numpy."""
+            # Get common frequencies
+            prepared_array = data_array.isel(f=freq_sel)
+            # Remove the normal dimension if it exists as a dimension (not just scalar coord)
+            if normal_dim in prepared_array.dims:
+                prepared_array = prepared_array.squeeze(dim=normal_dim, drop=True)
+            # Add mode index dimension if not present
+            if "mode_index" not in input_coords:
+                prepared_array = prepared_array.expand_dims(dim={"mode_index": [0]}, axis=1)
+            # Finally transpose input arrays so that they are all ordered to make use of numpy style broadcasting
+            prepared_array = prepared_array.transpose("f", "mode_index", *tangential_dims)
+            return prepared_array
+
+        # Reorder and add possible missing dim so data will be prepared for numpy broadcasting
+        prepped_fields_self = {
+            key: prepare_one_array(self_coords, normal_dim, tangential_dims, self_f_sel, field)
+            for key, field in fields_self.items()
+        }
+
+        prepped_fields_other = {
+            key: prepare_one_array(other_coords, normal_dim, tangential_dims, other_f_sel, field)
+            for key, field in fields_other.items()
+        }
+
+        # Prepare final coords
+        final_coords = {
+            "f": common_freqs.to_numpy(),
+        }
+
+        # Determine mode_index coordinate handling:
+        # - Neither has mode_index → just f → FreqDataArray
+        # - At least one has mode_index → f + mode_index_0 + mode_index_1 → MixedModeDataArray
+        #   (when only one has mode_index, the other gets a dummy [0] coordinate)
+        self_has_mode = "mode_index" in self_coords
+        other_has_mode = "mode_index" in other_coords
+
+        if self_has_mode and other_has_mode:
+            final_coords["mode_index_0"] = self_coords["mode_index"].to_numpy()
+            final_coords["mode_index_1"] = other_coords["mode_index"].to_numpy()
+        elif self_has_mode:
+            final_coords["mode_index_0"] = self_coords["mode_index"].to_numpy()
+        elif other_has_mode:
+            final_coords["mode_index_1"] = other_coords["mode_index"].to_numpy()
+
+        return final_coords, prepped_fields_self, prepped_fields_other
+
     def outer_dot(
-        self, field_data: Union[FieldData, ModeData], conjugate: bool = True
-    ) -> MixedModeDataArray:
+        self, field_data: FieldData | ModeData | ModeSolverData, conjugate: bool = True
+    ) -> FreqDataArray | MixedModeDataArray:
         r"""Dot product (modal overlap) with another :class:`.FieldData` object.
 
         The tangential fields from ``field_data`` are interpolated to this object's grid, so the
         data arrays don't need to have the same discretization.  The calculation is performed for
-        all common frequencies between data arrays.  In the output, ``mode_index_0`` and
-        ``mode_index_1`` are the mode indices from this object and ``field_data``, respectively, if
-        they are instances of ``ModeData``.
+        all common frequencies between data arrays.
 
         The dot product is defined as:
 
         .. math:
 
-           \frac{1}{4} \int \left( E_0 \times H_1^* + H_0^* \times E_1 \) \, {\rm d}S
+           \frac{1}{4} \int \left( E_0^* \times H_1 + H_0^* \times E_1 \right) \, {\rm d}S
 
         Parameters
         ----------
-        field_data : :class:`ElectromagneticFieldData`
+        field_data : :class:`.FieldData` | :class:`.ModeData` | :class:`.ModeSolverData`
             A data instance to compute the dot product with.
         conjugate : bool = True
             If ``True`` (default), the dot product is defined as above. If ``False``, the definition
-            is similar, but without the complex conjugation of the $H$ fields.
+            is similar, but without the complex conjugation of the fields.
 
         Returns
         -------
-        :class:`xarray.DataArray`
-            Data array with the complex-valued modal overlaps between the two mode data.
+        :class:`.FreqDataArray` | :class:`.MixedModeDataArray`
+            Data array with the complex-valued modal overlaps.
+
+            - If neither dataset has ``mode_index``: returns :class:`.FreqDataArray`.
+            - If only ``self`` has ``mode_index``: returns :class:`.MixedModeDataArray` with
+              ``mode_index_0`` coordinate.
+            - If only ``field_data`` has ``mode_index``: returns :class:`.MixedModeDataArray` with
+              ``mode_index_1`` coordinate.
+            - If both datasets have ``mode_index``: returns :class:`.MixedModeDataArray` with
+              ``mode_index_0`` and ``mode_index_1`` coordinates.
 
         See also
         --------
@@ -895,157 +1327,63 @@ def outer_dot(
         if not all(a == b for a, b in zip(tan_dims, field_data._tangential_dims)):
             raise DataError("Tangential dimensions must match between the two monitors.")
 
-        # Tangential fields for current
-        fields_self = self._colocated_tangential_fields
-        if conjugate:
-            fields_self = {component: field.conj() for component, field in fields_self.items()}
-
-        # Tangential fields for other data
-
-        fields_other = field_data._interpolated_tangential_fields(self._plane_grid_boundaries)
-
-        # Tangential field component names
-        dim1, dim2 = tan_dims
-        e_1 = "E" + dim1
-        e_2 = "E" + dim2
-        h_1 = "H" + dim1
-        h_2 = "H" + dim2
-
-        # Prepare array with proper dimensions for the dot product data
-        arrays = (fields_self[e_1], fields_other[e_1])
-        coords = (arrays[0].coords, arrays[1].coords)
+        dim1, dim2 = self._tangential_dims
 
-        # Common frequencies to both data arrays
-        freq_self = Index(coords[0]["f"].values)
-        freq_other = Index(coords[1]["f"].values)
-        common_freqs = freq_self.intersection(freq_other, sort=False)
-        f = common_freqs.to_numpy()
-        # Keep frequency order consistent with the current data while aligning the other dataset.
-        isel1 = freq_self.get_indexer(common_freqs)
-        isel2 = freq_other.get_indexer(common_freqs)
-
-        # Mode indices, if available
-        modes_in_self = "mode_index" in coords[0]
-        modes_in_other = "mode_index" in coords[1]
-
-        keys = (e_1, e_2, h_1, h_2)
-        for key in keys:
-            fields_self[key] = fields_self[key].isel(f=isel1)
-            if modes_in_self:
-                fields_self[key] = fields_self[key].rename(mode_index="mode_index_0")
-            else:
-                fields_self[key] = fields_self[key].expand_dims(
-                    dim={"mode_index_0": [0]}, axis=len(fields_self[key].shape)
-                )
-            fields_other[key] = fields_other[key].isel(f=isel2)
-            if modes_in_other:
-                fields_other[key] = fields_other[key].rename(mode_index="mode_index_1")
-            else:
-                fields_other[key] = fields_other[key].expand_dims(
-                    dim={"mode_index_1": [0]}, axis=len(fields_other[key].shape)
-                )
+        # Get the two sets of fields on the same grid, either colocated to boundaries or on the original
+        # uncolocated grid positions
+        use_interpolation = self.monitor.colocate or field_data.monitor.colocate
+        if use_interpolation:
+            fields_self = self._colocated_tangential_fields
+            fields_other = field_data._interpolated_tangential_fields(self._plane_grid_boundaries)
+            dS_EuHv, dS_EvHu, _, _ = self._diff_area_at_boundaries()
+            dS_numpy = (dS_EuHv.to_numpy(), dS_EvHu.to_numpy())
+        else:
+            fields_self = self._tangential_fields
+            fields_other = field_data._tangential_fields
+            # Sanity check
+            assert fields_self["E" + dim1].coords[dim1] == fields_other["E" + dim1].coords[dim1]
+            assert fields_self["E" + dim2].coords[dim2] == fields_other["E" + dim2].coords[dim2]
+            assert fields_self["E" + dim1].coords[dim1] == fields_self["H" + dim2].coords[dim1]
+            assert fields_self["E" + dim2].coords[dim2] == fields_self["H" + dim1].coords[dim2]
+            # Get boundary-adjusted differential areas for proper integration weights
+            dS_EuHv, dS_EvHu, _, _ = self._diff_area_at_yee_positions(
+                truncate_to_monitor_bounds=False
+            )
+            dS_numpy = (dS_EuHv.to_numpy(), dS_EvHu.to_numpy())
 
-        d_area = self._diff_area.expand_dims(dim={"f": f}, axis=2).to_numpy()
-
-        # function to apply at each pair of mode indices before integrating
-        def fn(fields_1, fields_2):
-            e_self_1 = fields_1[e_1]
-            e_self_2 = fields_1[e_2]
-            h_self_1 = fields_1[h_1]
-            h_self_2 = fields_1[h_2]
-            e_other_1 = fields_2[e_1]
-            e_other_2 = fields_2[e_2]
-            h_other_1 = fields_2[h_1]
-            h_other_2 = fields_2[h_2]
-
-            # Cross products of fields
-            e_self_x_h_other = e_self_1 * h_other_2 - e_self_2 * h_other_1
-            h_self_x_e_other = h_self_1 * e_other_2 - h_self_2 * e_other_1
-
-            summand = 0.25 * (e_self_x_h_other - h_self_x_e_other) * d_area
-            return summand
-
-        result = self._outer_fn_summation(
-            fields_1=fields_self,
-            fields_2=fields_other,
-            outer_dim_1="mode_index_0",
-            outer_dim_2="mode_index_1",
-            sum_dims=tan_dims,
-            fn=fn,
+        # Determine broadcast behavior and final dimensions
+        final_coords, prepped_fields_self, prepped_fields_other = (
+            self._prepare_data_arrays_for_outer_dot(fields_self, fields_other)
         )
 
-        # Remove mode index coordinate if the input did not have it
-        if not modes_in_self:
-            result = result.isel(mode_index_0=0, drop=True)
-        if not modes_in_other:
-            result = result.isel(mode_index_1=0, drop=True)
-
-        return result
-
-    @staticmethod
-    def _outer_fn_summation(
-        fields_1: dict[str, xr.DataArray],
-        fields_2: dict[str, xr.DataArray],
-        outer_dim_1: str,
-        outer_dim_2: str,
-        sum_dims: list[str],
-        fn: Callable,
-    ) -> DataArray:
-        """
-        Loop over ``outer_dim_1`` and ``outer_dim_2``, apply ``fn`` to ``fields_1`` and ``fields_2``, and sum over ``sum_dims``.
-        The resulting ``DataArray`` has has dimensions any dimensions in the fields which are not contained in sum_dims.
-        This can be more memory efficient than vectorizing over the ``outer_dims``, which can involve broadcasting and reshaping data.
-        It also converts to numpy arrays outside the loops to minimize xarray overhead.
-        """
-        # first, convert to numpy outside the loop to reduce xarray overhead
-        fields_1_numpy = {key: val.to_numpy() for key, val in fields_1.items()}
-        fields_2_numpy = {key: val.to_numpy() for key, val in fields_2.items()}
-
-        # get one of the data arrays to look at for indexing
-        # assuming all data arrays have the same structure
-        data_array_temp_1 = list(fields_1.values())[0]
-        data_array_temp_2 = list(fields_2.values())[0]
-        numpy_temp_1 = data_array_temp_1.to_numpy()
-        numpy_temp_2 = data_array_temp_2.to_numpy()
-
-        # find the numpy axes associated with the provided dimensions
-        outer_axis_1 = data_array_temp_1.get_axis_num(outer_dim_1)
-        outer_axis_2 = data_array_temp_2.get_axis_num(outer_dim_2)
-        sum_axes = [data_array_temp_1.get_axis_num(dim) for dim in sum_dims]
-
-        # coords and array for result of calculation
-        coords = {key: val.to_numpy() for key, val in data_array_temp_1.coords.items()}
-        for dim in sum_dims:
-            coords.pop(dim)
-        # last two inds are the outer_dims
-        coords.pop(outer_dim_1)
-        coords[outer_dim_1] = data_array_temp_1.coords[outer_dim_1].to_numpy()
-        coords[outer_dim_2] = data_array_temp_2.coords[outer_dim_2].to_numpy()
-        # drop scalar non-indexing dimensions
-        coords = {key: val for key, val in coords.items() if len(val.shape) != 0}
-        shape = [len(val) for val in coords.values()]
-        dtype = np.promote_types(numpy_temp_1.dtype, numpy_temp_2.dtype)
-        data = np.zeros(shape, dtype=dtype)
-
-        # indexing tuples
-        idx_1 = [slice(None)] * numpy_temp_1.ndim
-        idx_2 = [slice(None)] * numpy_temp_2.ndim
-        idx_data = [slice(None)] * data.ndim
-
-        # calculate the sums of products
-        for outer_1 in range(numpy_temp_1.shape[outer_axis_1]):
-            for outer_2 in range(numpy_temp_2.shape[outer_axis_2]):
-                idx_1[outer_axis_1] = outer_1
-                idx_2[outer_axis_2] = outer_2
-                idx_data[-2] = outer_1
-                idx_data[-1] = outer_2
-                fields_1_curr = {key: val[tuple(idx_1)] for key, val in fields_1_numpy.items()}
-                fields_2_curr = {key: val[tuple(idx_2)] for key, val in fields_2_numpy.items()}
-                summand_curr = fn(fields_1_curr, fields_2_curr)
-                data_curr = np.sum(summand_curr, axis=tuple(sum_axes))
-                data[tuple(idx_data)] = data_curr
-
-        return DataArray(data, coords=coords)
+        # Extract tangential components as tuples of numpy arrays
+        u, v = self._tangential_dims
+        E1 = (prepped_fields_self["E" + u].to_numpy(), prepped_fields_self["E" + v].to_numpy())
+        H1 = (prepped_fields_self["H" + u].to_numpy(), prepped_fields_self["H" + v].to_numpy())
+        E2 = (prepped_fields_other["E" + u].to_numpy(), prepped_fields_other["E" + v].to_numpy())
+        H2 = (prepped_fields_other["H" + u].to_numpy(), prepped_fields_other["H" + v].to_numpy())
+        numpy_result = _outer_dot_numpy(E1, H1, E2, H2, dS_numpy, conjugate=conjugate)
+
+        # Determine return type based on final_coords
+        # numpy_result shape is (n_freq, n_modes_0, n_modes_1)
+        has_mode_index_0 = "mode_index_0" in final_coords
+        has_mode_index_1 = "mode_index_1" in final_coords
+
+        if has_mode_index_0 and has_mode_index_1:
+            # Both have mode_index: return MixedModeDataArray (no squeeze)
+            return MixedModeDataArray(numpy_result, coords=final_coords)
+        elif has_mode_index_0:
+            # Only self has mode_index: squeeze axis 2 and return MixedModeDataArray
+            squeezed = numpy_result.squeeze(axis=2)
+            return MixedModeDataArray(squeezed, coords=final_coords)
+        elif has_mode_index_1:
+            # Only other has mode_index: squeeze axis 1 and return MixedModeDataArray
+            squeezed = numpy_result.squeeze(axis=1)
+            return MixedModeDataArray(squeezed, coords=final_coords)
+        else:
+            # Neither has mode_index: return FreqDataArray
+            squeezed = numpy_result.squeeze(axis=(1, 2))
+            return FreqDataArray(squeezed, coords=final_coords)
 
     @property
     def time_reversed_copy(self) -> FieldData:
@@ -1519,13 +1857,46 @@ def poynting(self) -> ScalarFieldTimeDataArray:
         e_x_h -= np.real(tan_fields["E" + dim2]) * np.real(tan_fields["H" + dim1])
         return e_x_h
 
+    @cached_property
+    def complex_flux(self) -> DataArray:
+        """Complex flux is not defined for time-domain data."""
+        raise DataError("Complex power flow is not defined for time-domain data.")
+
     @cached_property
     def flux(self) -> FluxTimeDataArray:
         """Flux for data corresponding to a 2D monitor."""
+        dim1, dim2 = self._tangential_dims
+        tangential_dims = self._tangential_dims
+        normal_dim = self._normal_dim
+
+        # Get fields based on colocate setting
+        if self.monitor.colocate:
+            fields = self._colocated_tangential_fields
+        else:
+            fields = self._tangential_fields
+
+        # Get differential area elements (handles both colocated and non-colocated)
+        dS_EuHv, dS_EvHu, _, _ = self._diff_area
+        dS_numpy = (dS_EuHv.to_numpy(), dS_EvHu.to_numpy())
 
-        # Compute flux by integrating Poynting vector in-plane
-        d_area = self._diff_area
-        return FluxTimeDataArray((self.poynting * d_area).sum(dim=d_area.dims))
+        # Prepare arrays: remove normal dim if present, put spatial dims last
+        def prepare(data_array: DataArray) -> DataArray:
+            prepared = data_array
+            if normal_dim in prepared.dims:
+                prepared = prepared.squeeze(dim=normal_dim, drop=True)
+            return prepared.transpose(..., *tangential_dims)
+
+        # Extract and prepare field components
+        Eu = np.real(prepare(fields["E" + dim1]).to_numpy())
+        Ev = np.real(prepare(fields["E" + dim2]).to_numpy())
+        Hu = np.real(prepare(fields["H" + dim1]).to_numpy())
+        Hv = np.real(prepare(fields["H" + dim2]).to_numpy())
+
+        flux_result = _instantaneous_power_flow_numpy((Eu, Ev), (Hu, Hv), dS_numpy)
+
+        return FluxTimeDataArray(
+            flux_result, coords={"t": fields["E" + dim1].coords["t"].to_numpy()}
+        )
 
     def dot(self, field_data: ElectromagneticFieldData, conjugate: bool = True) -> xr.DataArray:
         """Inner product is not defined for time-domain data."""
@@ -2027,7 +2398,7 @@ def _colocated_propagation_axes_field(self, field_name: Literal["E", "H"]) -> Da
         with normal axis along z to frame with propagation axis along z.
         """
         tan_dims = self._tangential_dims
-        normal_dim = "xyz"[self.monitor.zero_dims[0]]
+        normal_dim = self._normal_dim
         fields = self._colocated_fields
         fields = {key: val.squeeze(dim=normal_dim, drop=True) for key, val in fields.items()}
         mode_spec = self.monitor.mode_spec
@@ -2068,7 +2439,8 @@ def pol_fraction(self) -> xr.Dataset:
 
         tan_dims = self._tangential_dims
         e_field = self._colocated_propagation_axes_field("E")
-        diff_area = self._diff_area
+        # Use boundary areas since we're working with colocated fields
+        diff_area = self._diff_area_at_boundaries()[0]
         tm_int = (diff_area * np.abs(e_field.sel(component=1, drop=True)) ** 2).sum(dim=tan_dims)
         te_int = (diff_area * np.abs(e_field.sel(component=0, drop=True)) ** 2).sum(dim=tan_dims)
         te_frac = te_int / (te_int + tm_int)
@@ -2101,7 +2473,8 @@ def pol_fraction_waveguide(self) -> xr.Dataset:
         tan_dims = self._tangential_dims
         e_field = self._colocated_propagation_axes_field("E")
         h_field = self._colocated_propagation_axes_field("H")
-        diff_area = self._diff_area
+        # Use boundary areas since we're working with colocated fields
+        diff_area = self._diff_area_at_boundaries()[0]
 
         # te fraction
         field_int = [np.abs(e_field.sel(component=ind, drop=True)) ** 2 for ind in range(3)]
@@ -2506,6 +2879,7 @@ def normalize(self, source_spectrum_fn: Callable[[float], complex]) -> ModeSolve
     def _normalize_modes(self):
         """Normalize modes. Note: this modifies ``self`` in-place."""
         scaling = np.sqrt(np.abs(self.flux))
+        # scaling = np.sqrt(self.dot(self, conjugate=self.monitor.conjugated_dot_product))
         for field in self.field_components.values():
             field /= scaling
 
@@ -2687,7 +3061,7 @@ def interp_in_freq(
                     self.monitor.mode_spec,
                     update_dict["n_complex"],
                     self.monitor.direction,
-                    "xyz"[self.monitor._normal_axis],
+                    self._normal_dim,
                 )
             )
 
diff --git a/tidy3d/components/data/utils.py b/tidy3d/components/data/utils.py
index e2c96c0031..854bef9fdc 100644
--- a/tidy3d/components/data/utils.py
+++ b/tidy3d/components/data/utils.py
@@ -20,6 +20,230 @@
 CustomSpatialDataTypeAnnotated = Union[SpatialDataArray, annotate_type(UnstructuredGridDatasetType)]
 
 
+def _instantaneous_power_flow_numpy(
+    E: tuple[np.ndarray, np.ndarray],
+    H: tuple[np.ndarray, np.ndarray],
+    dS: tuple[np.ndarray, np.ndarray],
+) -> np.ndarray:
+    """Compute instantaneous power flow (Poynting vector integral) for real fields.
+
+    Computes (E x H) integrated over the plane. For time-domain (real) fields.
+
+    Parameters
+    ----------
+    E : tuple[np.ndarray, np.ndarray]
+        Tangential E-field components (Eu, Ev), each shape ``(..., nu, nv)``.
+    H : tuple[np.ndarray, np.ndarray]
+        Tangential H-field components (Hu, Hv), each shape ``(..., nu, nv)``.
+    dS : tuple[np.ndarray, np.ndarray]
+        Area elements at the two Yee grid locations, each shape ``(nu, nv)``.
+        ``dS[0]`` is for Eu/Hv location, ``dS[1]`` for Ev/Hu location.
+
+    Returns
+    -------
+    np.ndarray
+        Flux values with shape ``(...)``, spatial dimensions integrated out.
+    """
+    dS_EuHv, dS_EvHu = dS
+    Eu, Ev = E
+    Hu, Hv = H
+
+    # Instantaneous Poynting vector: S = E x H (normal component)
+    # For tangential fields (Eu, Ev) and (Hu, Hv):
+    # Sn = Eu * Hv - Ev * Hu
+    # Each term is at a different Yee grid location
+    term_EuHv = (Eu * Hv) * dS_EuHv
+    term_EvHu = (Ev * Hu) * dS_EvHu
+
+    # Sum over spatial dimensions (last two)
+    return np.sum(term_EuHv - term_EvHu, axis=(-2, -1))
+
+
+def _complex_power_flow_numpy(
+    E: tuple[np.ndarray, np.ndarray],
+    H: tuple[np.ndarray, np.ndarray],
+    dS: tuple[np.ndarray, np.ndarray],
+) -> np.ndarray:
+    """Compute power flow (Poynting vector integral) with Yee grid area elements.
+
+    Computes 0.5 * (E x H*) integrated over the plane.
+
+    Parameters
+    ----------
+    E : tuple[np.ndarray, np.ndarray]
+        Tangential E-field components (Eu, Ev), each shape ``(..., nu, nv)``.
+    H : tuple[np.ndarray, np.ndarray]
+        Tangential H-field components (Hu, Hv), each shape ``(..., nu, nv)``.
+    dS : tuple[np.ndarray, np.ndarray]
+        Area elements at the two Yee grid locations, each shape ``(nu, nv)``.
+        ``dS[0]`` is for Eu/Hv location, ``dS[1]`` for Ev/Hu location.
+
+    Returns
+    -------
+    np.ndarray
+        Flux values with shape ``(...)``, spatial dimensions integrated out.
+    """
+    dS_EuHv, dS_EvHu = dS
+    Eu, Ev = E
+    Hu, Hv = H
+
+    # Poynting vector: S = 0.5 E x H* (normal component)
+    # For tangential fields (Eu, Ev) and (Hu, Hv):
+    # Sn = Eu * Hv* - Ev * Hu*
+    # Each term is at a different Yee grid location
+    term_EuHv = (Eu * np.conj(Hv)) * dS_EuHv
+    term_EvHu = (Ev * np.conj(Hu)) * dS_EvHu
+
+    # Sum over spatial dimensions (last two)
+    return 0.5 * np.sum(term_EuHv - term_EvHu, axis=(-2, -1))
+
+
+def _dot_numpy(
+    E1: tuple[np.ndarray, np.ndarray],
+    H1: tuple[np.ndarray, np.ndarray],
+    E2: tuple[np.ndarray, np.ndarray],
+    H2: tuple[np.ndarray, np.ndarray],
+    dS: tuple[np.ndarray, np.ndarray],
+    conjugate: bool = False,
+    bidirectional: bool = True,
+) -> np.ndarray:
+    """Compute modal overlap integral.
+
+    By default computes the bidirectional overlap: 1/4 * integral(E1 x H2 + H1 x E2) dS.
+    With bidirectional=False, computes just: 1/2 * integral(E1 x H2) dS.
+
+    Parameters
+    ----------
+    E1 : tuple[np.ndarray, np.ndarray]
+        Tangential E-field components (Eu, Ev) from first dataset, each shape ``(..., nu, nv)``.
+    H1 : tuple[np.ndarray, np.ndarray]
+        Tangential H-field components (Hu, Hv) from first dataset, each shape ``(..., nu, nv)``.
+    E2 : tuple[np.ndarray, np.ndarray]
+        Tangential E-field components (Eu, Ev) from second dataset, each shape ``(..., nu, nv)``.
+    H2 : tuple[np.ndarray, np.ndarray]
+        Tangential H-field components (Hu, Hv) from second dataset, each shape ``(..., nu, nv)``.
+    dS : tuple[np.ndarray, np.ndarray]
+        Area elements at the two Yee grid locations, each shape ``(nu, nv)``.
+        ``dS[0]`` is for Eu/Hv location, ``dS[1]`` for Ev/Hu location.
+    conjugate : bool
+        If True, conjugate the first set of fields (E1, H1) before computing overlap.
+    bidirectional : bool
+        If True (default), computes symmetric overlap 1/4 * (E1 x H2 + H1 x E2).
+        If False, computes just 1/2 * (E1 x H2).
+
+    Returns
+    -------
+    np.ndarray
+        Overlap values with shape ``(...)``, spatial dimensions integrated out.
+    """
+    dS_EuHv, dS_EvHu = dS
+
+    E1u, E1v = E1
+    H1u, H1v = H1
+    E2u, E2v = E2
+    H2u, H2v = H2
+
+    assert E1u.shape == H1v.shape
+    assert E1v.shape == H1u.shape
+    assert E2u.shape == H2v.shape
+    assert E2v.shape == H2u.shape
+
+    if conjugate:
+        E1u, E1v = np.conj(E1u), np.conj(E1v)
+        H1u, H1v = np.conj(H1u), np.conj(H1v)
+
+    # Group terms by grid location for proper area weighting
+    if bidirectional:
+        # Bidirectional: 0.25 * integral(E1 x H2 + H1 x E2) dS
+        # At Eu/Hv location: E1u * H2v + H1v * E2u
+        # At Ev/Hu location: E1v * H2u + H1u * E2v
+        term_EuHv = (E1u * H2v + H1v * E2u) * dS_EuHv
+        term_EvHu = (E1v * H2u + H1u * E2v) * dS_EvHu
+        # Sum over spatial dimensions (last two)
+        return 0.25 * np.sum(term_EuHv - term_EvHu, axis=(-2, -1))
+    else:
+        # Non-bidirectional: 0.5 * integral(E1 x H2) dS
+        # At Eu/Hv location: E1u * H2v
+        # At Ev/Hu location: E1v * H2u
+        term_EuHv = (E1u * H2v) * dS_EuHv
+        term_EvHu = (E1v * H2u) * dS_EvHu
+        # Sum over spatial dimensions (last two)
+        return 0.5 * np.sum(term_EuHv - term_EvHu, axis=(-2, -1))
+
+
+def _outer_dot_numpy(
+    E1: tuple[np.ndarray, np.ndarray],
+    H1: tuple[np.ndarray, np.ndarray],
+    E2: tuple[np.ndarray, np.ndarray],
+    H2: tuple[np.ndarray, np.ndarray],
+    dS: tuple[np.ndarray, np.ndarray],
+    conjugate: bool = False,
+    bidirectional: bool = True,
+) -> np.ndarray:
+    """Compute pairwise modal overlap matrix.
+
+    Computes all elements of the overlap matrix S[i,j] = <mode_i | mode_j>.
+    By default computes the bidirectional overlap: 1/4 * integral(E1 x H2 + H1 x E2) dS.
+    With bidirectional=False, computes just: 1/2 * integral(E1 x H2) dS.
+
+    Parameters
+    ----------
+    E1 : tuple[np.ndarray, np.ndarray]
+        Tangential E-field components (Eu, Ev) from first dataset, each shape ``(..., n_modes_1, nu, nv)``.
+        Mode index is at -3, spatial dims at -2, -1.
+    H1 : tuple[np.ndarray, np.ndarray]
+        Tangential H-field components (Hu, Hv) from first dataset, each shape ``(..., n_modes_1, nu, nv)``.
+    E2 : tuple[np.ndarray, np.ndarray]
+        Tangential E-field components (Eu, Ev) from second dataset, each shape ``(..., n_modes_2, nu, nv)``.
+    H2 : tuple[np.ndarray, np.ndarray]
+        Tangential H-field components (Hu, Hv) from second dataset, each shape ``(..., n_modes_2, nu, nv)``.
+    dS : tuple[np.ndarray, np.ndarray]
+        Area elements at the two Yee grid locations, each shape ``(nu, nv)``.
+        ``dS[0]`` is for Eu/Hv location, ``dS[1]`` for Ev/Hu location.
+    conjugate : bool
+        If True, conjugate the first set of fields (E1, H1) before computing overlap.
+    bidirectional : bool
+        If True (default), computes symmetric overlap 1/4 * (E1 x H2 + H1 x E2).
+        If False, computes just 1/2 * (E1 x H2).
+
+    Returns
+    -------
+    np.ndarray
+        Overlap matrix with shape ``(..., n_modes_1, n_modes_2)``.
+    """
+    E1u, E1v = E1
+    H1u, H1v = H1
+    E2u, E2v = E2
+    H2u, H2v = H2
+
+    assert E1u.shape == H1v.shape
+    assert E1v.shape == H1u.shape
+    assert E2u.shape == H2v.shape
+    assert E2v.shape == H2u.shape
+
+    # Get number of modes and broadcast shape
+    n_modes_1 = E1u.shape[-3]
+    n_modes_2 = E2u.shape[-3]
+    broadcast_shape = E1u.shape[:-3]
+    dtype = E1u.dtype
+
+    # Initialize output matrix
+    S = np.zeros((*broadcast_shape, n_modes_1, n_modes_2), dtype=dtype)
+
+    # Compute all elements of the overlap matrix
+    for i in range(n_modes_1):
+        j_indices = np.arange(n_modes_2)
+
+        E1_i = (E1u[..., i : i + 1, :, :], E1v[..., i : i + 1, :, :])
+        H1_i = (H1u[..., i : i + 1, :, :], H1v[..., i : i + 1, :, :])
+        E2_j = (E2u[..., j_indices, :, :], E2v[..., j_indices, :, :])
+        H2_j = (H2u[..., j_indices, :, :], H2v[..., j_indices, :, :])
+
+        S[..., i, j_indices] = _dot_numpy(E1_i, H1_i, E2_j, H2_j, dS, conjugate, bidirectional)
+
+    return S
+
+
 def _get_numpy_array(data_array: Union[ArrayLike, DataArray, UnstructuredGridDataset]) -> ArrayLike:
     """Get numpy representation of dataarray/dataset values."""
     if isinstance(data_array, UnstructuredGridDataset):
diff --git a/tidy3d/components/eme/monitor.py b/tidy3d/components/eme/monitor.py
index 4b57259c06..bfdbb0fa3d 100644
--- a/tidy3d/components/eme/monitor.py
+++ b/tidy3d/components/eme/monitor.py
@@ -65,8 +65,8 @@ class EMEMonitor(AbstractMonitor, ABC):
         "Not all monitors support values different from 1.",
     )
 
-    colocate: Literal[True] = pd.Field(
-        True,
+    colocate: Literal[False] = pd.Field(
+        False,
         title="Colocate Fields",
         description="Defines whether fields are colocated to grid cell boundaries (i.e. to the "
         "primal grid) on-the-fly during a solver run. Can be toggled for field recording monitors "
@@ -148,7 +148,7 @@ class EMEModeSolverMonitor(EMEMonitor):
     )
 
     colocate: bool = pd.Field(
-        True,
+        False,
         title="Colocate Fields",
         description="Toggle whether fields should be colocated to grid cell boundaries (i.e. "
         "primal grid nodes). Default (False) is used internally in EME propagation.",
@@ -222,7 +222,7 @@ class EMEFieldMonitor(EMEMonitor, AbstractFieldMonitor):
     )
 
     colocate: bool = pd.Field(
-        True,
+        False,
         title="Colocate Fields",
         description="Toggle whether fields should be colocated to grid cell boundaries (i.e. "
         "primal grid nodes). Default (False) is used internally in EME propagation.",
diff --git a/tidy3d/components/microwave/data/monitor_data.py b/tidy3d/components/microwave/data/monitor_data.py
index cfccc9bcb9..f0b7c24d9f 100644
--- a/tidy3d/components/microwave/data/monitor_data.py
+++ b/tidy3d/components/microwave/data/monitor_data.py
@@ -374,15 +374,36 @@ def wave_impedance(self) -> ImpedanceFreqModeDataArray:
         """
         self._check_fields_stored(["Ex", "Ey", "Ez", "Hx", "Hy", "Hz"])
 
-        tan_fields = self._colocated_tangential_fields
         dim1, dim2 = self._tangential_dims
-        e1 = tan_fields["E" + dim1]
-        e2 = tan_fields["E" + dim2]
-        diff_area = self._diff_area
-        field_int = [np.abs(e_field) ** 2 for e_field in [e1, e2]]
-        tangential_intensity = (diff_area * (field_int[0] + field_int[1])).sum(
-            dim=self._tangential_dims
+
+        # Get fields based on colocate setting
+        if self.monitor.colocate:
+            fields = self._colocated_tangential_fields
+        else:
+            fields = self._tangential_fields
+
+        # Get differential area elements (handles both colocated and non-colocated)
+        dS_E1H2, dS_E2H1, _, _ = self._diff_area
+        # Add coordinates from field arrays for proper xarray broadcasting
+        e1_coords = fields["E" + dim1].coords
+        e2_coords = fields["E" + dim2].coords
+        dS_E1 = xr.DataArray(
+            dS_E1H2.values,
+            coords={dim1: e1_coords[dim1], dim2: e1_coords[dim2]},
+        )
+        dS_E2 = xr.DataArray(
+            dS_E2H1.values,
+            coords={dim1: e2_coords[dim1], dim2: e2_coords[dim2]},
         )
+
+        e1 = fields["E" + dim1]
+        e2 = fields["E" + dim2]
+
+        # Integrate |E|² with proper area elements at each Yee location
+        intensity_E1 = (np.abs(e1) ** 2 * dS_E1).sum(dim=self._tangential_dims)
+        intensity_E2 = (np.abs(e2) ** 2 * dS_E2).sum(dim=self._tangential_dims)
+        tangential_intensity = intensity_E1 + intensity_E2
+
         direction = self.monitor.store_fields_direction
         P = self.complex_flux if direction == "+" else -self.complex_flux
         Z_wave = tangential_intensity / P / 2
diff --git a/tidy3d/components/microwave/impedance_calculator.py b/tidy3d/components/microwave/impedance_calculator.py
index 8f80c57320..d7ac0fcf4d 100644
--- a/tidy3d/components/microwave/impedance_calculator.py
+++ b/tidy3d/components/microwave/impedance_calculator.py
@@ -131,16 +131,18 @@ def compute_impedance(
             flux_sign = 1 if em_field.monitor.store_fields_direction == "+" else -1
 
         if self.voltage_integral is None:
-            flux = flux_sign * em_field.complex_flux
             if isinstance(em_field, FieldTimeData):
+                flux = flux_sign * em_field.flux
                 impedance = flux / np.real(current) ** 2
             else:
+                flux = flux_sign * em_field.complex_flux
                 impedance = 2 * flux / (current * np.conj(current))
         elif self.current_integral is None:
-            flux = flux_sign * em_field.complex_flux
             if isinstance(em_field, FieldTimeData):
+                flux = flux_sign * em_field.flux
                 impedance = np.real(voltage) ** 2 / flux
             else:
+                flux = flux_sign * em_field.complex_flux
                 impedance = (voltage * np.conj(voltage)) / (2 * np.conj(flux))
         else:
             if isinstance(em_field, FieldTimeData):
diff --git a/tidy3d/components/mode/solver.py b/tidy3d/components/mode/solver.py
index e533c3d084..bd8fbe516b 100644
--- a/tidy3d/components/mode/solver.py
+++ b/tidy3d/components/mode/solver.py
@@ -7,6 +7,7 @@
 import numpy as np
 
 from tidy3d.components.base import Tidy3dBaseModel
+from tidy3d.components.data.utils import _dot_numpy, _outer_dot_numpy
 from tidy3d.components.types import EpsSpecType, ModeSolverType, Numpy
 from tidy3d.constants import C_0, ETA_0, fp_eps, pec_val
 
@@ -1214,27 +1215,22 @@ def _identify_degenerate_modes_with_dot(
     def _dot(
         E: np.ndarray,
         H: np.ndarray,
-        dl_primal: np.ndarray,
-        dl_dual: np.ndarray,
+        dl_primal: list[np.ndarray],
+        dl_dual: list[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])
+        dS = (np.outer(dl_primal[0], dl_dual[1]), 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
-        term2 = Ey[..., mode_1] * Hx[..., mode_2] + Ey[..., mode_2] * Hx[..., mode_1]
-        term2 *= Ey_Hx_dS
-        return (1 / 4) * np.sum(term1 - term2)
+        # Extract tangential field components for the two modes as tuples
+        E1 = (E[0, ..., mode_1], E[1, ..., mode_1])
+        H1 = (H[0, ..., mode_1], H[1, ..., mode_1])
+        E2 = (E[0, ..., mode_2], E[1, ..., mode_2])
+        H2 = (H[0, ..., mode_2], H[1, ..., mode_2])
+
+        return _dot_numpy(E1, H1, E2, H2, dS, conjugate=False)
 
     @staticmethod
     def _cauchy_schwarz_dot_bound(
@@ -1286,8 +1282,8 @@ def _cauchy_schwarz_dot_bound(
     def _outer_dot(
         E: np.ndarray,
         H: np.ndarray,
-        dl_primal: np.ndarray,
-        dl_dual: np.ndarray,
+        dl_primal: list[np.ndarray],
+        dl_dual: list[np.ndarray],
         mode_indices: set[int],
     ) -> np.ndarray:
         """Vectorized modal overlap matrix calculation for a set of modes.
@@ -1300,52 +1296,16 @@ def _outer_dot(
         """
         # Convert set to sorted list for consistent indexing
         mode_list = sorted(mode_indices)
-        n_modes = len(mode_list)
 
-        # 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]
+        # Extract tangential field components as tuples
+        # Fancy indexing with mode_list moves mode dimension to front: (n_modes, Nx, Ny)
+        E_tan = (E[0, ..., mode_list], E[1, ..., mode_list])
+        H_tan = (H[0, ..., mode_list], 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])
-
-        # Initialize output matrix
-        dtype = Ex_sel.dtype
-        S = np.zeros((n_modes, n_modes), dtype=dtype)
-
-        # Only compute upper triangle (including diagonal)
-        for i in range(n_modes):
-            # Vectorize over j >= i
-            j_indices = np.arange(i, n_modes)
-
-            # Extract mode i fields
-            Ex_i = Ex_sel[..., i : i + 1]  # (Nx, Ny, 1, 1)
-            Ey_i = Ey_sel[..., i : i + 1]
-            Hy_i = Hy_sel[..., i : i + 1]
-            Hx_i = Hx_sel[..., i : i + 1]
-
-            # Extract mode j fields (vectorized for j >= i)
-            Ex_j = Ex_sel[..., j_indices]  # (Nx, Ny, 1, n_j)
-            Ey_j = Ey_sel[..., j_indices]
-            Hy_j = Hy_sel[..., j_indices]
-            Hx_j = Hx_sel[..., j_indices]
-
-            # 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]
-
-            # Compute term2: (Ey[i] * Hx[j] + Ey[j] * Hx[i]) * dS
-            term2 = (Ey_i * Hx_j + Ey_j * Hx_i) * Ey_Hx_dS[..., np.newaxis]
-
-            # Sum over spatial dimensions to get S[i, j] for j >= i
-            S[i, j_indices] = (1 / 4) * np.sum(term1 - term2, axis=(0, 1))
-
-        # Fill lower triangle by symmetry
-        S = S + S.T - np.diag(np.diag(S))
+        dS = (np.outer(dl_primal[0], dl_dual[1]), np.outer(dl_dual[0], dl_primal[1]))
 
-        return S
+        return _outer_dot_numpy(E_tan, H_tan, E_tan, H_tan, dS, conjugate=False)
 
     @staticmethod
     def _normalize_modes(
diff --git a/tidy3d/components/monitor.py b/tidy3d/components/monitor.py
index df4ec33049..b3290e0d3b 100644
--- a/tidy3d/components/monitor.py
+++ b/tidy3d/components/monitor.py
@@ -658,6 +658,14 @@ class SurfaceIntegrationMonitor(Monitor, ABC):
         description="Surfaces to exclude in the integration, if a volume monitor.",
     )
 
+    colocate: bool = pydantic.Field(
+        True,
+        title="Colocate Fields",
+        description="Defines whether fields are colocated to grid cell boundaries (i.e. to the "
+        "primal grid). Can be toggled for field recording monitors and is hard-coded for other "
+        "monitors depending on their specific function.",
+    )
+
     @property
     def integration_surfaces(self):
         """Surfaces of the monitor where fields will be recorded for subsequent integration."""