diff --git a/CMakeLists.txt b/CMakeLists.txt
index e9bf3fef4..e8a92a47b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -12,7 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-cmake_minimum_required(VERSION 3.31)
+cmake_minimum_required(VERSION 3.28)
 project(qsim LANGUAGES CXX)
 
 include(CheckLanguage)
diff --git a/Makefile b/Makefile
index ec8adddb0..b4e16da06 100644
--- a/Makefile
+++ b/Makefile
@@ -171,9 +171,9 @@ run-tests tests: $(TESTS)
 
 .PHONY: check-cuquantum-root-set
 check-cuquantum-root-set:
-	@if [[ -z "$(CUQUANTUM_ROOT)" ]]; then \
-	    echo Error: '$$CUQUANTUM_ROOT must be set in order to use cuStateVec.' \
-	    exit 1 \
+	@if test -z "$(CUQUANTUM_ROOT)"; then \
+	    echo Error: '$$CUQUANTUM_ROOT must be set in order to use cuStateVec.'; \
+	    exit 1; \
 	fi
 
 eigen:
diff --git a/apps/qsim_qtrajectory_cuda.cu b/apps/qsim_qtrajectory_cuda.cu
index 0d513cfa9..006d3944f 100644
--- a/apps/qsim_qtrajectory_cuda.cu
+++ b/apps/qsim_qtrajectory_cuda.cu
@@ -30,6 +30,7 @@
 #include "../lib/gates_qsim.h"
 #include "../lib/io_file.h"
 #include "../lib/qtrajectory.h"
