diff --git a/benches/third_party/bench_common.jl b/benches/third_party/bench_common.jl
new file mode 100644
index 00000000..7dc6aac8
--- /dev/null
+++ b/benches/third_party/bench_common.jl
@@ -0,0 +1,60 @@
+# Copyright 2026 Algorithmiq
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+# Shared across the third-party Julia benchmark scripts in this project (pauli_prop,
+# majorana_prop). Included via `include(joinpath(@__DIR__, "..", "bench_common.jl"))`.
+
+if Threads.nthreads() == 1
+    @warn "Julia is running with only 1 thread: the peak-memory sampler runs as a background " *
+          "task, but a single-threaded Julia can only switch to it between steps, not while the " *
+          "propagation call itself is running, so the memory figures will undercount any " *
+          "transient spike freed before that call returns. Re-run with `julia --threads=auto` " *
+          "(or set JULIA_NUM_THREADS) for accurate peak-memory measurements."
+end
+
+# Linux-only: reads the kernel's live RSS for this process directly, in kB.
+function current_rss_bytes()::Int
+    for line in eachline("/proc/self/status")
+        if startswith(line, "VmRSS:")
+            return parse(Int, split(line)[2]) * 1024
+        end
+    end
+    return 0
+end
+
+# Tracks this process's peak RSS within a resettable window: the kernel's own high-water mark
+# (Sys.maxrss()) is monotonic for the whole process and can't be reset per step, so this instead
+# polls the *current* RSS from a background task and keeps the max seen since the last reset!.
+mutable struct RssPeakSampler
+    peak_bytes::Threads.Atomic{Int}
+    running::Threads.Atomic{Bool}
+end
+
+function start_sampler(interval_s::Float64=1e-3)
+    sampler = RssPeakSampler(Threads.Atomic{Int}(current_rss_bytes()), Threads.Atomic{Bool}(true))
+    Threads.@spawn begin
+        while sampler.running[]
+            rss = current_rss_bytes()
+            if rss > sampler.peak_bytes[]
+                sampler.peak_bytes[] = rss
+            end
+            sleep(interval_s)
+        end
+    end
+    return sampler
+end
+
+reset!(sampler::RssPeakSampler) = (sampler.peak_bytes[] = current_rss_bytes())
+peak_mb(sampler::RssPeakSampler) = sampler.peak_bytes[] / 1024^2
+stop!(sampler::RssPeakSampler) = (sampler.running[] = false)
diff --git a/benches/third_party/bench_common.py b/benches/third_party/bench_common.py
new file mode 100644
index 00000000..29098d45
--- /dev/null
+++ b/benches/third_party/bench_common.py
@@ -0,0 +1,65 @@
+# Copyright 2026 Algorithmiq
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Shared across the third-party benchmark scripts in this project (pauli_prop, majorana_prop)."""
+
+from __future__ import annotations
+
+import threading
+
+import psutil
+
+
+class RssPeakSampler:
+    """Tracks this process's peak resident set size within a resettable window.
+
+    The kernel-tracked high-water mark (``resource.getrusage().ru_maxrss``) is monotonic for the
+    whole process lifetime and can never be reset, so it cannot isolate a single step's peak from
+    the steps around it. This instead polls the *current* RSS from a background thread at a high
+    fixed frequency and keeps the max seen since the last :meth:`reset`, so each step gets its own
+    peak, independent of what earlier or later steps did.
+
+    Caveat: this can only sample while the calling thread doesn't hold the GIL — a C extension
+    call that never releases it (e.g. monoprop's ``propagate()``) blocks this thread out for its
+    whole duration, so the sampler only catches whatever RSS growth is already visible by the time
+    control returns to Python.
+    """
+
+    def __init__(self, interval_s: float = 1e-3) -> None:
+        self._process = psutil.Process()
+        self._interval_s = interval_s
+        self._peak_bytes = 0
+        self._stop_event = threading.Event()
+        self._thread = threading.Thread(target=self._poll_loop, daemon=True)
+
+    def _poll_loop(self) -> None:
+        while not self._stop_event.wait(self._interval_s):
+            rss = self._process.memory_info().rss
+            if rss > self._peak_bytes:
+                self._peak_bytes = rss
+
+    def __enter__(self) -> RssPeakSampler:
+        self._thread.start()
+        return self
+
+    def __exit__(self, *exc_info: object) -> None:
+        self._stop_event.set()
+        self._thread.join()
+
+    def reset(self) -> None:
+        """Start a new measurement window, floored at the current RSS."""
+        self._peak_bytes = self._process.memory_info().rss
+
+    def peak_mb(self) -> float:
+        return self._peak_bytes / 1024**2
diff --git a/benches/third_party/majorana_prop/julia_hubbard1d_benchmark.jl b/benches/third_party/majorana_prop/julia_hubbard1d_benchmark.jl
index 0b8bc649..811f590a 100644
--- a/benches/third_party/majorana_prop/julia_hubbard1d_benchmark.jl
+++ b/benches/third_party/majorana_prop/julia_hubbard1d_benchmark.jl
@@ -17,6 +17,8 @@ using BenchmarkTools
 using ArgParse
 using JSON
 
+include(joinpath(@__DIR__, "..", "bench_common.jl"))
+
 
 function experiment(N_spinful_sites, fock_state, circ_single, thetas_single, n_layers)
     site_index = N_spinful_sites ÷ 2
@@ -26,33 +28,64 @@ function experiment(N_spinful_sites, fock_state, circ_single, thetas_single, n_l
 
     obs = VectorMajoranaSum(MajoranaSum(N_spinful_sites, :nup, site_index))
 
-    res = zeros(n_layers + 1)
-    res[1] = overlapwithfock(obs, fock_state)
+    values = zeros(n_layers + 1)
+    term_counts = zeros(Int, n_layers + 1)
+    cumulative_runtimes = zeros(n_layers + 1)
+    memory_size = zeros(n_layers + 1)
+    native_memory_size = zeros(n_layers + 1)
+
+    sampler = start_sampler()
+
+    reset!(sampler)
+    cumulative_runtimes[1] = @elapsed (values[1] = overlapwithfock(obs, fock_state))
+    term_counts[1] = length(obs)
+    memory_size[1] = peak_mb(sampler)
+    native_memory_size[1] = Base.summarysize(obs) / 1024^2
+    for k = 1:n_layers
+        reset!(sampler)
+        step_runtime = @elapsed propagate!(circ_single, obs, thetas_single, min_abs_coeff=min_abs_coeff, max_unpaired=max_unpaired)
+        values[k+1] = overlapwithfock(obs, fock_state)
+        term_counts[k+1] = length(obs)
+        cumulative_runtimes[k+1] = cumulative_runtimes[k] + step_runtime
+        memory_size[k+1] = peak_mb(sampler)
+        native_memory_size[k+1] = Base.summarysize(obs) / 1024^2
+    end
 
+    stop!(sampler)
 
-    loop_elapsed = @elapsed for k = 1:n_layers
-        propagate!(circ_single, obs, thetas_single, min_abs_coeff=min_abs_coeff, max_unpaired=max_unpaired)
-        res[k+1] = overlapwithfock(obs, fock_state)
-    end
-    memory_size = Base.summarysize(obs) / 1024^2
-    return res, length(obs), loop_elapsed, memory_size
+    return values, term_counts, cumulative_runtimes, memory_size, native_memory_size
 
 
 end
-function save_result(output_path, N_spinful_sites, n_layers, obs_length, final_res, loop_elapsed, memory_size)
-    """Append one benchmark result as a JSON line, creating the parent directory if needed."""
-    record = Dict(
-        "n_spinful_sites" => N_spinful_sites,
-        "n_layers" => n_layers,
-        "num_terms" => obs_length,
-        "final_overlap" => final_res,
-        "runtime_seconds" => loop_elapsed,
-        "memory_MB" => memory_size,
-    )
-
-    open(output_path, "a") do io
-        JSON.print(io, record)
-        println(io)
+function save_result(output_path, source, N_spinful_sites, n_layers, term_counts, values, cumulative_runtimes, memory_size, native_memory_size, num_threads)
+    """Merge this run's per-step data into the shared results JSON file, keyed by source label."""
+    data = if isfile(output_path)
+        JSON.parsefile(output_path)
+    else
+        Dict(
+            "n_spinful_sites" => N_spinful_sites,
+            "n_layers" => n_layers,
+            "step_range" => collect(0:n_layers),
+            "num_threads" => Dict(),
+            "runtime_seconds" => Dict(),
+            "expectation_value" => Dict(),
+            "num_terms" => Dict(),
+            "memory_MB" => Dict(),
+            "native_memory_MB" => Dict(),
+        )
+    end
+    data["num_threads"] = get(data, "num_threads", Dict())
+    data["num_threads"][source] = num_threads
+    data["runtime_seconds"][source] = cumulative_runtimes
+    data["expectation_value"][source] = values
+    data["num_terms"][source] = term_counts
+    data["memory_MB"][source] = memory_size
+    data["native_memory_MB"] = get(data, "native_memory_MB", Dict())
+    data["native_memory_MB"][source] = native_memory_size
+
+    mkpath(dirname(output_path))
+    open(output_path, "w") do io
+        JSON.print(io, data, 4)
     end
 end
 
@@ -61,28 +94,27 @@ function main(args)
     s = ArgParseSettings(description="Arguments for the 1D Hubbard model benchmark.")
 
     @add_arg_table! s begin
-        "--case", "-c"
-        help = "Case pair to run."
+        "--n-spins", "-n"
+        help = "Number of spinful sites."
         arg_type = Int
-        default = 1
+        default = 60
+        dest_name = "n_spins"
+        "--max-layers", "-l"
+        help = "Number of Trotter layers."
+        arg_type = Int
+        default = 20
+        dest_name = "max_layers"
         "--output", "-o"
-        help = "Path to the JSONL file results are appended to."
+        help = "Path to the shared JSON file results are merged into."
         arg_type = String
-        default = joinpath(@__DIR__, "julia_hubbard1d_benchmark_results.jsonl")
+        default = joinpath(@__DIR__, "results.json")
 
     end
 
     parsed_args = parse_args(s)
 
-    spin_layers_pairs = []
-    for i in [20, 40, 60]
-        for j in 10:2:18
-            push!(spin_layers_pairs, (i, j))
-        end
-    end
-
-    case_pair = parsed_args["case"]
-    N_spinful_sites, n_layers = spin_layers_pairs[case_pair]
+    N_spinful_sites = parsed_args["n_spins"]
+    n_layers = parsed_args["max_layers"]
 
     t = 1.
     U = 1.5
@@ -117,11 +149,10 @@ function main(args)
 
     println("Number of threads: $(Threads.nthreads())")
 
-    res, obs_length, loop_elapsed, memory_size = experiment(N_spinful_sites, fock_state, circ_single, thetas_single, n_layers)
-    final_res = res[end]
-    println("$N_spinful_sites n_spin $n_layers layers $obs_length num_terms $final_res final overlap $loop_elapsed seconds")
+    values, term_counts, cumulative_runtimes, memory_size, native_memory_size = experiment(N_spinful_sites, fock_state, circ_single, thetas_single, n_layers)
+    println("$N_spinful_sites n_spin $n_layers layers $(term_counts[end]) num_terms $(values[end]) final overlap $(cumulative_runtimes[end]) seconds")
 
-    save_result(parsed_args["output"], N_spinful_sites, n_layers, obs_length, final_res, loop_elapsed, memory_size)
+    save_result(parsed_args["output"], "MajoranaPropagation.jl", N_spinful_sites, n_layers, term_counts, values, cumulative_runtimes, memory_size, native_memory_size, Threads.nthreads())
 
 end
 
diff --git a/benches/third_party/majorana_prop/julia_hubbard1d_benchmark_results.jsonl b/benches/third_party/majorana_prop/julia_hubbard1d_benchmark_results.jsonl
deleted file mode 100644
index c448107a..00000000
--- a/benches/third_party/majorana_prop/julia_hubbard1d_benchmark_results.jsonl
+++ /dev/null
@@ -1,15 +0,0 @@
-{"final_overlap":0.5540635634956823,"runtime_seconds":4.489840569,"n_spinful_sites":20,"memory_MB":31.664398193359375,"n_layers":10,"num_terms":597051}
-{"final_overlap":0.6265984784190225,"runtime_seconds":11.580550668,"n_spinful_sites":20,"memory_MB":98.66783905029297,"n_layers":12,"num_terms":1754896}
-{"final_overlap":0.5540635642561327,"runtime_seconds":8.734066081,"n_spinful_sites":40,"memory_MB":52.7220458984375,"n_layers":10,"num_terms":597311}
-{"final_overlap":0.5540635608915114,"runtime_seconds":11.348441772,"n_spinful_sites":60,"memory_MB":52.974884033203125,"n_layers":10,"num_terms":597601}
-{"final_overlap":0.6265984814098748,"runtime_seconds":34.35385869,"n_spinful_sites":60,"memory_MB":164.24754333496094,"n_layers":12,"num_terms":1754476}
-{"final_overlap":0.6265984871125446,"runtime_seconds":31.924598718,"n_spinful_sites":40,"memory_MB":164.2796630859375,"n_layers":12,"num_terms":1754526}
-{"final_overlap":0.646851959851176,"runtime_seconds":38.978723019,"n_spinful_sites":20,"memory_MB":277.8764343261719,"n_layers":14,"num_terms":4870330}
-{"final_overlap":0.6468519627999069,"runtime_seconds":78.097265415,"n_spinful_sites":40,"memory_MB":462.9164047241211,"n_layers":14,"num_terms":4867455}
-{"final_overlap":0.6468519475075789,"runtime_seconds":112.541124557,"n_spinful_sites":60,"memory_MB":462.87574005126953,"n_layers":14,"num_terms":4865089}
-{"final_overlap":0.6237908253934835,"runtime_seconds":117.691216266,"n_spinful_sites":20,"memory_MB":734.4124298095703,"n_layers":16,"num_terms":12875636}
-{"final_overlap":0.6237907501255578,"runtime_seconds":292.79322942,"n_spinful_sites":60,"memory_MB":1222.6302642822266,"n_layers":16,"num_terms":12859755}
-{"final_overlap":0.6237907699662404,"runtime_seconds":303.675066342,"n_spinful_sites":40,"memory_MB":1224.568588256836,"n_layers":16,"num_terms":12869611}
-{"final_overlap":0.5748922977223003,"runtime_seconds":348.358432812,"n_spinful_sites":20,"memory_MB":1899.7230987548828,"n_layers":18,"num_terms":32668247}
-{"final_overlap":0.5748922102989575,"runtime_seconds":602.853697251,"n_spinful_sites":40,"memory_MB":3167.8807373046875,"n_layers":18,"num_terms":32655012}
-{"final_overlap":0.5748921471542362,"runtime_seconds":981.613709723,"n_spinful_sites":60,"memory_MB":3162.911178588867,"n_layers":18,"num_terms":32625568}
diff --git a/benches/third_party/majorana_prop/monoprop_hubbard1d_benchmark.py b/benches/third_party/majorana_prop/monoprop_hubbard1d_benchmark.py
index eb8a290a..603a157e 100644
--- a/benches/third_party/majorana_prop/monoprop_hubbard1d_benchmark.py
+++ b/benches/third_party/majorana_prop/monoprop_hubbard1d_benchmark.py
@@ -17,6 +17,7 @@
 import argparse
 import json
 import os
+import sys
 from pathlib import Path
 from time import perf_counter
 
@@ -24,7 +25,9 @@
 from monoprop import Circuit, ExpGate, MajoranaPropagator
 from monoprop.fermi import FermiOperator
 
-os.environ["YAQS_LOG_LEVEL"] = "INFO"
+
+sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
+from bench_common import RssPeakSampler  # noqa: E402
 
 
 def mode(site, spin):
@@ -98,33 +101,67 @@ def number_operator_majorana(site, spin, num_qubits):
     )
 
 
