From d9a7b2c2d04cd297bb48581558774747497fdab8 Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Thu, 4 Jun 2026 20:27:14 -0400
Subject: [PATCH 01/21] Add closed-system dynamics propagator

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/__init__.py        |  1 +
 python/cudaq/contrib/propagator.py      | 86 +++++++++++++++++++++++++
 python/tests/contrib/test_propagator.py | 51 +++++++++++++++
 3 files changed, 138 insertions(+)
 create mode 100644 python/cudaq/contrib/propagator.py
 create mode 100644 python/tests/contrib/test_propagator.py

diff --git a/python/cudaq/contrib/__init__.py b/python/cudaq/contrib/__init__.py
index 968595cc7ac..1aa337465b8 100644
--- a/python/cudaq/contrib/__init__.py
+++ b/python/cudaq/contrib/__init__.py
@@ -7,3 +7,4 @@
 # ============================================================================ #
 
 from .qiskit_convert import from_qiskit, from_qasm
+from .propagator import propagator
diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
new file mode 100644
index 00000000000..d361fe2ee6a
--- /dev/null
+++ b/python/cudaq/contrib/propagator.py
@@ -0,0 +1,86 @@
+from __future__ import annotations
+
+from collections.abc import Mapping
+from typing import Optional
+
+import numpy as np
+
+
+def _total_dimension(dimensions: Mapping[int, int]) -> int:
+    dimension = 1
+    for local_dimension in dimensions.values():
+        dimension *= local_dimension
+    return dimension
+
+
+def _identity_state(dimension: int):
+    import cudaq
+
+    identity = np.eye(dimension, dtype=np.complex128).reshape(-1)
+    return cudaq.State.from_data(identity)
+
+
+def _state_to_matrix(state, dimension: int) -> np.ndarray:
+    data = np.array(state).reshape(-1)
+    expected_size = dimension * dimension
+
+    if data.size != expected_size:
+        raise RuntimeError(
+            "Expected final propagator state with size "
+            f"{expected_size}, got {data.size}.")
+
+    return data.reshape((dimension, dimension))
+
+
+def _generator(hamiltonian):
+    import cudaq
+
+    generator = cudaq.SuperOperator()
+    generator += cudaq.SuperOperator.left_multiply(-1j * hamiltonian)
+    return generator
+
+
+def propagator(
+    hamiltonian,
+    dimensions: Mapping[int, int],
+    schedule,
+    *,
+    integrator=None,
+    max_batch_size: Optional[int] = None,
+):
+    """Compute a closed-system dynamics propagator.
+
+    Computes the matrix U satisfying
+
+        dU/dt = -i H(t) U,  U(t_initial) = I,
+
+    by evolving the identity matrix under the left-multiplication
+    superoperator generated by ``-i H(t)``.
+
+    Args:
+        hamiltonian: CUDA-Q operator H(t).
+        dimensions: Mapping from degree-of-freedom index to local dimension.
+        schedule: CUDA-Q dynamics schedule.
+        integrator: Optional dynamics integrator.
+        max_batch_size: Optional maximum batch size for the dynamics backend.
+
+    Returns:
+        Dense complex NumPy array with shape ``(dim, dim)``.
+    """
+    import cudaq
+
+    dimension = _total_dimension(dimensions)
+
+    result = cudaq.evolve(
+        _generator(hamiltonian),
+        dimensions,
+        schedule,
+        _identity_state(dimension),
+        collapse_operators=[],
+        observables=[],
+        store_intermediate_results=cudaq.IntermediateResultSave.NONE,
+        integrator=integrator,
+        max_batch_size=max_batch_size,
+    )
+
+    return _state_to_matrix(result.final_state(), dimension)
diff --git a/python/tests/contrib/test_propagator.py b/python/tests/contrib/test_propagator.py
new file mode 100644
index 00000000000..1d30f542e72
--- /dev/null
+++ b/python/tests/contrib/test_propagator.py
@@ -0,0 +1,51 @@
+import numpy as np
+import pytest
+
+import cudaq
+import cudaq.contrib
+from cudaq import Schedule, spin
+
+scipy_linalg = pytest.importorskip("scipy.linalg")
+
+if cudaq.num_available_gpus() == 0:
+    pytest.skip("cudaq.contrib.propagator requires the dynamics backend",
+                allow_module_level=True)
+
+
+@pytest.fixture(autouse=True)
+def dynamics_target():
+    cudaq.set_target("dynamics")
+    yield
+    cudaq.reset_target()
+
+
+def test_propagator_constant_z_hamiltonian():
+    omega = 0.3
+    t_final = 1.2
+
+    hamiltonian = 0.5 * omega * spin.z(0)
+    schedule = Schedule(np.linspace(0.0, t_final, 21), ["t"])
+
+    computed = cudaq.contrib.propagator(hamiltonian, {0: 2}, schedule)
+
+    z_matrix = np.array([[1.0, 0.0], [0.0, -1.0]], dtype=np.complex128)
+    expected = scipy_linalg.expm(-1j * 0.5 * omega * z_matrix * t_final)
+
+    np.testing.assert_allclose(computed, expected, atol=1e-4)
+    np.testing.assert_allclose(computed.conj().T @ computed,
+                               np.eye(2),
+                               atol=1e-4)
+
+
+def test_propagator_constant_x_hamiltonian():
+    t_final = 0.4
+
+    hamiltonian = spin.x(0)
+    schedule = Schedule(np.linspace(0.0, t_final, 21), ["t"])
+
+    computed = cudaq.contrib.propagator(hamiltonian, {0: 2}, schedule)
+
+    x_matrix = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=np.complex128)
+    expected = scipy_linalg.expm(-1j * x_matrix * t_final)
+
+    np.testing.assert_allclose(computed, expected, atol=1e-4)

From 608c313c3747e6866a8b4eeeffc369f9fbc35fd6 Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Thu, 4 Jun 2026 20:44:29 -0400
Subject: [PATCH 02/21] Support batched and time-dependent propagators

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py      |  53 +++++++++--
 python/tests/contrib/test_propagator.py | 111 ++++++++++++++++++++++--
 2 files changed, 150 insertions(+), 14 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index d361fe2ee6a..697d73b5d28 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -1,6 +1,6 @@
 from __future__ import annotations
 
-from collections.abc import Mapping
+from collections.abc import Mapping, Sequence
 from typing import Optional
 
 import numpy as np
@@ -26,7 +26,7 @@ def _state_to_matrix(state, dimension: int) -> np.ndarray:
 
     if data.size != expected_size:
         raise RuntimeError(
-            "Expected final propagator state with size "
+            "Expected propagator state with size "
             f"{expected_size}, got {data.size}.")
 
     return data.reshape((dimension, dimension))
@@ -40,15 +40,27 @@ def _generator(hamiltonian):
     return generator
 
 
+def _extract_propagator(result, dimension: int,
+                        store_intermediate_results: bool):
+    if store_intermediate_results:
+        return [
+            _state_to_matrix(state, dimension)
+            for state in result.intermediate_states()
+        ]
+
+    return _state_to_matrix(result.final_state(), dimension)
+
+
 def propagator(
     hamiltonian,
     dimensions: Mapping[int, int],
     schedule,
     *,
+    store_intermediate_results: bool = False,
     integrator=None,
     max_batch_size: Optional[int] = None,
 ):
-    """Compute a closed-system dynamics propagator.
+    """Compute closed-system dynamics propagators.
 
     Computes the matrix U satisfying
 