+#include "../lib/run_qsim.h"
 #include "../lib/simulator_cuda.h"
 
 struct Options {
@@ -190,7 +191,7 @@ int main(int argc, char* argv[]) {
   using fp_type = float;
 
   struct Factory {
-    using Simulator = qsim::SimulatorCUDA<float>;
+    using Simulator = qsim::SimulatorCUDA<fp_type>;
     using StateSpace = Simulator::StateSpace;
 
     Factory(const StateSpace::Parameter& param) : param(param) {}
@@ -209,17 +210,18 @@ int main(int argc, char* argv[]) {
   using Simulator = Factory::Simulator;
   using StateSpace = Simulator::StateSpace;
   using State = StateSpace::State;
-  using Fuser = MultiQubitGateFuser<IO, GateQSim<fp_type>>;
-  using QTSimulator = QuantumTrajectorySimulator<IO, GateQSim<fp_type>,
-                                                 MultiQubitGateFuser,
-                                                 Simulator>;
+  using Gate = GateQSim<fp_type>;
+  using Fuser = MultiQubitGateFuser<IO, Gate>;
+  using FuserQT = MultiQubitGateFuser<IO, const Gate*>;
+  using RunnerQT = QSimRunner<IO, FuserQT, Factory>;
+  using QTSimulator = QuantumTrajectorySimulator<IO, Gate, RunnerQT>;
 
   auto opt = GetOptions(argc, argv);
   if (!ValidateOptions(opt)) {
     return 1;
   }
 
-  Circuit<GateQSim<fp_type>> circuit;
+  Circuit<Gate> circuit;
   unsigned maxtime = opt.times.back();
   if (!CircuitQsimParser<IOFile>::FromFile(maxtime, opt.circuit_file,
                                            circuit)) {
@@ -254,7 +256,7 @@ int main(int argc, char* argv[]) {
 
   auto noisy_circuits = AddNoise(circuit, opt.times, channel1, channel2);
 
-  auto observables = GetObservables<GateQSim<fp_type>>(circuit.num_qubits);
+  auto observables = GetObservables<Gate>(circuit.num_qubits);
 
   std::vector<std::vector<std::vector<std::complex<double>>>> results;
   results.reserve(opt.num_trajectories);
diff --git a/lib/BUILD b/lib/BUILD
index 6462d301d..5b2dd31aa 100644
--- a/lib/BUILD
+++ b/lib/BUILD
@@ -489,6 +489,7 @@ cc_library(
     name = "run_qsim",
     hdrs = ["run_qsim.h"],
     deps = [
+        ":circuit",
         ":gate",
         ":gate_appl",
         ":util",
@@ -650,7 +651,6 @@ cuda_library(
     name = "simulator_custatevec",
     hdrs = [
         "simulator_custatevec.h",
-        "simulator_custatevec_kernels.h",
     ],
     deps = [
         ":bits",
diff --git a/lib/circuit.h b/lib/circuit.h
index 59018eedf..7f0cbb503 100644
--- a/lib/circuit.h
+++ b/lib/circuit.h
@@ -31,6 +31,49 @@ struct Circuit {
   std::vector<Gate> gates;
 };
 
+namespace detail {
+
+/**
+ * An adapter for vectors of gates.
+ */
+template <typename Circuit>
+struct Gates;
+
+template <typename Gate>
+struct Gates<qsim::Circuit<Gate>> {
+  static const std::vector<Gate>& get(const qsim::Circuit<Gate>& circuit) {
+    return circuit.gates;
+  }
+
+  static const Gate& gate(const Gate& g) {
+    return g;
+  }
+};
+
+template <typename Gate>
+struct Gates<std::vector<Gate>> {
+  static const std::vector<Gate>& get(const std::vector<Gate>& gates) {
+    return gates;
+  }
+
+  static const Gate& gate(const Gate& g) {
+    return g;
+  }
+};
+
+template <typename Gate>
+struct Gates<std::vector<Gate*>> {
+  static const std::vector<Gate*>& get(const std::vector<Gate*>& gates) {
+    return gates;
+  }
+
+  static const Gate& gate(const Gate* g) {
+    return *g;
+  }
+};
+
+}  // namespace detail
+
 }  // namespace qsim
 
 #endif  // CIRCUIT_H_
diff --git a/lib/qtrajectory.h b/lib/qtrajectory.h
index 1da6692f6..3ea2a4fa5 100644
--- a/lib/qtrajectory.h
+++ b/lib/qtrajectory.h
@@ -31,19 +31,18 @@ namespace qsim {
  * Quantum trajectory simulator.
  */
 template <typename IO, typename Gate,
-          template <typename, typename> class FuserT, typename Simulator,
-          typename RGen = std::mt19937>
+          typename Runner, typename RGen = std::mt19937>
 class QuantumTrajectorySimulator {
  public:
-  using Fuser = FuserT<IO, const Gate*>;
+  using Simulator = typename Runner::Simulator;
   using StateSpace = typename Simulator::StateSpace;
-  using State = typename Simulator::State;
+  using State = typename StateSpace::State;
   using MeasurementResult = typename StateSpace::MeasurementResult;
 
   /**
    * User-specified parameters for the simulator.
    */
-  struct Parameter : public Fuser::Parameter {
+  struct Parameter : public Runner::Parameter {
     /**
      * If true, collect statistics of sampled Kraus operator indices.
      */
@@ -265,7 +264,8 @@ class QuantumTrajectorySimulator {
       if (channel[0].kind == gate::kMeasurement) {
         // Measurement channel.
 
-        if (!ApplyDeferredOps(param, num_qubits, simulator, gates, state)) {
+        if (!ApplyDeferredOps(
+                 param, num_qubits, state_space, simulator, gates, state)) {
           return false;
         }
 
@@ -309,7 +309,8 @@ class QuantumTrajectorySimulator {
 
       if (r < cp) continue;
 
-      if (!ApplyDeferredOps(param, num_qubits, simulator, gates, state)) {
+      if (!ApplyDeferredOps(
+               param, num_qubits, state_space, simulator, gates, state)) {
         return false;
       }
 
@@ -351,7 +352,8 @@ class QuantumTrajectorySimulator {
     }
 
     if (apply_last_deferred_ops || !stat.primary) {
-      if (!ApplyDeferredOps(param, num_qubits, simulator, gates, state)) {
+      if (!ApplyDeferredOps(
+               param, num_qubits, state_space, simulator, gates, state)) {
         return false;
       }
 
@@ -372,20 +374,13 @@ class QuantumTrajectorySimulator {
   }
 
   static bool ApplyDeferredOps(
-      const Parameter& param, unsigned num_qubits, const Simulator& simulator,
+      const Parameter& param, unsigned num_qubits,
+      const StateSpace& state_space, const Simulator& simulator,
       std::vector<const Gate*>& gates, State& state) {
     if (gates.size() > 0) {
-      auto fgates = Fuser::FuseGates(param, num_qubits, gates);
-
+      bool rc = Runner::Run(param, gates, state_space, simulator, state);
       gates.resize(0);
-
-      if (fgates.size() == 0) {
-        return false;
-      }
-
-      for (const auto& fgate : fgates) {
-        ApplyFusedGate(simulator, fgate, state);
-      }
+      return rc;
     }
 
     return true;
diff --git a/lib/run_qsim.h b/lib/run_qsim.h
index 375291527..d718c6b31 100644
--- a/lib/run_qsim.h
+++ b/lib/run_qsim.h
@@ -19,6 +19,7 @@
 #include <string>
 #include <vector>
 
+#include "circuit.h"
 #include "gate.h"
 #include "gate_appl.h"
 #include "util.h"
@@ -174,6 +175,56 @@ struct QSimRunner final {
   static bool Run(const Parameter& param, const Factory& factory,
                   const Circuit& circuit, State& state,
                   std::vector<MeasurementResult>& measure_results) {
+    StateSpace state_space = factory.CreateStateSpace();
+    Simulator simulator = factory.CreateSimulator();
+
+    return Run(param, circuit, state_space, simulator, state, measure_results);
+  }
+
+  /**
+   * Runs the given circuit and make the final state available to the caller,
+   * discarding the result of any intermediate measurements in the circuit.
+   * @param param Options for gate fusion, parallelism and logging.
+   * @param factory Object to create simulators and state spaces.
+   * @param circuit The circuit to be simulated.
+   * @param state As an input parameter, this should contain the initial state
+   *   of the system. After a successful run, it will be populated with the
+   *   final state of the system.
+   * @return True if the simulation completed successfully; false otherwise.
+   */
+  template <typename Circuit>
+  static bool Run(const Parameter& param, const Factory& factory,
+                  const Circuit& circuit, State& state) {
+    StateSpace state_space = factory.CreateStateSpace();
+    Simulator simulator = factory.CreateSimulator();
+
+    std::vector<MeasurementResult> discarded_results;
+
+    return Run(
+        param, circuit, state_space, simulator, state, discarded_results);
+  }
+
+  /**
+   * Runs the given circuit and make the final state available to the caller,
+   * recording the result of any intermediate measurements in the circuit.
+   * @param param Options for gate fusion, parallelism and logging.
+   * @param circuit The circuit to be simulated.
+   * @param state_space StateSpace object required to perform measurements.
+   * @param simulator Simulator object. Provides specific implementations for
+   *   applying gates.
+   * @param state As an input parameter, this should contain the initial state
+   *   of the system. After a successful run, it will be populated with the
+   *   final state of the system.
+   * @param measure_results As an input parameter, this should be empty.
+   *   After a successful run, this will contain all measurements results from
+   *   the run, ordered by time and qubit index.
+   * @return True if the simulation completed successfully; false otherwise.
+   */
+  template <typename Circuit>
+  static bool Run(const Parameter& param, const Circuit& circuit,
+                  const StateSpace& state_space, const Simulator& simulator,
+                  State& state,
+                  std::vector<MeasurementResult>& measure_results) {
     double t0 = 0.0;
     double t1 = 0.0;
 
@@ -183,19 +234,18 @@ struct QSimRunner final {
 
     RGen rgen(param.seed);
 
-    StateSpace state_space = factory.CreateStateSpace();
-    Simulator simulator = factory.CreateSimulator();
-
     if (param.verbosity > 1) {
       t1 = GetTime();
       IO::messagef("init time is %g seconds.\n", t1 - t0);
       t0 = GetTime();
     }
 
-    auto fused_gates = Fuser::FuseGates(param, circuit.num_qubits,
-                                        circuit.gates);
+    using Gates = detail::Gates<Circuit>;
+    const auto& gates = Gates::get(circuit);
 
-    if (fused_gates.size() == 0 && circuit.gates.size() > 0) {
+    auto fused_gates = Fuser::FuseGates(param, state.num_qubits(), gates);
+
+    if (fused_gates.size() == 0 && gates.size() > 0) {
       return false;
     }
 
@@ -242,18 +292,23 @@ struct QSimRunner final {
    * Runs the given circuit and make the final state available to the caller,
    * discarding the result of any intermediate measurements in the circuit.
    * @param param Options for gate fusion, parallelism and logging.
-   * @param factory Object to create simulators and state spaces.
    * @param circuit The circuit to be simulated.
+   * @param state_space StateSpace object required to perform measurements.
+   * @param simulator Simulator object. Provides specific implementations for
+   *   applying gates.
    * @param state As an input parameter, this should contain the initial state
    *   of the system. After a successful run, it will be populated with the
    *   final state of the system.
    * @return True if the simulation completed successfully; false otherwise.
    */
   template <typename Circuit>
-  static bool Run(const Parameter& param, const Factory& factory,
-                  const Circuit& circuit, State& state) {
+  static bool Run(const Parameter& param, const Circuit& circuit,
+                  const StateSpace& state_space, const Simulator& simulator,
+                  State& state) {
     std::vector<MeasurementResult> discarded_results;
-    return Run(param, factory, circuit, state, discarded_results);
+
+    return Run(
+        param, circuit, state_space, simulator, state, discarded_results);
   }
 };
 
diff --git a/lib/simulator_custatevec.h b/lib/simulator_custatevec.h
index 40d1902de..b3f3cb8fa 100644
--- a/lib/simulator_custatevec.h
+++ b/lib/simulator_custatevec.h
@@ -66,7 +66,7 @@ class SimulatorCuStateVec final {
       uint64_t size = uint64_t{1} << state.num_qubits();
 
       if (StateSpace::is_float) {
-        cuComplex a = {matrix[0], matrix[1]};
+        cuComplex a = {(float) matrix[0], (float) matrix[1]};
         auto p = (cuComplex*) state.get();
         ErrorCheck(cublasCscal(cublas_handle_, size, &a, p, 1));
       } else {
diff --git a/lib/vectorspace.h b/lib/vectorspace.h
index 7b33a538f..8040b58db 100644
--- a/lib/vectorspace.h
+++ b/lib/vectorspace.h
@@ -74,7 +74,7 @@ class VectorSpace {
       return num_qubits_;
     }
 
-    bool requires_copy_to_host() const {
+    static constexpr bool requires_copy_to_host() {
       return false;
     }
 
@@ -174,7 +174,7 @@ class VectorSpace {
     return true;
   }
 
-  void DeviceSync() {}
+  static void DeviceSync() {}
 
  protected:
   For for_;
diff --git a/lib/vectorspace_cuda.h b/lib/vectorspace_cuda.h
index fd91553d7..5cfd4e834 100644
--- a/lib/vectorspace_cuda.h
+++ b/lib/vectorspace_cuda.h
@@ -72,7 +72,7 @@ class VectorSpaceCUDA {
       return num_qubits_;
     }
 
-    bool requires_copy_to_host() const {
+    static constexpr bool requires_copy_to_host() {
       return true;
     }
 
@@ -160,7 +160,7 @@ class VectorSpaceCUDA {
     return true;
   }
 
-  void DeviceSync() {
+  static void DeviceSync() {
     ErrorCheck(cudaDeviceSynchronize());
   }
 
diff --git a/pybind_interface/avx2/CMakeLists.txt b/pybind_interface/avx2/CMakeLists.txt
index 16eed37d2..86a4e9a73 100644
--- a/pybind_interface/avx2/CMakeLists.txt
+++ b/pybind_interface/avx2/CMakeLists.txt
@@ -12,7 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-cmake_minimum_required(VERSION 3.31)
+cmake_minimum_required(VERSION 3.28)
 project(qsim)
 
 IF (WIN32)
diff --git a/pybind_interface/avx2/pybind_main_avx2.cpp b/pybind_interface/avx2/pybind_main_avx2.cpp
index 1f7dc0842..a70f1b0f2 100644
--- a/pybind_interface/avx2/pybind_main_avx2.cpp
+++ b/pybind_interface/avx2/pybind_main_avx2.cpp
@@ -15,6 +15,10 @@
 #include "pybind_main_avx2.h"
 
 #include "../../lib/formux.h"
+#include "../../lib/fuser_mqubit.h"
+#include "../../lib/gates_cirq.h"
+#include "../../lib/io.h"
+#include "../../lib/run_qsim.h"
 #include "../../lib/simulator_avx.h"
 #include "../../lib/util_cpu.h"
 
@@ -33,6 +37,14 @@ namespace qsim {
     using Simulator = qsim::Simulator<For>;
     using StateSpace = Simulator::StateSpace;
 
+    using Gate = Cirq::GateCirq<float>;
+    using Runner = QSimRunner<IO, MultiQubitGateFuser<IO, Gate>, Factory>;
+    using RunnerQT =
+        QSimRunner<IO, MultiQubitGateFuser<IO, const Gate*>, Factory>;
+    using RunnerParameter = Runner::Parameter;
+    using NoisyRunner = qsim::QuantumTrajectorySimulator<IO, Gate, RunnerQT>;
+    using NoisyRunnerParameter = NoisyRunner::Parameter;
+
     StateSpace CreateStateSpace() const {
       return StateSpace(num_threads);
     }
diff --git a/pybind_interface/avx512/CMakeLists.txt b/pybind_interface/avx512/CMakeLists.txt
index 5ca2da26e..892d14c24 100644
--- a/pybind_interface/avx512/CMakeLists.txt
+++ b/pybind_interface/avx512/CMakeLists.txt
@@ -12,7 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-cmake_minimum_required(VERSION 3.31)
+cmake_minimum_required(VERSION 3.28)
 project(qsim)
 
 IF (WIN32)
diff --git a/pybind_interface/avx512/pybind_main_avx512.cpp b/pybind_interface/avx512/pybind_main_avx512.cpp
index b87ececb6..548bf40bf 100644
--- a/pybind_interface/avx512/pybind_main_avx512.cpp
+++ b/pybind_interface/avx512/pybind_main_avx512.cpp
@@ -15,6 +15,10 @@
 #include "pybind_main_avx512.h"
 
 #include "../../lib/formux.h"
+#include "../../lib/fuser_mqubit.h"
+#include "../../lib/gates_cirq.h"
+#include "../../lib/io.h"
+#include "../../lib/run_qsim.h"
 #include "../../lib/simulator_avx512.h"
 #include "../../lib/util_cpu.h"
 
@@ -33,6 +37,14 @@ namespace qsim {
     using Simulator = qsim::Simulator<For>;
     using StateSpace = Simulator::StateSpace;
 
+    using Gate = Cirq::GateCirq<float>;
+    using Runner = QSimRunner<IO, MultiQubitGateFuser<IO, Gate>, Factory>;
+    using RunnerQT =
+        QSimRunner<IO, MultiQubitGateFuser<IO, const Gate*>, Factory>;
+    using RunnerParameter = Runner::Parameter;
+    using NoisyRunner = qsim::QuantumTrajectorySimulator<IO, Gate, RunnerQT>;
+    using NoisyRunnerParameter = NoisyRunner::Parameter;
+
     StateSpace CreateStateSpace() const {
       return StateSpace(num_threads);
     }
diff --git a/pybind_interface/basic/CMakeLists.txt b/pybind_interface/basic/CMakeLists.txt
index ca9e6c263..3d5750887 100644
--- a/pybind_interface/basic/CMakeLists.txt
+++ b/pybind_interface/basic/CMakeLists.txt
@@ -12,7 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-cmake_minimum_required(VERSION 3.31)
+cmake_minimum_required(VERSION 3.28)
 project(qsim)
 
 if(WIN32)
diff --git a/pybind_interface/basic/pybind_main_basic.cpp b/pybind_interface/basic/pybind_main_basic.cpp
index 51fdc6081..9b37191b5 100644
--- a/pybind_interface/basic/pybind_main_basic.cpp
+++ b/pybind_interface/basic/pybind_main_basic.cpp
@@ -15,6 +15,10 @@
 #include "pybind_main_basic.h"
 
 #include "../../lib/formux.h"
+#include "../../lib/fuser_mqubit.h"
+#include "../../lib/gates_cirq.h"
+#include "../../lib/io.h"
+#include "../../lib/run_qsim.h"
 #include "../../lib/simulator_basic.h"
 #include "../../lib/util_cpu.h"
 
@@ -33,6 +37,14 @@ namespace qsim {
     using Simulator = qsim::Simulator<For>;
     using StateSpace = Simulator::StateSpace;
 
+    using Gate = Cirq::GateCirq<float>;
+    using Runner = QSimRunner<IO, MultiQubitGateFuser<IO, Gate>, Factory>;
+    using RunnerQT =
+        QSimRunner<IO, MultiQubitGateFuser<IO, const Gate*>, Factory>;
+    using RunnerParameter = Runner::Parameter;
+    using NoisyRunner = qsim::QuantumTrajectorySimulator<IO, Gate, RunnerQT>;
+    using NoisyRunnerParameter = NoisyRunner::Parameter;
+
     StateSpace CreateStateSpace() const {
       return StateSpace(num_threads);
     }
diff --git a/pybind_interface/cuda/CMakeLists.txt b/pybind_interface/cuda/CMakeLists.txt
index 8c57e3b81..24ebbe8e5 100644
--- a/pybind_interface/cuda/CMakeLists.txt
+++ b/pybind_interface/cuda/CMakeLists.txt
@@ -12,7 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-cmake_minimum_required(VERSION 3.31)
+cmake_minimum_required(VERSION 3.28)
 project(qsim LANGUAGES CXX CUDA)
 
 if(WIN32)
diff --git a/pybind_interface/cuda/pybind_main_cuda.cpp b/pybind_interface/cuda/pybind_main_cuda.cpp
index 74801408d..df9a293ee 100644
--- a/pybind_interface/cuda/pybind_main_cuda.cpp
+++ b/pybind_interface/cuda/pybind_main_cuda.cpp
@@ -14,15 +14,16 @@
 
 #include "pybind_main_cuda.h"
 
+#include "../../lib/fuser_mqubit.h"
+#include "../../lib/gates_cirq.h"
+#include "../../lib/io.h"
+#include "../../lib/run_qsim.h"
 #include "../../lib/simulator_cuda.h"
 
 namespace qsim {
   using Simulator = SimulatorCUDA<float>;
 
   struct Factory {
-    using Simulator = qsim::Simulator;
-    using StateSpace = Simulator::StateSpace;
-
     Factory(
       unsigned num_sim_threads,
       unsigned num_state_threads,
@@ -31,6 +32,18 @@ namespace qsim {
       ss_params.num_threads = num_state_threads;
       ss_params.num_dblocks = num_dblocks;
     }
+
+    using Simulator = qsim::Simulator;
+    using StateSpace = Simulator::StateSpace;
+
+    using Gate = Cirq::GateCirq<float>;
+    using Runner = QSimRunner<IO, MultiQubitGateFuser<IO, Gate>, Factory>;
+    using RunnerQT =
+        QSimRunner<IO, MultiQubitGateFuser<IO, const Gate*>, Factory>;
+    using RunnerParameter = Runner::Parameter;
+    using NoisyRunner = qsim::QuantumTrajectorySimulator<IO, Gate, RunnerQT>;
+    using NoisyRunnerParameter = NoisyRunner::Parameter;
+
     StateSpace CreateStateSpace() const {
       return StateSpace(ss_params);
     }
diff --git a/pybind_interface/custatevec/CMakeLists.txt b/pybind_interface/custatevec/CMakeLists.txt
index d5b280e72..eac992c2d 100644
--- a/pybind_interface/custatevec/CMakeLists.txt
+++ b/pybind_interface/custatevec/CMakeLists.txt
@@ -12,7 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-cmake_minimum_required(VERSION 3.31)
+cmake_minimum_required(VERSION 3.28)
 project(qsim LANGUAGES CXX CUDA)
 
 if(WIN32)
diff --git a/pybind_interface/custatevec/pybind_main_custatevec.cpp b/pybind_interface/custatevec/pybind_main_custatevec.cpp
index 1fe1d70bf..b4523c8ff 100644
--- a/pybind_interface/custatevec/pybind_main_custatevec.cpp
+++ b/pybind_interface/custatevec/pybind_main_custatevec.cpp
@@ -17,45 +17,55 @@
 
 #include "pybind_main_custatevec.h"
 
+#include "../../lib/fuser_mqubit.h"
+#include "../../lib/gates_cirq.h"
+#include "../../lib/io.h"
+#include "../../lib/run_qsim.h"
 #include "../../lib/simulator_custatevec.h"
 
 namespace qsim {
+  using Simulator = SimulatorCuStateVec<float>;
 
-using Simulator = SimulatorCuStateVec<float>;
+  struct Factory {
+    // num_sim_threads, num_state_threads and num_dblocks are unused, but kept
+    // for consistency with other factories.
+    Factory(unsigned num_sim_threads,
+            unsigned num_state_threads,
+            unsigned num_dblocks) {
+      ErrorCheck(cublasCreate(&cublas_handle));
+      ErrorCheck(custatevecCreate(&custatevec_handle));
+    }
 
-struct Factory {
-  using Simulator = qsim::Simulator;
-  using StateSpace = Simulator::StateSpace;
+    ~Factory() {
+      ErrorCheck(cublasDestroy(cublas_handle));
+      ErrorCheck(custatevecDestroy(custatevec_handle));
+    }
 
-  // num_sim_threads, num_state_threads and num_dblocks are unused, but kept
-  // for consistency with other factories.
-  Factory(unsigned num_sim_threads,
-          unsigned num_state_threads,
-          unsigned num_dblocks) {
-    ErrorCheck(cublasCreate(&cublas_handle));
-    ErrorCheck(custatevecCreate(&custatevec_handle));
-  }
+    using Simulator = qsim::Simulator;
+    using StateSpace = Simulator::StateSpace;
 
-  ~Factory() {
-    ErrorCheck(cublasDestroy(cublas_handle));
-    ErrorCheck(custatevecDestroy(custatevec_handle));
-  }
+    using Gate = Cirq::GateCirq<float>;
+    using Runner = QSimRunner<IO, MultiQubitGateFuser<IO, Gate>, Factory>;
+    using RunnerQT =
+        QSimRunner<IO, MultiQubitGateFuser<IO, const Gate*>, Factory>;
+    using RunnerParameter = Runner::Parameter;
+    using NoisyRunner = qsim::QuantumTrajectorySimulator<IO, Gate, RunnerQT>;
+    using NoisyRunnerParameter = NoisyRunner::Parameter;
 
-  StateSpace CreateStateSpace() const {
-    return StateSpace(cublas_handle, custatevec_handle);
-  }
+    StateSpace CreateStateSpace() const {
+      return StateSpace(cublas_handle, custatevec_handle);
+    }
 
-  Simulator CreateSimulator() const {
-    return Simulator(cublas_handle, custatevec_handle);
-  }
+    Simulator CreateSimulator() const {
+      return Simulator(cublas_handle, custatevec_handle);
+    }
 
-  cublasHandle_t cublas_handle;
-  custatevecHandle_t custatevec_handle;
-};
-
-inline void SetFlushToZeroAndDenormalsAreZeros() {}
-inline void ClearFlushToZeroAndDenormalsAreZeros() {}
+    cublasHandle_t cublas_handle;
+    custatevecHandle_t custatevec_handle;
+  };
 
+  inline void SetFlushToZeroAndDenormalsAreZeros() {}
+  inline void ClearFlushToZeroAndDenormalsAreZeros() {}
 }
 
 #include "../pybind_main.cpp"
diff --git a/pybind_interface/decide/CMakeLists.txt b/pybind_interface/decide/CMakeLists.txt
index 8b6413821..dce80e51d 100644
--- a/pybind_interface/decide/CMakeLists.txt
+++ b/pybind_interface/decide/CMakeLists.txt
@@ -12,7 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-cmake_minimum_required(VERSION 3.31)
+cmake_minimum_required(VERSION 3.28)
 project(qsim LANGUAGES CXX)
 
 include(CheckLanguage)
diff --git a/pybind_interface/hip/CMakeLists.txt b/pybind_interface/hip/CMakeLists.txt
index 6956b01b9..4b3a7cd13 100644
--- a/pybind_interface/hip/CMakeLists.txt
+++ b/pybind_interface/hip/CMakeLists.txt
@@ -12,7 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-cmake_minimum_required(VERSION 3.31)
+cmake_minimum_required(VERSION 3.28)
 project(qsim LANGUAGES CXX HIP)
 
 if(WIN32)
diff --git a/pybind_interface/hip/pybind_main_hip.cpp b/pybind_interface/hip/pybind_main_hip.cpp
index d6ee181f9..98a9e3fcd 100644
--- a/pybind_interface/hip/pybind_main_hip.cpp
+++ b/pybind_interface/hip/pybind_main_hip.cpp
@@ -14,21 +14,33 @@
 
 #include "pybind_main_hip.h"
 
+#include "../../lib/fuser_mqubit.h"
+#include "../../lib/gates_cirq.h"
+#include "../../lib/io.h"
+#include "../../lib/run_qsim.h"
 #include "../../lib/simulator_cuda.h"
 
 namespace qsim {
   using Simulator = SimulatorCUDA<float>;
 
   struct Factory {
-    using Simulator = qsim::Simulator;
-    using StateSpace = Simulator::StateSpace;
-
     Factory(
       unsigned num_sim_threads,
       unsigned num_state_threads,
       unsigned num_dblocks
     ) : ss_params{num_state_threads, num_dblocks} {}
 
+    using Simulator = qsim::Simulator;
+    using StateSpace = Simulator::StateSpace;
+
+    using Gate = Cirq::GateCirq<float>;
+    using Runner = QSimRunner<IO, MultiQubitGateFuser<IO, Gate>, Factory>;
+    using RunnerQT =
+        QSimRunner<IO, MultiQubitGateFuser<IO, const Gate*>, Factory>;
+    using RunnerParameter = Runner::Parameter;
+    using NoisyRunner = qsim::QuantumTrajectorySimulator<IO, Gate, RunnerQT>;
+    using NoisyRunnerParameter = NoisyRunner::Parameter;
+
     StateSpace CreateStateSpace() const {
       return StateSpace(ss_params);
     }
diff --git a/pybind_interface/pybind_main.cpp b/pybind_interface/pybind_main.cpp
index b820ca190..5f6644469 100644
--- a/pybind_interface/pybind_main.cpp
+++ b/pybind_interface/pybind_main.cpp
@@ -26,10 +26,8 @@
 #include "../lib/expect.h"
 #include "../lib/formux.h"
 #include "../lib/fuser_mqubit.h"
-#include "../lib/gates_qsim.h"
 #include "../lib/io.h"
 #include "../lib/qtrajectory.h"
-#include "../lib/run_qsim.h"
 #include "../lib/run_qsimh.h"
 
 using namespace qsim;
@@ -46,17 +44,17 @@ T parseOptions(const py::dict &options, const char *key) {
   return value.cast<T>();
 }
 
-Circuit<Cirq::GateCirq<float>> getCircuit(const py::dict &options) {
+Circuit<Factory::Gate> getCircuit(const py::dict &options) {
   try {
-    return options["c\0"].cast<Circuit<Cirq::GateCirq<float>>>();
+    return options["c\0"].cast<Circuit<Factory::Gate>>();
   } catch (const std::invalid_argument &exp) {
     throw;
   }
 }
 
-NoisyCircuit<Cirq::GateCirq<float>> getNoisyCircuit(const py::dict &options) {
+NoisyCircuit<Factory::Gate> getNoisyCircuit(const py::dict &options) {
   try {
-    return options["c\0"].cast<NoisyCircuit<Cirq::GateCirq<float>>>();
+    return options["c\0"].cast<NoisyCircuit<Factory::Gate>>();
   } catch (const std::invalid_argument &exp) {
     throw;
   }
@@ -81,7 +79,7 @@ std::vector<Bitstring> getBitstrings(const py::dict &options, int num_qubits) {
 
 }  // namespace
 
-Cirq::GateCirq<float> create_gate(const qsim::Cirq::GateKind gate_kind,
+Factory::Gate create_gate(const qsim::Cirq::GateKind gate_kind,
                                   const unsigned time,
                                   const std::vector<unsigned>& qubits,
                                   const std::map<std::string, float>& params) {
@@ -203,7 +201,7 @@ Cirq::GateCirq<float> create_gate(const qsim::Cirq::GateKind gate_kind,
       return Cirq::CCX<float>::Create(time, qubits[0], qubits[1], qubits[2]);
     case Cirq::kMeasurement: {
       std::vector<unsigned> qubits_ = qubits;
-      return gate::Measurement<Cirq::GateCirq<float>>::Create(
+      return gate::Measurement<Factory::Gate>::Create(
         time, std::move(qubits_));
       }
     // Matrix gates are handled in the add_matrix methods below.
@@ -212,7 +210,7 @@ Cirq::GateCirq<float> create_gate(const qsim::Cirq::GateKind gate_kind,
   }
 }
 
-Cirq::GateCirq<float> create_diagonal_gate(const unsigned time,
+Factory::Gate create_diagonal_gate(const unsigned time,
                                            const std::vector<unsigned>& qubits,
                                            const std::vector<float>& angles) {
   switch (qubits.size()) {
@@ -228,7 +226,7 @@ Cirq::GateCirq<float> create_diagonal_gate(const unsigned time,
   }
 }
 
-Cirq::GateCirq<float> create_matrix_gate(const unsigned time,
+Factory::Gate create_matrix_gate(const unsigned time,
                                          const std::vector<unsigned>& qubits,
                                          const std::vector<float>& matrix) {
   switch (qubits.size()) {
@@ -253,26 +251,26 @@ Cirq::GateCirq<float> create_matrix_gate(const unsigned time,
 void add_gate(const qsim::Cirq::GateKind gate_kind, const unsigned time,
               const std::vector<unsigned>& qubits,
               const std::map<std::string, float>& params,
-              Circuit<Cirq::GateCirq<float>>* circuit) {
+              Circuit<Factory::Gate>* circuit) {
   circuit->gates.push_back(create_gate(gate_kind, time, qubits, params));
 }
 
 void add_diagonal_gate(const unsigned time, const std::vector<unsigned>& qubits,
                        const std::vector<float>& angles,
-                       Circuit<Cirq::GateCirq<float>>* circuit) {
+                       Circuit<Factory::Gate>* circuit) {
   circuit->gates.push_back(create_diagonal_gate(time, qubits, angles));
 }
 
 void add_matrix_gate(const unsigned time,
                      const std::vector<unsigned>& qubits,
                      const std::vector<float>& matrix,
-                     Circuit<Cirq::GateCirq<float>>* circuit) {
+                     Circuit<Factory::Gate>* circuit) {
   circuit->gates.push_back(create_matrix_gate(time, qubits, matrix));
 }
 
 void control_last_gate(const std::vector<unsigned>& qubits,
                        const std::vector<unsigned>& values,
-                       Circuit<Cirq::GateCirq<float>>* circuit) {
+                       Circuit<Factory::Gate>* circuit) {
   MakeControlledGate(qubits, values, circuit->gates.back());
 }
 
@@ -288,7 +286,7 @@ Channel<Gate> create_single_gate_channel(Gate gate) {
 void add_gate_channel(const qsim::Cirq::GateKind gate_kind, const unsigned time,
                       const std::vector<unsigned>& qubits,
                       const std::map<std::string, float>& params,
-                      NoisyCircuit<Cirq::GateCirq<float>>* ncircuit) {
+                      NoisyCircuit<Factory::Gate>* ncircuit) {
   ncircuit->channels.push_back(create_single_gate_channel(
     create_gate(gate_kind, time, qubits, params)));
 }
@@ -296,7 +294,7 @@ void add_gate_channel(const qsim::Cirq::GateKind gate_kind, const unsigned time,
 void add_diagonal_gate_channel(const unsigned time,
                                const std::vector<unsigned>& qubits,
                                const std::vector<float>& angles,
-                               NoisyCircuit<Cirq::GateCirq<float>>* ncircuit) {
+                               NoisyCircuit<Factory::Gate>* ncircuit) {
   ncircuit->channels.push_back(create_single_gate_channel(
     create_diagonal_gate(time, qubits, angles)));
 }
@@ -304,19 +302,19 @@ void add_diagonal_gate_channel(const unsigned time,
 void add_matrix_gate_channel(const unsigned time,
                              const std::vector<unsigned>& qubits,
                              const std::vector<float>& matrix,
-                             NoisyCircuit<Cirq::GateCirq<float>>* ncircuit) {
+                             NoisyCircuit<Factory::Gate>* ncircuit) {
   ncircuit->channels.push_back(create_single_gate_channel(
     create_matrix_gate(time, qubits, matrix)));
 }
 
 void control_last_gate_channel(const std::vector<unsigned>& qubits,
                                const std::vector<unsigned>& values,
-                               NoisyCircuit<Cirq::GateCirq<float>>* ncircuit) {
+                               NoisyCircuit<Factory::Gate>* ncircuit) {
   if (ncircuit->channels.back().size() > 1) {
     throw std::invalid_argument(
         "Control cannot be added to noisy channels.");
   }
-  for (Cirq::GateCirq<float>& op : ncircuit->channels.back()[0].ops) {
+  for (Factory::Gate& op : ncircuit->channels.back()[0].ops) {
     MakeControlledGate(qubits, values, op);
   }
 }
@@ -325,10 +323,9 @@ void add_channel(const unsigned time,
                  const std::vector<unsigned>& qubits,
                  const std::vector<std::tuple<float, std::vector<float>, bool>>&
                      prob_matrix_unitary_triples,
-                 NoisyCircuit<Cirq::GateCirq<float>>* ncircuit) {
+                 NoisyCircuit<Factory::Gate>* ncircuit) {
   // Adds a channel to the noisy circuit.
-  using Gate = Cirq::GateCirq<float>;
-  Channel<Gate> channel;
+  Channel<Factory::Gate> channel;
   // prob_matrix_unitary_triples contains triples with these elements:
   //   0. The lower-bound probability of applying the matrix.
   //   1. The matrix to be applied.
@@ -337,9 +334,9 @@ void add_channel(const unsigned time,
     const float prob = std::get<0>(triple);
     const std::vector<float>& mat = std::get<1>(triple);
     bool is_unitary = std::get<2>(triple);
-    Gate gate = create_matrix_gate(time, qubits, mat);
-    channel.emplace_back(KrausOperator<Gate>{
-      KrausOperator<Gate>::kNormal, is_unitary, prob, {gate}
+    Factory::Gate gate = create_matrix_gate(time, qubits, mat);
+    channel.emplace_back(KrausOperator<Factory::Gate>{
+      KrausOperator<Factory::Gate>::kNormal, is_unitary, prob, {gate}
     });
     if (!is_unitary) {
       channel.back().CalculateKdKMatrix();
@@ -350,21 +347,21 @@ void add_channel(const unsigned time,
 
 void add_gate_to_opstring(const Cirq::GateKind gate_kind,
                           const std::vector<unsigned>& qubits,
-                          OpString<Cirq::GateCirq<float>>* opstring) {
+                          OpString<Factory::Gate>* opstring) {
   static std::map<std::string, float> params;
   opstring->ops.push_back(create_gate(gate_kind, 0, qubits, params));
 }
 
 void add_matrix_gate_to_opstring(const std::vector<unsigned>& qubits,
                                  const std::vector<float>& matrix,
-                                 OpString<Cirq::GateCirq<float>>* opstring) {
+                                 OpString<Factory::Gate>* opstring) {
   opstring->ops.push_back(create_matrix_gate(0, qubits, matrix));
 }
 
 // Methods for simulating amplitudes.
 
 std::vector<std::complex<float>> qsim_simulate(const py::dict &options) {
-  Circuit<Cirq::GateCirq<float>> circuit;
+  Circuit<Factory::Gate> circuit;
   std::vector<Bitstring> bitstrings;
   try {
     circuit = getCircuit(options);
@@ -375,6 +372,8 @@ std::vector<std::complex<float>> qsim_simulate(const py::dict &options) {
   }
 
   using Simulator = Factory::Simulator;
+  using Runner = Factory::Runner;
+  using RunnerParameter = Factory::RunnerParameter;
   using StateSpace = Simulator::StateSpace;
   using State = StateSpace::State;
 
@@ -390,16 +389,13 @@ std::vector<std::complex<float>> qsim_simulate(const py::dict &options) {
     }
   };
 
-  using Runner = QSimRunner<IO, MultiQubitGateFuser<IO, Cirq::GateCirq<float>>,
-                            Factory>;
-
   bool use_gpu;
   bool denormals_are_zeros;
   unsigned gpu_mode;
   unsigned num_sim_threads = 0;
   unsigned num_state_threads = 0;
   unsigned num_dblocks = 0;
-  Runner::Parameter param;
+  RunnerParameter param;
   try {
     use_gpu = parseOptions<unsigned>(options, "g\0");
     gpu_mode = parseOptions<unsigned>(options, "gmode\0");
@@ -431,7 +427,7 @@ std::vector<std::complex<float>> qsim_simulate(const py::dict &options) {
 }
 
 std::vector<std::complex<float>> qtrajectory_simulate(const py::dict &options) {
-  NoisyCircuit<Cirq::GateCirq<float>> ncircuit;
+  NoisyCircuit<Factory::Gate> ncircuit;
   unsigned num_qubits;
   std::vector<Bitstring> bitstrings;
   try {
@@ -444,6 +440,8 @@ std::vector<std::complex<float>> qtrajectory_simulate(const py::dict &options) {
   }
 
   using Simulator = Factory::Simulator;
+  using NoisyRunner = Factory::NoisyRunner;
+  using NoisyRunnerParameter = Factory::NoisyRunnerParameter;
   using StateSpace = Simulator::StateSpace;
   using State = StateSpace::State;
 
@@ -451,11 +449,7 @@ std::vector<std::complex<float>> qtrajectory_simulate(const py::dict &options) {
   std::vector<std::complex<float>> amplitudes;
   amplitudes.reserve(bitstrings.size());
 
-  using Runner = qsim::QuantumTrajectorySimulator<IO, Cirq::GateCirq<float>,
-                                                  MultiQubitGateFuser,
-                                                  Simulator>;
-
-  Runner::Parameter param;
+  NoisyRunnerParameter param;
   bool use_gpu;
   bool denormals_are_zeros;
   unsigned gpu_mode;
@@ -487,7 +481,7 @@ std::vector<std::complex<float>> qtrajectory_simulate(const py::dict &options) {
   StateSpace state_space = factory.CreateStateSpace();
 
   auto measure = [&bitstrings, &amplitudes, &state_space](
-                  unsigned k, const State &state, Runner::Stat& stat) {
+                  unsigned k, const State &state, NoisyRunner::Stat& stat) {
     for (const auto &b : bitstrings) {
       amplitudes.push_back(state_space.GetAmpl(state, b));
     }
@@ -499,26 +493,54 @@ std::vector<std::complex<float>> qtrajectory_simulate(const py::dict &options) {
     ClearFlushToZeroAndDenormalsAreZeros();
   }
 
-  if (!Runner::RunBatch(param, ncircuit, seed, seed + 1, state_space,
-                        simulator, measure)) {
+  if (!NoisyRunner::RunBatch(param, ncircuit, seed, seed + 1, state_space,
+                             simulator, measure)) {
     IO::errorf("qtrajectory simulation of the circuit errored out.\n");
     return {};
   }
   return amplitudes;
 }
 
+template <bool requires_copy_to_host>
+struct Releaser;
+
+template <>
+struct Releaser<true> {
+  template <typename StateSpace, typename State>
+  static py::array_t<float> release(
+      const StateSpace& state_space, State& state, uint64_t fsv_size) {
+    auto* fsv = new float[state_space.MinSize(state.num_qubits())];
+    state_space.Copy(state, fsv);
+    // Cast on delete to silence warnings.
+    auto capsule = py::capsule(
+      fsv, [](void *data) { delete [] (float*)data; });
+    return py::array_t<float>(fsv_size, fsv, capsule);
+  }
+};
+
+template <>
+struct Releaser<false> {
+  template <typename StateSpace, typename State>
+  static py::array_t<float> release(
+      const StateSpace& state_space, State& state, uint64_t fsv_size) {
+    float* fsv = state.release();
+    auto capsule = py::capsule(
+        fsv, [](void *data) { detail::free(data); });
+    return py::array_t<float>(fsv_size, fsv, capsule);
+  }
+};
+
 // Helper class for simulating circuits of all types.
 class SimulatorHelper {
  public:
   using Simulator = Factory::Simulator;
+  using Runner = Factory::Runner;
+  using RunnerParameter = Factory::RunnerParameter;
+  using NoisyRunner = Factory::NoisyRunner;
+  using NoisyRunnerParameter = Factory::NoisyRunnerParameter;
   using StateSpace = Factory::StateSpace;
   using State = StateSpace::State;
 
-  using Gate = Cirq::GateCirq<float>;
-  using Runner = QSimRunner<IO, MultiQubitGateFuser<IO, Gate>, Factory>;
-  using NoisyRunner = qsim::QuantumTrajectorySimulator<
-      IO, Gate, MultiQubitGateFuser, Simulator>;
-
   SimulatorHelper() = delete;
 
   template <typename StateType>
@@ -544,7 +566,7 @@ class SimulatorHelper {
   static std::vector<std::complex<double>> simulate_expectation_values(
       const py::dict &options,
       const std::vector<std::tuple<
-                            std::vector<OpString<Cirq::GateCirq<float>>>,
+                            std::vector<OpString<Factory::Gate>>,
                             unsigned>>& opsums_and_qubit_counts,
       bool is_noisy, const StateType& input_state) {
     auto helper = SimulatorHelper(options, is_noisy);
@@ -582,8 +604,8 @@ class SimulatorHelper {
   simulate_moment_expectation_values(
       const py::dict &options,
       const std::vector<std::tuple<uint64_t, std::vector<
-        std::tuple<std::vector<OpString<Cirq::GateCirq<float>>>, unsigned>
-      >>>& opsums_and_qubit_counts,
+          std::tuple<std::vector<OpString<Factory::Gate>>,
+                     unsigned>>>>& opsums_and_qubit_counts,
       bool is_noisy, const StateType& input_state) {
     auto helper = SimulatorHelper(options, is_noisy);
     if (!helper.is_valid) {
@@ -707,16 +729,16 @@ class SimulatorHelper {
     state_space.NormalToInternalOrder(state);
   }
 
-  Runner::Parameter get_params() const {
-    Runner::Parameter params;
+  RunnerParameter get_params() const {
+    RunnerParameter params;
     params.max_fused_size = max_fused_size;
     params.verbosity = verbosity;
     params.seed = seed;
     return params;
   }
 
-  NoisyRunner::Parameter get_noisy_params() const {
-    NoisyRunner::Parameter params;
+  NoisyRunnerParameter get_noisy_params() const {
+    NoisyRunnerParameter params;
     params.max_fused_size = max_fused_size;
     params.verbosity = verbosity;
     return params;
@@ -759,9 +781,9 @@ class SimulatorHelper {
         seed, state_space, simulator, state, stat
       );
     } else {
-      Circuit<Gate> subcircuit;
+      Circuit<Factory::Gate> subcircuit;
       subcircuit.num_qubits = circuit.num_qubits;
-      subcircuit.gates = std::vector<Gate>(
+      subcircuit.gates = std::vector<Factory::Gate>(
         circuit.gates.begin() + begin,
         circuit.gates.begin() + end
       );
@@ -780,27 +802,17 @@ class SimulatorHelper {
     StateSpace state_space = factory.CreateStateSpace();
     state_space.InternalToNormalOrder(state);
     uint64_t fsv_size = 2 * (uint64_t{1} << num_qubits);
-    if (state.requires_copy_to_host()) {
-      auto* fsv = new float[state_space.MinSize(state.num_qubits())];
-      state_space.Copy(state, fsv);
-      // Cast on delete to silence warnings.
-      auto capsule = py::capsule(
-        fsv, [](void *data) { delete [] (float*)data; });
-      return py::array_t<float>(fsv_size, fsv, capsule);
-    } else {
-      float* fsv = state.release();
-      auto capsule = py::capsule(
-          fsv, [](void *data) { detail::free(data); });
-      return py::array_t<float>(fsv_size, fsv, capsule);
-    }
+
+    return Releaser<State::requires_copy_to_host()>::release(
+        state_space, state, fsv_size);
   }
 
   std::vector<std::complex<double>> get_expectation_value(
-      const std::vector<std::tuple<std::vector<OpString<Gate>>,
+      const std::vector<std::tuple<std::vector<OpString<Factory::Gate>>,
                                    unsigned>>& opsums_and_qubit_counts) {
     Simulator simulator = factory.CreateSimulator();
     StateSpace state_space = factory.CreateStateSpace();
-    using Fuser = MultiQubitGateFuser<IO, Gate>;
+    using Fuser = MultiQubitGateFuser<IO, Factory::Gate>;
 
     std::vector<std::complex<double>> results;
     results.reserve(opsums_and_qubit_counts.size());
@@ -829,8 +841,8 @@ class SimulatorHelper {
 
   bool is_noisy;
   // Only one of these will be populated, as specified by is_noisy.
-  Circuit<Gate> circuit;
-  NoisyCircuit<Gate> ncircuit;
+  Circuit<Factory::Gate> circuit;
+  NoisyCircuit<Factory::Gate> ncircuit;
 
   Factory factory;
   State state;
@@ -878,7 +890,7 @@ py::array_t<float> qtrajectory_simulate_fullstate(
 std::vector<std::complex<double>> qsim_simulate_expectation_values(
     const py::dict &options,
     const std::vector<std::tuple<
-                          std::vector<OpString<Cirq::GateCirq<float>>>,
+                          std::vector<OpString<Factory::Gate>>,
                           unsigned>>& opsums_and_qubit_counts,
     uint64_t input_state) {
   return SimulatorHelper::simulate_expectation_values(
@@ -888,7 +900,7 @@ std::vector<std::complex<double>> qsim_simulate_expectation_values(
 std::vector<std::complex<double>> qsim_simulate_expectation_values(
     const py::dict &options,
     const std::vector<std::tuple<
-                          std::vector<OpString<Cirq::GateCirq<float>>>,
+                          std::vector<OpString<Factory::Gate>>,
                           unsigned>>& opsums_and_qubit_counts,
     const py::array_t<float> &input_vector) {
   return SimulatorHelper::simulate_expectation_values(
@@ -899,7 +911,7 @@ std::vector<std::vector<std::complex<double>>>
 qsim_simulate_moment_expectation_values(
     const py::dict &options,
     const std::vector<std::tuple<uint64_t, std::vector<
-      std::tuple<std::vector<OpString<Cirq::GateCirq<float>>>, unsigned>
+      std::tuple<std::vector<OpString<Factory::Gate>>, unsigned>
     >>>& opsums_and_qubit_counts,
     uint64_t input_state) {
   return SimulatorHelper::simulate_moment_expectation_values(
@@ -910,7 +922,7 @@ std::vector<std::vector<std::complex<double>>>
 qsim_simulate_moment_expectation_values(
     const py::dict &options,
     const std::vector<std::tuple<uint64_t, std::vector<
-      std::tuple<std::vector<OpString<Cirq::GateCirq<float>>>, unsigned>
+      std::tuple<std::vector<OpString<Factory::Gate>>, unsigned>
     >>>& opsums_and_qubit_counts,
     const py::array_t<float> &input_vector) {
   return SimulatorHelper::simulate_moment_expectation_values(
@@ -920,7 +932,7 @@ qsim_simulate_moment_expectation_values(
 std::vector<std::complex<double>> qtrajectory_simulate_expectation_values(
     const py::dict &options,
     const std::vector<std::tuple<
-                          std::vector<OpString<Cirq::GateCirq<float>>>,
+                          std::vector<OpString<Factory::Gate>>,
                           unsigned>>& opsums_and_qubit_counts,
     uint64_t input_state) {
   return SimulatorHelper::simulate_expectation_values(
@@ -930,7 +942,7 @@ std::vector<std::complex<double>> qtrajectory_simulate_expectation_values(
 std::vector<std::complex<double>> qtrajectory_simulate_expectation_values(
     const py::dict &options,
     const std::vector<std::tuple<
-                          std::vector<OpString<Cirq::GateCirq<float>>>,
+                          std::vector<OpString<Factory::Gate>>,
                           unsigned>>& opsums_and_qubit_counts,
     const py::array_t<float> &input_vector) {
   return SimulatorHelper::simulate_expectation_values(
@@ -941,7 +953,7 @@ std::vector<std::vector<std::complex<double>>>
 qtrajectory_simulate_moment_expectation_values(
     const py::dict &options,
     const std::vector<std::tuple<uint64_t, std::vector<
-      std::tuple<std::vector<OpString<Cirq::GateCirq<float>>>, unsigned>
+      std::tuple<std::vector<OpString<Factory::Gate>>, unsigned>
     >>>& opsums_and_qubit_counts,
     uint64_t input_state) {
   return SimulatorHelper::simulate_moment_expectation_values(
@@ -952,7 +964,7 @@ std::vector<std::vector<std::complex<double>>>
 qtrajectory_simulate_moment_expectation_values(
     const py::dict &options,
     const std::vector<std::tuple<uint64_t, std::vector<
-      std::tuple<std::vector<OpString<Cirq::GateCirq<float>>>, unsigned>
+      std::tuple<std::vector<OpString<Factory::Gate>>, unsigned>
     >>>& opsums_and_qubit_counts,
     const py::array_t<float> &input_vector) {
   return SimulatorHelper::simulate_moment_expectation_values(
@@ -972,7 +984,7 @@ std::vector<uint64_t> qtrajectory_sample_final(
 }
 
 std::vector<unsigned> qsim_sample(const py::dict &options) {
-  Circuit<Cirq::GateCirq<float>> circuit;
+  Circuit<Factory::Gate> circuit;
   try {
     circuit = getCircuit(options);
   } catch (const std::invalid_argument &exp) {
@@ -981,11 +993,11 @@ std::vector<unsigned> qsim_sample(const py::dict &options) {
   }
 
   using Simulator = Factory::Simulator;
+  using Runner = Factory::Runner;
+  using RunnerParameter = Factory::RunnerParameter;
   using StateSpace = Simulator::StateSpace;
   using State = StateSpace::State;
   using MeasurementResult = StateSpace::MeasurementResult;
-  using Runner = QSimRunner<IO, MultiQubitGateFuser<IO, Cirq::GateCirq<float>>,
-                            Factory>;
 
   bool use_gpu;
   bool denormals_are_zeros;
@@ -993,7 +1005,7 @@ std::vector<unsigned> qsim_sample(const py::dict &options) {
   unsigned num_sim_threads = 0;
   unsigned num_state_threads = 0;
   unsigned num_dblocks = 0;
-  Runner::Parameter param;
+  RunnerParameter param;
   try {
     use_gpu = parseOptions<unsigned>(options, "g\0");
     gpu_mode = parseOptions<unsigned>(options, "gmode\0");
@@ -1038,7 +1050,7 @@ std::vector<unsigned> qsim_sample(const py::dict &options) {
 }
 
 std::vector<unsigned> qtrajectory_sample(const py::dict &options) {
-  NoisyCircuit<Cirq::GateCirq<float>> ncircuit;
+  NoisyCircuit<Factory::Gate> ncircuit;
   try {
     ncircuit = getNoisyCircuit(options);
   } catch (const std::invalid_argument &exp) {
@@ -1047,13 +1059,12 @@ std::vector<unsigned> qtrajectory_sample(const py::dict &options) {
   }
 
   using Simulator = Factory::Simulator;
+  using NoisyRunner = Factory::NoisyRunner;
+  using NoisyRunnerParameter = Factory::NoisyRunnerParameter;
   using StateSpace = Simulator::StateSpace;
   using State = StateSpace::State;
-  using Runner = qsim::QuantumTrajectorySimulator<IO, Cirq::GateCirq<float>,
-                                                  MultiQubitGateFuser,
-                                                  Simulator>;
 
-  Runner::Parameter param;
+  NoisyRunnerParameter param;
   bool use_gpu;
   bool denormals_are_zeros;
   unsigned gpu_mode;
@@ -1088,7 +1099,7 @@ std::vector<unsigned> qtrajectory_sample(const py::dict &options) {
   std::vector<std::vector<unsigned>> results;
 
   auto measure = [&results, &ncircuit](
-                  unsigned k, const State& state, Runner::Stat& stat) {
+                  unsigned k, const State& state, NoisyRunner::Stat& stat) {
     // Converts stat (which matches the MeasurementResult 'bits' field) into
     // bitstrings matching the MeasurementResult 'bitstring' field.
     unsigned idx = 0;
@@ -1116,8 +1127,8 @@ std::vector<unsigned> qtrajectory_sample(const py::dict &options) {
     ClearFlushToZeroAndDenormalsAreZeros();
   }
 
-  if (!Runner::RunBatch(param, ncircuit, seed, seed + 1,
-                        state_space, simulator, measure)) {
+  if (!NoisyRunner::RunBatch(param, ncircuit, seed, seed + 1,
+                             state_space, simulator, measure)) {
     IO::errorf("qtrajectory sampling of the circuit errored out.\n");
     return {};
   }
@@ -1132,11 +1143,11 @@ std::vector<unsigned> qtrajectory_sample(const py::dict &options) {
 // Method for running the hybrid simulator.
 
 std::vector<std::complex<float>> qsimh_simulate(const py::dict &options) {
-  using HybridSimulator = HybridSimulator<IO, Cirq::GateCirq<float>,
+  using HybridSimulator = HybridSimulator<IO, Factory::Gate,
                                           MultiQubitGateFuser, For>;
   using Runner = QSimHRunner<IO, HybridSimulator>;
 
-  Circuit<Cirq::GateCirq<float>> circuit;
+  Circuit<Factory::Gate> circuit;
   std::vector<Bitstring> bitstrings;
   Runner::Parameter param;
   py::list dense_parts;
diff --git a/pybind_interface/sse/CMakeLists.txt b/pybind_interface/sse/CMakeLists.txt
index e9ddea57a..e31360a2e 100644
--- a/pybind_interface/sse/CMakeLists.txt
+++ b/pybind_interface/sse/CMakeLists.txt
@@ -12,7 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-cmake_minimum_required(VERSION 3.31)
+cmake_minimum_required(VERSION 3.28)
 project(qsim)
 
 IF (WIN32)
diff --git a/pybind_interface/sse/pybind_main_sse.cpp b/pybind_interface/sse/pybind_main_sse.cpp
index 674aecad4..f04925135 100644
--- a/pybind_interface/sse/pybind_main_sse.cpp
+++ b/pybind_interface/sse/pybind_main_sse.cpp
@@ -15,6 +15,10 @@
 #include "pybind_main_sse.h"
 
 #include "../../lib/formux.h"
+#include "../../lib/fuser_mqubit.h"
+#include "../../lib/gates_cirq.h"
+#include "../../lib/io.h"
+#include "../../lib/run_qsim.h"
 #include "../../lib/simulator_sse.h"
 #include "../../lib/util_cpu.h"
 
@@ -33,6 +37,14 @@ namespace qsim {
     using Simulator = qsim::Simulator<For>;
     using StateSpace = Simulator::StateSpace;
 
+    using Gate = Cirq::GateCirq<float>;
+    using Runner = QSimRunner<IO, MultiQubitGateFuser<IO, Gate>, Factory>;
+    using RunnerQT =
+        QSimRunner<IO, MultiQubitGateFuser<IO, const Gate*>, Factory>;
+    using RunnerParameter = Runner::Parameter;
+    using NoisyRunner = qsim::QuantumTrajectorySimulator<IO, Gate, RunnerQT>;
+    using NoisyRunnerParameter = NoisyRunner::Parameter;
+
     StateSpace CreateStateSpace() const {
       return StateSpace(num_threads);
     }
diff --git a/tests/BUILD b/tests/BUILD
index edf5be81f..b42455884 100644
--- a/tests/BUILD
+++ b/tests/BUILD
@@ -77,6 +77,7 @@ cc_test(
         "//lib:gates_cirq",
         "//lib:io",
         "//lib:qtrajectory",
+        "//lib:run_qsim",
         "//lib:simulator",
     ],
     copts = select({
@@ -260,6 +261,10 @@ cc_test(
     deps = [
         ":qtrajectory_testfixture",
         "@com_google_googletest//:gtest_main",
+        "//lib:fuser_mqubit",
+        "//lib:gates_cirq",
+        "//lib:io",
+        "//lib:run_qsim",
         "//lib:seqfor",
         "//lib:simulator_avx",
     ],
diff --git a/tests/channels_cirq_test.cc b/tests/channels_cirq_test.cc
index dedb922ef..e82fcb3cf 100644
--- a/tests/channels_cirq_test.cc
+++ b/tests/channels_cirq_test.cc
@@ -25,18 +25,36 @@
 #include "../lib/gates_cirq.h"
 #include "../lib/io.h"
 #include "../lib/qtrajectory.h"
+#include "../lib/run_qsim.h"
 #include "../lib/simmux.h"
 
 namespace qsim {
 
 namespace {
 
-using StateSpace = Simulator<For>::StateSpace;
+template <typename For>
+struct Factory {
+  using Simulator = qsim::Simulator<For>;
+  using StateSpace = typename Simulator::StateSpace;
+  using fp_type = typename StateSpace::fp_type;
+
+  StateSpace CreateStateSpace() const {
+    return StateSpace(1);
+  }
+
+  Simulator CreateSimulator() const {
+    return Simulator(1);
+  }
+};
+
+using Simulator = Factory<For>::Simulator;
+using StateSpace = Simulator::StateSpace;
 using State = StateSpace::State;
 using fp_type = StateSpace::fp_type;
 using Gate = Cirq::GateCirq<fp_type>;
-using QTSimulator = QuantumTrajectorySimulator<IO, Gate, MultiQubitGateFuser,
-                                               Simulator<For>>;
+using Fuser = MultiQubitGateFuser<IO, const Gate*>;
+using Runner = QSimRunner<IO, Fuser, Factory<For>>;
+using QTSimulator = QuantumTrajectorySimulator<IO, Gate, Runner>;
 
 Circuit<Gate> CleanCircuit() {
   using Hd = Cirq::H<fp_type>;
@@ -79,8 +97,8 @@ void RunBatch(const NoisyCircuit<Gate>& ncircuit,
   QTSimulator::Parameter param;
   param.collect_mea_stat = true;
 
-  Simulator<For> simulator(num_threads);
-  Simulator<For>::StateSpace state_space(num_threads);
+  Simulator simulator(num_threads);
+  StateSpace state_space(num_threads);
 
   EXPECT_TRUE(QTSimulator::RunBatch(param, ncircuit, 0, num_reps, state_space,
                                     simulator, measure, histogram));
diff --git a/tests/qtrajectory_avx_test.cc b/tests/qtrajectory_avx_test.cc
index 2dc126d3e..983e899ff 100644
--- a/tests/qtrajectory_avx_test.cc
+++ b/tests/qtrajectory_avx_test.cc
@@ -16,6 +16,10 @@
 
 #include "gtest/gtest.h"
 
+#include "../lib/fuser_mqubit.h"
+#include "../lib/gates_cirq.h"
+#include "../lib/io.h"
+#include "../lib/run_qsim.h"
 #include "../lib/seqfor.h"
 #include "../lib/simulator_avx.h"
 
@@ -37,31 +41,45 @@ struct Factory {
 };
 
 TEST(QTrajectoryAVXTest, BitFlip) {
-  TestBitFlip(qsim::Factory<SequentialFor>());
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory<SequentialFor>>;
+  TestBitFlip<Runner>(Factory<SequentialFor>());
 }
 
 TEST(QTrajectoryAVXTest, GenDump) {
-  TestGenDump(qsim::Factory<SequentialFor>());
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory<SequentialFor>>;
+  TestGenDump<Runner>(Factory<SequentialFor>());
 }
 
 TEST(QTrajectoryAVXTest, ReusingResults) {
-  TestReusingResults(qsim::Factory<SequentialFor>());
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory<SequentialFor>>;
+  TestReusingResults<Runner>(Factory<SequentialFor>());
 }
 
 TEST(QTrajectoryAVXTest, CollectKopStat) {
-  TestCollectKopStat(qsim::Factory<SequentialFor>());
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory<SequentialFor>>;
+  TestCollectKopStat<Runner>(Factory<SequentialFor>());
 }
 
 TEST(QTrajectoryAVXTest, CleanCircuit) {
-  TestCleanCircuit(qsim::Factory<SequentialFor>());
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory<SequentialFor>>;
+  TestCleanCircuit<Runner>(Factory<SequentialFor>());
 }
 
 TEST(QTrajectoryAVXTest, InitialState) {
-  TestInitialState(qsim::Factory<SequentialFor>());
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory<SequentialFor>>;
+  TestInitialState<Runner>(Factory<SequentialFor>());
 }
 
 TEST(QTrajectoryAVXTest, UncomputeFinalState) {
-  TestUncomputeFinalState(qsim::Factory<SequentialFor>());
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory<SequentialFor>>;
+  TestUncomputeFinalState<Runner>(Factory<SequentialFor>());
 }
 
 }  // namespace qsim
diff --git a/tests/qtrajectory_cuda_test.cu b/tests/qtrajectory_cuda_test.cu
index 459a25f63..e113e0521 100644
--- a/tests/qtrajectory_cuda_test.cu
+++ b/tests/qtrajectory_cuda_test.cu
@@ -16,6 +16,10 @@
 
 #include "gtest/gtest.h"
 
+#include "../lib/fuser_mqubit.h"
+#include "../lib/gates_cirq.h"
+#include "../lib/io.h"
+#include "../lib/run_qsim.h"
 #include "../lib/simulator_cuda.h"
 
 namespace qsim {
@@ -41,51 +45,65 @@ struct Factory {
 
 TEST(QTrajectoryCUDATest, BitFlip) {
   using Factory = qsim::Factory<float>;
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory>;
   Factory::StateSpace::Parameter param;
   Factory factory(param);
-  TestBitFlip(factory);
+  TestBitFlip<Runner>(factory);
 }
 
 TEST(QTrajectoryCUDATest, GenDump) {
   using Factory = qsim::Factory<float>;
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory>;
   Factory::StateSpace::Parameter param;
   Factory factory(param);
-  TestGenDump(factory);
+  TestGenDump<Runner>(factory);
 }
 
 TEST(QTrajectoryCUDATest, ReusingResults) {
   using Factory = qsim::Factory<float>;
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory>;
   Factory::StateSpace::Parameter param;
   Factory factory(param);
-  TestReusingResults(factory);
+  TestReusingResults<Runner>(factory);
 }
 
 TEST(QTrajectoryCUDATest, CollectKopStat) {
   using Factory = qsim::Factory<float>;
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory>;
   Factory::StateSpace::Parameter param;
   Factory factory(param);
-  TestCollectKopStat(factory);
+  TestCollectKopStat<Runner>(factory);
 }
 
 TEST(QTrajectoryCUDATest, CleanCircuit) {
   using Factory = qsim::Factory<float>;
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory>;
   Factory::StateSpace::Parameter param;
   Factory factory(param);
-  TestCleanCircuit(factory);
+  TestCleanCircuit<Runner>(factory);
 }
 
 TEST(QTrajectoryCUDATest, InitialState) {
   using Factory = qsim::Factory<float>;
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory>;
   Factory::StateSpace::Parameter param;
   Factory factory(param);
-  TestInitialState(factory);
+  TestInitialState<Runner>(factory);
 }
 
 TEST(QTrajectoryCUDATest, UncomputeFinalState) {
   using Factory = qsim::Factory<float>;
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory>;
   Factory::StateSpace::Parameter param;
   Factory factory(param);
-  TestUncomputeFinalState(factory);
+  TestUncomputeFinalState<Runner>(factory);
 }
 
 }  // namespace qsim
diff --git a/tests/qtrajectory_custatevec_test.cu b/tests/qtrajectory_custatevec_test.cu
index 1001805e5..ee283bd6c 100644
--- a/tests/qtrajectory_custatevec_test.cu
+++ b/tests/qtrajectory_custatevec_test.cu
@@ -19,6 +19,10 @@
 
 #include "gtest/gtest.h"
 
+#include "../lib/fuser_mqubit.h"
+#include "../lib/gates_cirq.h"
+#include "../lib/io.h"
+#include "../lib/run_qsim.h"
 #include "../lib/simulator_custatevec.h"
 
 namespace qsim {
@@ -52,31 +56,45 @@ struct Factory {
 };
 
 TEST(QTrajectoryCuStateVecTest, BitFlip) {
-  TestBitFlip(qsim::Factory<float>());
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory<float>>;
+  TestBitFlip<Runner>(Factory<float>());
 }
 
 TEST(QTrajectoryCuStateVecTest, GenDump) {
-  TestGenDump(qsim::Factory<float>());
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory<float>>;
+  TestGenDump<Runner>(Factory<float>());
 }
 
 TEST(QTrajectoryCuStateVecTest, ReusingResults) {
-  TestReusingResults(qsim::Factory<float>());
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory<float>>;
+  TestReusingResults<Runner>(Factory<float>());
 }
 
 TEST(QTrajectoryCuStateVecTest, CollectKopStat) {
-  TestCollectKopStat(qsim::Factory<float>());
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory<float>>;
+  TestCollectKopStat<Runner>(Factory<float>());
 }
 
 TEST(QTrajectoryCuStateVecTest, CleanCircuit) {
-  TestCleanCircuit(qsim::Factory<float>());
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory<float>>;
+  TestCleanCircuit<Runner>(Factory<float>());
 }
 
 TEST(QTrajectoryCuStateVecTest, InitialState) {
-  TestInitialState(qsim::Factory<float>());
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory<float>>;
+  TestInitialState<Runner>(Factory<float>());
 }
 
 TEST(QTrajectoryCuStateVecTest, UncomputeFinalState) {
-  TestUncomputeFinalState(qsim::Factory<float>());
+  using Fuser = MultiQubitGateFuser<IO, const Cirq::GateCirq<float>*>;
+  using Runner = QSimRunner<IO, Fuser, Factory<float>>;
+  TestUncomputeFinalState<Runner>(Factory<float>());
 }
 
 }  // namespace qsim
diff --git a/tests/qtrajectory_testfixture.h b/tests/qtrajectory_testfixture.h
index 786d997c8..c0aa5bef4 100644
--- a/tests/qtrajectory_testfixture.h
+++ b/tests/qtrajectory_testfixture.h
@@ -341,14 +341,13 @@ NoisyCircuit<Gate> GenerateNoisyCircuit(
   return ncircuit;
 }
 
-template <typename Factory, typename Gate>
+template <typename Runner, typename Factory, typename Gate>
 void RunBatch(const Factory& factory, const NoisyCircuit<Gate>& ncircuit,
               const std::vector<double>& expected_results) {
   using Simulator = typename Factory::Simulator;
   using StateSpace = typename Simulator::StateSpace;
   using State = typename StateSpace::State;
-  using QTSimulator = QuantumTrajectorySimulator<IO, Gate, MultiQubitGateFuser,
-                                                 Simulator>;
+  using QTSimulator = QuantumTrajectorySimulator<IO, Gate, Runner>;
 
   unsigned num_qubits = 2;
   unsigned num_reps = 25000;
@@ -376,15 +375,14 @@ void RunBatch(const Factory& factory, const NoisyCircuit<Gate>& ncircuit,
   }
 }
 
-template <typename Factory, typename Gate>
+template <typename Runner, typename Factory, typename Gate>
 void RunOnceRepeatedly(const Factory& factory,
                        const NoisyCircuit<Gate>& ncircuit,
                        const std::vector<double>& expected_results) {
   using Simulator = typename Factory::Simulator;
   using StateSpace = typename Factory::StateSpace;
   using State = typename StateSpace::State;
-  using QTSimulator = QuantumTrajectorySimulator<IO, Gate, MultiQubitGateFuser,
-                                                 Simulator>;
+  using QTSimulator = QuantumTrajectorySimulator<IO, Gate, Runner>;
 
   unsigned num_qubits = 2;
   unsigned num_reps = 25000;
@@ -417,7 +415,7 @@ void RunOnceRepeatedly(const Factory& factory,
   }
 }
 
-template <typename Factory, typename Gate>
+template <typename Runner, typename Factory, typename Gate>
 std::vector<std::complex<double>> ExpValsRunBatch(
     const Factory& factory, const NoisyCircuit<Gate>& ncircuit,
     bool reuse_results) {
@@ -425,8 +423,7 @@ std::vector<std::complex<double>> ExpValsRunBatch(
   using StateSpace = typename Factory::StateSpace;
   using State = typename StateSpace::State;
   using Fuser = MultiQubitGateFuser<IO, Gate>;
-  using QTSimulator = QuantumTrajectorySimulator<IO, Gate, MultiQubitGateFuser,
-                                                 Simulator>;
+  using QTSimulator = QuantumTrajectorySimulator<IO, Gate, Runner>;
 
   unsigned num_qubits = 2;
   unsigned num_reps = 25000;
@@ -507,7 +504,7 @@ std::vector<std::complex<double>> ExpValsRunBatch(
   return results;
 }
 
-template <typename Factory, typename Gate>
+template <typename Runner, typename Factory, typename Gate>
 std::vector<std::complex<double>> ExpValsRunOnceRepeatedly(
     const Factory& factory, const NoisyCircuit<Gate>& ncircuit,
     bool reuse_results) {
@@ -515,8 +512,7 @@ std::vector<std::complex<double>> ExpValsRunOnceRepeatedly(
   using StateSpace = typename Factory::StateSpace;
   using State = typename StateSpace::State;
   using Fuser = MultiQubitGateFuser<IO, Gate>;
-  using QTSimulator = QuantumTrajectorySimulator<IO, Gate, MultiQubitGateFuser,
-                                                 Simulator>;
+  using QTSimulator = QuantumTrajectorySimulator<IO, Gate, Runner>;
 
   unsigned num_qubits = 2;
   unsigned num_reps = 25000;
@@ -592,7 +588,7 @@ std::vector<std::complex<double>> ExpValsRunOnceRepeatedly(
   return results;
 }
 
-template <typename Factory>
+template <typename Runner, typename Factory>
 void TestBitFlip(const Factory& factory) {
 /* The expected results are obtained with the following Cirq code.
 
@@ -633,10 +629,10 @@ for key, val in sorted(res.histogram(key='m').items()):
 
   auto ncircuit = GenerateNoisyCircuit<Gate>(0.01, AddBitFlipNoise1<Gate>,
                                              AddBitFlipNoise2<Gate>);
-  RunBatch(factory, ncircuit, expected_results);
+  RunBatch<Runner>(factory, ncircuit, expected_results);
 }
 
-template <typename Factory>
+template <typename Runner, typename Factory>
 void TestGenDump(const Factory& factory) {
 /* The expected results are obtained with the following Cirq code.
 
@@ -692,27 +688,27 @@ for key, val in sorted(res.histogram(key='m').items()):
   {
     auto ncircuit = GenerateNoisyCircuit<Gate>(0.1, AddGenAmplDumpNoise1<Gate>,
                                                AddGenAmplDumpNoise2<Gate>);
-    RunOnceRepeatedly(factory, ncircuit, expected_results);
+    RunOnceRepeatedly<Runner>(factory, ncircuit, expected_results);
   }
 
   {
     auto ncircuit = GenerateNoisyCircuit<Gate>(0.1, AddGenAmplDumpNoise1<Gate>,
                                                AddGenAmplDumpNoise2Alt<Gate>);
-    RunOnceRepeatedly(factory, ncircuit, expected_results);
+    RunOnceRepeatedly<Runner>(factory, ncircuit, expected_results);
   }
 }
 
-template <typename Factory>
+template <typename Runner, typename Factory>
 void TestReusingResults(const Factory& factory) {
   using Gate = Cirq::GateCirq<typename Factory::fp_type>;
 
   auto ncircuit = GenerateNoisyCircuit<Gate>(0.02, AddAmplDumpNoise1<Gate>,
                                              AddAmplDumpNoise2<Gate>, false);
 
-  auto results1 = ExpValsRunOnceRepeatedly(factory, ncircuit, false);
-  auto results2 = ExpValsRunOnceRepeatedly(factory, ncircuit, true);
-  auto results3 = ExpValsRunBatch(factory, ncircuit, false);
-  auto results4 = ExpValsRunBatch(factory, ncircuit, true);
+  auto results1 = ExpValsRunOnceRepeatedly<Runner>(factory, ncircuit, false);
+  auto results2 = ExpValsRunOnceRepeatedly<Runner>(factory, ncircuit, true);
+  auto results3 = ExpValsRunBatch<Runner>(factory, ncircuit, false);
+  auto results4 = ExpValsRunBatch<Runner>(factory, ncircuit, true);
 
   for (std::size_t k = 0; k < results1.size(); ++k) {
     EXPECT_NEAR(std::real(results1[k]), std::real(results2[k]), 1e-8);
@@ -724,16 +720,14 @@ void TestReusingResults(const Factory& factory) {
   }
 }
 
-template <typename Factory>
+template <typename Runner, typename Factory>
 void TestCollectKopStat(const Factory& factory) {
   using Simulator = typename Factory::Simulator;
   using StateSpace = typename Factory::StateSpace;
   using State = typename StateSpace::State;
   using fp_type = typename StateSpace::fp_type;
   using GateCirq = Cirq::GateCirq<fp_type>;
-  using QTSimulator = QuantumTrajectorySimulator<IO, GateCirq,
-                                                 MultiQubitGateFuser,
-                                                 Simulator>;
+  using QTSimulator = QuantumTrajectorySimulator<IO, GateCirq, Runner>;
 
   unsigned num_qubits = 4;
   unsigned num_reps = 20000;
@@ -799,16 +793,14 @@ void TestCollectKopStat(const Factory& factory) {
   }
 }
 
-template <typename Factory>
+template <typename Runner, typename Factory>
 void TestCleanCircuit(const Factory& factory) {
   using Simulator = typename Factory::Simulator;
   using StateSpace = typename Factory::StateSpace;
   using State = typename StateSpace::State;
   using fp_type = typename StateSpace::fp_type;
   using GateCirq = Cirq::GateCirq<fp_type>;
-  using QTSimulator = QuantumTrajectorySimulator<IO, GateCirq,
-                                                 MultiQubitGateFuser,
-                                                 Simulator>;
+  using QTSimulator = QuantumTrajectorySimulator<IO, GateCirq, Runner>;
 
   unsigned num_qubits = 4;
   auto size = uint64_t{1} << num_qubits;
@@ -892,16 +884,14 @@ void TestCleanCircuit(const Factory& factory) {
 }
 
 // Test that QTSimulator::Run does not overwrite initial states.
-template <typename Factory>
+template <typename Runner, typename Factory>
 void TestInitialState(const Factory& factory) {
   using Simulator = typename Factory::Simulator;
   using StateSpace = typename Factory::StateSpace;
   using State = typename StateSpace::State;
   using fp_type = typename StateSpace::fp_type;
   using GateCirq = Cirq::GateCirq<fp_type>;
-  using QTSimulator = QuantumTrajectorySimulator<IO, GateCirq,
-                                                 MultiQubitGateFuser,
-                                                 Simulator>;
+  using QTSimulator = QuantumTrajectorySimulator<IO, GateCirq, Runner>;
 
   unsigned num_qubits = 3;
 
@@ -942,16 +932,14 @@ void TestInitialState(const Factory& factory) {
   }
 }
 
-template <typename Factory>
+template <typename Runner, typename Factory>
 void TestUncomputeFinalState(const Factory& factory) {
   using Simulator = typename Factory::Simulator;
   using StateSpace = typename Factory::StateSpace;
   using State = typename StateSpace::State;
   using fp_type = typename StateSpace::fp_type;
   using GateCirq = Cirq::GateCirq<fp_type>;
-  using QTSimulator = QuantumTrajectorySimulator<IO, GateCirq,
-                                                 MultiQubitGateFuser,
-                                                 Simulator>;
+  using QTSimulator = QuantumTrajectorySimulator<IO, GateCirq, Runner>;
 
   unsigned num_qubits = 4;