-def save_result(output_path, record):
-    """Append one benchmark result as a JSON line, creating the parent directory if needed."""
+SOURCE_LABEL = "monoprop"
+
+
+def save_result(
+    output_path,
+    n_spinful_sites,
+    n_layers,
+    values,
+    term_counts,
+    cumulative_runtimes,
+    memory_size,
+    native_memory_size,
+):
+    """Merge this run's per-step data into the shared results JSON file, keyed by source label."""
     output_path = Path(output_path)
     output_path.parent.mkdir(parents=True, exist_ok=True)
-    with output_path.open("a") as f:
-        f.write(json.dumps(record) + "\n")
+    if output_path.exists():
+        with output_path.open() as f:
+            data = json.load(f)
+    else:
+        data = {
+            "n_spinful_sites": n_spinful_sites,
+            "n_layers": n_layers,
+            "step_range": list(range(n_layers + 1)),
+            "num_threads": {},
+            "runtime_seconds": {},
+            "expectation_value": {},
+            "num_terms": {},
+            "memory_MB": {},
+            "native_memory_MB": {},
+        }
+    data.setdefault("num_threads", {})[SOURCE_LABEL] = os.environ.get(
+        "monoprop_NUM_THREADS", "not set"
+    )
+    data["runtime_seconds"][SOURCE_LABEL] = cumulative_runtimes
+    data["expectation_value"][SOURCE_LABEL] = values
+    data["num_terms"][SOURCE_LABEL] = term_counts
+    data["memory_MB"][SOURCE_LABEL] = memory_size
+    data.setdefault("native_memory_MB", {})[SOURCE_LABEL] = native_memory_size
+    with output_path.open("w") as f:
+        json.dump(data, f, indent=4)
 
 
 def main():
     parser = argparse.ArgumentParser(description="Benchmark for 1D Hubbard model")
-    parser.add_argument("--case", "-c", help="Case pair to run.", type=int, default=0)
+    parser.add_argument(
+        "--n-spins", "-n", help="Number of spinful sites.", type=int, default=60
+    )
+    parser.add_argument(
+        "--max-layers", "-l", help="Number of Trotter layers.", type=int, default=20
+    )
     parser.add_argument(
         "--output",
         "-o",
-        help="Path to the JSONL file results are appended to.",
-        default=Path(__file__).with_name("monoprop_hubbard1d_benchmark_results.jsonl"),
+        help="Path to the shared JSON file results are merged into.",
+        default=Path(__file__).with_name("results.json"),
     )
 
     args = parser.parse_args()
 
-    spin_layer_cases = []
-    for i in [20, 40, 60]:
-        for j in range(10, 19, 2):
-            spin_layer_cases.append((i, j))
-
-    case_pair = args.case
-    n_spinful_sites, n_layers = spin_layer_cases[case_pair]
+    n_spinful_sites, n_layers = args.n_spins, args.max_layers
     trotter_steps = n_layers
     t = 1.0
     u = 1.5
@@ -159,31 +196,42 @@ def main():
 
     values = np.empty(trotter_steps + 1)
     term_counts = np.empty(trotter_steps + 1, dtype=int)
-
-    values[0] = simulator.expectation_value()
-    term_counts[0] = simulator.size()
-
-    t_start = perf_counter()
-    for step in range(trotter_steps):
-        simulator.propagate(fermi_circuit)
-        values[step + 1] = simulator.expectation_value()
-        term_counts[step + 1] = simulator.size()
-    t_total = perf_counter() - t_start
-    memory_size = simulator._simulator.operator_memory_bytes() / 1024**2
+    cumulative_runtimes = np.empty(trotter_steps + 1)
+    memory_size = np.empty(trotter_steps + 1)
+    native_memory_size = np.empty(trotter_steps + 1)
+
+    with RssPeakSampler() as sampler:
+        sampler.reset()
+        t_start = perf_counter()
+        values[0] = simulator.expectation_value()
+        cumulative_runtimes[0] = perf_counter() - t_start
+        term_counts[0] = simulator.size()
+        memory_size[0] = sampler.peak_mb()
+        native_memory_size[0] = simulator._simulator.operator_memory_bytes() / 1024**2
+        for step in range(trotter_steps):
+            sampler.reset()
+            step_start = perf_counter()
+            simulator.propagate(fermi_circuit)
+            step_runtime = perf_counter() - step_start
+            values[step + 1] = simulator.expectation_value()
+            term_counts[step + 1] = simulator.size()
+            cumulative_runtimes[step + 1] = cumulative_runtimes[step] + step_runtime
+            memory_size[step + 1] = sampler.peak_mb()
+            native_memory_size[step + 1] = (
+                simulator._simulator.operator_memory_bytes() / 1024**2
+            )
     print(
-        f"{n_spinful_sites} n_spin {n_layers} layers {term_counts[-1]} num_terms {values[-1]} final overlap  runtime {t_total:.3f} seconds"
+        f"{n_spinful_sites} n_spin {n_layers} layers {term_counts[-1]} num_terms {values[-1]} final overlap  runtime {cumulative_runtimes[-1]:.3f} seconds"
     )
     save_result(
         args.output,
-        {
-            "n_spinful_sites": n_spinful_sites,
-            "n_layers": n_layers,
-            "num_threads": os.environ.get("monoprop_NUM_THREADS", "not set"),
-            "runtime_seconds": t_total,
-            "final_overlap": values[-1],
-            "num_terms": term_counts[-1],
-            "memory_MB": memory_size,
-        },
+        n_spinful_sites,
+        n_layers,
+        values.tolist(),
+        term_counts.tolist(),
+        cumulative_runtimes.tolist(),
+        memory_size.tolist(),
+        native_memory_size.tolist(),
     )
 
 
diff --git a/benches/third_party/majorana_prop/monoprop_hubbard1d_benchmark_results.jsonl b/benches/third_party/majorana_prop/monoprop_hubbard1d_benchmark_results.jsonl
deleted file mode 100644
index b751e379..00000000
--- a/benches/third_party/majorana_prop/monoprop_hubbard1d_benchmark_results.jsonl
+++ /dev/null
@@ -1,15 +0,0 @@
-{"n_spinful_sites": 20, "n_layers": 10, "num_threads": "8", "runtime_seconds": 1.430209699086845, "memory_MB": 79.20969867706299, "final_overlap": 0.5540633497210324, "num_terms": 883516}
-{"n_spinful_sites": 60, "n_layers": 10, "num_threads": "8", "runtime_seconds": 3.816878085024655, "memory_MB": 125.21032428741455, "final_overlap": 0.5540633497210324, "num_terms": 883516}
-{"n_spinful_sites": 40, "n_layers": 10, "num_threads": "8", "runtime_seconds": 2.7269934061914682, "memory_MB": 102.21001148223877, "final_overlap": 0.5540633497210324, "num_terms": 883516}
-{"n_spinful_sites": 20, "n_layers": 12, "num_threads": "8", "runtime_seconds": 5.15137222991325, "memory_MB": 215.9552125930786, "final_overlap": 0.6265986462342832, "num_terms": 2649785}
-{"n_spinful_sites": 60, "n_layers": 12, "num_threads": "8", "runtime_seconds": 13.834477874450386, "memory_MB": 339.9558382034302, "final_overlap": 0.6265986462342832, "num_terms": 2649785}
-{"n_spinful_sites": 20, "n_layers": 14, "num_threads": "8", "runtime_seconds": 14.984117835760117, "memory_MB": 640.432110786438, "final_overlap": 0.646851559598512, "num_terms": 7425174}
-{"n_spinful_sites": 40, "n_layers": 12, "num_threads": "8", "runtime_seconds": 7.933134694583714, "memory_MB": 277.9555253982544, "final_overlap": 0.6265986462342832, "num_terms": 2649785}
-{"n_spinful_sites": 40, "n_layers": 14, "num_threads": "8", "runtime_seconds": 25.00145892612636, "memory_MB": 824.4324235916138, "final_overlap": 0.646851559598512, "num_terms": 7425178}
-{"n_spinful_sites": 60, "n_layers": 14, "num_threads": "8", "runtime_seconds": 48.1063445257023, "memory_MB": 1008.4327363967896, "final_overlap": 0.646851559598512, "num_terms": 7425178}
-{"n_spinful_sites": 20, "n_layers": 16, "num_threads": "8", "runtime_seconds": 56.15366280730814, "memory_MB": 1722.5340089797974, "final_overlap": 0.6237897634273656, "num_terms": 19568323}
-{"n_spinful_sites": 40, "n_layers": 16, "num_threads": "8", "runtime_seconds": 81.4321228240151, "memory_MB": 2218.5346269607544, "final_overlap": 0.6237897634273649, "num_terms": 19568345}
-{"n_spinful_sites": 60, "n_layers": 16, "num_threads": "8", "runtime_seconds": 120.94276295881718, "memory_MB": 2714.53493976593, "final_overlap": 0.6237897634273649, "num_terms": 19568345}
-{"n_spinful_sites": 20, "n_layers": 18, "num_threads": "8", "runtime_seconds": 171.3473529824987, "memory_MB": 3545.3169374465942, "final_overlap": 0.5748900593909845, "num_terms": 48479790}
-{"n_spinful_sites": 40, "n_layers": 18, "num_threads": "8", "runtime_seconds": 293.89168079406954, "memory_MB": 4537.319264411926, "final_overlap": 0.5748900593909397, "num_terms": 48480054}
-{"n_spinful_sites": 60, "n_layers": 18, "num_threads": "8", "runtime_seconds": 303.75193184800446, "memory_MB": 5529.319577217102, "final_overlap": 0.5748900593909397, "num_terms": 48480054}
diff --git a/benches/third_party/majorana_prop/plot_results.py b/benches/third_party/majorana_prop/plot_results.py
index 835b1b79..1a63d1b5 100644
--- a/benches/third_party/majorana_prop/plot_results.py
+++ b/benches/third_party/majorana_prop/plot_results.py
@@ -15,66 +15,55 @@
 from __future__ import annotations
 
 import argparse
