diff --git a/crates/ppvm-python-native/src/interface.rs b/crates/ppvm-python-native/src/interface.rs
index 88afc07e3..78f884a77 100644
--- a/crates/ppvm-python-native/src/interface.rs
+++ b/crates/ppvm-python-native/src/interface.rs
@@ -210,6 +210,17 @@ macro_rules! create_interface {
                 self.inner.truncate();
             }
 
+            // U(1)-conserving exchange/Heisenberg-style gates
+            pub fn exchange(&mut self, addr0: usize, addr1: usize, theta: f64) {
+                self.inner.exchange(addr0, addr1, theta);
+                self.inner.truncate();
+            }
+
+            pub fn xyzz(&mut self, addr0: usize, addr1: usize, theta_xy: f64, theta_zz: f64) {
+                self.inner.xyzz(addr0, addr1, theta_xy, theta_zz);
+                self.inner.truncate();
+            }
+
             // noise
             pub fn pauli_error(&mut self, addr0: usize, p: [f64; 3]) {
                 self.inner.pauli_error(addr0, p);
diff --git a/crates/ppvm-runtime/src/sum/mod.rs b/crates/ppvm-runtime/src/sum/mod.rs
index 280dd5a60..5e356ee7d 100644
--- a/crates/ppvm-runtime/src/sum/mod.rs
+++ b/crates/ppvm-runtime/src/sum/mod.rs
@@ -10,6 +10,7 @@ mod proj;
 mod rot1;
 mod rot2;
 mod trace;
+mod u1;
 
 #[cfg(feature = "approx")]
 mod approx;
diff --git a/crates/ppvm-runtime/src/sum/u1.rs b/crates/ppvm-runtime/src/sum/u1.rs
new file mode 100644
index 000000000..37ab8a4b1
--- /dev/null
+++ b/crates/ppvm-runtime/src/sum/u1.rs
@@ -0,0 +1,375 @@
+// SPDX-FileCopyrightText: 2026 The PPVM Authors
+// SPDX-License-Identifier: Apache-2.0
+
+use crate::traits::*;
+use crate::{config::Config, sum::PauliSum};
+
+/// Z-magnetization-conserving two-qubit propagation on a `PauliSum`.
+///
+/// The current implementation expresses each gate as a composition of
+/// the existing `rxx` / `ryy` / `rzz` primitives, which is mathematically
+/// equivalent to a fully fused single-pass implementation because the
+/// XX, YY, and ZZ generators all pairwise commute. A future revision is
+/// free to specialize this impl with a fused single `map_insert_multiple`
+/// pass — the public signature is stable.
+impl<T: Config> U1Conserving<T> for PauliSum<T>
+where
+    PauliSum<T>: RotationTwo<T>,
+    T::Coeff: Clone,
+{
+    fn exchange(&mut self, a: usize, b: usize, theta: impl Into<T::Coeff>) {
+        let theta: T::Coeff = theta.into();
+        self.rxx(a, b, theta.clone());
+        self.ryy(a, b, theta);
+    }
+
+    fn xyzz(
+        &mut self,
+        a: usize,
+        b: usize,
+        theta_xy: impl Into<T::Coeff>,
+        theta_zz: impl Into<T::Coeff>,
+    ) {
+        self.exchange(a, b, theta_xy);
+        self.rzz(a, b, theta_zz);
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::config::fxhash::ByteF64;
+    use crate::strategy::CoefficientThreshold;
+    use std::f64::consts::PI;
+
+    // Use CoefficientThreshold so the propagated maps drop the cos(θ)·0 and
+    // sin(θ)·0 ghost entries that rotate_2 leaves behind at θ = 0 or when a
+    // sub-branch happens to vanish. Without this, `assert_eq!` would compare
+    // maps that differ only in formally-zero terms.
+    type Cfg = ByteF64<2, CoefficientThreshold>;
+
+    fn ps(n: usize, term: &str, coeff: f64) -> PauliSum<Cfg> {
+        let mut s: PauliSum<Cfg> = PauliSum::builder()
+            .n_qubits(n)
+            .strategy(CoefficientThreshold(1e-12))
+            .build();
+        s += (term, coeff);
+        s
+    }
+
+    /// `exchange(a, b, 0)` is the identity.
+    #[test]
+    fn exchange_zero_angle_is_identity() {
+        for term in ["II", "IZ", "ZI", "ZZ", "XY", "YX", "XX", "YY"] {
+            let expect = ps(2, term, 1.0);
+            let mut got = ps(2, term, 1.0);
+            got.exchange(0, 1, 0.0);
+            got.truncate();
+            assert_eq!(got, expect, "exchange(0) altered {}", term);
+        }
+    }
+
+    /// `exchange(a, b, θ)` matches `rxx(θ)` followed by `ryy(θ)`.
+    #[test]
+    fn exchange_matches_rxx_then_ryy() {
+        let thetas = [0.0_f64, 0.123, 0.7, -0.4, PI / 3.0, PI / 2.0];
+        let terms = [
+            "II", "IZ", "ZI", "ZZ", "XY", "YX", "XX", "YY", "IX", "XI", "IY", "YI", "XZ", "ZX",
+            "YZ", "ZY",
+        ];
+        for theta in thetas {
+            for term in terms {
+                let mut fused = ps(2, term, 1.0);
+                let mut composed = ps(2, term, 1.0);
+                fused.exchange(0, 1, theta);
+                composed.rxx(0, 1, theta);
+                composed.ryy(0, 1, theta);
+                fused.truncate();
+                composed.truncate();
+                assert_eq!(
+                    fused, composed,
+                    "exchange disagrees on {} at θ={}",
+                    term, theta
+                );
+            }
+        }
+    }
+
+    /// `exchange` preserves the identity term: `I ↦ I` with the same coefficient.
+    #[test]
+    fn exchange_preserves_identity() {
+        let expect = ps(2, "II", 0.7);
+        let mut s = ps(2, "II", 0.7);
+        s.exchange(0, 1, 0.42);
+        s.truncate();
+        assert_eq!(s, expect);
+    }
+
+    /// `Z_a + Z_b` commutes with `X_a X_b + Y_a Y_b`, so total magnetization
+    /// must propagate trivially through `exchange`.
+    #[test]
+    fn exchange_preserves_total_z() {
+        let mut expect: PauliSum<Cfg> = PauliSum::builder()
+            .n_qubits(2)
+            .strategy(CoefficientThreshold(1e-12))
+            .build();
+        expect += ("IZ", 1.0);
+        expect += ("ZI", 1.0);
+        let mut s = expect.clone();
+        s.exchange(0, 1, 0.37);
+        s.exchange(0, 1, -1.1);
+        s.truncate();
+        assert_eq!(s, expect, "total Z was perturbed by exchange");
+    }
+
+    /// `xyzz(θ_xy, θ_zz)` matches `exchange(θ_xy)` followed by `rzz(θ_zz)`.
+    #[test]
+    fn xyzz_matches_exchange_then_rzz() {
+        let pairs = [
+            (0.0_f64, 0.0_f64),
+            (0.3, 0.1),
+            (-0.4, 0.7),
+            (PI / 5.0, -PI / 7.0),
+        ];
+        let terms = ["IZ", "ZI", "XY", "YX", "XX", "YY", "ZZ", "IX", "YI"];
+        for (theta_xy, theta_zz) in pairs {
+            for term in terms {
+                let mut combined = ps(2, term, 1.0);
+                let mut stepwise = ps(2, term, 1.0);
+                combined.xyzz(0, 1, theta_xy, theta_zz);
+                stepwise.exchange(0, 1, theta_xy);
+                stepwise.rzz(0, 1, theta_zz);
+                combined.truncate();
+                stepwise.truncate();
+                assert_eq!(
+                    combined, stepwise,
+                    "xyzz disagrees on {} at θ_xy={}, θ_zz={}",
+                    term, theta_xy, theta_zz
+                );
+            }
+        }
+    }
+
+    /// `xyzz` with both angles zero is the identity.
+    #[test]
+    fn xyzz_zero_angles_is_identity() {
+        for term in ["II", "IZ", "ZI", "XY", "YX", "ZZ"] {
+            let expect = ps(2, term, 1.0);
+            let mut got = ps(2, term, 1.0);
+            got.xyzz(0, 1, 0.0, 0.0);
+            got.truncate();
+            assert_eq!(got, expect, "xyzz(0,0) altered {}", term);
+        }
+    }
+
+    /// The magnetization difference and the "current" operator
+    /// `Y_a X_b − X_a Y_b` form a closed two-dimensional U(1) sector;
+    /// `exchange(θ)` rotates between them with the doubled-angle
+    /// `cos(2θ)` / `sin(2θ)` (the ppvm angle convention is
+    /// `exp(−i θ/2 G)`, so two stacked rotations contribute a 2θ).
+    #[test]
+    fn exchange_rotates_magnetization_difference_into_current() {
+        let theta = 0.3_f64;
+        let mut s: PauliSum<Cfg> = PauliSum::builder()
+            .n_qubits(2)
+            .strategy(CoefficientThreshold(1e-12))
+            .build();
+        s += ("ZI", 0.5);
+        s += ("IZ", -0.5);
+        s.exchange(0, 1, theta);
+        s.truncate();
+
+        let c = (2.0 * theta).cos();
+        let si = (2.0 * theta).sin();
+        let mut want: PauliSum<Cfg> = PauliSum::builder()
+            .n_qubits(2)
+            .strategy(CoefficientThreshold(1e-12))
+            .build();
+        want += ("ZI", 0.5 * c);
+        want += ("IZ", -0.5 * c);
+        want += ("YX", 0.5 * si);
+        want += ("XY", -0.5 * si);
+
+        // Compare via terms() under a tolerance — the maps store f64s, so
+        // structural equality is too strict here.
+        let mut got_terms: Vec<(String, f64)> =
+            s.data().iter().map(|(k, v)| (k.to_string(), *v)).collect();
+        let mut want_terms: Vec<(String, f64)> = want
+            .data()
+            .iter()
+            .map(|(k, v)| (k.to_string(), *v))
+            .collect();
+        got_terms.sort_by(|a, b| a.0.cmp(&b.0));
+        want_terms.sort_by(|a, b| a.0.cmp(&b.0));
+        assert_eq!(got_terms.len(), want_terms.len());
+        for ((kg, vg), (kw, vw)) in got_terms.iter().zip(want_terms.iter()) {
+            assert_eq!(kg, kw);
+            assert!(
+                (vg - vw).abs() < 1e-12,
+                "coeff for {} differs: got {}, want {}",
+                kg,
+                vg,
+                vw
+            );
+        }
+    }
+
+    // ===========================================================================
+    // Truncation-robust conservation hardening tests
+    // ===========================================================================
+    //
+    // These four tests pin down the realistic floating-point-precision guarantee
+    // that `exchange` / `xyzz` / `rzz` provide on observables built from `{I, Z}`
+    // Pauli strings (e.g. `Σ_i Z_i`, `Σ Z_i Z_j`). The conserved-sector
+    // coefficients recover to their starting values modulo per-gate ε (≈ 1e-15),
+    // so we compare under a `1e-10` tolerance — three orders of magnitude above
+    // accumulated drift for a short circuit, ten orders of magnitude below the
+    // conserved-coefficient magnitude (1.0).
+
+    use crate::strategy::{CombinedStrategy, MaxPauliWeight};
+
+    /// Sort-and-compare PauliSums under a per-term tolerance.
+    fn assert_close<T>(left: &PauliSum<T>, right: &PauliSum<T>, tol: f64)
+    where
+        T: crate::config::Config<Coeff = f64>,
+        T::Map: for<'a> crate::traits::ACMapIter<
+                'a,
+                Item = (&'a <T as crate::config::Config>::PauliWordType, &'a f64),
+            >,
+    {
+        let mut left_terms: Vec<(String, f64)> = left
+            .data()
+            .iter()
+            .map(|(k, v)| (k.to_string(), *v))
+            .collect();
+        let mut right_terms: Vec<(String, f64)> = right
+            .data()
+            .iter()
+            .map(|(k, v)| (k.to_string(), *v))
+            .collect();
+        left_terms.sort_by(|a, b| a.0.cmp(&b.0));
+        right_terms.sort_by(|a, b| a.0.cmp(&b.0));
+        assert_eq!(
+            left_terms.len(),
+            right_terms.len(),
+            "term count differs: left={:?}, right={:?}",
+            left_terms,
+            right_terms
+        );
+        for ((kl, vl), (kr, vr)) in left_terms.iter().zip(right_terms.iter()) {
+            assert_eq!(kl, kr, "key mismatch: {} vs {}", kl, kr);
+            assert!(
+                (vl - vr).abs() < tol,
+                "coeff for {} differs by more than {}: got {}, want {}",
+                kl,
+                tol,
+                vl,
+                vr
+            );
+        }
+    }
+
+    /// `Σ Z_i` is preserved across a chain of `exchange` calls when truncation
+    /// uses an aggressive `CoefficientThreshold` (well below the conserved
+    /// coefficient magnitude). The transient cross terms produced internally
+    /// by each gate cancel back to ε before truncation, so the cutoff (0.5)
+    /// drops only the ε residues and leaves `Σ Z_i` intact.
+    #[test]
+    fn exchange_preserves_total_z_under_coefficient_truncation() {
+        type C = ByteF64<1, CoefficientThreshold>;
+        let mut s: PauliSum<C> = PauliSum::builder()
+            .n_qubits(4)
+            .strategy(CoefficientThreshold(0.5))
+            .build();
+        let mut expect: PauliSum<C> = PauliSum::builder()
+            .n_qubits(4)
+            .strategy(CoefficientThreshold(0.5))
+            .build();
+        for term in ["ZIII", "IZII", "IIZI", "IIIZ"] {
+            s += (term, 1.0);
+            expect += (term, 1.0);
+        }
+        for (a, b) in [(0, 1), (1, 2), (2, 3)] {
+            s.exchange(a, b, 0.37);
+            s.truncate();
+        }
+        assert_close(&s, &expect, 1e-10);
+    }
+
+    /// Same conservation under `MaxPauliWeight(1)`: the cross terms have weight
+    /// 2 and are dropped by the discrete weight check, independent of any
+    /// floating-point sensitivity in the threshold comparison.
+    #[test]
+    fn exchange_preserves_total_z_under_weight_truncation() {
+        type C = ByteF64<1, CombinedStrategy<CoefficientThreshold, MaxPauliWeight>>;
+        let strat = CombinedStrategy(CoefficientThreshold(1e-12), MaxPauliWeight(1));
+        let mut s: PauliSum<C> = PauliSum::builder().n_qubits(4).strategy(strat).build();
+        let mut expect: PauliSum<C> = PauliSum::builder().n_qubits(4).strategy(strat).build();
+        for term in ["ZIII", "IZII", "IIZI", "IIIZ"] {
+            s += (term, 1.0);
+            expect += (term, 1.0);
+        }
+        for (a, b) in [(0, 1), (1, 2), (2, 3)] {
+            s.exchange(a, b, 0.41);
+            s.truncate();
+        }
+        assert_close(&s, &expect, 1e-10);
+    }
+
+    /// `Σ_{i<j} Z_i Z_j` (sum over *all* pairs, not just nearest neighbours)
+    /// commutes with `XX+YY+ZZ` on every edge, so `xyzz` on any edge
+    /// preserves the full sum. The nearest-neighbour sum is *not*
+    /// conserved: the cross terms that exchange generates from
+    /// `Z_aZ_c + Z_bZ_c` cancel only when the matching `Z_aZ_c` and
+    /// `Z_bZ_c` pair is present for every "third site" `c`. This
+    /// catches any regression in `rzz` semantics under aggressive
+    /// truncation and is a tighter check on `xyzz`'s ZZ piece than the
+    /// single-Z tests above.
+    #[test]
+    fn zz_correlation_preserves_under_truncation() {
+        type C = ByteF64<1, CombinedStrategy<CoefficientThreshold, MaxPauliWeight>>;
+        let strat = CombinedStrategy(CoefficientThreshold(0.5), MaxPauliWeight(2));
+        let mut s: PauliSum<C> = PauliSum::builder().n_qubits(4).strategy(strat).build();
+        let mut expect: PauliSum<C> = PauliSum::builder().n_qubits(4).strategy(strat).build();
+        // All C(4, 2) = 6 unordered pairs.
+        for term in ["ZZII", "ZIZI", "ZIIZ", "IZZI", "IZIZ", "IIZZ"] {
+            s += (term, 1.0);
+            expect += (term, 1.0);
+        }
+        for (a, b) in [(0, 1), (1, 2), (2, 3)] {
+            s.xyzz(a, b, 0.21, 0.07);
+            s.truncate();
+        }
+        assert_close(&s, &expect, 1e-10);
+    }
+
+    /// One full Trotter cycle (xyzz on every edge + rz on every site) on
+    /// `Σ Z_i`, repeated three times, under combined coefficient + weight
+    /// truncation. The Trotter helper itself lives in Python, so we
+    /// replicate its gate order manually here.
+    #[test]
+    fn u1_trotter_total_z_under_combined_truncation() {
+        type C = ByteF64<1, CombinedStrategy<CoefficientThreshold, MaxPauliWeight>>;
+        let strat = CombinedStrategy(CoefficientThreshold(0.1), MaxPauliWeight(1));
+        let mut s: PauliSum<C> = PauliSum::builder().n_qubits(5).strategy(strat).build();
+        let mut expect: PauliSum<C> = PauliSum::builder().n_qubits(5).strategy(strat).build();
+        for term in ["ZIIII", "IZIII", "IIZII", "IIIZI", "IIIIZ"] {
+            s += (term, 1.0);
+            expect += (term, 1.0);
+        }
+        let theta_xy = 0.15;
+        let theta_zz = 0.05;
+        let h = 0.03;
+        for _ in 0..3 {
+            for (a, b) in [(0, 1), (1, 2), (2, 3), (3, 4)] {
+                s.xyzz(a, b, theta_xy, theta_zz);
+                s.truncate();
+            }
+            for site in 0..5 {
+                s.rz(site, h);
+                s.truncate();
+            }
+        }
+        assert_close(&s, &expect, 1e-10);
+    }
+}
diff --git a/crates/ppvm-runtime/src/traits/branch/mod.rs b/crates/ppvm-runtime/src/traits/branch/mod.rs
index 0c50ddd43..9521b8f6b 100644
--- a/crates/ppvm-runtime/src/traits/branch/mod.rs
+++ b/crates/ppvm-runtime/src/traits/branch/mod.rs
@@ -6,6 +6,7 @@ mod proj;
 mod rot1;
 mod rot2;
 mod tgate;
+mod u1;
 mod u3;
 
 pub use crx::CRx;
@@ -13,4 +14,5 @@ pub use proj::Projection;
 pub use rot1::RotationOne;
 pub use rot2::RotationTwo;
 pub use tgate::TGate;
+pub use u1::U1Conserving;
 pub use u3::U3Gate;
diff --git a/crates/ppvm-runtime/src/traits/branch/u1.rs b/crates/ppvm-runtime/src/traits/branch/u1.rs
new file mode 100644
index 000000000..aef7f436d
--- /dev/null
+++ b/crates/ppvm-runtime/src/traits/branch/u1.rs
@@ -0,0 +1,77 @@
+// SPDX-FileCopyrightText: 2026 The PPVM Authors
+// SPDX-License-Identifier: Apache-2.0
+
+use crate::config::Config;
+
+/// Z-magnetization-conserving two-qubit gates.
+///
+/// These gates generate unitaries that commute with the total
+/// Z-magnetization `Σ_i Z_i`, so they preserve the U(1) symmetry sector
+/// of any observable that already commutes with `Σ_i Z_i`. They are the
+/// natural building blocks for XY, Heisenberg, and related spin-model
+/// dynamics.
+///
+/// Like every other Pauli-propagation gate in `ppvm-runtime`, these
+/// methods act in the **Heisenberg picture** — they conjugate the
+/// observable held in a [`PauliSum`](crate::sum::PauliSum) by the gate's
+/// unitary. Compose them with the rest of the circuit in reverse order,
+/// the same way `rxx` / `ryy` / `rzz` are composed.
+///
+/// The default implementations express each gate as a composition of
+/// existing `rxx` / `ryy` / `rzz` calls so that any backend implementing
+/// [`RotationTwo`](crate::traits::RotationTwo) automatically supports
+/// U(1)-conserving dynamics. Backends are free to override with a fused
+/// implementation when they can avoid the intermediate branching.
+///
+/// # Conservation under truncation
+///
+/// `exchange`, `xyzz`, and `rzz` all commute with the total Z
+/// magnetization `Σ_k Z_k`. As a consequence, an observable built from
+/// `{I, Z}`-only Pauli strings (`Σ_i Z_i`, `Σ_{i<j} Z_iZ_j`, products
+/// thereof) is conserved under propagation **up to per-gate
+/// floating-point precision** — typically `< 1e-14` per gate,
+/// accumulating linearly over a circuit. The off-diagonal cross terms
+/// that propagation introduces from individual `Z_i` summands cancel by
+/// coefficient when the contributions are summed, so they never grow
+/// large enough to be retained by a sensible truncation cutoff.
+///
+/// Standard truncation thresholds — anything well below the conserved
+/// coefficient magnitude (`~1`) and well above the per-gate ε floor —
+/// preserve conservation in practice. Setting
+/// [`CoefficientThreshold`](crate::strategy::CoefficientThreshold) close
+/// to 1, or running so many gates that the accumulated drift becomes
+/// comparable to the threshold, can break this guarantee. For more
+/// general U(1)-conserving Hermitian observables, conjugate-paired
+/// Pauli strings (`Y_iX_j` ↔ `X_iY_j`, etc.) carry coefficients of
+/// equal magnitude, so symmetric coefficient truncation keeps the
+/// observable inside the U(1) sub-algebra even when individual
+/// coefficients are perturbed.
+pub trait U1Conserving<T: Config> {
+    /// Fused XY exchange rotation:
+    /// `exp(-i θ/2 (X_a X_b + Y_a Y_b))`.
+    ///
+    /// `X_a X_b` and `Y_a Y_b` commute, so this is mathematically
+    /// equivalent to `rxx(a, b, theta)` followed by `ryy(a, b, theta)`.
+    /// Exposed as a single routine so backends can minimize intermediate
+    /// branching and so users can write U(1)-symmetric dynamics without
+    /// re-deriving the XX / YY decomposition every time.
+    fn exchange(&mut self, a: usize, b: usize, theta: impl Into<T::Coeff>);
+
+    /// Combined XY + ZZ Heisenberg-style interaction:
+    ///
+    /// ```text
+    /// exp(-i θ_xy/2 (X_a X_b + Y_a Y_b)) · exp(-i θ_zz/2 Z_a Z_b)
+    /// ```
+    ///
+    /// The two generators commute, so the factorization is exact. With
+    /// `theta_zz = 0` this reduces to [`exchange`](Self::exchange); with
+    /// `theta_xy = 0` it reduces to `rzz`. With both non-zero it
+    /// implements one Trotter slice of an XXZ-style Hamiltonian.
+    fn xyzz(
+        &mut self,
+        a: usize,
+        b: usize,
+        theta_xy: impl Into<T::Coeff>,
+        theta_zz: impl Into<T::Coeff>,
+    );
+}
diff --git a/crates/ppvm-runtime/src/traits/mod.rs b/crates/ppvm-runtime/src/traits/mod.rs
index 853e435e8..0bc3dd2fc 100644
--- a/crates/ppvm-runtime/src/traits/mod.rs
+++ b/crates/ppvm-runtime/src/traits/mod.rs
@@ -14,7 +14,7 @@ mod strategy;
 mod trace;
 mod word_trait;
 
-pub use branch::{CRx, Projection, RotationOne, RotationTwo, TGate, U3Gate};
+pub use branch::{CRx, Projection, RotationOne, RotationTwo, TGate, U1Conserving, U3Gate};
 pub use clifford::{Clifford, CliffordExtensions};
 pub use coefficient::{Coefficient, ComplexCoefficient};
 pub use map::{
diff --git a/docs/notebooks/u1_heisenberg.py b/docs/notebooks/u1_heisenberg.py
new file mode 100644
index 000000000..23c775ebf
--- /dev/null
+++ b/docs/notebooks/u1_heisenberg.py
@@ -0,0 +1,108 @@
+# ---
+# jupyter:
+#   jupytext:
+#     cell_metadata_filter: -all
+#     custom_cell_magics: kql
+#     text_representation:
+#       extension: .py
+#       format_name: percent
+#       format_version: '1.3'
+#       jupytext_version: 1.19.1
+#   kernelspec:
+#     display_name: ppvm (3.12.12)
+#     language: python
+#     name: python3
+# ---
+
+# %% [markdown]
+# # U(1)-Conserving Trotter Dynamics: XY / Heisenberg Chain
+#
+# This example shows how to use ppvm's U(1)-symmetric gate helpers to simulate
+# z-magnetization-conserving spin dynamics. The Hamiltonian we evolve is the
+# XXZ Heisenberg chain with an on-site Z field,
+#
+# $$ H = \sum_{(i,j)} J_{ij} \big( X_i X_j + Y_i Y_j \big) + \sum_{(i,j)} \Delta_{ij} Z_i Z_j + \sum_i h_i Z_i . $$
+#
+# Each term commutes with the total magnetization $\sum_k Z_k$, so the
+# dynamics preserve the U(1) symmetry sector of any starting observable that
+# already lives in that sector. ppvm exposes three pieces of API for this:
+#
+# 1. `PauliSum.exchange(a, b, theta)` — fused $\mathrm{exp}\!\big({-i\,\theta/2\,(X_a X_b + Y_a Y_b)}\big)$.
+# 2. `PauliSum.xyzz(a, b, theta_xy, theta_zz)` — combined XY + ZZ interaction.
+# 3. `PauliSum.apply_u1_trotter_step(edges, theta_xy, theta_zz, fields_z)` —
+#    a single Trotter slice expressed in terms of the above.
+#
+# All three act in the **Heisenberg picture**, like every other PauliSum gate.
+
+# %%
+from ppvm import PauliSum
+
+# %% [markdown]
+# ## Parameters
+
+# %%
+n = 4
+edges = [(i, i + 1) for i in range(n - 1)]
+J_xy = 0.5
+Delta = 0.2
+field = 0.1
+dt = 0.1
+n_steps = 5
+
+# %% [markdown]
+# ## Observable: total magnetization
+#
+# We initialise the observable $O = \sum_i Z_i$. Because the XXZ Hamiltonian
+# commutes with $\sum_i Z_i$, the propagated observable must stay numerically
+# identical to its starting form (up to truncation noise) — a sharp built-in
+# check of U(1) symmetry.
+
+# %%
+state = PauliSum.new(
+    n_qubits=n,
+    terms=[f"Z{i}" for i in range(n)],
+    min_abs_coeff=1e-12,
+)
+initial_terms = dict(state.terms)
+
+# %% [markdown]
+# ## Trotterized evolution
+#
+# Each `apply_u1_trotter_step` call applies, for every edge $(i,j)$:
+# `xyzz(i, j, J_xy * dt, Delta * dt)`, then a uniform $R_z(field \cdot dt)$ on
+# every site. The full circuit is one Heisenberg-picture sweep per call.
+
+# %%
+for _ in range(n_steps):
+    state.apply_u1_trotter_step(
+        edges=edges,
+        theta_xy=J_xy * dt,
+        theta_zz=Delta * dt,
+        fields_z=[field * dt] * n,
+    )
+
+# %% [markdown]
+# ## Symmetry check
+#
+# Total magnetization is conserved, so the final and initial PauliSums are
+# the same map of Pauli strings to coefficients.
+
+# %%
+final_terms = dict(state.terms)
+for term, coeff in initial_terms.items():
+    assert abs(final_terms.get(term, 0.0) - coeff) < 1e-10, term
+print("Total Z preserved across", n_steps, "Trotter steps.")
+
+# %% [markdown]
+# ## A starting observable outside the conserved sector
+#
+# If we instead start from $Z_0 - Z_1$ — the magnetization *difference* — the
+# dynamics rotate it inside the two-dimensional U(1) sector spanned by
+# $Z_0 - Z_1$ and the current-like operator $Y_0 X_1 - X_0 Y_1$. This is the
+# operator-level manifestation of "spin currents flowing across an XY bond".
+
+# %%
+diff = PauliSum.new(2, [("Z0", 0.5), ("Z1", -0.5)], min_abs_coeff=1e-12)
+diff.exchange(0, 1, theta=0.4)
+for name, coeff in diff.terms:
+    print(f"  {name}: {coeff:+.4f}")
diff --git a/ppvm-python/src/ppvm/mixins.py b/ppvm-python/src/ppvm/mixins.py
index 3e798746b..2f85c1005 100644
--- a/ppvm-python/src/ppvm/mixins.py
+++ b/ppvm-python/src/ppvm/mixins.py
@@ -155,6 +155,72 @@ def rzz(self, addr0: int, addr1: int, theta: float) -> None:
         """
         self._interface.rzz(addr0, addr1, theta)
 
+    # U(1)-conserving two-qubit gates for XY / Heisenberg-style dynamics.
+    def exchange(self, addr0: int, addr1: int, theta: float) -> None:
+        """Apply the fused XY exchange rotation to two qubits.
+
+        ```math
+        \\mathrm{exchange}(\\theta) = e^{-i \\frac{\\theta}{2} (X \\otimes X + Y \\otimes Y)}
+        ```
+
+        The `XX` and `YY` generators commute, so this is mathematically
+        equivalent to ``rxx(addr0, addr1, theta)`` followed by
+        ``ryy(addr0, addr1, theta)``, but exposed as a single routine for
+        building U(1)-symmetric (z-magnetization-conserving) dynamics such
+        as XY or Heisenberg spin models without manually composing the
+        XX / YY decomposition. Like every other PauliSum gate this acts in
+        the Heisenberg picture.
+
+        Conservation under truncation:
+            ``exchange`` commutes with the total Z magnetization
+            :math:`\\sum_k Z_k`. Observables built from `{I, Z}`-only Pauli
+            strings (e.g. :math:`\\sum_i Z_i`, :math:`\\sum_{i<j} Z_i Z_j`)
+            are conserved through repeated ``exchange`` calls up to
+            per-gate floating-point precision (~1e-15 per gate). Standard
+            truncation settings — `min_abs_coeff` well below 1 and any
+            `max_pauli_weight` that accommodates the starting observable —
+            preserve this guarantee. See :py:meth:`apply_u1_trotter_step`
+            for the broader story when composing several U(1) gates.
+
+        Args:
+            addr0: The index of the first qubit.
+            addr1: The index of the second qubit.
+            theta: The rotation angle in radians.
+        """
+        self._interface.exchange(addr0, addr1, theta)
+
+    def xyzz(self, addr0: int, addr1: int, theta_xy: float, theta_zz: float) -> None:
+        """Apply a combined XY exchange + ZZ interaction to two qubits.
+
+        ```math
+        e^{-i \\frac{\\theta_{xy}}{2} (X \\otimes X + Y \\otimes Y)}
+        \\cdot e^{-i \\frac{\\theta_{zz}}{2} Z \\otimes Z}
+        ```
+
+        The two factors commute, so the ordering of the (Heisenberg-picture)
+        XX+YY and ZZ pieces does not matter mathematically. This is the
+        convenience method to use for one site-pair of an XXZ-style
+        z-magnetization-conserving Hamiltonian.
+
+        Conservation under truncation:
+            Like :py:meth:`exchange`, ``xyzz`` commutes with the total Z
+            magnetization :math:`\\sum_k Z_k`. The same truncation-robust
+            guarantee applies: `{I, Z}`-polynomial observables (e.g.
+            :math:`\\sum_i Z_i`, :math:`\\sum_{i<j} Z_i Z_j` over all
+            pairs) are conserved through repeated ``xyzz`` calls up to
+            per-gate floating-point precision, as long as the truncation
+            cutoff stays comfortably below the conserved coefficient
+            magnitude. See :py:meth:`apply_u1_trotter_step` for the
+            broader discussion.
+
+        Args:
+            addr0: The index of the first qubit.
+            addr1: The index of the second qubit.
+            theta_xy: Rotation angle for the `XX + YY` (exchange) generator.
+            theta_zz: Rotation angle for the `ZZ` generator.
+        """
+        self._interface.xyzz(addr0, addr1, theta_xy, theta_zz)
+
 
 class CliffordExtensionMixin:
     _interface: Any
diff --git a/ppvm-python/src/ppvm/paulisum.py b/ppvm-python/src/ppvm/paulisum.py
index 4ad16248b..21d399220 100644
--- a/ppvm-python/src/ppvm/paulisum.py
+++ b/ppvm-python/src/ppvm/paulisum.py
@@ -129,6 +129,11 @@ class PauliSum(
         # Expectation value of ZZZ for GHZ state is 0
         result = ps.overlap_with_zero()
         ```
+
+        For z-magnetization-conserving (XY / Heisenberg-style) dynamics, use
+        the U(1) helpers — `exchange`, `xyzz`, and `apply_u1_trotter_step` —
+        which propagate the appropriate XX+YY and ZZ generators in
+        Heisenberg-picture order without manual rxx/ryy/rzz composition.
     """
 
     initial_terms: Sequence[str]
@@ -379,6 +384,134 @@ def amplitude_damping(self, addr0: int, gamma: float):
 
         self._interface.amplitude_damping(addr0, gamma)
 
+    def apply_u1_trotter_step(
+        self,
+        edges: "Sequence[tuple[int, int] | tuple[int, int, float]]",
+        theta_xy: "float | Sequence[float]",
+        theta_zz: "float | Sequence[float] | None" = None,
+        fields_z: "float | Sequence[float] | None" = None,
+    ) -> None:
+        """Apply one Trotter slice of a U(1)-symmetric (z-magnetization-conserving)
+        Hamiltonian.
+
+        For each edge `(i, j)` (or `(i, j, J)`) in `edges` this applies the
+        Heisenberg-picture image of
+
+        ```math
+        e^{-i \\frac{\\theta_{xy}}{2} (X_i X_j + Y_i Y_j)}
+        \\cdot e^{-i \\frac{\\theta_{zz}}{2} Z_i Z_j}
+        ```
+
+        followed by per-site `rz(\\theta_h)` rotations for any non-zero entries
+        in `fields_z`. The combined unitary commutes with `\\sum_k Z_k`, so any
+        observable in the total-Z sector stays in that sector. The angle
+        convention matches the rest of PauliSum: each rotation is
+        `exp(-i θ/2 G)`.
+
+        Args:
+            edges: A sequence of `(i, j)` pairs or `(i, j, J_ij)` triples
+                identifying the two-site interactions. The optional third
+                entry scales the XX+YY and ZZ angles for that edge.
+            theta_xy: The base `XX + YY` angle for every edge, or a per-edge
+                sequence of the same length as `edges`. Set to zero to skip
+                the exchange step.
+            theta_zz: Optional `ZZ` angle. Same shape rules as `theta_xy`.
+                If `None` or zero, no `ZZ` rotation is applied.
+            fields_z: Optional per-site Z field. Either a scalar applied to
+                every site, a sequence of length `n_qubits`, or `None`. If
+                `None` or zero on a site, the `Z` rotation for that site is
+                skipped.
+
+        Notes:
+            Gates are applied in the order:
+                exchange(edge), rzz(edge), then rz(site) for each site.
+            All terms generated by the unitary commute with the total `Z`
+            operator, so the ordering is mathematically irrelevant — only
+            the per-edge product matters. The fixed call order is chosen
+            for predictable test output.
+
+        Conservation under truncation:
+            Every gate emitted by this helper commutes with the total Z
+            magnetization :math:`\\sum_k Z_k`, so an observable built from
+            `{I, Z}`-only Pauli strings — :math:`\\sum_i Z_i`,
+            :math:`\\sum_{i<j} Z_i Z_j` (the sum over *all* pairs, not
+            just nearest neighbours), or any polynomial in `Z_i` — is
+            conserved through repeated calls to this helper up to
+            per-gate floating-point precision (~1e-15 per gate, growing
+            linearly with the number of gates). Standard truncation
+            settings preserve this guarantee: pick `min_abs_coeff` well
+            below the conserved coefficient magnitude (usually 1) and
+            well above the per-gate ε floor (≈ 1e-12 for short circuits,
+            looser for long ones). Setting `min_abs_coeff` close to 1, or
+            running enough gates for the accumulated drift to approach
+            the cutoff, can break this guarantee. For more general
+            U(1)-conserving Hermitian observables (spin currents,
+            higher-weight correlators), conjugate-paired Pauli strings
+            carry coefficients of equal magnitude, so symmetric truncation
+            keeps the observable inside the U(1) sub-algebra even when
+            individual coefficient values drift.
+
+        Example:
+            ```python
+            ps = PauliSum.new(4, ["Z0", "Z1", "Z2", "Z3"])
+            ps.apply_u1_trotter_step(
+                edges=[(0, 1), (1, 2), (2, 3)],
+                theta_xy=0.2,
+                theta_zz=0.05,
+                fields_z=[0.1, 0.0, -0.1, 0.0],
+            )
+            ```
+        """
+        n_edges = len(edges)
+        n = self.n_qubits
+
+        def _resolve_per_edge(value):
+            if value is None:
+                return [0.0] * n_edges
+            if isinstance(value, (int, float)):
+                return [float(value)] * n_edges
+            value = list(value)
+            if len(value) != n_edges:
+                raise ValueError(f"Expected {n_edges} entries for per-edge angle, got {len(value)}")
+            return [float(v) for v in value]
+
+        theta_xy_list = _resolve_per_edge(theta_xy)
+        theta_zz_list = _resolve_per_edge(theta_zz)
+
+        if fields_z is None:
+            fields_list: list[float] = [0.0] * n
+        elif isinstance(fields_z, (int, float)):
+            fields_list = [float(fields_z)] * n
+        else:
+            fields_list = [float(v) for v in fields_z]
+            if len(fields_list) != n:
+                raise ValueError(
+                    f"fields_z must have length n_qubits ({n}), got {len(fields_list)}"
+                )
+
+        for k, edge in enumerate(edges):
+            if len(edge) == 2:
+                i, j = edge
+                coupling = 1.0
+            elif len(edge) == 3:
+                i, j, coupling = edge
+                coupling = float(coupling)
+            else:
+                raise ValueError(f"Edge {edge!r} must be a 2- or 3-tuple (i, j[, coupling])")
+
+            theta_xy_eff = theta_xy_list[k] * coupling
+            theta_zz_eff = theta_zz_list[k] * coupling
+            if theta_xy_eff != 0.0 and theta_zz_eff != 0.0:
+                self.xyzz(i, j, theta_xy_eff, theta_zz_eff)
+            elif theta_xy_eff != 0.0:
+                self.exchange(i, j, theta_xy_eff)
+            elif theta_zz_eff != 0.0:
+                self.rzz(i, j, theta_zz_eff)
+
+        for site, h in enumerate(fields_list):
+            if h != 0.0:
+                self.rz(site, h)
+
 
 @dataclass(frozen=True)
 class LossyPauliSum(PauliSum):
diff --git a/ppvm-python/test/test_u1.py b/ppvm-python/test/test_u1.py
new file mode 100644
index 000000000..574cc2bb4
--- /dev/null
+++ b/ppvm-python/test/test_u1.py
@@ -0,0 +1,288 @@
+# SPDX-FileCopyrightText: 2026 The PPVM Authors
+# SPDX-License-Identifier: Apache-2.0
+
+"""Tests for U(1)-symmetric (z-magnetization-conserving) gate helpers."""
+
+import math
+
+import pytest
+
+from ppvm import PauliSum
+
+
+def _terms_dict(state: PauliSum) -> dict[str, float]:
+    return dict(state.terms)
+
+
+def _approx_equal(a: dict[str, float], b: dict[str, float], tol: float = 1e-10) -> bool:
+    keys = set(a) | set(b)
+    return all(abs(a.get(k, 0.0) - b.get(k, 0.0)) < tol for k in keys)
+
+
+def test_exchange_exists_and_runs():
+    """``PauliSum.exchange`` is callable and propagates a simple observable."""
+    ps = PauliSum.new(2, "IZ")
+    ps.exchange(0, 1, 0.3)
+    # IZ is in the (1,1) sector, so the exchange must produce a non-trivial spread.
+    assert len(ps) >= 1
+
+
+def test_xyzz_exists_and_runs():
+    """``PauliSum.xyzz`` is callable and propagates a simple observable."""
+    ps = PauliSum.new(2, "IZ")
+    ps.xyzz(0, 1, theta_xy=0.2, theta_zz=0.05)
+    assert len(ps) >= 1
+
+
+def test_exchange_matches_rxx_then_ryy():
+    """``exchange(θ)`` is equivalent to ``rxx(θ)`` then ``ryy(θ)``."""
+    terms = ["IZ", "ZI", "ZZ", "XY", "YX", "IX", "XI", "YZ"]
+    angles = [0.0, 0.3, -0.7, math.pi / 3]
+    for term in terms:
+        for theta in angles:
+            fused = PauliSum.new(2, term)
+            composed = PauliSum.new(2, term)
+            fused.exchange(0, 1, theta)
+            composed.rxx(0, 1, theta)
+            composed.ryy(0, 1, theta)
+            assert _approx_equal(_terms_dict(fused), _terms_dict(composed)), (
+                f"exchange disagrees on {term} at θ={theta}: "
+                f"{_terms_dict(fused)} vs {_terms_dict(composed)}"
+            )
+
+
+def test_xyzz_matches_exchange_then_rzz():
+    """``xyzz`` is equivalent to ``exchange`` followed by ``rzz``."""
+    terms = ["IZ", "ZI", "ZZ", "XY", "YX", "XX", "YY"]
+    pairs = [(0.0, 0.0), (0.3, 0.1), (-0.4, 0.7), (math.pi / 5, -math.pi / 7)]
+    for term in terms:
+        for theta_xy, theta_zz in pairs:
+            combined = PauliSum.new(2, term)
+            stepwise = PauliSum.new(2, term)
+            combined.xyzz(0, 1, theta_xy=theta_xy, theta_zz=theta_zz)
+            stepwise.exchange(0, 1, theta_xy)
+            stepwise.rzz(0, 1, theta_zz)
+            assert _approx_equal(_terms_dict(combined), _terms_dict(stepwise)), (
+                f"xyzz disagrees on {term} at θ_xy={theta_xy}, θ_zz={theta_zz}: "
+                f"{_terms_dict(combined)} vs {_terms_dict(stepwise)}"
+            )
+
+
+def test_exchange_zero_angle_is_identity():
+    """``exchange(a, b, 0)`` leaves any observable unchanged."""
+    for term in ["IZ", "ZI", "XY", "YX", "II", "ZZ"]:
+        ps = PauliSum.new(2, term)
+        before = _terms_dict(ps)
+        ps.exchange(0, 1, 0.0)
+        assert _approx_equal(_terms_dict(ps), before), (
+            f"exchange(0) altered {term}: got {_terms_dict(ps)}"
+        )
+
+
+def test_exchange_preserves_total_z():
+    """``Σ_i Z_i`` commutes with `XX + YY`; exchange propagation is a no-op."""
+    ps = PauliSum.new(3, [("Z0", 1.0), ("Z1", 1.0), ("Z2", 1.0)])
+    before = _terms_dict(ps)
+    ps.exchange(0, 1, 0.41)
+    ps.exchange(1, 2, -1.2)
+    assert _approx_equal(_terms_dict(ps), before), (
+        f"exchange perturbed total Z: {_terms_dict(ps)} vs {before}"
+    )
+
+
+def test_xyzz_preserves_total_z():
+    """`XX + YY + ZZ` also commutes with `Σ_i Z_i`."""
+    ps = PauliSum.new(3, [("Z0", 1.0), ("Z1", 1.0), ("Z2", 1.0)])
+    before = _terms_dict(ps)
+    ps.xyzz(0, 1, theta_xy=0.5, theta_zz=0.3)
+    ps.xyzz(1, 2, theta_xy=-0.2, theta_zz=0.7)
+    assert _approx_equal(_terms_dict(ps), before)
+
+
+def test_apply_u1_trotter_step_matches_manual_application():
+    """The high-level helper must reproduce the documented manual gate order."""
+    edges = [(0, 1), (1, 2), (2, 3)]
+    theta_xy = 0.21
+    theta_zz = 0.07
+    fields_z = [0.1, 0.0, -0.1, 0.05]
+
+    auto = PauliSum.new(4, [("Z0", 1.0), ("Z1", -1.0), ("Z2", 0.5)])
+    manual = PauliSum.new(4, [("Z0", 1.0), ("Z1", -1.0), ("Z2", 0.5)])
+
+    auto.apply_u1_trotter_step(
+        edges=edges,
+        theta_xy=theta_xy,
+        theta_zz=theta_zz,
+        fields_z=fields_z,
+    )
+    for i, j in edges:
+        manual.xyzz(i, j, theta_xy=theta_xy, theta_zz=theta_zz)
+    for site, h in enumerate(fields_z):
+        if h != 0.0:
+            manual.rz(site, h)
+
+    assert _approx_equal(_terms_dict(auto), _terms_dict(manual)), (
+        f"u1 trotter helper disagrees with manual application:\n"
+        f"auto={_terms_dict(auto)}\nmanual={_terms_dict(manual)}"
+    )
+
+
+def test_apply_u1_trotter_step_with_per_edge_couplings():
+    """Per-edge `(i, j, J)` triples must scale each edge's angles by `J`."""
+    base_xy = 0.3
+    couplings = [1.0, 0.5, -2.0]
+    edges = [(0, 1, couplings[0]), (1, 2, couplings[1]), (2, 3, couplings[2])]
+
+    auto = PauliSum.new(4, "Z0")
+    manual = PauliSum.new(4, "Z0")
+
+    auto.apply_u1_trotter_step(edges=edges, theta_xy=base_xy)
+    for i, j, J in edges:
+        manual.exchange(i, j, base_xy * J)
+
+    assert _approx_equal(_terms_dict(auto), _terms_dict(manual))
+
+
+def test_apply_u1_trotter_step_per_edge_angles():
+    """`theta_xy` may be a list of per-edge angles."""
+    edges = [(0, 1), (1, 2)]
+    angles = [0.4, -0.25]
+
+    auto = PauliSum.new(3, "Z0")
+    manual = PauliSum.new(3, "Z0")
+
+    auto.apply_u1_trotter_step(edges=edges, theta_xy=angles)
+    for k, (i, j) in enumerate(edges):
+        manual.exchange(i, j, angles[k])
+
+    assert _approx_equal(_terms_dict(auto), _terms_dict(manual))
+
+
+def test_apply_u1_trotter_step_preserves_total_z():
+    """A full Trotter slice still preserves any observable in the total-Z sector."""
+    ps = PauliSum.new(4, [("Z0", 1.0), ("Z1", 1.0), ("Z2", 1.0), ("Z3", 1.0)])
+    before = _terms_dict(ps)
+    ps.apply_u1_trotter_step(
+        edges=[(0, 1), (1, 2), (2, 3)],
+        theta_xy=0.3,
+        theta_zz=0.15,
+        fields_z=[0.1, -0.2, 0.05, 0.0],
+    )
+    # `\\sum Z_k` is invariant under the full Hamiltonian, including the
+    # diagonal Z field (which is itself in the total-Z algebra).
+    assert _approx_equal(_terms_dict(ps), before), (
+        f"trotter step perturbed total Z: {_terms_dict(ps)} vs {before}"
+    )
+
+
+def test_apply_u1_trotter_step_rejects_bad_lengths():
+    """The helper validates per-edge / per-site list lengths."""
+    ps = PauliSum.new(3, "Z0")
+    with pytest.raises(ValueError, match="per-edge"):
+        ps.apply_u1_trotter_step(edges=[(0, 1), (1, 2)], theta_xy=[0.1])
+    with pytest.raises(ValueError, match="fields_z"):
+        ps.apply_u1_trotter_step(edges=[(0, 1)], theta_xy=0.1, fields_z=[0.1, 0.0])
+
+
+def test_xy_dynamics_against_xxz_evolution():
+    """Two consecutive `xyzz` calls on the same edge with `θ_xy = 0`, `θ_zz ≠ 0`
+    must commute with everything in the Z-diagonal sector — a direct check
+    that the ZZ piece of `xyzz` matches a pure rzz."""
+    ps_xyzz = PauliSum.new(2, "ZZ")
+    ps_rzz = PauliSum.new(2, "ZZ")
+    ps_xyzz.xyzz(0, 1, theta_xy=0.0, theta_zz=0.42)
+    ps_rzz.rzz(0, 1, 0.42)
+    assert _approx_equal(_terms_dict(ps_xyzz), _terms_dict(ps_rzz))
+
+
+# =============================================================================
+# Truncation-robust conservation hardening
+# =============================================================================
+#
+# `apply_u1_trotter_step` propagates {I, Z}-polynomial observables exactly
+# in the conserved sector, modulo per-gate floating-point ε. These tests
+# pin that guarantee against aggressive truncation policies: as long as
+# the cutoff (`min_abs_coeff` or `max_pauli_weight`) leaves headroom above
+# the per-gate ε and below the conserved-coefficient magnitude (~1),
+# conservation is preserved within the `1e-10` tolerance used by
+# `_approx_equal`.
+
+
+def test_apply_u1_trotter_step_total_z_under_aggressive_coefficient_truncation():
+    """`Σ Z_i` survives many Trotter sweeps with `min_abs_coeff = 0.5`.
+
+    The conserved coefficients are O(1) — well above the cutoff. The
+    transient cross terms cancel to O(ε) ≪ 0.5 before truncation runs,
+    so the cutoff drops only the ε residues and leaves `Σ Z_i` intact.
+    """
+    n = 5
+    edges = [(i, i + 1) for i in range(n - 1)]
+    ps = PauliSum.new(
+        n,
+        [(f"Z{i}", 1.0) for i in range(n)],
+        min_abs_coeff=0.5,
+    )
+    expected = _terms_dict(ps)
+    for _ in range(4):
+        ps.apply_u1_trotter_step(
+            edges=edges,
+            theta_xy=0.3,
+            theta_zz=0.1,
+            fields_z=[0.05] * n,
+        )
+    assert _approx_equal(_terms_dict(ps), expected), (
+        f"total Z drifted under aggressive coefficient truncation:\n"
+        f"got={_terms_dict(ps)}\nwant={expected}"
+    )
+
+
+def test_apply_u1_trotter_step_total_z_under_max_pauli_weight_one():
+    """`Σ Z_i` survives Trotter sweeps with `max_pauli_weight = 1`.
+
+    Cross terms have weight 2 and are dropped by the discrete weight
+    check, independent of any floating-point sensitivity in the cutoff.
+    """
+    n = 5
+    edges = [(i, i + 1) for i in range(n - 1)]
+    ps = PauliSum.new(
+        n,
+        [(f"Z{i}", 1.0) for i in range(n)],
+        max_pauli_weight=1,
+    )
+    expected = _terms_dict(ps)
+    for _ in range(3):
+        ps.apply_u1_trotter_step(
+            edges=edges,
+            theta_xy=0.25,
+            theta_zz=0.07,
+            fields_z=[0.04] * n,
+        )
+    assert _approx_equal(_terms_dict(ps), expected), (
+        f"total Z drifted under MaxPauliWeight=1 truncation:\n"
+        f"got={_terms_dict(ps)}\nwant={expected}"
+    )
+
+
+def test_full_zz_correlator_under_aggressive_truncation():
+    """`Σ_{i<j} Z_iZ_j` (the sum over *all* pairs, not just nearest
+    neighbours) commutes with `XX+YY+ZZ` on every edge, so the full
+    correlator survives Trotter sweeps under aggressive truncation. The
+    nearest-neighbour-only sum is not conserved, so the full pair
+    decomposition is what locks down `xyzz`'s ZZ semantics."""
+    n = 4
+    # All C(n, 2) = 6 unordered pairs: (0,1),(0,2),(0,3),(1,2),(1,3),(2,3).
+    all_pairs = [(i, j) for i in range(n) for j in range(i + 1, n)]
+    terms = [(f"Z{i}Z{j}", 1.0) for i, j in all_pairs]
+    ps = PauliSum.new(n, terms, min_abs_coeff=0.5, max_pauli_weight=2)
+    expected = _terms_dict(ps)
+    edges = [(i, i + 1) for i in range(n - 1)]
+    for _ in range(2):
+        ps.apply_u1_trotter_step(
+            edges=edges,
+            theta_xy=0.21,
+            theta_zz=0.07,
+        )
+    assert _approx_equal(_terms_dict(ps), expected), (
+        f"full ZZ correlator drifted under aggressive truncation:\n"
+        f"got={_terms_dict(ps)}\nwant={expected}"
+    )