diff --git a/meanas/fdfd/functional.py b/meanas/fdfd/functional.py
index ddd074b..1637c51 100644
--- a/meanas/fdfd/functional.py
+++ b/meanas/fdfd/functional.py
@@ -226,6 +226,8 @@ def poynting_e_cross_h(dxes: dx_lists_t) -> Callable[[cfdfield, cfdfield], cfdfi
         For time-average power, call it as `f(E, H.conj())` and take `Re(...) / 2`.
     """
     def exh(e: cfdfield, h: cfdfield) -> cfdfield_t:
+        dtype = numpy.result_type(e, h, *dxes[0], *dxes[1])
+        e, h = numpy.asarray(e, dtype=dtype), numpy.asarray(h, dtype=dtype)
         s = numpy.empty_like(e)
         ex = e[0] * dxes[0][0][:, None, None]
         ey = e[1] * dxes[0][1][None, :, None]
diff --git a/meanas/fdfd/waveguide_2d.py b/meanas/fdfd/waveguide_2d.py
index c3c9053..958d1d0 100644
--- a/meanas/fdfd/waveguide_2d.py
+++ b/meanas/fdfd/waveguide_2d.py
@@ -836,7 +836,10 @@ def sensitivity(
         = H_{yx}^\star \partial_{\epsilon_i}(A_E) E_{xy} + H_{yx}^\star A_E \partial_{\epsilon_i} E_{xy}
     $$
 
-    However, $H_{yx}^\star$ is actually a left-eigenvector of $A_E$. This can be verified by inspecting
+    On a graded Yee grid, the left projection includes the quadrature weights:
+    $H_y^\star \Delta x_e \Delta y_h$ and $-H_x^\star \Delta x_h \Delta y_e$.
+    These weights multiply the left vector, before it acts on any derivative operators.
+    With these weights understood, $H_{yx}^\star$ is a left-eigenvector of $A_E$. This can be verified by inspecting
     the form of `operator_h` ($A_H$) and comparing its conjugate transpose to `operator_e` ($A_E$). Also, note
     $H_{yx}^\star \cdot E_{xy} = H^\star \times E$ recalls the mode orthogonality relation. See doi:10.5194/ars-9-85-201
     for a similar approach. Therefore,
@@ -883,11 +886,12 @@ def sensitivity(
     mu_x, mu_y, _mu_z = numpy.split(mu, 3)
     mu_yx = sparse.diags_array(numpy.hstack((mu_y, mu_x)))
 
-    da_exxhyy = vec(dxes[1][0][:, None] * dxes[0][1][None, :])
-    da_eyyhxx = vec(dxes[1][1][None, :] * dxes[0][0][:, None])
-    ev_xy = numpy.concatenate(numpy.split(e_norm, 3)[:2]) * numpy.concatenate([da_exxhyy, da_eyyhxx])
-    hx, hy, hz = numpy.split(h_norm, 3)
-    hv_yx_conj = numpy.conj(numpy.concatenate([hy, -hx]))
+    da_exxhyy = vec(dxes[0][0][:, None] * dxes[1][1][None, :])
+    da_eyyhxx = vec(dxes[1][0][:, None] * dxes[0][1][None, :])
+    ev_xy = numpy.concatenate(numpy.split(e_norm, 3)[:2])
+    hx, hy, _hz = numpy.split(h_norm, 3)
+    # The Yee-weighted H projection is the left eigenvector on a graded grid.
+    hv_yx_conj = numpy.conj(numpy.concatenate([hy, -hx])) * numpy.concatenate([da_exxhyy, da_eyyhxx])
 
     sens_xy1 = (hv_yx_conj @ (omega * omega * mu_yx)) * ev_xy
     sens_xy2 = (hv_yx_conj @ sparse.vstack((Dfx, Dfy)) @ eps_z_inv @ sparse.hstack((Dbx, Dby))) * ev_xy
diff --git a/meanas/fdfd/waveguide_3d.py b/meanas/fdfd/waveguide_3d.py
index 019a7c7..8955a1a 100644
--- a/meanas/fdfd/waveguide_3d.py
+++ b/meanas/fdfd/waveguide_3d.py
@@ -161,8 +161,9 @@ def mode_from_2d(
     e = unvec(ve, shape)
     h = unvec(vh, shape)
 
-    # Adjust for propagation direction
-    h *= polarity           # type: ignore          # mypy issue with numpy
+    # Reversing beta changes Ez and transverse H, while Hz is unchanged.
+    e[2] *= polarity
+    h[:2] *= polarity
 
     # Apply phase shift to H-field
     h[:2] *= numpy.exp(-1j * polarity * 0.5 * wavenumber * dx_prop)
diff --git a/meanas/test/test_fdfd_functional.py b/meanas/test/test_fdfd_functional.py
index 9c1f771..75d9d44 100644
--- a/meanas/test/test_fdfd_functional.py
+++ b/meanas/test/test_fdfd_functional.py
@@ -1,4 +1,5 @@
 import numpy
+import pytest
 
 from ..fdmath import unvec, vec
 from ..fdfd import functional, operators
@@ -120,3 +121,17 @@ def test_poynting_e_cross_h_matches_sparse_operator() -> None:
     functional_result = functional.poynting_e_cross_h(DXES)(E_FIELD, H_FIELD)
 
     assert_fields_match(functional_result, matrix_result)
+
+
+@pytest.mark.parametrize('wide_input', ['electric', 'magnetic', 'widths'])
+def test_poynting_uses_promoted_precision_before_products(wide_input) -> None:
+    shape = (2, 2, 2)
+    e = numpy.zeros((3, *shape), dtype=numpy.complex128 if wide_input == 'electric' else numpy.complex64)
+    h = numpy.zeros((3, *shape), dtype=numpy.complex128 if wide_input == 'magnetic' else numpy.complex64)
+    e[1], h[2] = 1e20, 1e20j
+    widths = [[numpy.ones(n, dtype=numpy.float64 if wide_input == 'widths' else numpy.float32) for n in shape] for _ in range(2)]
+    with numpy.errstate(over='raise', invalid='raise'):
+        result = functional.poynting_e_cross_h(widths)(e, h)
+    expected = unvec(operators.poynting_e_cross(vec(e.astype(complex)), widths) @ vec(h.astype(complex)), shape)
+    assert result.dtype == numpy.complex128
+    assert_fields_close(result, expected, atol=0, rtol=1e-14)
diff --git a/meanas/test/test_waveguide_2d_numerics.py b/meanas/test/test_waveguide_2d_numerics.py
index a7be517..d590be6 100644
--- a/meanas/test/test_waveguide_2d_numerics.py
+++ b/meanas/test/test_waveguide_2d_numerics.py
@@ -1,6 +1,7 @@
 import numpy
 from numpy.linalg import norm
 from numpy.testing import assert_allclose
+from scipy import linalg
 
 from ..fdmath import vec
 from ..fdfd import waveguide_2d
@@ -120,6 +121,33 @@ def test_waveguide_2d_sensitivity_matches_finite_difference() -> None:
     assert 0.4 < ratio < 1.8
 
 
+def test_sensitivity_matches_central_differences_on_graded_yee_mesh() -> None:
+    rng = numpy.random.default_rng(4)
+    shape = (4, 3)
+    omega = 0.7
+    epsilon = vec(1 + rng.uniform(0, 3, (3, *shape)))
+    magnetic_widths = [rng.uniform(0.2, 1.2, n) for n in shape]
+    electric_widths = [0.5 * (dx + numpy.roll(dx, -1)) for dx in magnetic_widths]
+    dxes = [electric_widths, magnetic_widths]
+    values, vectors = linalg.eig(waveguide_2d.operator_e(omega, dxes, epsilon).toarray())
+    selected = numpy.argmax(values.real)
+    beta = numpy.sqrt(values[selected])
+    e, h = waveguide_2d.normalized_fields_e(vectors[:, selected], beta, omega, dxes, epsilon)
+    sensitivity = waveguide_2d.sensitivity(e, h, beta, omega, dxes, epsilon)
+    delta = 1e-5
+    finite_difference = numpy.empty(epsilon.size, dtype=complex)
+    for index in range(epsilon.size):
+        roots = []
+        for direction in (-1, 1):
+            perturbed = epsilon.copy()
+            perturbed[index] += direction * delta
+            roots_squared = linalg.eigvals(waveguide_2d.operator_e(omega, dxes, perturbed).toarray())
+            candidates = numpy.sqrt(roots_squared)
+            roots.append(candidates[numpy.argmin(numpy.abs(candidates - beta))])
+        finite_difference[index] = (roots[1] - roots[0]) / (2 * delta)
+    assert_allclose(sensitivity, finite_difference, atol=1e-8, rtol=1e-5)
+
+
 def test_waveguide_2d_normalized_fields_h_are_finite_and_unit_normalized_with_mu() -> None:
     epsilon = build_asymmetric_epsilon()
     mu = build_mu_profile()
diff --git a/meanas/test/test_waveguide_mode_helpers.py b/meanas/test/test_waveguide_mode_helpers.py
index c375073..0531a59 100644
--- a/meanas/test/test_waveguide_mode_helpers.py
+++ b/meanas/test/test_waveguide_mode_helpers.py
@@ -2,17 +2,48 @@ import contextlib
 import io
 import numpy
 from numpy.linalg import norm
+from scipy import linalg
 import pytest
 import warnings
 
 from ..fdmath import vec, unvec
 from ..fdmath.operators import deriv_forward, deriv_back
+from ..fdmath.functional import curl_forward, curl_back
 from ..fdfd import waveguide_2d, waveguide_3d, waveguide_cyl
 
 
 OMEGA = 1 / 1500
 
 
+@pytest.mark.parametrize('axis', [0, 1, 2])
+@pytest.mark.parametrize('polarity', [-1, 1])
+def test_embedded_hybrid_mode_obeys_both_maxwell_curls(axis, polarity) -> None:
+    omega = 0.7
+    # Cyclic rotation preserves the orientation of the vector cross product.
+    permutation = numpy.roll(numpy.arange(3), axis + 1)
+    canonical = numpy.ones((3, 4, 3, 9))
+    canonical[:, 1:3, :, :] = 4
+    epsilon = canonical[permutation].transpose(0, *(permutation + 1))
+    shape = epsilon.shape[1:]
+    dxes = [[numpy.full(n, 0.4) for n in shape] for _ in range(2)]
+    cross_dxes = [[numpy.full(4, 0.4), numpy.full(3, 0.4)] for _ in range(2)]
+    values, vectors = linalg.eig(waveguide_2d.operator_e(omega, cross_dxes, vec(canonical[:, :, :, 0])).toarray())
+    selected = numpy.argmax(values.real)
+    slices = [slice(None)] * 3
+    slices[axis] = slice(4, 5)
+    mode = waveguide_3d.mode_from_2d(
+        vectors[:, selected], numpy.sqrt(values[selected]), omega=omega, dxes=dxes,
+        axis=axis, polarity=polarity, slices=slices, epsilon=epsilon,
+    )
+    e, h = [waveguide_3d.expand_e(mode[field], mode['wavenumber'], dxes, axis, polarity, slices) for field in ('E', 'H')]
+    assert norm(h[axis]) > 0.01 * norm(h)
+    region = [slice(None)] * 4
+    region[axis + 1] = slice(1, -1)
+    ce, ch = curl_forward(dxes[0])(e), curl_back(dxes[1])(h)
+    assert relative_residual(ce[tuple(region)], (-1j * omega * h)[tuple(region)]) < 1e-11
+    assert relative_residual(ch[tuple(region)], (1j * omega * epsilon * e)[tuple(region)]) < 1e-11
+
+
 def build_waveguide_3d_mode(
         *,
         slice_start: int,