+import json
 from pathlib import Path
 
 import matplotlib.pyplot as plt
-import pandas as pd
-
-
-def load_benchmark(path: Path, source: str) -> pd.DataFrame:
-    """Load a benchmark JSONL file into a DataFrame of runtime/term-count rows."""
-    df = pd.read_json(path, lines=True)
-    df = df.rename(
-        columns={
-            "n_spinful_sites": "n_spin",
-            "n_layers": "layers",
-            "runtime_seconds": "seconds",
-            "memory_MB": "memory",
-            "final_overlap": "overlap",
-        }
-    )
-    df["source"] = source
-    return df[
-        ["n_spin", "layers", "num_terms", "seconds", "memory", "overlap", "source"]
-    ].sort_values(["n_spin", "layers"])
-
-
-def plot_metric(ax, data: pd.DataFrame, metric: str, ylabel: str) -> None:
-    """Plot ``metric`` vs. layers for each n_spin/source combination onto ``ax``."""
-    styles = {"monoprop": "-o", "MajoranaPropagation.jl": "--x"}
-    colors = plt.cm.tab10.colors
-    for i, n_spin in enumerate(sorted(data["n_spin"].unique())):
-        color = colors[i % len(colors)]
-        for source, style in styles.items():
-            subset = data[(data["n_spin"] == n_spin) & (data["source"] == source)]
-            if subset.empty:
-                continue
+
+STYLES = {"monoprop": "-o", "MajoranaPropagation.jl": "--x"}
+
+
+def plot_metric(
+    ax,
+    step_range: list[int],
+    metric_dict: dict[str, list[float]],
+    ylabel: str,
+    secondary_dict: dict[str, list[float]] | None = None,
+) -> None:
+    """Plot ``metric_dict[source]`` vs. ``step_range`` for each source onto ``ax``.
+
+    ``secondary_dict``, where given, is drawn as a faint unlabeled line reusing each source's own
+    color and linestyle (only reduced alpha, no marker, distinguishes it) — a different
+    measurement of the same source, not a new series.
+    """
+    colors_by_source = {}
+    for source, values in metric_dict.items():
+        (line,) = ax.plot(step_range, values, STYLES.get(source, "-o"), label=source)
+        colors_by_source[source] = line.get_color()
+    if secondary_dict:
+        for source, values in secondary_dict.items():
+            linestyle = "--" if STYLES.get(source, "-o").startswith("--") else "-"
             ax.plot(
-                subset["layers"],
-                subset[metric],
-                style,
-                color=color,
-                label=f"n={n_spin} ({source})",
+                step_range,
+                values,
+                linestyle=linestyle,
+                color=colors_by_source.get(source),
+                alpha=0.4,
+                linewidth=1,
             )
     ax.set_xlabel("layers")
     ax.set_ylabel(ylabel)
-    ax.legend(fontsize="small", ncol=2)
+    ax.legend(fontsize="small")
     ax.grid(True, alpha=0.3)
 
 
 def main() -> None:
     parser = argparse.ArgumentParser(description=__doc__)
     parser.add_argument(
-        "--monoprop-results",
-        type=Path,
-        default=Path("monoprop_hubbard1d_benchmark_results.jsonl"),
-        help="Path to the monoprop benchmark results JSONL file.",
-    )
-    parser.add_argument(
-        "--julia-results",
+        "--results",
         type=Path,
-        default=Path("julia_hubbard1d_benchmark_results.jsonl"),
-        help="Path to the julia benchmark results JSONL file.",
+        default=Path(__file__).with_name("results.json"),
+        help="Path to the shared benchmark results JSON file.",
     )
     parser.add_argument(
         "--output-dir",
@@ -89,27 +78,43 @@ def main() -> None:
     )
     args = parser.parse_args()
 
-    monoprop_df = load_benchmark(args.monoprop_results, "monoprop")
-    julia_df = load_benchmark(args.julia_results, "MajoranaPropagation.jl")
-    df = pd.concat([monoprop_df, julia_df], ignore_index=True)
+    with args.results.open() as file:
+        data = json.load(file)
 
+    step_range = data["step_range"]
     args.output_dir.mkdir(parents=True, exist_ok=True)
 
     fig, axes = plt.subplots(2, 2, figsize=(14, 10))
-    plot_metric(axes[0, 0], df, "seconds", "time (seconds)")
-    axes[0, 0].set_title("Runtime vs layers")
+    plot_metric(axes[0, 0], step_range, data["runtime_seconds"], "time (seconds)")
+    axes[0, 0].set_title(f"Runtime vs layers (n_spin={data['n_spinful_sites']})")
 
-    plot_metric(axes[0, 1], df, "num_terms", "number of terms")
+    plot_metric(axes[0, 1], step_range, data["num_terms"], "number of terms")
     axes[0, 1].set_title("Number of terms vs layers")
 
-    plot_metric(axes[1, 0], df, "memory", "memory (MB)")
+    native_memory_dict = data.get("native_memory_MB", {})
+    plot_metric(
+        axes[1, 0],
+        step_range,
+        data["memory_MB"],
+        "memory (MB)",
+        secondary_dict=native_memory_dict,
+    )
     axes[1, 0].set_title("Memory vs layers")
 
-    plot_metric(axes[1, 1], df, "overlap", "final overlap")
-    axes[1, 1].set_title("Final overlap vs layers")
+    plot_metric(axes[1, 1], step_range, data["expectation_value"], "expectation value")
+    axes[1, 1].set_title("Expectation value vs layers")
 
     fig.tight_layout()
-    fig.savefig(args.output_dir / "majorana_results.png")
+    if native_memory_dict:
+        fig.text(
+            0.5,
+            0.005,
+            "Faint lines: each engine's own native memory accounting (reference only, not the plotted peak)",
+            ha="center",
+            fontsize=8,
+            color="gray",
+        )
+    fig.savefig(args.output_dir / "majorana_results.png", bbox_inches="tight")
 
     if args.show:
         plt.show()
diff --git a/benches/third_party/majorana_prop/results.json b/benches/third_party/majorana_prop/results.json
new file mode 100644
index 00000000..8f53f37d
--- /dev/null
+++ b/benches/third_party/majorana_prop/results.json
@@ -0,0 +1,271 @@
+{
+    "runtime_seconds": {
+        "monoprop": [
+            0.0002520989983167965,
+            0.008602955000242218,
+            0.023853869999584276,
+            0.03253005099759321,
+            0.04263711499879719,
+            0.054794172996480484,
+            0.0691948399944522,
+            0.08690634199228953,
+            0.12709497199466568,
+            0.16048369099371484,
+            0.21669280699279625,
+            0.30345294799190015,
+            0.44508027899064473,
+            0.7188304689916549,
+            1.3468647859917837,
+            2.438142914990749,
+            4.438301431990112,
+            7.850639691991091,
+            13.51275252499181,
+            23.016640379992168,
+            37.62995160999344
+        ],
+        "MajoranaPropagation.jl": [
+            2.4566e-5,
+            0.06570643400000001,
+            1.066572892,
+            1.5099781319999999,
+            1.9726161699999998,
+            2.837843683,
+            3.549592132,
+            4.467205515,
+            5.8747669469999995,
+            7.879885148,
+            11.006226405,
+            16.263505637,
+            25.217318649,
+            40.227838588,
+            64.894778695,
+            103.74468536399999,
+            165.595709796,
+            265.030177113,
+            425.03952550199995,
+            680.1839463509999,
+            1082.0488366109998
+        ]
+    },
+    "n_spinful_sites": 60,
+    "native_memory_MB": {
+        "monoprop": [
+            0.01990985870361328,
+            0.020415306091308594,
+            0.08391475677490234,
+            0.30586719512939453,
+            0.7510213851928711,
+            1.8966608047485352,
+            4.013327598571777,
+            8.111088752746582,
+            12.2099027633667,
+            24.486376762390137,
+            44.62830066680908,
+            81.94080448150635,
+            149.2356767654419,
+            223.22717571258545,
+            397.6318521499634,
+            722.331093788147,
+            1286.3422632217407,
+            1910.7156629562378,
+            3150.523093223572,
+            5411.497309684753,
+            7801.361777305603
+        ],
+        "MajoranaPropagation.jl": [
+            0.00018310546875,
+            0.00518035888671875,
+            0.11733245849609375,
+            0.5146713256835938,
+            1.0477371215820312,
+            2.1161575317382812,
+            4.2556304931640625,
+            8.202163696289062,
+            15.763961791992188,
+            28.70745849609375,
+            52.974884033203125,
+            92.66452026367188,
+            164.24754333496094,
+            304.7784729003906,
+            462.87574005126953,
+            760.7351608276367,
+            1222.6302642822266,
+            1976.4950866699219,
+            3162.911178588867,
+            5007.087661743164,
+            7825.958740234375
+        ]
+    },
+    "memory_MB": {
+        "monoprop": [
+            38.58984375,
+            38.58984375,
+            38.58984375,
+            38.58984375,
+            38.58984375,
+            43.62109375,
+            44.28515625,
+            52.13671875,
+            54.19921875,
+            64.51171875,
+            86.52734375,
+            122.26953125,
+            176.3203125,
+            285.68359375,
+            433.17578125,
+            726.55859375,
+            1208.3515625,
+            2059.21875,
+            3048.30078125,
+            4426.33984375,
+            10636.85546875
+        ],
+        "MajoranaPropagation.jl": [
+            534.9140625,
+            551.109375,
+            557.359375,
+            557.3125,
+            563.76171875,
+            577.4765625,
+            597.67578125,
+            612.3359375,
+            640.3359375,
+            730.30078125,
+            768.875,
+            898.4765625,
+            1081.39453125,
+            1229.63671875,
+            1765.05859375,
+            2044.484375,
+            2949.34375,
+            4322.48046875,
+            6477.3671875,
+            10206.3359375,
+            15694.859375
+        ]
+    },
+    "step_range": [
+        0,
+        1,
+        2,
+        3,
+        4,
+        5,
+        6,
+        7,
+        8,
+        9,
+        10,
+        11,
+        12,
+        13,
+        14,
+        15,
+        16,
+        17,
+        18,
+        19,
+        20
+    ],
+    "num_threads": {
+        "monoprop": "24",
+        "MajoranaPropagation.jl": 24
+    },
+    "expectation_value": {
+        "monoprop": [
+            0.0,
+            0.009736140422307249,
+            0.03827070394062437,
+            0.08364455094749212,
+            0.1427925477622065,
+            0.21181820471470336,
+            0.2863324191238933,
+            0.36181749808823005,
+            0.43398016578994325,
+            0.49905876932925186,
+            0.5540634920911527,
+            0.5969321865128213,
+            0.6265985937915881,
+            0.6429730022355784,
+            0.646851751012522,
+            0.6397661922659064,
+            0.6237899786874557,
+            0.6013251208948387,
+            0.5748905367524423,
+            0.5469245051110052,
+            0.5196156652467371
+        ],
+        "MajoranaPropagation.jl": [
+            0.0,
+            0.009736140422307242,
+            0.03827070330991876,
+            0.08364453430830862,
+            0.1427924873717963,
+            0.21181811458903516,
+            0.2863323455315943,
+            0.3618176213290471,
+            0.4339803099823046,
+            0.49905894198645734,
+            0.5540635608915114,
+            0.5969324107189943,
+            0.6265984814098748,
+            0.6429728162492787,
+            0.6468519475075789,
+            0.6397668861944836,
+            0.6237907501255578,
+            0.6013260020506138,
+            0.5748921471542362,
+            0.5469264564869873,
+            0.51961452777913
+        ]
+    },
+    "n_layers": 20,
+    "num_terms": {
+        "monoprop": [
+            1,
+            16,
+            1344,
+            5794,
+            15637,
+            35883,
+            73597,
+            141845,
+            267386,
+            488053,
+            872870,
+            1513509,
+            2606766,
+            4385183,
+            7276415,
+            11882982,
+            19078830,
+            30227167,
+            47122113,
+            72352237,
+            109381056
+        ],
+        "MajoranaPropagation.jl": [
+            2,
+            17,
+            1162,
+            4351,
+            11380,
+            25649,
+            52784,
+            101198,
+            187217,
+            338317,
+            597601,
+            1025597,
+            1754476,
+            2933083,
+            4865089,
+            7957480,
+            12859755,
+            20582020,
+            32625568,
+            51103291,
+            79406422
+        ]
+    }
+}
diff --git a/benches/third_party/majorana_prop/run_benchmarks.sh b/benches/third_party/majorana_prop/run_benchmarks.sh
index 5e6a2898..02e8b860 100755
--- a/benches/third_party/majorana_prop/run_benchmarks.sh
+++ b/benches/third_party/majorana_prop/run_benchmarks.sh
@@ -4,18 +4,19 @@
 set -euo pipefail
 
 