@@ -58,29 +70,52 @@ def propagator(
     superoperator generated by ``-i H(t)``.
 
     Args:
-        hamiltonian: CUDA-Q operator H(t).
+        hamiltonian: CUDA-Q operator H(t), or a sequence of operators for
+            batched propagator computation.
         dimensions: Mapping from degree-of-freedom index to local dimension.
         schedule: CUDA-Q dynamics schedule.
+        store_intermediate_results: If True, return propagators at the
+            intermediate schedule points saved by the dynamics backend.
         integrator: Optional dynamics integrator.
         max_batch_size: Optional maximum batch size for the dynamics backend.
 
     Returns:
-        Dense complex NumPy array with shape ``(dim, dim)``.
+        For a single Hamiltonian, returns a dense complex NumPy array with shape
+        ``(dim, dim)``. If ``store_intermediate_results`` is True, returns a
+        list of dense matrices.
+
+        For a sequence of Hamiltonians, returns one such result per Hamiltonian.
     """
     import cudaq
 
     dimension = _total_dimension(dimensions)
+    is_batched = isinstance(hamiltonian, Sequence)
+
+    hamiltonians = list(hamiltonian) if is_batched else [hamiltonian]
+    generators = [_generator(h) for h in hamiltonians]
+    initial_states = [_identity_state(dimension) for _ in generators]
+
+    save_mode = (cudaq.IntermediateResultSave.ALL
+                 if store_intermediate_results else
+                 cudaq.IntermediateResultSave.NONE)
 
     result = cudaq.evolve(
-        _generator(hamiltonian),
+        generators if is_batched else generators[0],
         dimensions,
         schedule,
-        _identity_state(dimension),
+        initial_states if is_batched else initial_states[0],
         collapse_operators=[],
         observables=[],
-        store_intermediate_results=cudaq.IntermediateResultSave.NONE,
+        store_intermediate_results=save_mode,
         integrator=integrator,
         max_batch_size=max_batch_size,
     )
 
-    return _state_to_matrix(result.final_state(), dimension)
+    if is_batched:
+        return [
+            _extract_propagator(single_result, dimension,
+                                store_intermediate_results)
+            for single_result in result
+        ]
+
+    return _extract_propagator(result, dimension, store_intermediate_results)
diff --git a/python/tests/contrib/test_propagator.py b/python/tests/contrib/test_propagator.py
index 1d30f542e72..d04b9050315 100644
--- a/python/tests/contrib/test_propagator.py
+++ b/python/tests/contrib/test_propagator.py
@@ -3,8 +3,9 @@
 
 import cudaq
 import cudaq.contrib
-from cudaq import Schedule, spin
+from cudaq import ScalarOperator, Schedule, spin
 
+scipy_integrate = pytest.importorskip("scipy.integrate")
 scipy_linalg = pytest.importorskip("scipy.linalg")
 
 if cudaq.num_available_gpus() == 0:
@@ -19,6 +20,31 @@ def dynamics_target():
     cudaq.reset_target()
 
 
+_X = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=np.complex128)
+_Z = np.array([[1.0, 0.0], [0.0, -1.0]], dtype=np.complex128)
+
+
+def _solve_propagator_reference(hamiltonian_matrix, t_final):
+    dimension = hamiltonian_matrix(0.0).shape[0]
+
+    def rhs(t, flattened_u):
+        u = flattened_u.reshape((dimension, dimension))
+        du = -1j * hamiltonian_matrix(t) @ u
+        return du.reshape(-1)
+
+    initial = np.eye(dimension, dtype=np.complex128).reshape(-1)
+    solution = scipy_integrate.solve_ivp(
+        rhs,
+        (0.0, t_final),
+        initial,
+        rtol=1e-10,
+        atol=1e-12,
+    )
+
+    assert solution.success
+    return solution.y[:, -1].reshape((dimension, dimension))
+
+
 def test_propagator_constant_z_hamiltonian():
     omega = 0.3
     t_final = 1.2
@@ -28,8 +54,7 @@ def test_propagator_constant_z_hamiltonian():
 
     computed = cudaq.contrib.propagator(hamiltonian, {0: 2}, schedule)
 
-    z_matrix = np.array([[1.0, 0.0], [0.0, -1.0]], dtype=np.complex128)
-    expected = scipy_linalg.expm(-1j * 0.5 * omega * z_matrix * t_final)
+    expected = scipy_linalg.expm(-1j * 0.5 * omega * _Z * t_final)
 
     np.testing.assert_allclose(computed, expected, atol=1e-4)
     np.testing.assert_allclose(computed.conj().T @ computed,
@@ -44,8 +69,84 @@ def test_propagator_constant_x_hamiltonian():
     schedule = Schedule(np.linspace(0.0, t_final, 21), ["t"])
 
     computed = cudaq.contrib.propagator(hamiltonian, {0: 2}, schedule)
+    expected = scipy_linalg.expm(-1j * _X * t_final)
+
+    np.testing.assert_allclose(computed, expected, atol=1e-4)
+
+
+def test_propagator_batched_hamiltonians():
+    t_final = 0.7
+    omegas = [0.2, 0.5]
+
+    hamiltonians = [0.5 * omega * spin.z(0) for omega in omegas]
+    schedule = Schedule(np.linspace(0.0, t_final, 21), ["t"])
+
+    computed = cudaq.contrib.propagator(
+        hamiltonians,
+        {0: 2},
+        schedule,
+        max_batch_size=2,
+    )
+
+    assert len(computed) == len(omegas)
+
+    for propagator, omega in zip(computed, omegas):
+        expected = scipy_linalg.expm(-1j * 0.5 * omega * _Z * t_final)
+        np.testing.assert_allclose(propagator, expected, atol=1e-4)
+
+
+def test_propagator_intermediate_results():
+    omega = 0.3
+    steps = np.linspace(0.0, 1.2, 7)
+
+    hamiltonian = 0.5 * omega * spin.z(0)
+    schedule = Schedule(steps, ["t"])
+
+    computed = cudaq.contrib.propagator(
+        hamiltonian,
+        {0: 2},
+        schedule,
+        store_intermediate_results=True,
+    )
+
+    assert len(computed) == len(steps)
+
+    for propagator, t in zip(computed, steps):
+        expected = scipy_linalg.expm(-1j * 0.5 * omega * _Z * t)
+        np.testing.assert_allclose(propagator, expected, atol=1e-4)
+
+
+def test_propagator_commuting_time_dependent_hamiltonian():
+    omega0 = 0.2
+    omega1 = 0.4
+    t_final = 1.1
+
+    hamiltonian = (0.5 * omega0 * spin.z(0) +
+                   0.5 * omega1 * ScalarOperator(lambda t: t) * spin.z(0))
+    schedule = Schedule(np.linspace(0.0, t_final, 31), ["t"])
+
+    computed = cudaq.contrib.propagator(hamiltonian, {0: 2}, schedule)
+
+    phase_integral = omega0 * t_final + 0.5 * omega1 * t_final**2
+    expected = scipy_linalg.expm(-1j * 0.5 * phase_integral * _Z)
+
+    np.testing.assert_allclose(computed, expected, atol=1e-4)
+
+
+def test_propagator_noncommuting_time_dependent_hamiltonian():
+    ax = 0.2
+    bz = 0.35
+    t_final = 0.8
+
+    hamiltonian = (ax * ScalarOperator(lambda t: t) * spin.x(0) +
+                   bz * ScalarOperator(lambda t: 1.0 - t) * spin.z(0))
+    schedule = Schedule(np.linspace(0.0, t_final, 41), ["t"])
+
+    computed = cudaq.contrib.propagator(hamiltonian, {0: 2}, schedule)
+
+    def hamiltonian_matrix(t):
+        return ax * t * _X + bz * (1.0 - t) * _Z
 
-    x_matrix = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=np.complex128)
-    expected = scipy_linalg.expm(-1j * x_matrix * t_final)
+    expected = _solve_propagator_reference(hamiltonian_matrix, t_final)
 
     np.testing.assert_allclose(computed, expected, atol=1e-4)

From f527c3dff72b5da11ce4ae042b4c4ad47160b1cc Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Thu, 4 Jun 2026 21:42:11 -0400
Subject: [PATCH 03/21] Support open-system propagators

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py      | 77 ++++++++++++++++---------
 python/tests/contrib/test_propagator.py | 53 +++++++++++++++++
 2 files changed, 103 insertions(+), 27 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index 697d73b5d28..044e3994133 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -32,7 +32,7 @@ def _state_to_matrix(state, dimension: int) -> np.ndarray:
     return data.reshape((dimension, dimension))
 
 
-def _generator(hamiltonian):
+def _closed_system_generator(hamiltonian):
     import cudaq
 
     generator = cudaq.SuperOperator()
@@ -40,6 +40,27 @@ def _generator(hamiltonian):
     return generator
 
 
+def _open_system_generator(hamiltonian, collapse_operators):
+    import cudaq
+
+    generator = cudaq.SuperOperator()
+    generator += cudaq.SuperOperator.left_multiply(-1j * hamiltonian)
+    generator += cudaq.SuperOperator.right_multiply(1j * hamiltonian)
+
+    for collapse_operator in collapse_operators:
+        collapse_operator_dagger = collapse_operator.dagger()
+        collapse_operator_product = collapse_operator_dagger * collapse_operator
+
+        generator += cudaq.SuperOperator.left_right_multiply(
+            collapse_operator, collapse_operator_dagger)
+        generator += cudaq.SuperOperator.left_multiply(
+            -0.5 * collapse_operator_product)
+        generator += cudaq.SuperOperator.right_multiply(
+            -0.5 * collapse_operator_product)
+
+    return generator
+
+
 def _extract_propagator(result, dimension: int,
                         store_intermediate_results: bool):
     if store_intermediate_results:
@@ -56,44 +77,45 @@ def propagator(
     dimensions: Mapping[int, int],
     schedule,
     *,
+    collapse_operators=None,
     store_intermediate_results: bool = False,
     integrator=None,
     max_batch_size: Optional[int] = None,
 ):
-    """Compute closed-system dynamics propagators.
-
-    Computes the matrix U satisfying
+    """Compute dynamics propagators.
 
-        dU/dt = -i H(t) U,  U(t_initial) = I,
+    For closed-system dynamics, computes the matrix U satisfying
 
-    by evolving the identity matrix under the left-multiplication
-    superoperator generated by ``-i H(t)``.
+        dU/dt = -i H(t) U,  U(t_initial) = I.
 
-    Args:
-        hamiltonian: CUDA-Q operator H(t), or a sequence of operators for
-            batched propagator computation.
-        dimensions: Mapping from degree-of-freedom index to local dimension.
-        schedule: CUDA-Q dynamics schedule.
-        store_intermediate_results: If True, return propagators at the
-            intermediate schedule points saved by the dynamics backend.
-        integrator: Optional dynamics integrator.
-        max_batch_size: Optional maximum batch size for the dynamics backend.
+    For open-system dynamics with collapse operators, computes the
+    superoperator propagator S satisfying
 
-    Returns:
-        For a single Hamiltonian, returns a dense complex NumPy array with shape
-        ``(dim, dim)``. If ``store_intermediate_results`` is True, returns a
-        list of dense matrices.
+        d vec(rho)/dt = L(t) vec(rho),  S(t_initial) = I,
 
-        For a sequence of Hamiltonians, returns one such result per Hamiltonian.
+    where L(t) is the Lindblad generator.
     """
     import cudaq
 
-    dimension = _total_dimension(dimensions)
-    is_batched = isinstance(hamiltonian, Sequence)
+    collapse_operators = [] if collapse_operators is None else list(
+        collapse_operators)
+    open_system = len(collapse_operators) > 0
+
+    system_dimension = _total_dimension(dimensions)
+    propagator_dimension = (system_dimension * system_dimension
+                            if open_system else system_dimension)
 
+    is_batched = isinstance(hamiltonian, Sequence)
     hamiltonians = list(hamiltonian) if is_batched else [hamiltonian]
-    generators = [_generator(h) for h in hamiltonians]
-    initial_states = [_identity_state(dimension) for _ in generators]
+
+    if open_system:
+        generators = [
+            _open_system_generator(h, collapse_operators) for h in hamiltonians
+        ]
+    else:
+        generators = [_closed_system_generator(h) for h in hamiltonians]
+
+    initial_states = [_identity_state(propagator_dimension) for _ in generators]
 
     save_mode = (cudaq.IntermediateResultSave.ALL
                  if store_intermediate_results else
@@ -113,9 +135,10 @@ def propagator(
 
     if is_batched:
         return [
-            _extract_propagator(single_result, dimension,
+            _extract_propagator(single_result, propagator_dimension,
                                 store_intermediate_results)
             for single_result in result
         ]
 
-    return _extract_propagator(result, dimension, store_intermediate_results)
+    return _extract_propagator(result, propagator_dimension,
+                               store_intermediate_results)
diff --git a/python/tests/contrib/test_propagator.py b/python/tests/contrib/test_propagator.py
index d04b9050315..029675a7d14 100644
--- a/python/tests/contrib/test_propagator.py
+++ b/python/tests/contrib/test_propagator.py
@@ -150,3 +150,56 @@ def hamiltonian_matrix(t):
     expected = _solve_propagator_reference(hamiltonian_matrix, t_final)
 
     np.testing.assert_allclose(computed, expected, atol=1e-4)
+
+
+def _amplitude_damping_liouvillian(gamma):
+    sigma_minus = np.array([[0.0, 0.0], [1.0, 0.0]], dtype=np.complex128)
+    sigma_plus = sigma_minus.conj().T
+    identity = np.eye(2, dtype=np.complex128)
+
+    collapse_product = sigma_plus @ sigma_minus
+
+    return (
+        gamma * np.kron(sigma_minus.conj(), sigma_minus)
+        - 0.5 * gamma * np.kron(identity, collapse_product)
+        - 0.5 * gamma * np.kron(collapse_product.T, identity)
+    )
+
+
+def test_open_system_propagator_amplitude_damping():
+    gamma = 0.2
+    t_final = 0.7
+
+    hamiltonian = 0.0 * spin.z(0)
+    collapse_operators = [np.sqrt(gamma) * spin.minus(0)]
+    schedule = Schedule(np.linspace(0.0, t_final, 21), ["t"])
+
+    computed = cudaq.contrib.propagator(
+        hamiltonian,
+        {0: 2},
+        schedule,
+        collapse_operators=collapse_operators,
+    )
+
+    expected = scipy_linalg.expm(_amplitude_damping_liouvillian(gamma) *
+                                 t_final)
+
+    np.testing.assert_allclose(computed, expected, atol=1e-4)
+
+
+def test_open_system_propagator_shape():
+    gamma = 0.1
+    t_final = 0.3
+
+    hamiltonian = 0.0 * spin.z(0)
+    collapse_operators = [np.sqrt(gamma) * spin.minus(0)]
+    schedule = Schedule(np.linspace(0.0, t_final, 11), ["t"])
+
+    computed = cudaq.contrib.propagator(
+        hamiltonian,
+        {0: 2},
+        schedule,
+        collapse_operators=collapse_operators,
+    )
+
+    assert computed.shape == (4, 4)

From c5a247ecd0b9ac8d4c90d2edc1c6ace1c8e27249 Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Thu, 4 Jun 2026 22:24:09 -0400
Subject: [PATCH 04/21] Fix propagator docstring spelling

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py | 49 ++++++++++++++++++++----------
 1 file changed, 33 insertions(+), 16 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index 044e3994133..34f6c7cf807 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -25,9 +25,8 @@ def _state_to_matrix(state, dimension: int) -> np.ndarray:
     expected_size = dimension * dimension
 
     if data.size != expected_size:
-        raise RuntimeError(
-            "Expected propagator state with size "
-            f"{expected_size}, got {data.size}.")
+        raise RuntimeError("Expected propagator state with size "
+                           f"{expected_size}, got {data.size}.")
 
     return data.reshape((dimension, dimension))
 
@@ -84,16 +83,35 @@ def propagator(
 ):
     """Compute dynamics propagators.
 
-    For closed-system dynamics, computes the matrix U satisfying
-
-        dU/dt = -i H(t) U,  U(t_initial) = I.
-
-    For open-system dynamics with collapse operators, computes the
-    superoperator propagator S satisfying
-
-        d vec(rho)/dt = L(t) vec(rho),  S(t_initial) = I,
-
-    where L(t) is the Lindblad generator.
+    For closed-system dynamics, computes the matrix U satisfying the
+    Schrodinger-picture propagator equation with initial condition U(t0) = I.
+
+    For open-system dynamics with collapse operators, computes the Lindblad
+    map S with initial condition S(t0) = I. This map acts on vectorized density
+    matrices and propagates rho(t0) to rho(t).
+
+    Args:
+        hamiltonian: CUDA-Q operator H(t), or a sequence of operators for
+            batched propagator computation.
+        dimensions: Mapping from degree-of-freedom index to local dimension.
+        schedule: CUDA-Q dynamics schedule.
+        collapse_operators: Optional sequence of Lindblad collapse operators.
+            If provided, the helper returns the Lindblad map.
+        store_intermediate_results: If True, return propagators at the
+            intermediate schedule points saved by the dynamics backend.
+        integrator: Optional dynamics integrator.
+        max_batch_size: Optional maximum batch size for the dynamics backend.
+
+    Returns:
+        For closed-system dynamics, returns a dense complex NumPy array with
+        shape ``(dim, dim)``.
+
+        For open-system dynamics, returns a dense complex NumPy array with
+        shape ``(dim**2, dim**2)``.
+
+        If ``store_intermediate_results`` is True, returns a list of dense
+        matrices. For a sequence of Hamiltonians, returns one such result per
+        Hamiltonian.
     """
     import cudaq
 
@@ -117,9 +135,8 @@ def propagator(
 
     initial_states = [_identity_state(propagator_dimension) for _ in generators]
 
-    save_mode = (cudaq.IntermediateResultSave.ALL
-                 if store_intermediate_results else
-                 cudaq.IntermediateResultSave.NONE)
+    save_mode = (cudaq.IntermediateResultSave.ALL if store_intermediate_results
+                 else cudaq.IntermediateResultSave.NONE)
 
     result = cudaq.evolve(
         generators if is_batched else generators[0],

From 5beddd50a6892b77c1972cf9f1a60d5516db322b Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Thu, 4 Jun 2026 22:40:55 -0400
Subject: [PATCH 05/21] Fix propagator spelling and formatting

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py      |  5 +++--
 python/tests/contrib/test_propagator.py | 12 +++++-------
 2 files changed, 8 insertions(+), 9 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index 34f6c7cf807..43e14ecbb58 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -87,8 +87,9 @@ def propagator(
     Schrodinger-picture propagator equation with initial condition U(t0) = I.
 
     For open-system dynamics with collapse operators, computes the Lindblad
-    map S with initial condition S(t0) = I. This map acts on vectorized density
-    matrices and propagates rho(t0) to rho(t).
+    map S with initial condition S(t0) = I. This map acts on density
+    matrices after matrix-to-vector reshaping and propagates rho(t0) to
+    rho(t).
 
     Args:
         hamiltonian: CUDA-Q operator H(t), or a sequence of operators for
diff --git a/python/tests/contrib/test_propagator.py b/python/tests/contrib/test_propagator.py
index 029675a7d14..f2be078a84a 100644
--- a/python/tests/contrib/test_propagator.py
+++ b/python/tests/contrib/test_propagator.py
@@ -159,11 +159,9 @@ def _amplitude_damping_liouvillian(gamma):
 
     collapse_product = sigma_plus @ sigma_minus
 
-    return (
-        gamma * np.kron(sigma_minus.conj(), sigma_minus)
-        - 0.5 * gamma * np.kron(identity, collapse_product)
-        - 0.5 * gamma * np.kron(collapse_product.T, identity)
-    )
+    return (gamma * np.kron(sigma_minus.conj(), sigma_minus) -
+            0.5 * gamma * np.kron(identity, collapse_product) -
+            0.5 * gamma * np.kron(collapse_product.T, identity))
 
 
 def test_open_system_propagator_amplitude_damping():
@@ -181,8 +179,8 @@ def test_open_system_propagator_amplitude_damping():
         collapse_operators=collapse_operators,
     )
 
-    expected = scipy_linalg.expm(_amplitude_damping_liouvillian(gamma) *
-                                 t_final)
+    expected = scipy_linalg.expm(
+        _amplitude_damping_liouvillian(gamma) * t_final)
 
     np.testing.assert_allclose(computed, expected, atol=1e-4)
 

From 1a406fb98672dffcdc18aff89626ee87a1f8c9e9 Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Thu, 4 Jun 2026 22:49:31 -0400
Subject: [PATCH 06/21] Fix propagator matrix convention and collapse adjoint
 operators

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py      | 27 ++++++++++++++++++++-----
 python/tests/contrib/test_propagator.py |  2 ++
 2 files changed, 24 insertions(+), 5 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index 43e14ecbb58..037abdc7da3 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -28,7 +28,7 @@ def _state_to_matrix(state, dimension: int) -> np.ndarray:
         raise RuntimeError("Expected propagator state with size "
                            f"{expected_size}, got {data.size}.")
 
-    return data.reshape((dimension, dimension))
+    return data.reshape((dimension, dimension)).T
 
 
 def _closed_system_generator(hamiltonian):
@@ -39,15 +39,16 @@ def _closed_system_generator(hamiltonian):
     return generator
 
 
-def _open_system_generator(hamiltonian, collapse_operators):
+def _open_system_generator(hamiltonian, collapse_operators,
+                           collapse_operator_adjoint_ops):
     import cudaq
 
     generator = cudaq.SuperOperator()
     generator += cudaq.SuperOperator.left_multiply(-1j * hamiltonian)
     generator += cudaq.SuperOperator.right_multiply(1j * hamiltonian)
 
-    for collapse_operator in collapse_operators:
-        collapse_operator_dagger = collapse_operator.dagger()
+    for collapse_operator, collapse_operator_dagger in zip(
+            collapse_operators, collapse_operator_adjoint_ops):
         collapse_operator_product = collapse_operator_dagger * collapse_operator
 
         generator += cudaq.SuperOperator.left_right_multiply(
@@ -77,6 +78,7 @@ def propagator(
     schedule,
     *,
     collapse_operators=None,
+    collapse_operator_adjoint_ops=None,
     store_intermediate_results: bool = False,
     integrator=None,
     max_batch_size: Optional[int] = None,
@@ -98,6 +100,8 @@ def propagator(
         schedule: CUDA-Q dynamics schedule.
         collapse_operators: Optional sequence of Lindblad collapse operators.
             If provided, the helper returns the Lindblad map.
+        collapse_operator_adjoint_ops: Optional sequence containing the adjoint
+            of each collapse operator.
         store_intermediate_results: If True, return propagators at the
             intermediate schedule points saved by the dynamics backend.
         integrator: Optional dynamics integrator.
@@ -118,7 +122,18 @@ def propagator(
 
     collapse_operators = [] if collapse_operators is None else list(
         collapse_operators)
+    collapse_operator_adjoint_ops = ([] if collapse_operator_adjoint_ops is None
+                                     else list(collapse_operator_adjoint_ops))
+
+    if collapse_operator_adjoint_ops and len(
+            collapse_operator_adjoint_ops) != len(collapse_operators):
+        raise ValueError("collapse_operator_adjoint_ops must have the same "
+                         "length as collapse_operators.")
+
     open_system = len(collapse_operators) > 0
+    if open_system and not collapse_operator_adjoint_ops:
+        raise ValueError("collapse_operator_adjoint_ops must be provided for "
+                         "open-system propagators.")
 
     system_dimension = _total_dimension(dimensions)
     propagator_dimension = (system_dimension * system_dimension
@@ -129,7 +144,9 @@ def propagator(
 
     if open_system:
         generators = [
-            _open_system_generator(h, collapse_operators) for h in hamiltonians
+            _open_system_generator(h, collapse_operators,
+                                   collapse_operator_adjoint_ops)
+            for h in hamiltonians
         ]
     else:
         generators = [_closed_system_generator(h) for h in hamiltonians]
diff --git a/python/tests/contrib/test_propagator.py b/python/tests/contrib/test_propagator.py
index f2be078a84a..c7227d78ca9 100644
--- a/python/tests/contrib/test_propagator.py
+++ b/python/tests/contrib/test_propagator.py
@@ -177,6 +177,7 @@ def test_open_system_propagator_amplitude_damping():
         {0: 2},
         schedule,
         collapse_operators=collapse_operators,
+        collapse_operator_adjoints=[np.sqrt(gamma) * spin.plus(0)],
     )
 
     expected = scipy_linalg.expm(
@@ -198,6 +199,7 @@ def test_open_system_propagator_shape():
         {0: 2},
         schedule,
         collapse_operators=collapse_operators,
+        collapse_operator_adjoints=[np.sqrt(gamma) * spin.plus(0)],
     )
 
     assert computed.shape == (4, 4)

From 336311ba5973425fc3f77ad2acb952236fcb623e Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Thu, 4 Jun 2026 22:53:47 -0400
Subject: [PATCH 07/21] Add propagator license headers

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py      | 8 ++++++++
 python/tests/contrib/test_propagator.py | 8 ++++++++
 2 files changed, 16 insertions(+)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index 037abdc7da3..800f84ecfb0 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -1,3 +1,11 @@
+# ============================================================================ #
+# Copyright (c) 2022 - 2026 NVIDIA Corporation & Affiliates.                   #
+# All rights reserved.                                                         #
+#                                                                              #
+# This source code and the accompanying materials are made available under     #
+# the terms of the Apache License 2.0 which accompanies this distribution.     #
+# ============================================================================ #
+
 from __future__ import annotations
 
 from collections.abc import Mapping, Sequence
diff --git a/python/tests/contrib/test_propagator.py b/python/tests/contrib/test_propagator.py
index c7227d78ca9..93c36da3df6 100644
--- a/python/tests/contrib/test_propagator.py
+++ b/python/tests/contrib/test_propagator.py
@@ -1,3 +1,11 @@
+# ============================================================================ #
+# Copyright (c) 2022 - 2026 NVIDIA Corporation & Affiliates.                   #
+# All rights reserved.                                                         #
+#                                                                              #
+# This source code and the accompanying materials are made available under     #
+# the terms of the Apache License 2.0 which accompanies this distribution.     #
+# ============================================================================ #
+
 import numpy as np
 import pytest
 

From 83573621abb517c63108fb2ffa4fc23e57240725 Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Thu, 4 Jun 2026 23:15:08 -0400
Subject: [PATCH 08/21] Assemble open-system propagators from basis evolution

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py      | 110 ++++++++++++++++++++----
 python/tests/contrib/test_propagator.py |   4 +-
 2 files changed, 96 insertions(+), 18 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index 800f84ecfb0..a73fa61c373 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -28,8 +28,20 @@ def _identity_state(dimension: int):
     return cudaq.State.from_data(identity)
 
 
+def _basis_state(index: int, dimension: int):
+    import cudaq
+
+    data = np.zeros(dimension, dtype=np.complex128)
+    data[index] = 1.0
+    return cudaq.State.from_data(data)
+
+
+def _state_data(state) -> np.ndarray:
+    return np.array(state).reshape(-1)
+
+
 def _state_to_matrix(state, dimension: int) -> np.ndarray:
-    data = np.array(state).reshape(-1)
+    data = _state_data(state)
     expected_size = dimension * dimension
 
     if data.size != expected_size:
@@ -55,12 +67,13 @@ def _open_system_generator(hamiltonian, collapse_operators,
     generator += cudaq.SuperOperator.left_multiply(-1j * hamiltonian)
     generator += cudaq.SuperOperator.right_multiply(1j * hamiltonian)
 
-    for collapse_operator, collapse_operator_dagger in zip(
+    for collapse_operator, collapse_operator_adjoint in zip(
             collapse_operators, collapse_operator_adjoint_ops):
-        collapse_operator_product = collapse_operator_dagger * collapse_operator
+        collapse_operator_product = (collapse_operator_adjoint *
+                                     collapse_operator)
 
         generator += cudaq.SuperOperator.left_right_multiply(
-            collapse_operator, collapse_operator_dagger)
+            collapse_operator, collapse_operator_adjoint)
         generator += cudaq.SuperOperator.left_multiply(
             -0.5 * collapse_operator_product)
         generator += cudaq.SuperOperator.right_multiply(
@@ -69,8 +82,8 @@ def _open_system_generator(hamiltonian, collapse_operators,
     return generator
 
 
-def _extract_propagator(result, dimension: int,
-                        store_intermediate_results: bool):
+def _extract_closed_system_propagator(result, dimension: int,
+                                      store_intermediate_results: bool):
     if store_intermediate_results:
         return [
             _state_to_matrix(state, dimension)
@@ -80,6 +93,61 @@ def _extract_propagator(result, dimension: int,
     return _state_to_matrix(result.final_state(), dimension)
 
 
+def _stack_columns(states) -> np.ndarray:
+    return np.column_stack([_state_data(state) for state in states])
+
+
+def _extract_open_system_propagator(results, store_intermediate_results: bool):
+    if store_intermediate_results:
+        intermediate_states = [
+            result.intermediate_states() for result in results
+        ]
+        if not intermediate_states:
+            return []
+
+        num_intermediate_states = len(intermediate_states[0])
+        return [
+            _stack_columns(
+                [column_states[index]
+                 for column_states in intermediate_states])
+            for index in range(num_intermediate_states)
+        ]
+
+    return _stack_columns([result.final_state() for result in results])
+
+
+def _evolve_open_system_propagator(
+    generator,
+    dimensions,
+    schedule,
+    liouville_dimension: int,
+    store_intermediate_results: bool,
+    integrator,
+    max_batch_size: Optional[int],
+):
+    import cudaq
+
+    save_mode = (cudaq.IntermediateResultSave.ALL if store_intermediate_results
+                 else cudaq.IntermediateResultSave.NONE)
+
+    results = []
+    for column_index in range(liouville_dimension):
+        results.append(
+            cudaq.evolve(
+                generator,
+                dimensions,
+                schedule,
+                _basis_state(column_index, liouville_dimension),
+                collapse_operators=[],
+                observables=[],
+                store_intermediate_results=save_mode,
+                integrator=integrator,
+                max_batch_size=max_batch_size,
+            ))
+
+    return _extract_open_system_propagator(results, store_intermediate_results)
+
+
 def propagator(
     hamiltonian,
     dimensions: Mapping[int, int],
@@ -109,7 +177,7 @@ def propagator(
         collapse_operators: Optional sequence of Lindblad collapse operators.
             If provided, the helper returns the Lindblad map.
         collapse_operator_adjoint_ops: Optional sequence containing the adjoint
-            of each collapse operator.
+            operator for each collapse operator.
         store_intermediate_results: If True, return propagators at the
             intermediate schedule points saved by the dynamics backend.
         integrator: Optional dynamics integrator.
@@ -144,22 +212,32 @@ def propagator(
                          "open-system propagators.")
 
     system_dimension = _total_dimension(dimensions)
-    propagator_dimension = (system_dimension * system_dimension
-                            if open_system else system_dimension)
 
     is_batched = isinstance(hamiltonian, Sequence)
     hamiltonians = list(hamiltonian) if is_batched else [hamiltonian]
 
     if open_system:
+        liouville_dimension = system_dimension * system_dimension
         generators = [
             _open_system_generator(h, collapse_operators,
                                    collapse_operator_adjoint_ops)
             for h in hamiltonians
         ]
-    else:
-        generators = [_closed_system_generator(h) for h in hamiltonians]
+        results = [
+            _evolve_open_system_propagator(
+                generator,
+                dimensions,
+                schedule,
+                liouville_dimension,
+                store_intermediate_results,
+                integrator,
+                max_batch_size,
+            ) for generator in generators
+        ]
+        return results if is_batched else results[0]
 
-    initial_states = [_identity_state(propagator_dimension) for _ in generators]
+    generators = [_closed_system_generator(h) for h in hamiltonians]
+    initial_states = [_identity_state(system_dimension) for _ in generators]
 
     save_mode = (cudaq.IntermediateResultSave.ALL if store_intermediate_results
                  else cudaq.IntermediateResultSave.NONE)
@@ -178,10 +256,10 @@ def propagator(
 
     if is_batched:
         return [
-            _extract_propagator(single_result, propagator_dimension,
-                                store_intermediate_results)
+            _extract_closed_system_propagator(single_result, system_dimension,
+                                              store_intermediate_results)
             for single_result in result
         ]
 
-    return _extract_propagator(result, propagator_dimension,
-                               store_intermediate_results)
+    return _extract_closed_system_propagator(result, system_dimension,
+                                             store_intermediate_results)
diff --git a/python/tests/contrib/test_propagator.py b/python/tests/contrib/test_propagator.py
index 93c36da3df6..2f0513ade48 100644
--- a/python/tests/contrib/test_propagator.py
+++ b/python/tests/contrib/test_propagator.py
@@ -185,7 +185,7 @@ def test_open_system_propagator_amplitude_damping():
         {0: 2},
         schedule,
         collapse_operators=collapse_operators,
-        collapse_operator_adjoints=[np.sqrt(gamma) * spin.plus(0)],
+        collapse_operator_adjoint_ops=[np.sqrt(gamma) * spin.plus(0)],
     )
 
     expected = scipy_linalg.expm(
@@ -207,7 +207,7 @@ def test_open_system_propagator_shape():
         {0: 2},
         schedule,
         collapse_operators=collapse_operators,
-        collapse_operator_adjoints=[np.sqrt(gamma) * spin.plus(0)],
+        collapse_operator_adjoint_ops=[np.sqrt(gamma) * spin.plus(0)],
     )
 
     assert computed.shape == (4, 4)

From 5386a725da883e6feb9b8dd7e306ae91bdc46009 Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Fri, 5 Jun 2026 00:05:30 -0400
Subject: [PATCH 09/21] Use direct Liouvillian for open-system propagators

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py | 179 ++++++++++++-----------------
 1 file changed, 72 insertions(+), 107 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index a73fa61c373..4242be02c4d 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -10,6 +10,7 @@
 
 from collections.abc import Mapping, Sequence
 from typing import Optional
+import uuid
 
 import numpy as np
 
@@ -28,20 +29,8 @@ def _identity_state(dimension: int):
     return cudaq.State.from_data(identity)
 
 
-def _basis_state(index: int, dimension: int):
-    import cudaq
-
-    data = np.zeros(dimension, dtype=np.complex128)
-    data[index] = 1.0
-    return cudaq.State.from_data(data)
-
-
-def _state_data(state) -> np.ndarray:
-    return np.array(state).reshape(-1)
-
-
 def _state_to_matrix(state, dimension: int) -> np.ndarray:
-    data = _state_data(state)
+    data = np.array(state).reshape(-1)
     expected_size = dimension * dimension
 
     if data.size != expected_size:
@@ -59,93 +48,75 @@ def _closed_system_generator(hamiltonian):
     return generator
 
 
-def _open_system_generator(hamiltonian, collapse_operators,
-                           collapse_operator_adjoint_ops):
-    import cudaq
+def _operator_matrix(operator, dimensions: Mapping[int, int], **parameters):
+    if parameters:
+        return np.asarray(operator.to_matrix(dimensions, **parameters),
+                          dtype=np.complex128)
+    return np.asarray(operator.to_matrix(dimensions), dtype=np.complex128)
 
-    generator = cudaq.SuperOperator()
-    generator += cudaq.SuperOperator.left_multiply(-1j * hamiltonian)
-    generator += cudaq.SuperOperator.right_multiply(1j * hamiltonian)
+
+def _open_system_matrix(hamiltonian, dimensions, collapse_operators,
+                        collapse_operator_adjoint_ops, **parameters):
+    system_dimension = _total_dimension(dimensions)
+    identity = np.eye(system_dimension, dtype=np.complex128)
+
+    hamiltonian_matrix = _operator_matrix(hamiltonian, dimensions, **parameters)
+
+    matrix = -1j * np.kron(identity, hamiltonian_matrix)
+    matrix += 1j * np.kron(hamiltonian_matrix.T, identity)
 
     for collapse_operator, collapse_operator_adjoint in zip(
             collapse_operators, collapse_operator_adjoint_ops):
-        collapse_operator_product = (collapse_operator_adjoint *
-                                     collapse_operator)
+        collapse_matrix = _operator_matrix(collapse_operator, dimensions,
+                                           **parameters)
+        collapse_adjoint_matrix = _operator_matrix(collapse_operator_adjoint,
+                                                   dimensions, **parameters)
+        collapse_product = collapse_adjoint_matrix @ collapse_matrix
 
-        generator += cudaq.SuperOperator.left_right_multiply(
-            collapse_operator, collapse_operator_adjoint)
-        generator += cudaq.SuperOperator.left_multiply(
-            -0.5 * collapse_operator_product)
-        generator += cudaq.SuperOperator.right_multiply(
-            -0.5 * collapse_operator_product)
+        matrix += np.kron(collapse_matrix.conj(), collapse_matrix)
+        matrix += -0.5 * np.kron(identity, collapse_product)
+        matrix += -0.5 * np.kron(collapse_product.T, identity)
 
-    return generator
+    return matrix
 
 
-def _extract_closed_system_propagator(result, dimension: int,
-                                      store_intermediate_results: bool):
-    if store_intermediate_results:
-        return [
-            _state_to_matrix(state, dimension)
-            for state in result.intermediate_states()
-        ]
-
-    return _state_to_matrix(result.final_state(), dimension)
+def _open_system_generator(hamiltonian, dimensions, collapse_operators,
+                           collapse_operator_adjoint_ops):
+    import cudaq
+    from cudaq import operators
 
+    system_dimension = _total_dimension(dimensions)
+    propagator_dimension = system_dimension * system_dimension
+    operator_id = f"propagator_open_system_{uuid.uuid4().hex}"
+
+    def create(t=0.0, dimension=propagator_dimension):
+        if dimension != propagator_dimension:
+            raise ValueError("Unexpected open-system propagator dimension.")
+        return _open_system_matrix(
+            hamiltonian,
+            dimensions,
+            collapse_operators,
+            collapse_operator_adjoint_ops,
+            t=t,
+        )
+
+    operators.define(operator_id, [propagator_dimension], create)
+    generator_matrix = operators.instantiate(operator_id, 0)
 
-def _stack_columns(states) -> np.ndarray:
-    return np.column_stack([_state_data(state) for state in states])
+    generator = cudaq.SuperOperator()
+    generator += cudaq.SuperOperator.left_multiply(generator_matrix)
+    return generator
 
 
-def _extract_open_system_propagator(results, store_intermediate_results: bool):
+def _extract_propagator(result, dimension: int,
+                        store_intermediate_results: bool):
     if store_intermediate_results:
-        intermediate_states = [
-            result.intermediate_states() for result in results
-        ]
-        if not intermediate_states:
-            return []
-
-        num_intermediate_states = len(intermediate_states[0])
         return [
-            _stack_columns(
-                [column_states[index]
-                 for column_states in intermediate_states])
-            for index in range(num_intermediate_states)
+            _state_to_matrix(state, dimension)
+            for state in result.intermediate_states()
         ]
 
-    return _stack_columns([result.final_state() for result in results])
-
-
-def _evolve_open_system_propagator(
-    generator,
-    dimensions,
-    schedule,
-    liouville_dimension: int,
-    store_intermediate_results: bool,
-    integrator,
-    max_batch_size: Optional[int],
-):
-    import cudaq
-
-    save_mode = (cudaq.IntermediateResultSave.ALL if store_intermediate_results
-                 else cudaq.IntermediateResultSave.NONE)
-
-    results = []
-    for column_index in range(liouville_dimension):
-        results.append(
-            cudaq.evolve(
-                generator,
-                dimensions,
-                schedule,
-                _basis_state(column_index, liouville_dimension),
-                collapse_operators=[],
-                observables=[],
-                store_intermediate_results=save_mode,
-                integrator=integrator,
-                max_batch_size=max_batch_size,
-            ))
-
-    return _extract_open_system_propagator(results, store_intermediate_results)
+    return _state_to_matrix(result.final_state(), dimension)
 
 
 def propagator(
@@ -212,39 +183,33 @@ def propagator(
                          "open-system propagators.")
 
     system_dimension = _total_dimension(dimensions)
-
-    is_batched = isinstance(hamiltonian, Sequence)
+    propagator_dimension = (system_dimension * system_dimension
+                            if open_system else system_dimension)
+    evolution_dimensions = ({
+        0: propagator_dimension
+    } if open_system else dimensions)
+
+    is_batched = isinstance(hamiltonian,
+                            Sequence) and not hasattr(hamiltonian, "to_matrix")
     hamiltonians = list(hamiltonian) if is_batched else [hamiltonian]
 
     if open_system:
-        liouville_dimension = system_dimension * system_dimension
         generators = [
-            _open_system_generator(h, collapse_operators,
+            _open_system_generator(h, dimensions, collapse_operators,
                                    collapse_operator_adjoint_ops)
             for h in hamiltonians
         ]
-        results = [
-            _evolve_open_system_propagator(
-                generator,
-                dimensions,
-                schedule,
-                liouville_dimension,
-                store_intermediate_results,
-                integrator,
-                max_batch_size,
-            ) for generator in generators
-        ]
-        return results if is_batched else results[0]
+    else:
+        generators = [_closed_system_generator(h) for h in hamiltonians]
 
-    generators = [_closed_system_generator(h) for h in hamiltonians]
-    initial_states = [_identity_state(system_dimension) for _ in generators]
+    initial_states = [_identity_state(propagator_dimension) for _ in generators]
 
     save_mode = (cudaq.IntermediateResultSave.ALL if store_intermediate_results
                  else cudaq.IntermediateResultSave.NONE)
 
     result = cudaq.evolve(
         generators if is_batched else generators[0],
-        dimensions,
+        evolution_dimensions,
         schedule,
         initial_states if is_batched else initial_states[0],
         collapse_operators=[],
@@ -256,10 +221,10 @@ def propagator(
 
     if is_batched:
         return [
-            _extract_closed_system_propagator(single_result, system_dimension,
-                                              store_intermediate_results)
+            _extract_propagator(single_result, propagator_dimension,
+                                store_intermediate_results)
             for single_result in result
         ]
 
-    return _extract_closed_system_propagator(result, system_dimension,
-                                             store_intermediate_results)
+    return _extract_propagator(result, propagator_dimension,
+                               store_intermediate_results)

From b8a8a1905b54a7fbf2d3aaaf49560d5d0b2d361a Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Fri, 5 Jun 2026 00:57:24 -0400
Subject: [PATCH 10/21] Move cudaq import to module scope

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py | 8 +-------
 1 file changed, 1 insertion(+), 7 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index 4242be02c4d..7796df331dc 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -12,6 +12,7 @@
 from typing import Optional
 import uuid
 
+import cudaq
 import numpy as np
 
 
@@ -23,8 +24,6 @@ def _total_dimension(dimensions: Mapping[int, int]) -> int:
 
 
 def _identity_state(dimension: int):
-    import cudaq
-
     identity = np.eye(dimension, dtype=np.complex128).reshape(-1)
     return cudaq.State.from_data(identity)
 
@@ -41,8 +40,6 @@ def _state_to_matrix(state, dimension: int) -> np.ndarray:
 
 
 def _closed_system_generator(hamiltonian):
-    import cudaq
-
     generator = cudaq.SuperOperator()
     generator += cudaq.SuperOperator.left_multiply(-1j * hamiltonian)
     return generator
@@ -82,7 +79,6 @@ def _open_system_matrix(hamiltonian, dimensions, collapse_operators,
 
 def _open_system_generator(hamiltonian, dimensions, collapse_operators,
                            collapse_operator_adjoint_ops):
-    import cudaq
     from cudaq import operators
 
     system_dimension = _total_dimension(dimensions)
@@ -165,8 +161,6 @@ def propagator(
         matrices. For a sequence of Hamiltonians, returns one such result per
         Hamiltonian.
     """
-    import cudaq
-
     collapse_operators = [] if collapse_operators is None else list(
         collapse_operators)
     collapse_operator_adjoint_ops = ([] if collapse_operator_adjoint_ops is None

From a99d137bb0dc925ca8195121c2df345ffb4df261 Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Fri, 5 Jun 2026 01:02:10 -0400
Subject: [PATCH 11/21] Add multi-degree propagator test

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/tests/contrib/test_propagator.py | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/python/tests/contrib/test_propagator.py b/python/tests/contrib/test_propagator.py
index 2f0513ade48..f145b4820ab 100644
--- a/python/tests/contrib/test_propagator.py
+++ b/python/tests/contrib/test_propagator.py
@@ -211,3 +211,18 @@ def test_open_system_propagator_shape():
     )
 
     assert computed.shape == (4, 4)
+
+
+@pytest.mark.skipif(
+    not _has_dynamics_backend(),
+    reason="cudaq.contrib.propagator requires the dynamics backend")
+def test_propagator_multi_degree_of_freedom():
+    hamiltonian = 0.25 * spin.z(0) + 0.5 * spin.x(1) + 0.75 * spin.z(2)
+    dimensions = {0: 2, 1: 2, 2: 2}
+    schedule = Schedule([0.0, 0.4], ["t"])
+
+    actual = cudaq.contrib.propagator(hamiltonian, dimensions, schedule)
+
+    expected = scipy.linalg.expm(-1j * hamiltonian.to_matrix(dimensions) * 0.4)
+
+    np.testing.assert_allclose(actual, expected, atol=1e-5)

From 5913a4d69727590576fa3833a5783e498ce0acda Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Fri, 5 Jun 2026 01:24:06 -0400
Subject: [PATCH 12/21] Pass schedule parameters through propagator callback

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py      |  4 ++--
 python/tests/contrib/test_propagator.py | 20 ++++++++++++++++++++
 2 files changed, 22 insertions(+), 2 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index 7796df331dc..e56c0e5990b 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -85,7 +85,7 @@ def _open_system_generator(hamiltonian, dimensions, collapse_operators,
     propagator_dimension = system_dimension * system_dimension
     operator_id = f"propagator_open_system_{uuid.uuid4().hex}"
 
-    def create(t=0.0, dimension=propagator_dimension):
+    def create(dimension=propagator_dimension, **parameters):
         if dimension != propagator_dimension:
             raise ValueError("Unexpected open-system propagator dimension.")
         return _open_system_matrix(
@@ -93,7 +93,7 @@ def create(t=0.0, dimension=propagator_dimension):
             dimensions,
             collapse_operators,
             collapse_operator_adjoint_ops,
-            t=t,
+            **parameters,
         )
 
     operators.define(operator_id, [propagator_dimension], create)
diff --git a/python/tests/contrib/test_propagator.py b/python/tests/contrib/test_propagator.py
index f145b4820ab..8bf1a75db97 100644
--- a/python/tests/contrib/test_propagator.py
+++ b/python/tests/contrib/test_propagator.py
@@ -226,3 +226,23 @@ def test_propagator_multi_degree_of_freedom():
     expected = scipy.linalg.expm(-1j * hamiltonian.to_matrix(dimensions) * 0.4)
 
     np.testing.assert_allclose(actual, expected, atol=1e-5)
+
+
+@pytest.mark.skipif(
+    not _has_dynamics_backend(),
+    reason="cudaq.contrib.propagator requires the dynamics backend")
+def test_open_system_propagator_schedule_parameter_name():
+    gamma = 0.2
+    hamiltonian = 0.1 * spin.z(0)
+    collapse_operator = np.sqrt(gamma) * spin.minus(0)
+    collapse_operator_adjoint = np.sqrt(gamma) * spin.plus(0)
+
+    computed = cudaq.contrib.propagator(
+        hamiltonian,
+        {0: 2},
+        Schedule([0.0, 0.1], ["time"]),
+        collapse_operators=[collapse_operator],
+        collapse_operator_adjoint_ops=[collapse_operator_adjoint],
+    )
+
+    assert computed.shape == (4, 4)

From 95b1864c67f7715094735714eec44272a77d1a13 Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Fri, 5 Jun 2026 17:18:08 -0400
Subject: [PATCH 13/21] Use structured SuperOperator for open-system
 propagators

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py      | 69 +++++--------------------
 python/tests/contrib/test_propagator.py | 20 +++++++
 2 files changed, 32 insertions(+), 57 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index e56c0e5990b..fceefe8b904 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -10,7 +10,6 @@
 
 from collections.abc import Mapping, Sequence
 from typing import Optional
-import uuid
 
 import cudaq
 import numpy as np
@@ -45,62 +44,20 @@ def _closed_system_generator(hamiltonian):
     return generator
 
 
-def _operator_matrix(operator, dimensions: Mapping[int, int], **parameters):
-    if parameters:
-        return np.asarray(operator.to_matrix(dimensions, **parameters),
-                          dtype=np.complex128)
-    return np.asarray(operator.to_matrix(dimensions), dtype=np.complex128)
-
-
-def _open_system_matrix(hamiltonian, dimensions, collapse_operators,
-                        collapse_operator_adjoint_ops, **parameters):
-    system_dimension = _total_dimension(dimensions)
-    identity = np.eye(system_dimension, dtype=np.complex128)
-
-    hamiltonian_matrix = _operator_matrix(hamiltonian, dimensions, **parameters)
-
-    matrix = -1j * np.kron(identity, hamiltonian_matrix)
-    matrix += 1j * np.kron(hamiltonian_matrix.T, identity)
+def _open_system_generator(hamiltonian, collapse_operators,
+                           collapse_operator_adjoint_ops):
+    generator = cudaq.SuperOperator()
+    generator += cudaq.SuperOperator.left_multiply(-1j * hamiltonian)
+    generator += cudaq.SuperOperator.right_multiply(1j * hamiltonian)
 
     for collapse_operator, collapse_operator_adjoint in zip(
             collapse_operators, collapse_operator_adjoint_ops):
-        collapse_matrix = _operator_matrix(collapse_operator, dimensions,
-                                           **parameters)
-        collapse_adjoint_matrix = _operator_matrix(collapse_operator_adjoint,
-                                                   dimensions, **parameters)
-        collapse_product = collapse_adjoint_matrix @ collapse_matrix
+        collapse_product = collapse_operator_adjoint * collapse_operator
+        generator += cudaq.SuperOperator.left_right_multiply(
+            collapse_operator, collapse_operator_adjoint)
+        generator += cudaq.SuperOperator.left_multiply(-0.5 * collapse_product)
+        generator += cudaq.SuperOperator.right_multiply(-0.5 * collapse_product)
 
-        matrix += np.kron(collapse_matrix.conj(), collapse_matrix)
-        matrix += -0.5 * np.kron(identity, collapse_product)
-        matrix += -0.5 * np.kron(collapse_product.T, identity)
-
-    return matrix
-
-
-def _open_system_generator(hamiltonian, dimensions, collapse_operators,
-                           collapse_operator_adjoint_ops):
-    from cudaq import operators
-
-    system_dimension = _total_dimension(dimensions)
-    propagator_dimension = system_dimension * system_dimension
-    operator_id = f"propagator_open_system_{uuid.uuid4().hex}"
-
-    def create(dimension=propagator_dimension, **parameters):
-        if dimension != propagator_dimension:
-            raise ValueError("Unexpected open-system propagator dimension.")
-        return _open_system_matrix(
-            hamiltonian,
-            dimensions,
-            collapse_operators,
-            collapse_operator_adjoint_ops,
-            **parameters,
-        )
-
-    operators.define(operator_id, [propagator_dimension], create)
-    generator_matrix = operators.instantiate(operator_id, 0)
-
-    generator = cudaq.SuperOperator()
-    generator += cudaq.SuperOperator.left_multiply(generator_matrix)
     return generator
 
 
@@ -179,9 +136,7 @@ def propagator(
     system_dimension = _total_dimension(dimensions)
     propagator_dimension = (system_dimension * system_dimension
                             if open_system else system_dimension)
-    evolution_dimensions = ({
-        0: propagator_dimension
-    } if open_system else dimensions)
+    evolution_dimensions = dimensions
 
     is_batched = isinstance(hamiltonian,
                             Sequence) and not hasattr(hamiltonian, "to_matrix")
@@ -189,7 +144,7 @@ def propagator(
 
     if open_system:
         generators = [
-            _open_system_generator(h, dimensions, collapse_operators,
+            _open_system_generator(h, collapse_operators,
                                    collapse_operator_adjoint_ops)
             for h in hamiltonians
         ]
diff --git a/python/tests/contrib/test_propagator.py b/python/tests/contrib/test_propagator.py
index 8bf1a75db97..6b563f35c4a 100644
--- a/python/tests/contrib/test_propagator.py
+++ b/python/tests/contrib/test_propagator.py
@@ -246,3 +246,23 @@ def test_open_system_propagator_schedule_parameter_name():
     )
 
     assert computed.shape == (4, 4)
+
+
+@pytest.mark.skipif(
+    not _has_dynamics_backend(),
+    reason="cudaq.contrib.propagator requires the dynamics backend")
+def test_open_system_propagator_product_collapse_operator():
+    hamiltonian = 0.1 * spin.z(0) + 0.2 * spin.z(1)
+    collapse_operator = 0.3 * spin.minus(0) * spin.minus(1)
+    collapse_operator_adjoint = 0.3 * spin.plus(0) * spin.plus(1)
+    dimensions = {0: 2, 1: 2}
+
+    computed = cudaq.contrib.propagator(
+        hamiltonian,
+        dimensions,
+        Schedule([0.0, 0.1], ["t"]),
+        collapse_operators=[collapse_operator],
+        collapse_operator_adjoint_ops=[collapse_operator_adjoint],
+    )
+
+    assert computed.shape == (16, 16)

From 4cdb7a17d926cce306cc1e287884e626eb9dfcb6 Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Fri, 5 Jun 2026 17:52:31 -0400
Subject: [PATCH 14/21] Use structured SuperOperator for open-system
 propagators

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py      | 34 ++++++++++++++++++++++++-
 python/tests/contrib/test_propagator.py | 26 -------------------
 2 files changed, 33 insertions(+), 27 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index fceefe8b904..c223ff45f2d 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -27,6 +27,15 @@ def _identity_state(dimension: int):
     return cudaq.State.from_data(identity)
 
 
+def _basis_states(dimension: int):
+    states = []
+    for index in range(dimension):
+        data = np.zeros(dimension, dtype=np.complex128)
+        data[index] = 1.0
+        states.append(cudaq.State.from_data(data))
+    return states
+
+
 def _state_to_matrix(state, dimension: int) -> np.ndarray:
     data = np.array(state).reshape(-1)
     expected_size = dimension * dimension
@@ -72,6 +81,14 @@ def _extract_propagator(result, dimension: int,
     return _state_to_matrix(result.final_state(), dimension)
 
 
+def _extract_batched_basis_propagator(results):
+    columns = [
+        np.array(single_result.final_state()).reshape(-1)
+        for single_result in results
+    ]
+    return np.column_stack(columns)
+
+
 def propagator(
     hamiltonian,
     dimensions: Mapping[int, int],
@@ -151,7 +168,14 @@ def propagator(
     else:
         generators = [_closed_system_generator(h) for h in hamiltonians]
 
-    initial_states = [_identity_state(propagator_dimension) for _ in generators]
+    if open_system:
+        initial_states = [
+            _basis_states(propagator_dimension) for _ in generators
+        ]
+    else:
+        initial_states = [
+            _identity_state(propagator_dimension) for _ in generators
+        ]
 
     save_mode = (cudaq.IntermediateResultSave.ALL if store_intermediate_results
                  else cudaq.IntermediateResultSave.NONE)
@@ -168,6 +192,14 @@ def propagator(
         max_batch_size=max_batch_size,
     )
 
+    if open_system:
+        if is_batched:
+            return [
+                _extract_batched_basis_propagator(single_result)
+                for single_result in result
+            ]
+        return _extract_batched_basis_propagator(result)
+
     if is_batched:
         return [
             _extract_propagator(single_result, propagator_dimension,
diff --git a/python/tests/contrib/test_propagator.py b/python/tests/contrib/test_propagator.py
index 6b563f35c4a..86a68159785 100644
--- a/python/tests/contrib/test_propagator.py
+++ b/python/tests/contrib/test_propagator.py
@@ -213,9 +213,6 @@ def test_open_system_propagator_shape():
     assert computed.shape == (4, 4)
 
 
-@pytest.mark.skipif(
-    not _has_dynamics_backend(),
-    reason="cudaq.contrib.propagator requires the dynamics backend")
 def test_propagator_multi_degree_of_freedom():
     hamiltonian = 0.25 * spin.z(0) + 0.5 * spin.x(1) + 0.75 * spin.z(2)
     dimensions = {0: 2, 1: 2, 2: 2}
@@ -228,9 +225,6 @@ def test_propagator_multi_degree_of_freedom():
     np.testing.assert_allclose(actual, expected, atol=1e-5)
 
 
-@pytest.mark.skipif(
-    not _has_dynamics_backend(),
-    reason="cudaq.contrib.propagator requires the dynamics backend")
 def test_open_system_propagator_schedule_parameter_name():
     gamma = 0.2
     hamiltonian = 0.1 * spin.z(0)
@@ -246,23 +240,3 @@ def test_open_system_propagator_schedule_parameter_name():
     )
 
     assert computed.shape == (4, 4)
-
-
-@pytest.mark.skipif(
-    not _has_dynamics_backend(),
-    reason="cudaq.contrib.propagator requires the dynamics backend")
-def test_open_system_propagator_product_collapse_operator():
-    hamiltonian = 0.1 * spin.z(0) + 0.2 * spin.z(1)
-    collapse_operator = 0.3 * spin.minus(0) * spin.minus(1)
-    collapse_operator_adjoint = 0.3 * spin.plus(0) * spin.plus(1)
-    dimensions = {0: 2, 1: 2}
-
-    computed = cudaq.contrib.propagator(
-        hamiltonian,
-        dimensions,
-        Schedule([0.0, 0.1], ["t"]),
-        collapse_operators=[collapse_operator],
-        collapse_operator_adjoint_ops=[collapse_operator_adjoint],
-    )
-
-    assert computed.shape == (16, 16)

From 88d0e6ebf10380b9ed323fb5d8ddae8e0fc354ec Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Sun, 7 Jun 2026 01:37:42 -0400
Subject: [PATCH 15/21] Use evolve collapse operators for open-system
 propagators

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py      | 37 +------------------------
 python/tests/contrib/test_propagator.py |  5 ----
 2 files changed, 1 insertion(+), 41 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index c223ff45f2d..f89615fa574 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -53,23 +53,6 @@ def _closed_system_generator(hamiltonian):
     return generator
 
 
-def _open_system_generator(hamiltonian, collapse_operators,
-                           collapse_operator_adjoint_ops):
-    generator = cudaq.SuperOperator()
-    generator += cudaq.SuperOperator.left_multiply(-1j * hamiltonian)
-    generator += cudaq.SuperOperator.right_multiply(1j * hamiltonian)
-
-    for collapse_operator, collapse_operator_adjoint in zip(
-            collapse_operators, collapse_operator_adjoint_ops):
-        collapse_product = collapse_operator_adjoint * collapse_operator
-        generator += cudaq.SuperOperator.left_right_multiply(
-            collapse_operator, collapse_operator_adjoint)
-        generator += cudaq.SuperOperator.left_multiply(-0.5 * collapse_product)
-        generator += cudaq.SuperOperator.right_multiply(-0.5 * collapse_product)
-
-    return generator
-
-
 def _extract_propagator(result, dimension: int,
                         store_intermediate_results: bool):
     if store_intermediate_results:
@@ -95,7 +78,6 @@ def propagator(
     schedule,
     *,
     collapse_operators=None,
-    collapse_operator_adjoint_ops=None,
     store_intermediate_results: bool = False,
     integrator=None,
     max_batch_size: Optional[int] = None,
@@ -117,8 +99,6 @@ def propagator(
         schedule: CUDA-Q dynamics schedule.
         collapse_operators: Optional sequence of Lindblad collapse operators.
             If provided, the helper returns the Lindblad map.
-        collapse_operator_adjoint_ops: Optional sequence containing the adjoint
-            operator for each collapse operator.
         store_intermediate_results: If True, return propagators at the
             intermediate schedule points saved by the dynamics backend.
         integrator: Optional dynamics integrator.
@@ -137,18 +117,7 @@ def propagator(
     """
     collapse_operators = [] if collapse_operators is None else list(
         collapse_operators)
-    collapse_operator_adjoint_ops = ([] if collapse_operator_adjoint_ops is None
-                                     else list(collapse_operator_adjoint_ops))
-
-    if collapse_operator_adjoint_ops and len(
-            collapse_operator_adjoint_ops) != len(collapse_operators):
-        raise ValueError("collapse_operator_adjoint_ops must have the same "
-                         "length as collapse_operators.")
-
     open_system = len(collapse_operators) > 0
-    if open_system and not collapse_operator_adjoint_ops:
-        raise ValueError("collapse_operator_adjoint_ops must be provided for "
-                         "open-system propagators.")
 
     system_dimension = _total_dimension(dimensions)
     propagator_dimension = (system_dimension * system_dimension
@@ -160,11 +129,7 @@ def propagator(
     hamiltonians = list(hamiltonian) if is_batched else [hamiltonian]
 
     if open_system:
-        generators = [
-            _open_system_generator(h, collapse_operators,
-                                   collapse_operator_adjoint_ops)
-            for h in hamiltonians
-        ]
+        generators = hamiltonians
     else:
         generators = [_closed_system_generator(h) for h in hamiltonians]
 
diff --git a/python/tests/contrib/test_propagator.py b/python/tests/contrib/test_propagator.py
index 86a68159785..a8b68f20d47 100644
--- a/python/tests/contrib/test_propagator.py
+++ b/python/tests/contrib/test_propagator.py
@@ -185,7 +185,6 @@ def test_open_system_propagator_amplitude_damping():
         {0: 2},
         schedule,
         collapse_operators=collapse_operators,
-        collapse_operator_adjoint_ops=[np.sqrt(gamma) * spin.plus(0)],
     )
 
     expected = scipy_linalg.expm(
@@ -207,7 +206,6 @@ def test_open_system_propagator_shape():
         {0: 2},
         schedule,
         collapse_operators=collapse_operators,
-        collapse_operator_adjoint_ops=[np.sqrt(gamma) * spin.plus(0)],
     )
 
     assert computed.shape == (4, 4)
@@ -229,14 +227,11 @@ def test_open_system_propagator_schedule_parameter_name():
     gamma = 0.2
     hamiltonian = 0.1 * spin.z(0)
     collapse_operator = np.sqrt(gamma) * spin.minus(0)
-    collapse_operator_adjoint = np.sqrt(gamma) * spin.plus(0)
-
     computed = cudaq.contrib.propagator(
         hamiltonian,
         {0: 2},
         Schedule([0.0, 0.1], ["time"]),
         collapse_operators=[collapse_operator],
-        collapse_operator_adjoint_ops=[collapse_operator_adjoint],
     )
 
     assert computed.shape == (4, 4)

From a1c0ecf9aa39acee784e7663450bb2afe3aa6415 Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Sun, 7 Jun 2026 01:52:27 -0400
Subject: [PATCH 16/21] Pass collapse operators through evolve

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py      | 2 +-
 python/tests/contrib/test_propagator.py | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index f89615fa574..0a827f10331 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -150,7 +150,7 @@ def propagator(
         evolution_dimensions,
         schedule,
         initial_states if is_batched else initial_states[0],
-        collapse_operators=[],
+        collapse_operators=collapse_operators,
         observables=[],
         store_intermediate_results=save_mode,
         integrator=integrator,
diff --git a/python/tests/contrib/test_propagator.py b/python/tests/contrib/test_propagator.py
index a8b68f20d47..08f20235dcd 100644
--- a/python/tests/contrib/test_propagator.py
+++ b/python/tests/contrib/test_propagator.py
@@ -214,11 +214,11 @@ def test_open_system_propagator_shape():
 def test_propagator_multi_degree_of_freedom():
     hamiltonian = 0.25 * spin.z(0) + 0.5 * spin.x(1) + 0.75 * spin.z(2)
     dimensions = {0: 2, 1: 2, 2: 2}
-    schedule = Schedule([0.0, 0.4], ["t"])
+    schedule = Schedule(np.linspace(0.0, 0.4, 21), ["t"])
 
     actual = cudaq.contrib.propagator(hamiltonian, dimensions, schedule)
 
-    expected = scipy.linalg.expm(-1j * hamiltonian.to_matrix(dimensions) * 0.4)
+    expected = scipy_linalg.expm(-1j * hamiltonian.to_matrix(dimensions) * 0.4)
 
     np.testing.assert_allclose(actual, expected, atol=1e-5)
 

From 2fd6184a9abf77ee2cc06ce3bdeefceba0db7b8a Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Tue, 9 Jun 2026 19:52:48 -0400
Subject: [PATCH 17/21] Support batched collapse operators

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py      | 51 +++++++++++++++--
 python/tests/contrib/test_propagator.py | 73 +++++++++++++++++++++++++
 2 files changed, 118 insertions(+), 6 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index 0a827f10331..a257481c0dc 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -16,6 +16,7 @@
 
 
 def _total_dimension(dimensions: Mapping[int, int]) -> int:
+    """Return the product of all local Hilbert-space dimensions."""
     dimension = 1
     for local_dimension in dimensions.values():
         dimension *= local_dimension
@@ -23,11 +24,13 @@ def _total_dimension(dimensions: Mapping[int, int]) -> int:
 
 
 def _identity_state(dimension: int):
+    """Return |I> used to evolve closed-system propagators directly."""
     identity = np.eye(dimension, dtype=np.complex128).reshape(-1)
     return cudaq.State.from_data(identity)
 
 
 def _basis_states(dimension: int):
+    """Return Liouville basis states used to reconstruct Lindblad maps."""
     states = []
     for index in range(dimension):
         data = np.zeros(dimension, dtype=np.complex128)
@@ -37,6 +40,7 @@ def _basis_states(dimension: int):
 
 
 def _state_to_matrix(state, dimension: int) -> np.ndarray:
+    """Convert a vectorized propagated identity state back to a matrix."""
     data = np.array(state).reshape(-1)
     expected_size = dimension * dimension
 
@@ -48,6 +52,7 @@ def _state_to_matrix(state, dimension: int) -> np.ndarray:
 
 
 def _closed_system_generator(hamiltonian):
+    """Represent dU/dt = -i H U as left multiplication by -i H."""
     generator = cudaq.SuperOperator()
     generator += cudaq.SuperOperator.left_multiply(-1j * hamiltonian)
     return generator
@@ -55,6 +60,7 @@ def _closed_system_generator(hamiltonian):
 
 def _extract_propagator(result, dimension: int,
                         store_intermediate_results: bool):
+    """Extract closed-system propagators from a CUDA-Q evolve result."""
     if store_intermediate_results:
         return [
             _state_to_matrix(state, dimension)
@@ -65,6 +71,7 @@ def _extract_propagator(result, dimension: int,
 
 
 def _extract_batched_basis_propagator(results):
+    """Stack evolved Liouville basis states into a dense Lindblad map."""
     columns = [
         np.array(single_result.final_state()).reshape(-1)
         for single_result in results
@@ -72,6 +79,31 @@ def _extract_batched_basis_propagator(results):
     return np.column_stack(columns)
 
 
+def _is_operator_like(value) -> bool:
+    """Return True for CUDA-Q operators accepted by the dynamics backend."""
+    return hasattr(value, "to_matrix")
+
+
+def _is_collapse_operator_batch(collapse_operators) -> bool:
+    """Return True when collapse operators are grouped per Hamiltonian."""
+    return bool(collapse_operators) and not _is_operator_like(
+        collapse_operators[0])
+
+
+def _collapse_operator_batches(collapse_operators, batch_size: int):
+    """Broadcast collapse operators to match the Hamiltonian batch size."""
+    if not collapse_operators:
+        return [[] for _ in range(batch_size)]
+
+    if _is_collapse_operator_batch(collapse_operators):
+        if len(collapse_operators) != batch_size:
+            raise ValueError("Batched collapse_operators must have the same "
+                             "length as the Hamiltonian batch.")
+        return [list(ops) for ops in collapse_operators]
+
+    return [collapse_operators for _ in range(batch_size)]
+
+
 def propagator(
     hamiltonian,
     dimensions: Mapping[int, int],
@@ -117,17 +149,19 @@ def propagator(
     """
     collapse_operators = [] if collapse_operators is None else list(
         collapse_operators)
-    open_system = len(collapse_operators) > 0
+
+    is_batched = isinstance(hamiltonian,
+                            Sequence) and not hasattr(hamiltonian, "to_matrix")
+    hamiltonians = list(hamiltonian) if is_batched else [hamiltonian]
+    collapse_operator_batches = _collapse_operator_batches(
+        collapse_operators, len(hamiltonians))
+    open_system = any(collapse_operator_batches)
 
     system_dimension = _total_dimension(dimensions)
     propagator_dimension = (system_dimension * system_dimension
                             if open_system else system_dimension)
     evolution_dimensions = dimensions
 
-    is_batched = isinstance(hamiltonian,
-                            Sequence) and not hasattr(hamiltonian, "to_matrix")
-    hamiltonians = list(hamiltonian) if is_batched else [hamiltonian]
-
     if open_system:
         generators = hamiltonians
     else:
@@ -145,12 +179,17 @@ def propagator(
     save_mode = (cudaq.IntermediateResultSave.ALL if store_intermediate_results
                  else cudaq.IntermediateResultSave.NONE)
 
+    evolve_collapse_operators = []
+    if open_system:
+        evolve_collapse_operators = (collapse_operator_batches if is_batched
+                                     else collapse_operator_batches[0])
+
     result = cudaq.evolve(
         generators if is_batched else generators[0],
         evolution_dimensions,
         schedule,
         initial_states if is_batched else initial_states[0],
-        collapse_operators=collapse_operators,
+        collapse_operators=evolve_collapse_operators,
         observables=[],
         store_intermediate_results=save_mode,
         integrator=integrator,
diff --git a/python/tests/contrib/test_propagator.py b/python/tests/contrib/test_propagator.py
index 08f20235dcd..95a99df2359 100644
--- a/python/tests/contrib/test_propagator.py
+++ b/python/tests/contrib/test_propagator.py
@@ -235,3 +235,76 @@ def test_open_system_propagator_schedule_parameter_name():
     )
 
     assert computed.shape == (4, 4)
+
+
+def test_open_system_propagator_nonzero_hamiltonian():
+    gamma = 0.2
+    omega = 0.5
+    t_final = 0.7
+
+    hamiltonian = 0.5 * omega * spin.z(0)
+    collapse_operators = [np.sqrt(gamma) * spin.minus(0)]
+    schedule = Schedule(np.linspace(0.0, t_final, 21), ["t"])
+
+    computed = cudaq.contrib.propagator(
+        hamiltonian,
+        {0: 2},
+        schedule,
+        collapse_operators=collapse_operators,
+    )
+
+    liouvillian = _amplitude_damping_liouvillian(gamma)
+    hamiltonian_matrix = hamiltonian.to_matrix({0: 2})
+    coherent = (-1j * np.kron(np.eye(2), hamiltonian_matrix) +
+                1j * np.kron(hamiltonian_matrix.T, np.eye(2)))
+    expected = scipy_linalg.expm((coherent + liouvillian) * t_final)
+
+    np.testing.assert_allclose(computed, expected, atol=1e-4)
+
+
+def test_open_system_propagator_batched_shared_collapse_operators():
+    gamma = 0.2
+    t_final = 0.7
+    omegas = [0.2, 0.5]
+
+    hamiltonians = [0.5 * omega * spin.z(0) for omega in omegas]
+    collapse_operators = [np.sqrt(gamma) * spin.minus(0)]
+    schedule = Schedule(np.linspace(0.0, t_final, 21), ["t"])
+
+    computed = cudaq.contrib.propagator(
+        hamiltonians,
+        {0: 2},
+        schedule,
+        collapse_operators=collapse_operators,
+    )
+
+    assert len(computed) == len(hamiltonians)
+    for actual, hamiltonian in zip(computed, hamiltonians):
+        hamiltonian_matrix = hamiltonian.to_matrix({0: 2})
+        coherent = (-1j * np.kron(np.eye(2), hamiltonian_matrix) +
+                    1j * np.kron(hamiltonian_matrix.T, np.eye(2)))
+        expected = scipy_linalg.expm(
+            (coherent + _amplitude_damping_liouvillian(gamma)) * t_final)
+        np.testing.assert_allclose(actual, expected, atol=1e-4)
+
+
+def test_open_system_propagator_batched_collapse_operator_lists():
+    gammas = [0.1, 0.2]
+    t_final = 0.7
+
+    hamiltonians = [0.0 * spin.z(0), 0.0 * spin.z(0)]
+    collapse_operators = [[np.sqrt(gamma) * spin.minus(0)] for gamma in gammas]
+    schedule = Schedule(np.linspace(0.0, t_final, 21), ["t"])
+
+    computed = cudaq.contrib.propagator(
+        hamiltonians,
+        {0: 2},
+        schedule,
+        collapse_operators=collapse_operators,
+    )
+
+    assert len(computed) == len(gammas)
+    for actual, gamma in zip(computed, gammas):
+        expected = scipy_linalg.expm(
+            _amplitude_damping_liouvillian(gamma) * t_final)
+        np.testing.assert_allclose(actual, expected, atol=1e-4)

From e3ee81ca247a639eff5caadac5e953cd555d0c7e Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Tue, 9 Jun 2026 20:08:52 -0400
Subject: [PATCH 18/21] Flatten batched open-system propagator states

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py | 24 ++++++++++++++++++++----
 1 file changed, 20 insertions(+), 4 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index a257481c0dc..c5dfd0e0d37 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -184,11 +184,25 @@ def propagator(
         evolve_collapse_operators = (collapse_operator_batches if is_batched
                                      else collapse_operator_batches[0])
 
+    evolve_generators = generators if is_batched else generators[0]
+    evolve_initial_states = initial_states if is_batched else initial_states[0]
+
+    if open_system and is_batched:
+        evolve_generators = []
+        evolve_initial_states = []
+        evolve_collapse_operators = []
+        for generator, basis_states, collapse_batch in zip(
+                generators, initial_states, collapse_operator_batches):
+            for basis_state in basis_states:
+                evolve_generators.append(generator)
+                evolve_initial_states.append(basis_state)
+                evolve_collapse_operators.append(collapse_batch)
+
     result = cudaq.evolve(
-        generators if is_batched else generators[0],
+        evolve_generators,
         evolution_dimensions,
         schedule,
-        initial_states if is_batched else initial_states[0],
+        evolve_initial_states,
         collapse_operators=evolve_collapse_operators,
         observables=[],
         store_intermediate_results=save_mode,
@@ -199,8 +213,10 @@ def propagator(
     if open_system:
         if is_batched:
             return [
-                _extract_batched_basis_propagator(single_result)
-                for single_result in result
+                _extract_batched_basis_propagator(
+                    result[index * propagator_dimension:(index + 1) *
+                           propagator_dimension])
+                for index in range(len(generators))
             ]
         return _extract_batched_basis_propagator(result)
 

From 807749876f39aacb4384bd06d7b9b1a75d51bf68 Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Tue, 9 Jun 2026 20:40:31 -0400
Subject: [PATCH 19/21] Escape technical terms in propagator comments

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index c5dfd0e0d37..46ff2edffd9 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -30,7 +30,7 @@ def _identity_state(dimension: int):
 
 
 def _basis_states(dimension: int):
-    """Return Liouville basis states used to reconstruct Lindblad maps."""
+    """Return `Liouville` basis states used to reconstruct Lindblad maps."""
     states = []
     for index in range(dimension):
         data = np.zeros(dimension, dtype=np.complex128)
@@ -40,7 +40,7 @@ def _basis_states(dimension: int):
 
 
 def _state_to_matrix(state, dimension: int) -> np.ndarray:
-    """Convert a vectorized propagated identity state back to a matrix."""
+    """Convert a `vectorized` propagated identity state back to a matrix."""
     data = np.array(state).reshape(-1)
     expected_size = dimension * dimension
 
@@ -52,7 +52,7 @@ def _state_to_matrix(state, dimension: int) -> np.ndarray:
 
 
 def _closed_system_generator(hamiltonian):
-    """Represent dU/dt = -i H U as left multiplication by -i H."""
+    """Represent `dU/dt = -i H U` as left multiplication by `-i H`."""
     generator = cudaq.SuperOperator()
     generator += cudaq.SuperOperator.left_multiply(-1j * hamiltonian)
     return generator

From ba05bf2a31802b5154d8746075d3bb53dd33b8eb Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Tue, 9 Jun 2026 20:50:06 -0400
Subject: [PATCH 20/21] Support open-system intermediate propagators

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py      | 22 ++++++++++++++++++----
 python/tests/contrib/test_propagator.py | 23 +++++++++++++++++++++++
 2 files changed, 41 insertions(+), 4 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index 46ff2edffd9..eb584a7a64c 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -70,8 +70,21 @@ def _extract_propagator(result, dimension: int,
     return _state_to_matrix(result.final_state(), dimension)
 
 
-def _extract_batched_basis_propagator(results):
-    """Stack evolved Liouville basis states into a dense Lindblad map."""
+def _extract_batched_basis_propagator(results,
+                                      store_intermediate_results: bool):
+    """Stack evolved `Liouville` basis states into dense Lindblad maps."""
+    if store_intermediate_results:
+        intermediate_states = [
+            list(single_result.intermediate_states())
+            for single_result in results
+        ]
+        return [
+            np.column_stack(
+                np.array(states[time_index]).reshape(-1)
+                for states in intermediate_states)
+            for time_index in range(len(intermediate_states[0]))
+        ]
+
     columns = [
         np.array(single_result.final_state()).reshape(-1)
         for single_result in results
@@ -215,10 +228,11 @@ def propagator(
             return [
                 _extract_batched_basis_propagator(
                     result[index * propagator_dimension:(index + 1) *
-                           propagator_dimension])
+                           propagator_dimension], store_intermediate_results)
                 for index in range(len(generators))
             ]
-        return _extract_batched_basis_propagator(result)
+        return _extract_batched_basis_propagator(result,
+                                                 store_intermediate_results)
 
     if is_batched:
         return [
diff --git a/python/tests/contrib/test_propagator.py b/python/tests/contrib/test_propagator.py
index 95a99df2359..3570f70508c 100644
--- a/python/tests/contrib/test_propagator.py
+++ b/python/tests/contrib/test_propagator.py
@@ -308,3 +308,26 @@ def test_open_system_propagator_batched_collapse_operator_lists():
         expected = scipy_linalg.expm(
             _amplitude_damping_liouvillian(gamma) * t_final)
         np.testing.assert_allclose(actual, expected, atol=1e-4)
+
+
+def test_open_system_propagator_intermediate_results():
+    gamma = 0.2
+    times = np.linspace(0.0, 0.7, 6)
+
+    hamiltonian = 0.0 * spin.z(0)
+    collapse_operators = [np.sqrt(gamma) * spin.minus(0)]
+    schedule = Schedule(times, ["t"])
+
+    computed = cudaq.contrib.propagator(
+        hamiltonian,
+        {0: 2},
+        schedule,
+        collapse_operators=collapse_operators,
+        store_intermediate_results=True,
+    )
+
+    assert len(computed) == len(times)
+    for actual, time in zip(computed, times):
+        expected = scipy_linalg.expm(
+            _amplitude_damping_liouvillian(gamma) * time)
+        np.testing.assert_allclose(actual, expected, atol=1e-4)

From 7fa78d61151203715e1b04cec1fe06fa703e1fbb Mon Sep 17 00:00:00 2001
From: Claudia Friedsam <claudia.friedsam@gmail.com>
Date: Tue, 9 Jun 2026 21:39:22 -0400
Subject: [PATCH 21/21] Use list for intermediate propagator stacking

Signed-off-by: Claudia Friedsam <claudia.friedsam@gmail.com>
---
 python/cudaq/contrib/propagator.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/python/cudaq/contrib/propagator.py b/python/cudaq/contrib/propagator.py
index eb584a7a64c..d432b9e13c9 100644
--- a/python/cudaq/contrib/propagator.py
+++ b/python/cudaq/contrib/propagator.py
@@ -79,9 +79,10 @@ def _extract_batched_basis_propagator(results,
             for single_result in results
         ]
         return [
-            np.column_stack(
+            np.column_stack([
                 np.array(states[time_index]).reshape(-1)
-                for states in intermediate_states)
+                for states in intermediate_states
+            ])
             for time_index in range(len(intermediate_states[0]))
         ]