-export JULIA_NUM_THREADS=8
-export monoprop_NUM_THREADS=8 # noqa: SIM112
+export JULIA_NUM_THREADS=24
+export monoprop_NUM_THREADS=24
+
+
+echo "Running monoprop benchmark (monoprop_NUM_THREADS=${monoprop_NUM_THREADS})"
+uv run python monoprop_hubbard1d_benchmark.py
 
 julia --project=@. -e 'using Pkg; Pkg.instantiate()'
 julia --project=@. -e 'using Pkg; Pkg.precompile()'
 
-echo "Running Julia benchmark (cases 1-15, JULIA_NUM_THREADS=${JULIA_NUM_THREADS})"
-for case in $(seq 1 15); do
-    julia --project=@. julia_hubbard1d_benchmark.jl --case "$case"
-done
+echo "Running Julia benchmark (JULIA_NUM_THREADS=${JULIA_NUM_THREADS})"
+julia --project=@. julia_hubbard1d_benchmark.jl
+
+
 
-echo "Running monoprop benchmark (cases 0-14, monoprop_NUM_THREADS=${monoprop_NUM_THREADS})"
-for case in $(seq 0 14); do
-    uv run python monoprop_hubbard1d_benchmark.py --case "$case"
-done
+uv python plot_results.py
diff --git a/benches/third_party/pauli_prop/_common.py b/benches/third_party/pauli_prop/_common.py
new file mode 100644
index 00000000..69c416fa
--- /dev/null
+++ b/benches/third_party/pauli_prop/_common.py
@@ -0,0 +1,159 @@
+# Copyright 2026 Algorithmiq
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from __future__ import annotations
+
+import json
+import sys
+from dataclasses import dataclass
+from pathlib import Path
+
+BENCH_DIR = Path(__file__).parent
+SETTINGS_FILE = BENCH_DIR / "settings.json"
+RESULTS_FILE = BENCH_DIR / "results.json"
+
+sys.path.insert(0, str(BENCH_DIR.parent))
+from bench_common import RssPeakSampler  # noqa: E402, F401
+
+
+@dataclass(frozen=True)
+class Settings:
+    """Simulation parameters shared by every engine, loaded from settings.json."""
+
+    nx: int
+    ny: int
+    nq: int
+    hx: float
+    hz: float
+    j: float
+    dt: float
+    theta_x: float
+    theta_z: float
+    theta_zz: float
+    step_range: range
+    lower_atol: float
+    max_pauli_weight: int
+    obs_qubits: tuple[int, int]
+    grid_edges: list[tuple[int, int]]
+
+
+def _grid_edges(nx: int, ny: int) -> list[tuple[int, int]]:
+    """Nearest-neighbor edges of an nx-by-ny grid, row-major qubit indexing."""
+    edges = []
+    for row in range(ny):
+        for col in range(nx):
+            idx = row * nx + col
+            if col + 1 < nx:
+                edges.append((idx, idx + 1))
+            if row + 1 < ny:
+                edges.append((idx, idx + nx))
+    return edges
+
+
+def load_settings() -> Settings:
+    with open(SETTINGS_FILE) as file:
+        raw = json.load(file)
+    nx, ny = raw["nx"], raw["ny"]
+    nq = nx * ny
+    return Settings(
+        nx=nx,
+        ny=ny,
+        nq=nq,
+        hx=raw["hx"],
+        hz=raw["hz"],
+        j=raw["j"],
+        dt=raw["dt"],
+        theta_x=raw["dt"] * raw["hx"],
+        theta_z=raw["dt"] * raw["hz"],
+        theta_zz=raw["dt"] * raw["j"],
+        step_range=range(raw["step_min"], raw["step_max"] + 1, raw["step_size"]),
+        lower_atol=raw["lower_atol"],
+        max_pauli_weight=nq if raw["cutoff"] is None else raw["cutoff"],
+        obs_qubits=tuple(raw["obs_qubits"]),
+        grid_edges=_grid_edges(nx, ny),
+    )
+
+
+_RESULTS_KEYS = (
+    "step_range",
+    "num_terms",
+    "runtime",
+    "memory",
+    "native_memory",
+    "expvals",
+)
+
+
+def init_results(settings: Settings) -> None:
+    """(Re)create results.json's skeleton from settings.json.
+
+    Called unconditionally by run_monoprop.py, since it's meant to start every full run from
+    scratch. Every other engine script instead calls ensure_results_file(), which only creates
+    the skeleton if results.json is missing/invalid, so it doesn't clobber earlier engines'
+    results when run after them.
+    """
+    RESULTS_FILE.write_text(
+        json.dumps(
+            dict.fromkeys(_RESULTS_KEYS[1:], {})
+            | {"step_range": list(settings.step_range)},
+            indent=4,
+        )
+    )
+
+
+def ensure_results_file(settings: Settings) -> None:
+    """Make sure results.json exists and has the expected shape before this engine starts.
+
+    Lets every engine script be run standalone, in any order, without results.json already
+    existing: without this, a missing/invalid/stale results.json would only surface as a
+    confusing JSONDecodeError from update_results() at the very end — after the engine already
+    ran its full (possibly expensive) simulation. Unlike init_results(), this leaves an
+    already-valid results.json (e.g. with other engines' results already in it) untouched.
+    """
+    try:
+        with open(RESULTS_FILE) as file:
+            data = json.load(file)
+        if not all(key in data for key in _RESULTS_KEYS):
+            raise ValueError("results.json is missing expected keys")
+    except (FileNotFoundError, json.JSONDecodeError, ValueError):
+        init_results(settings)
+
+
+def update_results(
+    label: str,
+    *,
+    runtime: list[float],
+    memory: list[float],
+    expvals: list[float],
+    num_terms: list[int],
+    native_memory: list[float] | None = None,
+) -> None:
+    """Merge one engine's results into the shared results.json (read-modify-write).
+
+    ``memory`` is the OS/driver-level peak (host RSS or GPU device memory, depending on the
+    engine) and is what gets plotted. ``native_memory``, where available, is that engine's own
+    internal accounting (e.g. monoprop's operator-memory API, cuPauliProp's cupy pool) kept for
+    reference/cross-checking — it is not necessarily the true peak, since it can miss memory the
+    engine allocates outside of what it tracks itself.
+    """
+    with open(RESULTS_FILE) as file:
+        data = json.load(file)
+    data["runtime"][label] = runtime
+    data["memory"][label] = memory
+    data["expvals"][label] = expvals
+    data["num_terms"][label] = num_terms
+    if native_memory is not None:
+        data.setdefault("native_memory", {})[label] = native_memory
+    with open(RESULTS_FILE, "w") as file:
+        json.dump(data, file, indent=4)
diff --git a/benches/third_party/pauli_prop/plot_results.py b/benches/third_party/pauli_prop/plot_results.py
index a1005392..266cc478 100644
--- a/benches/third_party/pauli_prop/plot_results.py
+++ b/benches/third_party/pauli_prop/plot_results.py
@@ -35,6 +35,7 @@
 step_range = data["step_range"]
 runtime_dict = data["runtime"]
 memory_dict = data["memory"]
+native_memory_dict = data.get("native_memory", {})
 expvals_dict = data["expvals"]
 
 
@@ -63,13 +64,25 @@ def _style_axes(ax: plt.Axes, ylabel: str) -> None:
 
 for label, runtime in runtime_dict.items():
     steps, values = _filter_from_min_step(step_range, runtime)
-    ax1.plot(steps, values, color=colors[label], label=label)
+    ax1.plot(steps, values, color=colors[label], label=label, marker="o", markersize=4)
 _style_axes(ax1, "Time per step [s]")
 
 for label, memory in memory_dict.items():
     steps, values = _filter_from_min_step(step_range, memory)
-    ax2.plot(steps, values, color=colors[label], label=label)
+    ax2.plot(steps, values, color=colors[label], label=label, marker="o", markersize=4)
+for label, native_memory in native_memory_dict.items():
+    steps, values = _filter_from_min_step(step_range, native_memory)
+    ax2.plot(steps, values, color=colors[label], linestyle="--", alpha=0.5)
 _style_axes(ax2, "Memory per step [MB]")
 
 fig.tight_layout()
-fig.savefig(Path(__file__).parent / "pauli_results.png", dpi=150)
+if native_memory_dict:
+    fig.text(
+        0.5,
+        -0.03,
+        "Dashed: each engine's own native memory accounting (reference only, not the plotted peak)",
+        ha="center",
+        fontsize=8,
+        color="gray",
+    )
+fig.savefig(Path(__file__).parent / "pauli_results.png", dpi=150, bbox_inches="tight")
diff --git a/benches/third_party/pauli_prop/results.json b/benches/third_party/pauli_prop/results.json
index 6e8c97c8..cec07aae 100644
--- a/benches/third_party/pauli_prop/results.json
+++ b/benches/third_party/pauli_prop/results.json
@@ -1,231 +1,301 @@
 {
     "expvals": {
-        "cuPauliProp (GPU)": [
-            0.997502082639013,
-            0.9903354444404697,
-            0.9794361767419918,
-            0.9661809649164091,
-            0.9522199951963057,
-            0.9392112109453805,
-            0.9286632604749514,
-            0.9217215256190198,
-            0.9192581471551303,
-            0.9213746654172797,
-            0.9275425648672762,
-            0.9368014183084314,
-            0.9479586047326162,
-            0.959807462945837,
-            0.9713666467866726,
-            0.9817625295735538,
-            0.9897022306685668,
-            0.9940854516308357,
-            0.9943906072586334,
-            0.9904540465925349,
-            0.9826744215990753
-        ],
         "monoprop": [
             0.997502082639013,
-            0.9903354444404695,
-            0.9794361767419921,
-            0.9661809649070634,
-            0.9522199947834874,
-            0.9392112074970561,
-            0.9286632441788407,
-            0.9217296258792056,
-            0.9192734598537711,
-            0.9213900001024592,
-            0.9275541488151822,
-            0.9368118309137169,
-            0.9479633902932676,
-            0.9598075713501336,
-            0.9713646524559834,
-            0.9817709140803824,
-            0.9897243499921878,
-            0.9941371201900285,
-            0.9944276721765963,
-            0.9904727245014008,
-            0.9826714907890624
+            0.9903354444404696,
+            0.9794363795243637,
+            0.9661812273018453,
+            0.9522201534337065,
+            0.9392107973137432,
+            0.928658478653306,
+            0.9217007651974013,
+            0.9191560692122138,
+            0.9211446372117517,
+            0.9271845753190968,
+            0.936262879961247,
+            0.9471923300179088,
+            0.9587376552397145,
+            0.969822741864718,
+            0.9796969929662139,
+            0.9872955993601161,
+            0.9916867408908515,
+            0.9923197727253984,
+            0.9889797246327313,
+            0.9820133275023492,
+            0.9725123646210144,
+            0.9615984971418958,
+            0.9506665763222667,
+            0.9410126650301269,
+            0.9335833847204386,
+            0.9293465294442301,
+            0.9286933588139175
         ],
         "PauliPropagation.jl": [
             0.997502082639013,
             0.9903354444404696,
-            0.9794361767419826,
-            0.9661809649157846,
-            0.9522199952275356,
-            0.9392112114322317,
-            0.9286632469647392,
-            0.9217215611595609,
-            0.9192583164556223,
-            0.9213751727716332,
-            0.9275437882758828,
-            0.9368026677065466,
-            0.947960191050668,
-            0.959808776235663,
-            0.9713620730260529,
-            0.9817573708646599,
-            0.9896923125323471,
-            0.9940732143365996,
-            0.9943759464284816,
-            0.9904334566154421,
-            0.9826514575855047
+            0.979436175612008,
+            0.9661809422038069,
+            0.9522198487388678,
+            0.939210514948294,
+            0.9286608751710758,
+            0.9217112490651745,
+            0.9191948496469634,
+            0.9211910356841434,
+            0.9271784314013242,
+            0.93621958436388,
+            0.9471566926984755,
+            0.9587790697230953,
+            0.9700810609412963,
+            0.980162477442118,
+            0.9878248525459249,
+            0.9920739814443645,
+            0.9924615969482351,
+            0.988932945770035,
+            0.9818791500670394,
+            0.9724983537122938,
+            0.9617580208511031,
+            0.9510289865354113,
+            0.9414094249138573,
+            0.9339856557135312,
+            0.9296925398484073,
+            0.9290590727386197
         ],
         "QuEra ppvm": [
             0.997502082639013,
             0.9903354444404695,
-            0.9794361709556907,
-            0.9661809639597803,
-            0.9522199930461324,
-            0.939211173609548,
-            0.9286630881228823,
-            0.9217301572595826,
-            0.9192646317194918,
-            0.9213799216222511,
-            0.9275449829326113,
-            0.9368045165263471,
-            0.947960505529967,
-            0.9598101077658275,
-            0.9713725242948469,
-            0.9817783051928322,
-            0.9897214879232376,
-            0.9941040455760576,
-            0.994408812093782,
-            0.9904678303543354,
-            0.9826766600224272
+            0.9794361695706763,
+            0.9661809419183106,
+            0.952219851108817,
+            0.9392104831471718,
+            0.9286607106018888,
+            0.9217238529580906,
+            0.9191949081786188,
+            0.9211849242571408,
+            0.9271566494169908,
+            0.9361971175113651,
+            0.9471358062576694,
+            0.9587629296265585,
+            0.9700533017012288,
+            0.9801074408254876,
+            0.9877817362610698,
+            0.9920378901252231,
+            0.9923792792918043,
+            0.9888079988294428,
+            0.9818048059131084,
+            0.9724237538531765,
+            0.9617113062712275,
+            0.9509518686115404,
+            0.9412854838029944,
+            0.9338317435872012,
+            0.9295056908756767,
+            0.9288856203834183
         ],
         "Qiskit pauli-prop": [
             0.997502082639013,
-            0.9902986523563847,
-            0.979410268438022,
-            0.9662328011860434,
-            0.95233216052613,
-            0.9393290880734805,
-            0.9286509501231658,
-            0.921402878684645,
-            0.918416768030144,
-            0.9200459604215262,
-            0.926371295132605,
-            0.9364729500095121,
-            0.9485473738161867,
-            0.9609797327252398,
-            0.9722174479178898,
-            0.9808890727517166,
-            0.987202501099926,
-            0.9913511314008273,
-            0.9925835343091947,
-            0.990196575726684,
-            0.9839317508567623
+            0.9902616416411947,
+            0.9793735672305858,
+            0.9662045535658592,
+            0.9523044065600075,
+            0.9393045616009582,
+            0.9286284339439106,
+            0.9213619296914645,
+            0.9183461621464206,
+            0.9197422831890821,
+            0.9257819331883357,
+            0.9355874815011092,
+            0.9473519646429599,
+            0.9595351288671997,
+            0.9706431446987511,
+            0.9793536253319963,
+            0.985549673916859,
+            0.989473983293285,
+            0.9906124875227108,
+            0.9885381901105181,
+            0.9829097988824232,
+            0.9739113012525142,
+            0.9622651253609089,
+            0.9500765782737872,
+            0.9390899221282094,
+            0.9307655543344951,
+            0.9263461450046511,
+            0.9258467634203383
+        ],
+        "cuPauliProp (GPU)": [
+            0.997502082639013,
+            0.9903354444404697,
+            0.9794361756120173,
+            0.966180942204426,
+            0.9522198487076173,
+            0.9392105144690244,
+            0.9286608888777017,
+            0.9217122128671503,
+            0.9191946799798976,
+            0.9211905093829498,
+            0.9271772949331979,
+            0.9362137063468555,
+            0.9471504289968796,
+            0.9587710132467064,
+            0.9700774110665801,
+            0.9801616099109822,
+            0.9878274775819245,
+            0.9920739099660344,
+            0.9924591024470226,
+            0.9889358100230531,
+            0.9818742136637044,
+            0.9724789876050834,
+            0.9617425377637874,
+            0.9510140541676235,
+            0.9413815283280191,
+            0.9339539276043874,
+            0.9296816450847989,
+            0.9290742115849888
         ]
     },
     "runtime": {
-        "cuPauliProp (GPU)": [
-            0.020176051184535027,
-            0.019714422058314085,
-            0.02136979578062892,
-            0.021861208137124777,
-            0.022774800192564726,
-            0.02481368323788047,
-            0.02854348300024867,
-            0.029668635223060846,
-            0.030517147853970528,
-            0.033659269101917744,
-            0.03531042719259858,
-            0.03751846496015787,
-            0.039435874205082655,
-            0.04148514196276665,
-            0.04600577801465988,
-            0.05237875320017338,
-            0.06031936779618263,
-            0.1818047002889216,
-            0.2003558687865734,
-            0.24681029003113508
-        ],
         "monoprop": [
-            0.0022588460706174374,
-            0.002730349078774452,
-            0.003645222634077072,
-            0.004602230153977871,
-            0.005807321984320879,
-            0.007761758286505938,
-            0.010728122666478157,
-            0.015644437167793512,
-            0.022155003156512976,
-            0.030679493211209774,
-            0.04855358274653554,
-            0.07472044182941318,
-            0.10821856698021293,
-            0.1829305151477456,
-            0.2592461039312184,
-            0.34765971498563886,
-            0.5380650600418448,
-            0.753288147971034,
-            1.133904397021979,
-            1.6045584678649902
+            0.0097709740002756,
+            0.009900220000417903,
+            0.01062245899811387,
+            0.010208111998508684,
+            0.010301656999217812,
+            0.011290021997410804,
+            0.0119530299998587,
+            0.012820420997741167,
+            0.014460218000749592,
+            0.01662767300149426,
+            0.021685593004804105,
+            0.025799656999879517,
+            0.03668873199785594,
+            0.049491324003611226,
+            0.06257617600203957,
+            0.08155569299560739,
+            0.10271872900193557,
+            0.15972542300005443,
+            0.26723438299814006,
+            0.4563296580017777,
+            0.5873996790032834,
+            0.8560777880047681,
+            1.2323196020006435,
+            1.717144444999576,
+            2.43780124100158,
+            3.313076704995183,
+            4.471609448999516
         ],
         "PauliPropagation.jl": [
-            0.000487617,
-            0.001201636,
-            0.002411776,
-            0.005773478,
-            0.010534435,
-            0.02769105,
-            0.049536102,
-            0.088493142,
-            0.15495462,
-            0.280121239,
-            0.528610672,
-            0.762137845,
-            1.214638218,
-            2.272540931,
-            3.176273919,
-            5.286394518,
-            7.245942986,
-            11.462868432,
-            15.114584049,
-            24.351451212
+            0.937622339,
+            0.632157357,
+            0.703072464,
+            0.627356295,
+            0.796526843,
+            0.659159013,
+            0.701493095,
+            0.833921075,
+            0.917434561,
+            1.408576301,
+            1.80445128,
+            2.311063839,
+            2.965923955,
+            4.003973926,
+            5.832220713,
+            8.279620075,
+            11.971029624,
+            16.877259221,
+            24.415467044,
+            32.886627487,
+            42.673318014,
+            57.788362696,
+            77.248532745,
+            103.575524758,
+            139.279354003,
+            185.886735186,
+            246.66351665
         ],
         "QuEra ppvm": [
-            0.00018265610560774803,
-            0.000364821869879961,
-            0.0006602783687412739,
-            0.0016023130156099796,
-            0.0034242477267980576,
-            0.007562015671283007,
-            0.0130642163567245,
-            0.021341342013329268,
-            0.04229155322536826,
-            0.07402261719107628,
-            0.128699810244143,
-            0.23124448582530022,
-            0.38588487124070525,
-            0.6836787946522236,
-            1.3010696759447455,
-            2.5424343938939273,
-            4.085273690987378,
-            6.896651620976627,
-            9.923423228785396,
-            14.719230208080262
+            0.0009772880002856255,
+            0.0019692240020958707,
+            0.003414564002014231,
+            0.006883168003696483,
+            0.01336483699560631,
+            0.026874936003878247,
+            0.049350055000104476,
+            0.08147332700173138,
+            0.1398485649988288,
+            0.24163390600006096,
+            0.41795454300154233,
+            0.6846707399963634,
+            1.1364982200029772,
+            2.14537202000065,
+            4.182404636005231,
+            6.776013484995929,
+            10.121395098001813,
+            14.356937502001529,
+            20.45524097100133,
+            29.319063470997207,
+            41.06935781400534,
+            58.04860244700103,
+            84.5264299709961,
+            119.09453202800069,
+            166.27936962900276,
+            240.10627245400246,
+            323.1610831039943
         ],
         "Qiskit pauli-prop": [
-            0.007913357112556696,
-            0.008004344068467617,
-            0.009662941563874483,
-            0.012174050323665142,
-            0.017478680703788996,
-            0.027277078945189714,
-            0.04581389995291829,
-            0.07399411406368017,
-            0.12357027316465974,
-            0.21096500102430582,
-            0.34313367400318384,
-            0.5556078860536218,
-            0.8972947858273983,
-            1.4494283101521432,
-            2.3114430508576334,
-            3.8287112680263817,
-            5.7641011090017855,
-            8.530939413700253,
-            12.57545457687229,
-            18.853173348121345
+            0.0910057479995885,
+            0.09350110300147207,
+            0.09966104299383005,
+            0.11175302699848544,
+            0.1350917930030846,
+            0.17215780800324865,
+            0.23493830300139962,
+            0.3418267029992421,
+            0.5214368059969274,
+            0.8170689979961026,
+            1.2857159330014838,
+            2.0053665779996663,
+            3.158475438001915,
+            4.8786017739985255,
+            7.400178270996548,
+            11.198802080994938,
+            16.48318638800265,
+            23.556144617999962,
+            33.42640477000532,
+            46.69200970899692,
+            65.27846676199988,
+            91.56253991199628,
+            125.43558071399457,
+            172.17275571200298,
+            235.64995327400538,
+            320.03587235399755,
+            437.8359219260019
+        ],
+        "cuPauliProp (GPU)": [
+            0.18526526999630732,
+            0.18972675999975763,
+            0.1926419649971649,
+            0.19678318400110584,
+            0.2010178460041061,
+            0.20950794600503286,
+            0.21115881600417197,
+            0.21557678499812027,
+            0.2187041910001426,
+            0.22579057099937927,
+            0.23606160299823387,
+            0.25517424199642846,
+            0.2714905280008679,
+            0.30626787999790395,
+            0.35749711799871875,
+            0.419371971001965,
+            0.49677398100175196,
+            0.5983771039973362,
+            0.890931271998852,
+            1.250609467002505,
+            1.8410984980000649,
+            2.3888087309969706,
+            3.464189799000451,
+            4.2792873379949015,
+            6.69145935200504,
+            8.090263483994931,
+            12.02248638200399
         ]
     },
     "step_range": [
@@ -249,240 +319,409 @@
         34,
         36,
         38,
-        40
+        40,
+        42,
+        44,
+        46,
+        48,
+        50,
+        52,
+        54
     ],
     "memory": {
-        "cuPauliProp (GPU)": [
-            0.0068359375,
-            0.02099609375,
-            0.0390625,
-            0.08251953125,
-            0.162109375,
-            0.2978515625,
-            0.5458984375,
-            0.958984375,
-            1.68505859375,
-            2.8603515625,
-            4.74658203125,
-            7.6689453125,
-            12.1708984375,
-            18.57080078125,
-            28.3544921875,
-            41.5859375,
-            60.4814453125,
-            86.59228515625,
-            122.46875,
-            173.025390625,
-            242.6591796875
-        ],
         "monoprop": [
-            0.02881336212158203,
-            0.04080677032470703,
-            0.06046581268310547,
-            0.14527225494384766,
-            0.2773103713989258,
-            0.5006303787231445,
-            0.8862504959106445,
-            1.7319231033325195,
-            3.103184700012207,
-            5.1228837966918945,
-            7.2314958572387695,
-            13.951741218566895,
-            24.81438159942627,
-            28.878422737121582,
-            50.83928394317627,
-            82.69517993927002,
-            113.43906116485596,
-            167.34633350372314,
-            228.9335069656372,
-            400.28574085235596,
-            458.79584217071533
+            99.24609375,
+            99.24609375,
+            99.24609375,
+            99.24609375,
+            99.24609375,
+            99.24609375,
+            99.24609375,
+            99.24609375,
+            104.984375,
+            106.0,
+            112.1875,
+            117.6015625,
+            129.71875,
+            140.33203125,
+            152.49609375,
+            177.30078125,
+            215.34375,
+            269.90625,
+            323.6875,
+            408.66796875,
+            545.92578125,
+            695.8046875,
+            939.47265625,
+            1238.2265625,
+            1702.04296875,
+            2192.26953125,
+            3068.6796875,
+            3920.0234375
         ],
         "PauliPropagation.jl": [
-            0.0062255859375,
-            0.0245361328125,
-            0.0977783203125,
-            0.0977783203125,
-            0.3907470703125,
-            0.3907470703125,
-            1.5626220703125,
-            1.5626220703125,
-            3.1251220703125,
-            6.2501220703125,
-            6.2501220703125,
-            12.5001220703125,
-            25.0001220703125,
-            25.0001220703125,
-            50.0001220703125,
-            50.0001220703125,
-            100.0001220703125,
-            200.0001220703125,
-            200.0001220703125,
-            200.0001220703125,
-            400.0001220703125
+            554.5078125,
+            562.78125,
+            563.6015625,
+            563.6015625,
+            563.6015625,
+            565.1171875,
+            559.171875,
+            578.06640625,
+            593.17578125,
+            614.23828125,
+            638.3359375,
+            669.609375,
+            714.609375,
+            756.48046875,
+            860.4921875,
+            912.06640625,
+            1006.12890625,
+            1198.44921875,
+            1331.4609375,
+            1656.81640625,
+            2051.0,
+            2623.4375,
+            3161.921875,
+            4197.296875,
+            5691.19921875,
+            7193.875,
+            8734.47265625,
+            11237.98046875
         ],
         "QuEra ppvm": [
-            0.99609375,
-            1.09765625,
-            1.09765625,
-            1.1328125,
-            1.3125,
-            1.66015625,
-            2.46875,
-            2.9609375,
-            4.42578125,
-            7.16796875,
-            12.328125,
-            17.77734375,
-            23.87109375,
-            32.74609375,
-            45.66015625,
-            64.7578125,
-            91.03125,
-            129.25,
-            180.68359375,
-            250.69921875,
-            348.30859375
+            150.1796875,
+            150.1796875,
+            150.1796875,
+            150.1796875,
+            151.12109375,
+            151.37890625,
+            151.63671875,
+            153.69921875,
+            158.2734375,
+            161.2890625,
+            168.16796875,
+            174.16796875,
+            183.9453125,
+            196.8984375,
+            217.1328125,
+            250.6328125,
+            299.640625,
+            356.64453125,
+            448.66015625,
+            548.64453125,
+            720.6640625,
+            900.65625,
+            1138.68359375,
+            1530.67578125,
+            1960.6875,
+            2674.6875,
+            3434.6953125,
+            4710.6953125
         ],
         "Qiskit pauli-prop": [
-            0.1875,
-            0.45703125,
-            0.6328125,
-            3.01171875,
-            3.546875,
-            6.04296875,
-            9.453125,
-            14.8046875,
-            20.88671875,
-            35.22265625,
-            46.765625,
-            65.19921875,
-            97.421875,
-            153.75390625,
-            242.32421875,
-            398.3203125,
-            508.140625,
-            667.93359375,
-            949.3125,
-            1063.67578125,
-            1388.71875
+            105.2734375,
+            106.00390625,
+            106.31640625,
+            107.0859375,
+            108.11328125,
+            110.171875,
+            116.05859375,
+            127.546875,
+            137.8125,
+            158.859375,
+            182.53515625,
+            220.7109375,
+            288.05859375,
+            391.26171875,
+            579.98828125,
+            816.6015625,
+            1063.75,
+            1450.38671875,
+            1951.3984375,
+            2738.16015625,
+            3688.72265625,
+            5139.078125,
+            6925.5,
+            9614.63671875,
+            12773.3515625,
+            17074.453125,
+            23311.59765625,
+            30644.8203125
+        ],
+        "cuPauliProp (GPU)": [
+            448.25,
+            448.25,
+            452.25,
+            458.25,
+            478.25,
+            530.25,
+            658.25,
+            906.25,
+            1230.25,
+            1754.25,
+            2652.25,
+            4046.25,
+            6220.25,
+            9414.25,
+            14098.25,
+            21972.25,
+            34032.25,
+            51172.25,
+            74560.25,
+            80806.25,
+            73070.25,
+            80468.25,
+            80306.25,
+            80012.25,
+            80066.25,
+            80654.25,
+            79984.25,
+            80494.25
         ]
     },
-    "num_terms": {
-        "cuPauliProp (GPU)": [
-            120,
-            430,
-            832,
-            1767,
-            3507,
-            6469,
-            11892,
-            20924,
-            36784,
-            62461,
-            103654,
-            167490,
-            265849,
-            405649,
-            619391,
-            908426,
-            1321207,
-            1891606,
-            2675337,
-            3779776,
-            5300936
+    "native_memory": {
+        "monoprop": [
+            0.021811485290527344,
+            0.053984642028808594,
+            0.09384822845458984,
+            0.19596195220947266,
+            0.36931705474853516,
+            0.6798334121704102,
+            1.1329317092895508,
+            2.244696617126465,
+            3.527798652648926,
+            5.652615547180176,
+            8.914322853088379,
+            16.2201509475708,
+            25.746647834777832,
+            39.50025653839111,
+            57.21512317657471,
+            87.02963733673096,
+            125.40126514434814,
+            151.31500720977783,
+            215.3877305984497,
+            325.40620136260986,
+            467.08703327178955,
+            672.3922700881958,
+            962.2750978469849,
+            1153.1309022903442,
+            1635.0377168655396,
+            2350.2291383743286,
+            3327.9866762161255,
+            3630.8824434280396
+        ],
+        "PauliPropagation.jl": [
+            0.02861785888671875,
+            0.07315826416015625,
+            0.1643218994140625,
+            0.1643218994140625,
+            0.3488922119140625,
+            0.7205963134765625,
+            1.466888427734375,
+            2.9627304077148438,
+            5.958045959472656,
+            9.327789306640625,
+            9.327789306640625,
+            28.211692810058594,
+            28.211692810058594,
+            52.738128662109375,
+            80.33037567138672,
+            80.33037567138672,
+            132.37162017822266,
+            132.37162017822266,
+            232.91806030273438,
+            430.03279876708984,
+            430.03279876708984,
+            651.7868728637695,
+            1069.2601928710938,
+            1069.2601928710938,
+            1538.9176712036133,
+            2067.282341003418,
+            2997.692581176758,
+            2997.692581176758
         ],
+        "cuPauliProp (GPU)": [
+            0.04736328125,
+            0.1572265625,
+            0.294921875,
+            0.61767578125,
+            1.21630859375,
+            2.18701171875,
+            3.9521484375,
+            6.82568359375,
+            11.75439453125,
+            19.64892578125,
+            31.93310546875,
+            50.68310546875,
+            79.15869140625,
+            119.01025390625,
+            178.7099609375,
+            258.55224609375,
+            370.36279296875,
+            521.47412109375,
+            727.70068359375,
+            1012.0078125,
+            1400.3779296875,
+            1931.33984375,
+            2657.013671875,
+            3634.3740234375,
+            4924.080078125,
+            6620.19482421875,
+            8834.146484375,
+            11651.52783203125
+        ]
+    },
+    "num_terms": {
         "monoprop": [
             120,
-            430,
-            832,
-            1768,
-            3522,
-            6512,
-            11979,
-            21078,
-            36994,
-            63049,
-            105073,
-            170823,
-            274222,
-            421161,
-            649218,
-            969166,
-            1430033,
-            2075465,
-            2966901,
-            4233793,
-            6002213
+            414,
+            768,
+            1630,
+            3255,
+            5818,
+            10454,
+            17934,
+            30953,
+            51534,
+            83801,
+            133713,
+            209965,
+            319328,
+            482566,
+            710792,
+            1032242,
+            1465836,
+            2056508,
+            2871315,
+            3982863,
+            5509645,
+            7585864,
+            10386236,
+            14099731,
+            18948708,
+            25241196,
+            33309328
         ],
         "PauliPropagation.jl": [
             120,
-            430,
-            831,
-            1765,
-            3506,
-            6462,
-            11886,
-            20914,
-            36810,
-            62448,
-            103589,
-            167411,
-            265753,
-            405091,
-            618061,
-            906588,
-            1318324,
-            1887517,
-            2670435,
-            3774991,
-            5295608
+            414,
+            779,
+            1641,
+            3242,
+            5837,
+            10554,
+            18235,
+            31461,
+            52522,
+            85320,
+            135506,
+            211626,
+            317843,
+            476925,
+            689973,
+            988078,
+            1390964,
+            1942287,
+            2701851,
+            3740573,
+            5160256,
+            7101940,
+            9718162,
+            13169143,
+            17708646,
+            23635470,
+            31177867
         ],
         "QuEra ppvm": [
             4,
             203,
-            483,
-            1035,
-            2335,
-            4582,
-            8182,
-            14694,
-            26055,
-            45190,
-            76374,
-            124092,
-            203066,
-            315571,
-            479217,
-            722068,
-            1052355,
-            1541676,
-            2189107,
-            3080583,
-            4332743
+            469,
+            1009,
+            2219,
+            4221,
+            7488,
+            13129,
+            22826,
+            38886,
+            64629,
+            103409,
+            166019,
+            253342,
+            378146,
+            560772,
+            805621,
+            1161443,
+            1624827,
+            2256599,
+            3122764,
+            4258841,
+            5765752,
+            7771081,
+            10461746,
+            14030516,
+            18743094,
+            24878659
         ],
         "Qiskit pauli-prop": [
             120,
-            430,
-            832,
-            1768,
-            3522,
-            6512,
-            11979,
-            21078,
-            36994,
-            63049,
-            105073,
-            170823,
-            274222,
-            421161,
-            649218,
-            969166,
-            1430033,
-            2075465,
-            2966901,
-            4233793,
-            6002213
+            414,
+            768,
+            1630,
+            3255,
+            5818,
+            10454,
+            17934,
+            30953,
+            51534,
+            83801,
+            133713,
+            209965,
+            319328,
+            482566,
+            710792,
+            1032242,
+            1465836,
+            2056508,
+            2871315,
+            3982863,
+            5509645,
+            7585864,
+            10386236,
+            14099731,
+            18948708,
+            25241196,
+            33309328
+        ],
+        "cuPauliProp (GPU)": [
+            120,
+            414,
+            780,
+            1643,
+            3244,
+            5843,
+            10562,
+            18249,
+            31430,
+            52548,
+            85406,
+            135554,
+            211720,
+            318312,
+            477990,
+            691552,
+            990613,
+            1394792,
+            1946394,
+            2706837,
+            3745619,
+            5165796,
+            7106774,
+            9720944,
+            13170557,
+            17707199,
+            23628909,
+            31164633
         ]
     }
-}
+}
\ No newline at end of file
diff --git a/benches/third_party/pauli_prop/run_cupauliprop.py b/benches/third_party/pauli_prop/run_cupauliprop.py
new file mode 100644
index 00000000..44610099
--- /dev/null
+++ b/benches/third_party/pauli_prop/run_cupauliprop.py
@@ -0,0 +1,194 @@
+# Copyright 2026 Algorithmiq
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from __future__ import annotations
+
+import threading
+import time
+
+import cupy as cp
+import numpy as np
+from cuquantum.pauliprop.experimental import (
+    LibraryHandle,
+    PauliExpansion,
+    PauliExpansionOptions,
+    PauliRotationGate,
+    Truncation,
+    get_num_packed_integers,
+)
+from tqdm import tqdm
+
+from _common import ensure_results_file, load_settings, update_results
+
+LABEL = "cuPauliProp (GPU)"
+
+
+def _pauli_string_to_packed_integers(
+    paulis: list[str], qubits: list[int], num_qubits: int
+) -> np.ndarray:
+    """Pack a Pauli string into the (x, z) bitfield layout expected by cuPauliProp."""
+    num_packed_ints = get_num_packed_integers(num_qubits)
+    out = np.zeros(num_packed_ints * 2, dtype=np.uint64)
+    x_ptr = out[:num_packed_ints]
+    z_ptr = out[num_packed_ints:]
+    for pauli, qubit in zip(paulis, qubits):
+        int_ind = qubit // 64
+        bit_ind = qubit % 64
+        if pauli in ("X", "Y"):
+            x_ptr[int_ind] |= 1 << bit_ind
+        if pauli in ("Z", "Y"):
+            z_ptr[int_ind] |= 1 << bit_ind
+    return out
+
+
+class GpuMemPeakSampler:
+    """Tracks the GPU's peak used device memory (driver-level, via cudaMemGetInfo) within a
+    resettable window.
+
+    This is the GPU analogue of RssPeakSampler: rather than trusting one library's own pool
+    accounting (which can miss allocations that bypass that pool, e.g. CUDA context or
+    library-handle overhead), it polls the CUDA driver's own free/total memory report from a
+    background thread, so it reflects everything actually resident on the device. Caveat: if
+    another process shares this GPU concurrently, its usage is included too.
+    """
+
+    def __init__(self, interval_s: float = 1e-3) -> None:
+        self._interval_s = interval_s
+        self._peak_bytes = 0
+        self._stop_event = threading.Event()
+        self._thread = threading.Thread(target=self._poll_loop, daemon=True)
+
+    @staticmethod
+    def _used_bytes() -> int:
+        free_bytes, total_bytes = cp.cuda.runtime.memGetInfo()
+        return total_bytes - free_bytes
+
+    def _poll_loop(self) -> None:
+        while not self._stop_event.wait(self._interval_s):
+            used = self._used_bytes()
+            if used > self._peak_bytes:
+                self._peak_bytes = used
+
+    def __enter__(self) -> GpuMemPeakSampler:
+        self._thread.start()
+        return self
+
+    def __exit__(self, *exc_info: object) -> None:
+        self._stop_event.set()
+        self._thread.join()
+
+    def reset(self) -> None:
+        self._peak_bytes = self._used_bytes()
+
+    def peak_mb(self) -> float:
+        return self._peak_bytes / 1024**2
+
+
+class _PoolPeakHook(cp.cuda.MemoryHook):
+    """Event-driven peak of cupy's default memory pool, for reference alongside the GPU-wide
+    peak. Fires on every allocation, so — unlike sampling on a timer — it cannot miss a spike
+    that is allocated and freed within a single step. Only sees allocations routed through
+    cupy's own pool (which cuQuantum's scratch workspace uses by default, but CUDA
+    context/library-handle overhead does not)."""
+
+    name = "PoolPeakHook"
+
+    def __init__(self) -> None:
+        self.peak_bytes = 0
+
+    def malloc_postprocess(self, **kwargs: object) -> None:
+        used = cp.get_default_memory_pool().used_bytes()
+        if used > self.peak_bytes:
+            self.peak_bytes = used
+
+    def reset(self) -> None:
+        self.peak_bytes = cp.get_default_memory_pool().used_bytes()
+
+
+settings = load_settings()
+ensure_results_file(settings)
+
+cupp_handle = LibraryHandle()
+
+num_packed = get_num_packed_integers(settings.nq)
+cupp_xz = cp.zeros((1, 2 * num_packed), dtype=cp.uint64)
+cupp_coefs = cp.ones((1,), dtype=cp.float64)
+cupp_xz[0] = cp.asarray(
+    _pauli_string_to_packed_integers(["Z", "Z"], list(settings.obs_qubits), settings.nq)
+)
+
+cupp_expansion = PauliExpansion(
+    library_handle=cupp_handle,
+    num_qubits=settings.nq,
+    num_terms=1,
+    xz_bits=cupp_xz,
+    coeffs=cupp_coefs,
+    options=PauliExpansionOptions(memory_limit="80%", blocking=True),
+)
+cupp_truncation = Truncation(
+    pauli_coeff_cutoff=settings.lower_atol,
+    pauli_weight_cutoff=settings.max_pauli_weight,
+)
+
+cupp_step_gates = [
+    PauliRotationGate(settings.theta_zz, ["Z", "Z"], [i, k])
+    for i, k in settings.grid_edges
+]
+cupp_step_gates += [
+    PauliRotationGate(settings.theta_z, ["Z"], [i]) for i in range(settings.nq)
+]
+cupp_step_gates += [
+    PauliRotationGate(settings.theta_x, ["X"], [i]) for i in range(settings.nq)
+]
+
+runtime: list[float] = []
+memory: list[float] = []
+expvals: list[float] = []
+num_terms: list[int] = []
+native_memory: list[float] = []
+
+pool_hook = _PoolPeakHook()
+
+with GpuMemPeakSampler() as gpu_sampler, pool_hook:
+    for step_idx, _ in enumerate(tqdm(settings.step_range, desc=LABEL)):
+        gpu_sampler.reset()
+        pool_hook.reset()
+        t1 = time.perf_counter()
+        for gate in reversed(cupp_step_gates):
+            cupp_expansion = cupp_expansion.apply_gate(
+                gate,
+                truncation=cupp_truncation,
+                adjoint=True,
+                sort_order=None,
+                keep_duplicates=False,
+            )
+        trace_significand, trace_exponent = cupp_expansion.trace_with_zero_state()
+        cupp_expval = float(trace_significand * np.exp2(trace_exponent))
+        t2 = time.perf_counter()
+
+        if step_idx > 0:
+            runtime.append(t2 - t1)
+        expvals.append(cupp_expval)
+        num_terms.append(cupp_expansion.num_terms)
+        memory.append(gpu_sampler.peak_mb())
+        native_memory.append(pool_hook.peak_bytes / 1024**2)
+
+update_results(
+    LABEL,
+    runtime=runtime,
+    memory=memory,
+    expvals=expvals,
+    num_terms=num_terms,
+    native_memory=native_memory,
+)
diff --git a/benches/third_party/pauli_prop/run_model.jl b/benches/third_party/pauli_prop/run_model.jl
index 2c9e20fe..4dc50cdd 100644
--- a/benches/third_party/pauli_prop/run_model.jl
+++ b/benches/third_party/pauli_prop/run_model.jl
@@ -16,6 +16,8 @@ using PauliPropagation
 using JSON
 using ProgressMeter
 
+include(joinpath(@__DIR__, "..", "bench_common.jl"))
+
 settings = JSON.parsefile(joinpath(@__DIR__, "settings.json"))
 
 nx, ny = settings["nx"], settings["ny"]
@@ -41,15 +43,19 @@ append!(step_parameters, fill(theta_zz, length(topology)))
 append!(step_parameters, fill(theta_z, nq))
 append!(step_parameters, fill(theta_x, nq))
 
-pauli_sum = PauliSum(nq)
+pauli_sum = VectorPauliSum(PauliSum(nq))
 add!(pauli_sum, [:Z, :Z], collect(obs_qubits), 1.0)
 
 num_terms = Int[]
 runtime = Float64[]
 memory = Float64[]
+native_memory = Float64[]
 expvals = Float64[]
 
+sampler = start_sampler()
+
 @showprogress for (step_idx, num_steps) in enumerate(step_range)
+    reset!(sampler)
     t1 = time_ns()
     global pauli_sum = propagate(
         step_circuit, pauli_sum, step_parameters;
@@ -62,15 +68,20 @@ expvals = Float64[]
         push!(runtime, (t2 - t1) / 1e9)
     end
     push!(num_terms, length(pauli_sum))
-    push!(memory, Base.summarysize(pauli_sum) / 1024^2)
+    push!(memory, peak_mb(sampler))
+    push!(native_memory, Base.summarysize(pauli_sum) / 1024^2)
     push!(expvals, expval)
 end
 
+stop!(sampler)
+
 results_file = joinpath(@__DIR__, "results.json")
 data = JSON.parsefile(results_file)
 data["num_terms"]["PauliPropagation.jl"] = num_terms
 data["runtime"]["PauliPropagation.jl"] = runtime
 data["memory"]["PauliPropagation.jl"] = memory
+haskey(data, "native_memory") || (data["native_memory"] = Dict())
+data["native_memory"]["PauliPropagation.jl"] = native_memory
 data["expvals"]["PauliPropagation.jl"] = expvals
 
 open(results_file, "w") do file
diff --git a/benches/third_party/pauli_prop/run_model.py b/benches/third_party/pauli_prop/run_model.py
deleted file mode 100644
index 78323884..00000000
--- a/benches/third_party/pauli_prop/run_model.py
+++ /dev/null
@@ -1,245 +0,0 @@
-# Copyright 2026 Algorithmiq
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#     http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-from __future__ import annotations
-
-import json
-import time
-from pathlib import Path
-
-import cupy as cp
-import numpy as np
-import psutil
-from cuquantum.pauliprop.experimental import (
-    LibraryHandle,
-    PauliExpansion,
-    PauliExpansionOptions,
-    PauliRotationGate,
-    Truncation,
-    get_num_packed_integers,
-)
-from monoprop import PauliPropagator
-from monoprop.qiskit_conversion import from_qiskit_circuit, from_qiskit_operator
-from ppvm import PauliSum
-from qiskit.circuit import QuantumCircuit
-from qiskit.quantum_info import SparsePauliOp
-from tqdm import tqdm
-
-from pauli_prop import propagate_through_circuit
-
-
-def _pauli_string_to_packed_integers(
-    paulis: list[str], qubits: list[int], num_qubits: int
-) -> np.ndarray:
-    """Pack a Pauli string into the (x, z) bitfield layout expected by cuPauliProp."""
-    num_packed_ints = get_num_packed_integers(num_qubits)
-    out = np.zeros(num_packed_ints * 2, dtype=np.uint64)
-    x_ptr = out[:num_packed_ints]
-    z_ptr = out[num_packed_ints:]
-    for pauli, qubit in zip(paulis, qubits):
-        int_ind = qubit // 64
-        bit_ind = qubit % 64
-        if pauli in ("X", "Y"):
-            x_ptr[int_ind] |= 1 << bit_ind
-        if pauli in ("Z", "Y"):
-            z_ptr[int_ind] |= 1 << bit_ind
-    return out
-
-
-def _grid_edges(nx: int, ny: int) -> list[tuple[int, int]]:
-    """Nearest-neighbor edges of an nx-by-ny grid, row-major qubit indexing."""
-    edges = []
-    for row in range(ny):
-        for col in range(nx):
-            idx = row * nx + col
-            if col + 1 < nx:
-                edges.append((idx, idx + 1))
-            if row + 1 < ny:
-                edges.append((idx, idx + nx))
-    return edges
-
-
-with open(Path(__file__).parent / "settings.json") as settings_file:
-    settings = json.load(settings_file)
-
-nx, ny = settings["nx"], settings["ny"]
-nq = nx * ny
-hx = settings["hx"]
-hz = settings["hz"]
-j = settings["j"]
-dt = settings["dt"]
-
-step_range = range(
-    settings["step_min"], settings["step_max"] + 1, settings["step_size"]
-)
-lower_atol = settings["lower_atol"]
-max_pauli_weight = nq if settings["cutoff"] is None else settings["cutoff"]
-obs_qubits = tuple(settings["obs_qubits"])
-
-labels = [
-    "monoprop",
-    "QuEra ppvm",
-    "Qiskit pauli-prop",
-    "cuPauliProp (GPU)",
-]
-runtime_dict = {label: [] for label in labels}
-expvals_dict = {label: [] for label in labels}
-num_terms_dict = {label: [] for label in labels}
-memory_dict = {label: [] for label in labels}
-
-process = psutil.Process()
-ppvm_mem_bytes = 0
-qiskit_mem_bytes = 0
-
-theta_x = dt * hx
-theta_z = dt * hz
-theta_zz = dt * j
-grid_edges = _grid_edges(nx, ny)
-
-step_circ = QuantumCircuit(nq)
-for i, k in grid_edges:
-    step_circ.rzz(theta_zz, i, k)
-for i in range(nq):
-    step_circ.rz(theta_z, i)
-for i in range(nq):
-    step_circ.rx(theta_x, i)
-
-obs = SparsePauliOp.from_sparse_list([("ZZ", list(obs_qubits), 1.0)], num_qubits=nq)
-
-mp_circ = from_qiskit_circuit(step_circ, initial_state=[])
-mp_obs = from_qiskit_operator(obs)
-mp = PauliPropagator(
-    initial_operator=mp_obs,
-    initial_state=mp_circ.initial_state,
-    cutoff=max_pauli_weight,
-    lower_atol=lower_atol,
-)
-
-ppvm_obs = PauliSum.new(
-    n_qubits=nq,
-    terms=[f"Z{obs_qubits[0]}Z{obs_qubits[1]}"],
-    min_abs_coeff=lower_atol,
-    max_pauli_weight=max_pauli_weight,
-)
-
-qiskit_obs = obs
-
-cupp_handle = LibraryHandle()
-
-num_packed = get_num_packed_integers(nq)
-cupp_xz = cp.zeros((1, 2 * num_packed), dtype=cp.uint64)
-cupp_coefs = cp.ones((1,), dtype=cp.float64)
-cupp_xz[0] = cp.asarray(
-    _pauli_string_to_packed_integers(["Z", "Z"], list(obs_qubits), nq)
-)
-
-cupp_expansion = PauliExpansion(
-    library_handle=cupp_handle,
-    num_qubits=nq,
-    num_terms=1,
-    xz_bits=cupp_xz,
-    coeffs=cupp_coefs,
-    options=PauliExpansionOptions(memory_limit="80%", blocking=True),
-)
-cupp_truncation = Truncation(
-    pauli_coeff_cutoff=lower_atol,
-    pauli_weight_cutoff=max_pauli_weight,
-)
-
-cupp_step_gates = [
-    PauliRotationGate(theta_zz, ["Z", "Z"], [i, k]) for i, k in grid_edges
-]
-cupp_step_gates += [PauliRotationGate(theta_z, ["Z"], [i]) for i in range(nq)]
-cupp_step_gates += [PauliRotationGate(theta_x, ["X"], [i]) for i in range(nq)]
-
-for step_idx, _ in enumerate(tqdm(step_range, desc="Running simulations")):
-    t1 = time.perf_counter()
-    mp.propagate(mp_circ)
-    expval = mp.expectation_value()
-    t2 = time.perf_counter()
-    if step_idx > 0:
-        runtime_dict["monoprop"].append(t2 - t1)
-    expvals_dict["monoprop"].append(expval)
-    num_terms_dict["monoprop"].append(mp.size())
-    memory_dict["monoprop"].append(mp._simulator.operator_memory_bytes() / 1024**2)
-
-    # ppvm exposes no memory accounting: approximate it via RSS growth over its own step.
-    mem_before = process.memory_info().rss
-    t1 = time.perf_counter()
-    for i, k in grid_edges:
-        ppvm_obs.rzz(i, k, theta_zz)
-    for i in range(nq):
-        ppvm_obs.rz(i, theta_z)
-    for i in range(nq):
-        ppvm_obs.rx(i, theta_x)
-    ppvm_expval = ppvm_obs.overlap_with_zero()
-    t2 = time.perf_counter()
-    ppvm_mem_bytes += max(0, process.memory_info().rss - mem_before)
-
-    if step_idx > 0:
-        runtime_dict["QuEra ppvm"].append(t2 - t1)
-    expvals_dict["QuEra ppvm"].append(ppvm_expval)
-    num_terms_dict["QuEra ppvm"].append(len(ppvm_obs))
-    memory_dict["QuEra ppvm"].append(ppvm_mem_bytes / 1024**2)
-
-    # max_terms tracks monoprop's own term count at this step: this API has no weight-based
-    # cutoff, only a mandatory positive term budget.
-    mem_before = process.memory_info().rss
-    max_terms = mp.size()
-    t1 = time.perf_counter()
-    qiskit_obs, _ = propagate_through_circuit(
-        qiskit_obs, step_circ, max_terms=max_terms, atol=lower_atol, frame="h"
-    )
-    qiskit_expval = float(qiskit_obs.coeffs[~qiskit_obs.paulis.x.any(axis=1)].sum())
-    t2 = time.perf_counter()
-    qiskit_mem_bytes += max(0, process.memory_info().rss - mem_before)
-
-    if step_idx > 0:
-        runtime_dict["Qiskit pauli-prop"].append(t2 - t1)
-    expvals_dict["Qiskit pauli-prop"].append(qiskit_expval)
-    num_terms_dict["Qiskit pauli-prop"].append(len(qiskit_obs))
-    memory_dict["Qiskit pauli-prop"].append(qiskit_mem_bytes / 1024**2)
-
-    t1 = time.perf_counter()
-    for gate in reversed(cupp_step_gates):
-        cupp_expansion = cupp_expansion.apply_gate(
-            gate,
-            truncation=cupp_truncation,
-            adjoint=True,
-            sort_order=None,
-            keep_duplicates=False,
-        )
-    trace_significand, trace_exponent = cupp_expansion.trace_with_zero_state()
-    cupp_expval = float(trace_significand * np.exp2(trace_exponent))
-    t2 = time.perf_counter()
-    if step_idx > 0:
-        runtime_dict["cuPauliProp (GPU)"].append(t2 - t1)
-    expvals_dict["cuPauliProp (GPU)"].append(cupp_expval)
-    num_terms_dict["cuPauliProp (GPU)"].append(cupp_expansion.num_terms)
-    memory_dict["cuPauliProp (GPU)"].append(
-        cp.get_default_memory_pool().used_bytes() / 1024**2
-    )
-
-with open(Path(__file__).parent / "results.json", "w") as file:
-    json.dump(
-        {
-            "step_range": list(step_range),
-            "num_terms": num_terms_dict,
-            "runtime": runtime_dict,
-            "memory": memory_dict,
-            "expvals": expvals_dict,
-        },
-        file,
-        indent=4,
-    )
diff --git a/benches/third_party/pauli_prop/run_model.sh b/benches/third_party/pauli_prop/run_model.sh
new file mode 100755
index 00000000..26d8317c
--- /dev/null
+++ b/benches/third_party/pauli_prop/run_model.sh
@@ -0,0 +1,18 @@
+#!/usr/bin/env bash
+# Run the Python Pauli-propagation engines (monoprop, QuEra ppvm, Qiskit pauli-prop,
+# cuPauliProp) back to back, each in its own process, writing results.json.
+#
+# Each engine runs to completion before the next starts: this keeps the peak-memory measurements
+# (see _common.RssPeakSampler / run_cupauliprop.py's GpuMemPeakSampler) uncontaminated by another
+# engine's allocations, and stops the engines from contending for the same CPU/GPU at the same
+# time, which would otherwise skew both the runtime and memory numbers.
+#
+# Order matters: run_monoprop.py (re)creates results.json from settings.json, and run_qiskit.py
+# reads run_monoprop.py's per-step term counts back out of it to set its own term budget.
+set -euo pipefail
+cd "$(dirname "${BASH_SOURCE[0]}")"
+
+uv run python run_monoprop.py
+uv run python run_ppvm.py
+uv run python run_qiskit.py
+uv run python run_cupauliprop.py
diff --git a/benches/third_party/pauli_prop/run_monoprop.py b/benches/third_party/pauli_prop/run_monoprop.py
new file mode 100644
index 00000000..26ab8385
--- /dev/null
+++ b/benches/third_party/pauli_prop/run_monoprop.py
@@ -0,0 +1,84 @@
+# Copyright 2026 Algorithmiq
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from __future__ import annotations
+
+import time
+
+from monoprop import PauliPropagator
+from monoprop.qiskit_conversion import from_qiskit_circuit, from_qiskit_operator
+from qiskit.circuit import QuantumCircuit
+from qiskit.quantum_info import SparsePauliOp
+from tqdm import tqdm
+
+from _common import RssPeakSampler, init_results, load_settings, update_results
+
+LABEL = "monoprop"
+
+settings = load_settings()
+init_results(settings)
+
+step_circ = QuantumCircuit(settings.nq)
+for i, k in settings.grid_edges:
+    step_circ.rzz(settings.theta_zz, i, k)
+for i in range(settings.nq):
+    step_circ.rz(settings.theta_z, i)
+for i in range(settings.nq):
+    step_circ.rx(settings.theta_x, i)
+
+obs = SparsePauliOp.from_sparse_list(
+    [("ZZ", list(settings.obs_qubits), 1.0)], num_qubits=settings.nq
+)
+
+mp_circ = from_qiskit_circuit(step_circ, initial_state=[])
+mp_obs = from_qiskit_operator(obs)
+mp = PauliPropagator(
+    initial_operator=mp_obs,
+    initial_state=mp_circ.initial_state,
+    cutoff=settings.max_pauli_weight,
+    lower_atol=settings.lower_atol,
+)
+
+runtime: list[float] = []
+memory: list[float] = []
+expvals: list[float] = []
+num_terms: list[int] = []
+native_memory: list[float] = []
+
+with RssPeakSampler() as sampler:
+    for step_idx, _ in enumerate(tqdm(settings.step_range, desc=LABEL)):
+        sampler.reset()
+        t1 = time.perf_counter()
+        mp.propagate(mp_circ)
+        expval = mp.expectation_value()
+        t2 = time.perf_counter()
+
+        if step_idx > 0:
+            runtime.append(t2 - t1)
+        expvals.append(expval)
+        num_terms.append(mp.size())
+        memory.append(sampler.peak_mb())
+        native_memory.append(
+            (mp._simulator.operator_memory_bytes() + mp._simulator.graph_memory_bytes())
+            / 1024**2
+        )
+
+update_results(
+    LABEL,
+    runtime=runtime,
+    memory=memory,
+    expvals=expvals,
+    num_terms=num_terms,
+    native_memory=native_memory,
+)
diff --git a/benches/third_party/pauli_prop/run_ppvm.py b/benches/third_party/pauli_prop/run_ppvm.py
new file mode 100644
index 00000000..bc3dd225
--- /dev/null
+++ b/benches/third_party/pauli_prop/run_ppvm.py
@@ -0,0 +1,62 @@
+# Copyright 2026 Algorithmiq
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from __future__ import annotations
+
+import time
+
+from ppvm import PauliSum
+from tqdm import tqdm
+
+from _common import RssPeakSampler, ensure_results_file, load_settings, update_results
+
+LABEL = "QuEra ppvm"
+
+settings = load_settings()
+ensure_results_file(settings)
+
+ppvm_obs = PauliSum.new(
+    n_qubits=settings.nq,
+    terms=[f"Z{settings.obs_qubits[0]}Z{settings.obs_qubits[1]}"],
+    min_abs_coeff=settings.lower_atol,
+    max_pauli_weight=settings.max_pauli_weight,
+)
+
+runtime: list[float] = []
+memory: list[float] = []
+expvals: list[float] = []
+num_terms: list[int] = []
+
+with RssPeakSampler() as sampler:
+    for step_idx, _ in enumerate(tqdm(settings.step_range, desc=LABEL)):
+        sampler.reset()
+        t1 = time.perf_counter()
+        for i, k in settings.grid_edges:
+            ppvm_obs.rzz(i, k, settings.theta_zz)
+        for i in range(settings.nq):
+            ppvm_obs.rz(i, settings.theta_z)
+        for i in range(settings.nq):
+            ppvm_obs.rx(i, settings.theta_x)
+        ppvm_expval = ppvm_obs.overlap_with_zero()
+        t2 = time.perf_counter()
+
+        if step_idx > 0:
+            runtime.append(t2 - t1)
+        expvals.append(ppvm_expval)
+        num_terms.append(len(ppvm_obs))
+        memory.append(sampler.peak_mb())
+
+update_results(
+    LABEL, runtime=runtime, memory=memory, expvals=expvals, num_terms=num_terms
+)
diff --git a/benches/third_party/pauli_prop/run_qiskit.py b/benches/third_party/pauli_prop/run_qiskit.py
new file mode 100644
index 00000000..1a4090c3
--- /dev/null
+++ b/benches/third_party/pauli_prop/run_qiskit.py
@@ -0,0 +1,89 @@
+# Copyright 2026 Algorithmiq
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from __future__ import annotations
+
+import json
+import time
+
+from pauli_prop import propagate_through_circuit
+from qiskit.circuit import QuantumCircuit
+from qiskit.quantum_info import SparsePauliOp
+from tqdm import tqdm
+
+from _common import (
+    RESULTS_FILE,
+    RssPeakSampler,
+    ensure_results_file,
+    load_settings,
+    update_results,
+)
+
+LABEL = "Qiskit pauli-prop"
+
+settings = load_settings()
+ensure_results_file(settings)
+
+step_circ = QuantumCircuit(settings.nq)
+for i, k in settings.grid_edges:
+    step_circ.rzz(settings.theta_zz, i, k)
+for i in range(settings.nq):
+    step_circ.rz(settings.theta_z, i)
+for i in range(settings.nq):
+    step_circ.rx(settings.theta_x, i)
+
+qiskit_obs = SparsePauliOp.from_sparse_list(
+    [("ZZ", list(settings.obs_qubits), 1.0)], num_qubits=settings.nq
+)
+
+# Qiskit's propagate_through_circuit API has no weight-based cutoff, only a mandatory positive
+# term budget, so we reuse monoprop's own term count at each step (already in results.json,
+# since run_model.sh always runs run_monoprop.py first) to keep the comparison apples-to-apples.
+with open(RESULTS_FILE) as file:
+    monoprop_num_terms = json.load(file)["num_terms"].get("monoprop")
+if monoprop_num_terms is None:
+    raise RuntimeError(
+        "run_qiskit.py needs monoprop's per-step term counts, which aren't in results.json yet — "
+        "run run_monoprop.py (or run_model.sh) first."
+    )
+
+runtime: list[float] = []
+memory: list[float] = []
+expvals: list[float] = []
+num_terms: list[int] = []
+
+with RssPeakSampler() as sampler:
+    for step_idx, _ in enumerate(tqdm(settings.step_range, desc=LABEL)):
+        max_terms = monoprop_num_terms[step_idx]
+        sampler.reset()
+        t1 = time.perf_counter()
+        qiskit_obs, _ = propagate_through_circuit(
+            qiskit_obs,
+            step_circ,
+            max_terms=max_terms,
+            atol=settings.lower_atol,
+            frame="h",
+        )
+        qiskit_expval = float(qiskit_obs.coeffs[~qiskit_obs.paulis.x.any(axis=1)].sum())
+        t2 = time.perf_counter()
+
+        if step_idx > 0:
+            runtime.append(t2 - t1)
+        expvals.append(qiskit_expval)
+        num_terms.append(len(qiskit_obs))
+        memory.append(sampler.peak_mb())
+
+update_results(
+    LABEL, runtime=runtime, memory=memory, expvals=expvals, num_terms=num_terms
+)
diff --git a/benches/third_party/pauli_prop/settings.json b/benches/third_party/pauli_prop/settings.json
index b847f655..5e46f86d 100644
--- a/benches/third_party/pauli_prop/settings.json
+++ b/benches/third_party/pauli_prop/settings.json
@@ -1,12 +1,12 @@
 {
-    "nx": 6,
-    "ny": 6,
+    "nx": 12,
+    "ny": 12,
     "hx": 1.0,
     "hz": 1.0,
     "j": 1.5,
     "dt": 0.05,
     "step_min": 0,
-    "step_max": 40,
+    "step_max": 55,
     "step_size": 2,
     "lower_atol": 1e-6,
     "cutoff": null,
diff --git a/docs/public/benchmarks/majorana_results.png b/docs/public/benchmarks/majorana_results.png
index f02fa45a..56ea268f 100644
Binary files a/docs/public/benchmarks/majorana_results.png and b/docs/public/benchmarks/majorana_results.png differ
diff --git a/docs/public/benchmarks/pauli_results.png b/docs/public/benchmarks/pauli_results.png
index 54f0e665..b78f7fae 100644
Binary files a/docs/public/benchmarks/pauli_results.png and b/docs/public/benchmarks/pauli_results.png differ