diff --git a/pytket/binders/passes.cpp b/pytket/binders/passes.cpp
index 9df5f70d0e..ef398f8a0d 100644
--- a/pytket/binders/passes.cpp
+++ b/pytket/binders/passes.cpp
@@ -970,14 +970,19 @@ PYBIND11_MODULE(passes, m) {
       "\n:param allow_zzphase: If set to True, allows the algorithm to "
       "implement 2-qubit rotations using ZZPhase gates when deemed "
       "optimal. Defaults to False."
-      "\n:param timeout: Sets maximum out of time spent finding solution."
+      "\n:param thread_timeout: Sets maximum out of time spent finding a "
+      "single solution in one thread."
       "\n:param only_reduce: Only returns modified circuit if it has "
       "fewer two-qubit gates."
+      "\n:param trials: Sets maximum number of found solutions. The "
+      "smallest circuit is returned, prioritising the number of 2qb-gates, "
+      "then the number of gates, then the depth."
       "\n:return: a pass to perform the simplification",
       py::arg("discount_rate") = 0.7, py::arg("depth_weight") = 0.3,
       py::arg("max_lookahead") = 500, py::arg("max_tqe_candidates") = 500,
       py::arg("seed") = 0, py::arg("allow_zzphase") = false,
-      py::arg("timeout") = 100, py::arg("only_reduce") = false);
+      py::arg("thread_timeout") = 100, py::arg("only_reduce") = false,
+      py::arg("trials") = 1);
   m.def(
       "PauliSquash", &PauliSquash,
       "Applies :py:meth:`PauliSimp` followed by "
diff --git a/pytket/binders/transform.cpp b/pytket/binders/transform.cpp
index 2aaf153972..286fc6062f 100644
--- a/pytket/binders/transform.cpp
+++ b/pytket/binders/transform.cpp
@@ -451,12 +451,16 @@ PYBIND11_MODULE(transform, m) {
           "\n:param allow_zzphase: If set to True, allows the algorithm to "
           "implement 2-qubit rotations using ZZPhase gates when deemed "
           "optimal. Defaults to False."
-          "\n:param timeout: Sets maximum out of time spent finding solution."
+          "\n:param thread_timeout: Sets maximum out of time spent finding a "
+          "single solution in one thread."
+          "\n:param trials: Sets maximum number of found solutions. The "
+          "smallest circuit is returned, prioritising the number of 2qb-gates, "
+          "then the number of gates, then the depth."
           "\n:return: a pass to perform the simplification",
           py::arg("discount_rate") = 0.7, py::arg("depth_weight") = 0.3,
           py::arg("max_tqe_candidates") = 500, py::arg("max_lookahead") = 500,
           py::arg("seed") = 0, py::arg("allow_zzphase") = false,
-          py::arg("timeout") = 100)
+          py::arg("thread_timeout") = 100, py::arg("trials") = 1)
       .def_static(
           "ZZPhaseToRz", &Transforms::ZZPhase_to_Rz,
           "Fixes all ZZPhase gate angles to [-1, 1) half turns.")
diff --git a/pytket/conanfile.py b/pytket/conanfile.py
index 9b76b1cb1a..123ac8339a 100644
--- a/pytket/conanfile.py
+++ b/pytket/conanfile.py
@@ -38,7 +38,7 @@ def requirements(self):
         self.requires("pybind11_json/0.2.14")
         self.requires("symengine/0.13.0")
         self.requires("tkassert/0.3.4@tket/stable")
-        self.requires("tket/1.3.47@tket/stable")
+        self.requires("tket/1.3.48@tket/stable")
         self.requires("tklog/0.3.3@tket/stable")
         self.requires("tkrng/0.3.3@tket/stable")
         self.requires("tktokenswap/0.3.9@tket/stable")
diff --git a/pytket/docs/changelog.rst b/pytket/docs/changelog.rst
index 368432fbb6..388211ff4a 100644
--- a/pytket/docs/changelog.rst
+++ b/pytket/docs/changelog.rst
@@ -12,8 +12,8 @@ Features:
   conditions.
 * Add `custom_deserialisation` argument to `BasePass` and `SequencePass` 
   `from_dict` method to support construction of `CustomPass` from json.
-* Add `timeout` argument to `GreedyPauliSimp`.
-* Add `only_reduce` argument to `GreedyPauliSimp`.
+* Add `thread_timeout`, `only_reduce`, and `trials` arguments 
+  to `GreedyPauliSimp`.
 * Add option to not relabel `ClassicalExpBox` when calling `rename_units`
   and `flatten_registers`
 * Implement `dagger()` and `transpose()` for `CustomGate`.
diff --git a/pytket/pytket/_tket/passes.pyi b/pytket/pytket/_tket/passes.pyi
index 6bfc1477e5..9a59fae388 100644
--- a/pytket/pytket/_tket/passes.pyi
+++ b/pytket/pytket/_tket/passes.pyi
@@ -447,7 +447,7 @@ def GlobalisePhasedX(squash: bool = True) -> BasePass:
     
     It is not recommended to use this pass with symbolic expressions, as in certain cases a blow-up in symbolic expression sizes may occur.
     """
-def GreedyPauliSimp(discount_rate: float = 0.7, depth_weight: float = 0.3, max_lookahead: int = 500, max_tqe_candidates: int = 500, seed: int = 0, allow_zzphase: bool = False, timeout: int = 100, only_reduce: bool = False) -> BasePass:
+def GreedyPauliSimp(discount_rate: float = 0.7, depth_weight: float = 0.3, max_lookahead: int = 500, max_tqe_candidates: int = 500, seed: int = 0, allow_zzphase: bool = False, thread_timeout: int = 100, only_reduce: bool = False, trials: int = 1) -> BasePass:
     """
     Construct a pass that converts a circuit into a graph of Pauli gadgets to account for commutation and phase folding, and resynthesises them using a greedy algorithm adapted from arxiv.org/abs/2103.08602. The method for synthesising the final Clifford operator is adapted from arxiv.org/abs/2305.10966.
     
@@ -457,8 +457,9 @@ def GreedyPauliSimp(discount_rate: float = 0.7, depth_weight: float = 0.3, max_l
     :param max_lookahead:  Maximum lookahead when evaluating each Clifford gate candidate. Default to 500.
     :param seed:  Unsigned integer seed used for sampling candidates and tie breaking. Default to 0.
     :param allow_zzphase: If set to True, allows the algorithm to implement 2-qubit rotations using ZZPhase gates when deemed optimal. Defaults to False.
-    :param timeout: Sets maximum out of time spent finding solution.
+    :param thread_timeout: Sets maximum out of time spent finding a single solution in one thread.
     :param only_reduce: Only returns modified circuit if it has fewer two-qubit gates.
+    :param trials: Sets maximum number of found solutions. The smallest circuit is returned, prioritising the number of 2qb-gates, then the number of gates, then the depth.
     :return: a pass to perform the simplification
     """
 def GuidedPauliSimp(strat: pytket._tket.transform.PauliSynthStrat = pytket._tket.transform.PauliSynthStrat.Sets, cx_config: pytket._tket.circuit.CXConfigType = pytket._tket.circuit.CXConfigType.Snake) -> BasePass:
diff --git a/pytket/pytket/_tket/transform.pyi b/pytket/pytket/_tket/transform.pyi
index 3648de26da..4b137ff45f 100644
--- a/pytket/pytket/_tket/transform.pyi
+++ b/pytket/pytket/_tket/transform.pyi
@@ -167,7 +167,7 @@ class Transform:
         It is not recommended to use this transformation with symbolic expressions, as in certain cases a blow-up in symbolic expression sizes may occur.
         """
     @staticmethod
-    def GreedyPauliSimp(discount_rate: float = 0.7, depth_weight: float = 0.3, max_tqe_candidates: int = 500, max_lookahead: int = 500, seed: int = 0, allow_zzphase: bool = False, timeout: int = 100) -> Transform:
+    def GreedyPauliSimp(discount_rate: float = 0.7, depth_weight: float = 0.3, max_tqe_candidates: int = 500, max_lookahead: int = 500, seed: int = 0, allow_zzphase: bool = False, thread_timeout: int = 100, trials: int = 1) -> Transform:
         """
         Convert a circuit into a graph of Pauli gadgets to account for commutation and phase folding, and resynthesises them using a greedy algorithm adapted from arxiv.org/abs/2103.08602. The method for synthesising the final Clifford operator is adapted from arxiv.org/abs/2305.10966.
         
@@ -177,7 +177,8 @@ class Transform:
         :param max_lookahead:  Maximum lookahead when evaluating each Clifford gate candidate. Default to 500.
         :param seed:  Unsigned integer seed used for sampling candidates and tie breaking. Default to 0.
         :param allow_zzphase: If set to True, allows the algorithm to implement 2-qubit rotations using ZZPhase gates when deemed optimal. Defaults to False.
-        :param timeout: Sets maximum out of time spent finding solution.
+        :param thread_timeout: Sets maximum out of time spent finding a single solution in one thread.
+        :param trials: Sets maximum number of found solutions. The smallest circuit is returned, prioritising the number of 2qb-gates, then the number of gates, then the depth.
         :return: a pass to perform the simplification
         """
     @staticmethod
diff --git a/pytket/tests/passes_serialisation_test.py b/pytket/tests/passes_serialisation_test.py
index 9c145df6ac..59db4ade06 100644
--- a/pytket/tests/passes_serialisation_test.py
+++ b/pytket/tests/passes_serialisation_test.py
@@ -296,8 +296,9 @@ def nonparam_predicate_dict(name: str) -> Dict[str, Any]:
             "max_tqe_candidates": 100,
             "seed": 2,
             "allow_zzphase": True,
-            "timeout": 5000,
+            "thread_timeout": 5000,
             "only_reduce": False,
+            "trials": 1,
         }
     ),
     # lists must be sorted by OpType value
diff --git a/pytket/tests/predicates_test.py b/pytket/tests/predicates_test.py
index 5de371e01a..1bda5404cc 100644
--- a/pytket/tests/predicates_test.py
+++ b/pytket/tests/predicates_test.py
@@ -1035,9 +1035,8 @@ def test_greedy_pauli_synth() -> None:
         rega[0], regb[0]
     ).SWAP(regb[1], rega[0])
     d = circ.copy()
-    pss = GreedyPauliSimp(0.5, 0.5)
-    assert not GreedyPauliSimp(0.5, 0.5, timeout=0, only_reduce=False).apply(d)
-    assert pss.apply(d)
+    assert GreedyPauliSimp(0.5, 0.5, thread_timeout=10, trials=5).apply(d)
+
     assert np.allclose(circ.get_unitary(), d.get_unitary())
     assert d.name == "test"
     # test gateset
@@ -1059,7 +1058,6 @@ def test_greedy_pauli_synth() -> None:
     circ.measure_all()
     circ.Reset(0)
     circ.add_pauliexpbox(pg1, [2, 3])
-    assert not GreedyPauliSimp(0.5, 0.5, 100, 100, 0, True, 0).apply(circ)
     assert GreedyPauliSimp(0.5, 0.5, 100, 100, 0, True, 100).apply(circ)
     # PauliExpBoxes implemented using ZZPhase
     d = Circuit(4, 4, name="test")
@@ -1084,6 +1082,16 @@ def test_greedy_pauli_synth() -> None:
         GreedyPauliSimp().apply(circ)
     err_msg = "Predicate requirements are not satisfied"
     assert err_msg in str(e.value)
+    # large circuit that doesn't complete within thread_timeout argument
+    c = Circuit(13)
+    for _ in range(20):
+        for i in range(13):
+            for j in range(i + 1, 13):
+                c.CX(i, j)
+                c.Rz(0.23, j)
+            c.H(i)
+    assert not GreedyPauliSimp(thread_timeout=1).apply(c)
+    assert GreedyPauliSimp().apply(c)
 
 
 def test_auto_rebase_deprecation(recwarn: Any) -> None:
diff --git a/schemas/compiler_pass_v1.json b/schemas/compiler_pass_v1.json
index 8a5e1a3b8c..3a3dd7d85f 100644
--- a/schemas/compiler_pass_v1.json
+++ b/schemas/compiler_pass_v1.json
@@ -374,13 +374,17 @@
           "type": "boolean",
           "definition": "parameter controlling the use of ZZPhase gates in \"GreedyPauliSimp\""
         },
-        "timeout": {
+        "thread_timeout": {
           "type": "number",
-          "definition": "parameter controlling the maximum runtime of \"GreedyPauliSimp\""
+          "definition": "parameter controlling the maximum runtime of a single thread in \"GreedyPauliSimp\""
         },
         "only_reduce": {
           "type": "boolean",
           "definition": "parameter controlling whether \"GreedyPauliSimp\" can return circuits with more two qubit gates"
+        },
+        "trials": {
+          "type": "number",
+          "definition": "parameter controlling the number of random solutions found when calling \"GreedyPauliSimp\""
         }
       },
       "required": [
@@ -913,10 +917,11 @@
               "max_tqe_candidates",
               "seed",
               "allow_zzphase",
-              "timeout",
-              "only_reduce"
+              "thread_timeout",
+              "only_reduce",
+              "trials"
             ],
-            "maxProperties": 9
+            "maxProperties": 10
           }
         },
         {
diff --git a/tket/conanfile.py b/tket/conanfile.py
index c479bca357..fa46ab76f8 100644
--- a/tket/conanfile.py
+++ b/tket/conanfile.py
@@ -23,7 +23,7 @@
 
 class TketConan(ConanFile):
     name = "tket"
-    version = "1.3.47"
+    version = "1.3.48"
     package_type = "library"
     license = "Apache 2"
     homepage = "https://github.com/CQCL/tket"
diff --git a/tket/include/tket/Predicates/PassGenerators.hpp b/tket/include/tket/Predicates/PassGenerators.hpp
index 4d3e8b1a31..65861a4659 100644
--- a/tket/include/tket/Predicates/PassGenerators.hpp
+++ b/tket/include/tket/Predicates/PassGenerators.hpp
@@ -355,15 +355,17 @@ PassPtr gen_special_UCC_synthesis(
  * @param max_tqe_candidates
  * @param seed
  * @param allow_zzphase
- * @param timeout
+ * @param thread_timeout
  * @param only_reduce
+ * @param trials
  * @return PassPtr
  */
 PassPtr gen_greedy_pauli_simp(
     double discount_rate = 0.7, double depth_weight = 0.3,
     unsigned max_lookahead = 500, unsigned max_tqe_candidates = 500,
-    unsigned seed = 0, bool allow_zzphase = false, unsigned timeout = 100,
-    bool only_reduce = false);
+    unsigned seed = 0, bool allow_zzphase = false,
+    unsigned thread_timeout = 100, bool only_reduce = false,
+    unsigned trials = 1);
 
 /**
  * Generate a pass to simplify the circuit where it acts on known basis states.
diff --git a/tket/include/tket/Transformations/GreedyPauliOptimisation.hpp b/tket/include/tket/Transformations/GreedyPauliOptimisation.hpp
index 3cc6348582..8f53b08bde 100644
--- a/tket/include/tket/Transformations/GreedyPauliOptimisation.hpp
+++ b/tket/include/tket/Transformations/GreedyPauliOptimisation.hpp
@@ -14,6 +14,8 @@
 
 #pragma once
 
+#include <atomic>
+
 #include "Transform.hpp"
 #include "tket/Circuit/Circuit.hpp"
 #include "tket/Clifford/UnitaryTableau.hpp"
@@ -600,6 +602,29 @@ class GPGraph {
 std::tuple<std::vector<PauliNode_ptr>, std::vector<PauliNode_ptr>>
 gpg_from_unordered_set(const std::vector<SymPauliTensor>& unordered_set);
 
+/**
+ * @brief Converts the given circuit into a GPGraph and conjugates each node
+ * by greedily applying 2-qubit Clifford gates until the node can be realised
+ * as a single-qubit gate, a measurement, or a reset. The final Clifford
+ * operator is synthesized in a similar fashion. Allows early termination
+ * from a thread via a stop_flag.
+ *
+ * @param circ
+ * @param stop_flag
+ * @param discount_rate
+ * @param depth_weight
+ * @param max_lookahead
+ * @param max_tqe_candidates
+ * @param seed
+ * @param allow_zzphase
+ * @return Circuit
+ */
+Circuit greedy_pauli_graph_synthesis_flag(
+    Circuit circ, std::atomic<bool>& stop_flag, double discount_rate = 0.7,
+    double depth_weight = 0.3, unsigned max_lookahead = 500,
+    unsigned max_tqe_candidates = 500, unsigned seed = 0,
+    bool allow_zzphase = false);
+
 /**
  * @brief Converts the given circuit into a GPGraph and conjugates each node
  * by greedily applying 2-qubit Clifford gates until the node can be realised
@@ -643,7 +668,8 @@ Circuit greedy_pauli_set_synthesis(
 Transform greedy_pauli_optimisation(
     double discount_rate = 0.7, double depth_weight = 0.3,
     unsigned max_lookahead = 500, unsigned max_tqe_candidates = 500,
-    unsigned seed = 0, bool allow_zzphase = false, unsigned timeout = 100);
+    unsigned seed = 0, bool allow_zzphase = false,
+    unsigned thread_timeout = 100, unsigned trials = 1);
 
 }  // namespace Transforms
 
diff --git a/tket/src/Predicates/CompilerPass.cpp b/tket/src/Predicates/CompilerPass.cpp
index 516e13cb53..e09262fcc7 100644
--- a/tket/src/Predicates/CompilerPass.cpp
+++ b/tket/src/Predicates/CompilerPass.cpp
@@ -522,11 +522,12 @@ PassPtr deserialise(
       unsigned max_lookahead = content.at("max_lookahead").get<unsigned>();
       unsigned seed = content.at("seed").get<unsigned>();
       bool allow_zzphase = content.at("allow_zzphase").get<bool>();
-      unsigned timeout = content.at("timeout").get<unsigned>();
+      unsigned timeout = content.at("thread_timeout").get<unsigned>();
       bool only_reduce = content.at("only_reduce").get<bool>();
+      unsigned trials = content.at("trials").get<unsigned>();
       pp = gen_greedy_pauli_simp(
           discount_rate, depth_weight, max_lookahead, max_tqe_candidates, seed,
-          allow_zzphase, timeout, only_reduce);
+          allow_zzphase, timeout, only_reduce, trials);
 
     } else if (passname == "PauliSimp") {
       // SEQUENCE PASS - DESERIALIZABLE ONLY
diff --git a/tket/src/Predicates/PassGenerators.cpp b/tket/src/Predicates/PassGenerators.cpp
index 4d99463d0b..88766bbb2e 100644
--- a/tket/src/Predicates/PassGenerators.cpp
+++ b/tket/src/Predicates/PassGenerators.cpp
@@ -1017,24 +1017,33 @@ PassPtr gen_synthesise_pauli_graph(
 PassPtr gen_greedy_pauli_simp(
     double discount_rate, double depth_weight, unsigned max_lookahead,
     unsigned max_tqe_candidates, unsigned seed, bool allow_zzphase,
-    unsigned timeout, bool only_reduce) {
-  Transform t =
-      Transform([discount_rate, depth_weight, max_lookahead, max_tqe_candidates,
-                 seed, allow_zzphase, timeout, only_reduce](Circuit& circ) {
-        Transform gpo = Transforms::greedy_pauli_optimisation(
-            discount_rate, depth_weight, max_lookahead, max_tqe_candidates,
-            seed, allow_zzphase, timeout);
-        if (only_reduce) {
-          Circuit gpo_circ = circ;
-          if (gpo.apply(gpo_circ) &&
-              gpo_circ.count_n_qubit_gates(2) < circ.count_n_qubit_gates(2)) {
-            circ = gpo_circ;
-            return true;
-          }
-          return false;
-        }
-        return gpo.apply(circ);
-      });
+    unsigned thread_timeout, bool only_reduce, unsigned trials) {
+  Transform t = Transform([discount_rate, depth_weight, max_lookahead,
+                           max_tqe_candidates, seed, allow_zzphase,
+                           thread_timeout, only_reduce, trials](Circuit& circ) {
+    Transform gpo = Transforms::greedy_pauli_optimisation(
+        discount_rate, depth_weight, max_lookahead, max_tqe_candidates, seed,
+        allow_zzphase, thread_timeout, trials);
+    if (only_reduce) {
+      Circuit gpo_circ = circ;
+      // comparison will be inaccurate if circuit has PauliExpBox
+      gpo_circ.decompose_boxes_recursively();
+      unsigned original_n_2qb_gates = gpo_circ.count_n_qubit_gates(2);
+      unsigned original_n_gates = gpo_circ.n_gates();
+      unsigned original_depth = gpo_circ.depth();
+      if (gpo.apply(gpo_circ) &&
+          std::make_tuple(
+              gpo_circ.count_n_qubit_gates(2), gpo_circ.n_gates(),
+              gpo_circ.depth()) <
+              std::make_tuple(
+                  original_n_2qb_gates, original_n_gates, original_depth)) {
+        circ = gpo_circ;
+        return true;
+      }
+      return false;
+    }
+    return gpo.apply(circ);
+  });
   OpTypeSet ins = {
       OpType::Z,
       OpType::X,
@@ -1083,8 +1092,9 @@ PassPtr gen_greedy_pauli_simp(
   j["max_tqe_candidates"] = max_tqe_candidates;
   j["seed"] = seed;
   j["allow_zzphase"] = allow_zzphase;
-  j["timeout"] = timeout;
+  j["thread_timeout"] = thread_timeout;
   j["only_reduce"] = only_reduce;
+  j["trials"] = trials;
 
   return std::make_shared<StandardPass>(precons, t, postcon, j);
 }
diff --git a/tket/src/Transformations/GreedyPauliOptimisation.cpp b/tket/src/Transformations/GreedyPauliOptimisation.cpp
index bf88e4a265..3c80ee1eed 100644
--- a/tket/src/Transformations/GreedyPauliOptimisation.cpp
+++ b/tket/src/Transformations/GreedyPauliOptimisation.cpp
@@ -17,6 +17,7 @@
 #include <algorithm>
 #include <chrono>
 #include <future>
+#include <queue>
 #include <random>
 
 #include "tket/Circuit/PauliExpBoxes.hpp"
@@ -327,7 +328,8 @@ struct DepthTracker {
 static void tableau_row_nodes_synthesis(
     std::vector<PauliNode_ptr>& rows, Circuit& circ,
     DepthTracker& depth_tracker, double depth_weight, unsigned max_lookahead,
-    unsigned max_tqe_candidates, unsigned seed) {
+    unsigned max_tqe_candidates, unsigned seed,
+    std::shared_ptr<std::atomic<bool>> stop_flag) {
   // only consider nodes with a non-zero cost
   std::vector<unsigned> remaining_indices;
   for (unsigned i = 0; i < rows.size(); i++) {
@@ -336,6 +338,10 @@ static void tableau_row_nodes_synthesis(
     }
   }
   while (remaining_indices.size() != 0) {
+    // check early termination
+    if (stop_flag.get()->load()) {
+      return;
+    };
     // get nodes with min cost
     std::vector<unsigned> min_nodes_indices = {remaining_indices[0]};
     unsigned min_cost = rows[remaining_indices[0]]->tqe_cost();
@@ -607,11 +613,18 @@ static void pauli_exps_synthesis(
     std::vector<PauliNode_ptr>& rows, Circuit& circ,
     DepthTracker& depth_tracker, double discount_rate, double depth_weight,
     unsigned max_lookahead, unsigned max_tqe_candidates, unsigned seed,
-    bool allow_zzphase) {
+    bool allow_zzphase, std::shared_ptr<std::atomic<bool>> stop_flag) {
   while (true) {
+    // check timeout
+    if (stop_flag.get()->load()) {
+      return;
+    };
+
     consume_nodes(
         rotation_sets, circ, depth_tracker, discount_rate, depth_weight);
+
     if (rotation_sets.empty()) break;
+
     std::vector<PauliNode_ptr>& first_set = rotation_sets[0];
     // get nodes with min cost
     std::vector<unsigned> min_nodes_indices = {0};
@@ -634,6 +647,7 @@ static void pauli_exps_synthesis(
     // sample
     std::vector<TQE> sampled_tqes =
         sample_tqes(tqe_candidates, max_tqe_candidates, seed);
+
     // for each tqe we compute costs which might subject to normalisation
     std::map<TQE, std::vector<double>> tqe_candidates_cost;
     for (const TQE& tqe : sampled_tqes) {
@@ -721,21 +735,23 @@ Circuit greedy_pauli_set_synthesis(
   std::vector<std::vector<PauliNode_ptr>> rotation_sets{rotation_set};
   DepthTracker depth_tracker(n_qubits);
   // synthesise Pauli exps
+  std::shared_ptr<std::atomic<bool>> dummy_stop_flag =
+      std::make_shared<std::atomic<bool>>(false);
   pauli_exps_synthesis(
       rotation_sets, rows, c, depth_tracker, 0, depth_weight, max_lookahead,
-      max_tqe_candidates, seed, allow_zzphase);
+      max_tqe_candidates, seed, allow_zzphase, dummy_stop_flag);
   // synthesise the tableau
   tableau_row_nodes_synthesis(
       rows, c, depth_tracker, depth_weight, max_lookahead, max_tqe_candidates,
-      seed);
+      seed, dummy_stop_flag);
   c.replace_SWAPs();
   return c;
 }
 
-Circuit greedy_pauli_graph_synthesis(
-    const Circuit& circ, double discount_rate, double depth_weight,
-    unsigned max_lookahead, unsigned max_tqe_candidates, unsigned seed,
-    bool allow_zzphase) {
+Circuit greedy_pauli_graph_synthesis_flag(
+    Circuit circ, std::shared_ptr<std::atomic<bool>> stop_flag,
+    double discount_rate, double depth_weight, unsigned max_lookahead,
+    unsigned max_tqe_candidates, unsigned seed, bool allow_zzphase) {
   if (max_lookahead == 0) {
     throw GreedyPauliSimpError("max_lookahead must be greater than 0.");
   }
@@ -757,16 +773,36 @@ Circuit greedy_pauli_graph_synthesis(
     rev_unit_map.insert({pair.second, pair.first});
   }
   GPGraph gpg(circ_flat);
+
+  // We regularly check whether the timeout has ocurred
+  if (stop_flag.get()->load()) {
+    return Circuit();
+  }
+
   auto [rotation_sets, rows, measures] = gpg.get_sequence();
+
+  if (stop_flag.get()->load()) {
+    return Circuit();
+  }
+
   DepthTracker depth_tracker(n_qubits);
   // synthesise Pauli exps
   pauli_exps_synthesis(
       rotation_sets, rows, new_circ, depth_tracker, discount_rate, depth_weight,
-      max_lookahead, max_tqe_candidates, seed, allow_zzphase);
+      max_lookahead, max_tqe_candidates, seed, allow_zzphase, stop_flag);
+
+  if (stop_flag.get()->load()) {
+    return Circuit();
+  }
   // synthesise the tableau
   tableau_row_nodes_synthesis(
       rows, new_circ, depth_tracker, depth_weight, max_lookahead,
-      max_tqe_candidates, seed);
+      max_tqe_candidates, seed, stop_flag);
+
+  if (stop_flag.get()->load()) {
+    return Circuit();
+  }
+
   for (auto it = measures.begin(); it != measures.end(); ++it) {
     new_circ.add_measure(it->left, it->right);
   }
@@ -775,26 +811,69 @@ Circuit greedy_pauli_graph_synthesis(
   return new_circ;
 }
 
+Circuit greedy_pauli_graph_synthesis(
+    const Circuit& circ, double discount_rate, double depth_weight,
+    unsigned max_lookahead, unsigned max_tqe_candidates, unsigned seed,
+    bool allow_zzphase) {
+  std::shared_ptr<std::atomic<bool>> dummy_stop_flag =
+      std::make_shared<std::atomic<bool>>(false);
+  return greedy_pauli_graph_synthesis_flag(
+      circ, dummy_stop_flag, discount_rate, depth_weight, max_lookahead,
+      max_tqe_candidates, seed, allow_zzphase);
+}
+
 }  // namespace GreedyPauliSimp
 
 Transform greedy_pauli_optimisation(
     double discount_rate, double depth_weight, unsigned max_lookahead,
     unsigned max_tqe_candidates, unsigned seed, bool allow_zzphase,
-    unsigned timeout) {
+    unsigned thread_timeout, unsigned trials) {
   return Transform([discount_rate, depth_weight, max_lookahead,
-                    max_tqe_candidates, seed, allow_zzphase,
-                    timeout](Circuit& circ) {
-    std::future<Circuit> future = std::async(
-        std::launch::async, GreedyPauliSimp::greedy_pauli_graph_synthesis, circ,
-        discount_rate, depth_weight, max_lookahead, max_tqe_candidates, seed,
-        allow_zzphase);
-    if (future.wait_for(std::chrono::seconds(timeout)) ==
-        std::future_status::ready) {
-      circ = future.get();
-      circ.decompose_boxes_recursively();
-      return true;
+                    max_tqe_candidates, seed, allow_zzphase, thread_timeout,
+                    trials](Circuit& circ) {
+    std::mt19937 seed_gen(seed);
+    std::vector<Circuit> circuits;
+    unsigned threads_started = 0;
+
+    while (threads_started < trials) {
+      std::shared_ptr<std::atomic<bool>> stop_flag =
+          std::make_shared<std::atomic<bool>>(false);
+      std::future<Circuit> future = std::async(
+          std::launch::async,
+          [&, stop_flag]() {  // Capture `stop_flag` explicitly in the lambda
+            return GreedyPauliSimp::greedy_pauli_graph_synthesis_flag(
+                circ, stop_flag, discount_rate, depth_weight, max_lookahead,
+                max_tqe_candidates, seed_gen(), allow_zzphase);
+          });
+      threads_started++;
+
+      if (future.wait_for(std::chrono::seconds(thread_timeout)) ==
+          std::future_status::ready) {
+        Circuit c = future.get();
+        c.decompose_boxes_recursively();
+        circuits.push_back(c);
+      } else {
+        // If the thread isn't complete within time, prompt cancelling the
+        // optimisation and break from while loop
+        *stop_flag = true;
+        break;
+      }
     }
-    return false;
+
+    // Return the smallest circuit if any were found within the single
+    // thread_timeout
+    // If none are found then return false
+    if (circuits.empty()) return false;
+    auto min = std::min_element(
+        circuits.begin(), circuits.end(),
+        [](const Circuit& a, const Circuit& b) {
+          return std::make_tuple(
+                     a.count_n_qubit_gates(2), a.n_gates(), a.depth()) <
+                 std::make_tuple(
+                     b.count_n_qubit_gates(2), b.n_gates(), b.depth());
+        });
+    circ = *min;
+    return true;
   });
 }
 
diff --git a/tket/test/src/test_GreedyPauli.cpp b/tket/test/src/test_GreedyPauli.cpp
index c803514183..76e0c0cfc6 100644
--- a/tket/test/src/test_GreedyPauli.cpp
+++ b/tket/test/src/test_GreedyPauli.cpp
@@ -760,5 +760,74 @@ SCENARIO("Test GreedyPauliSimp pass construction") {
                  ->apply(cu));
   }
 }
+
+SCENARIO("Test GreedyPauliSimp with multiple trials and threads") {
+  GIVEN("Large circuit with ZZPhase") {
+    Circuit circ(6);
+    circ.add_box(
+        PauliExpBox(SymPauliTensor({Pauli::X, Pauli::X}, 0.3)), {0, 1});
+    circ.add_box(
+        PauliExpBox(SymPauliTensor({Pauli::Z, Pauli::Y, Pauli::X}, 0.2)),
+        {0, 1, 2});
+    circ.add_box(
+        PauliExpCommutingSetBox({
+            {{Pauli::I, Pauli::Y, Pauli::I, Pauli::Z}, 1.2},
+            {{Pauli::X, Pauli::Y, Pauli::Z, Pauli::I}, 0.8},
+            {{Pauli::I, Pauli::I, Pauli::I, Pauli::Z}, 1.25},
+        }),
+        {1, 2, 3, 4});
+    circ.add_box(
+        PauliExpBox(SymPauliTensor({Pauli::Y, Pauli::X}, 0.1)), {2, 3});
+    circ.add_box(
+        PauliExpBox(SymPauliTensor({Pauli::Z, Pauli::Y, Pauli::X}, 0.11)),
+        {1, 3, 4});
+    circ.add_box(
+        PauliExpBox(SymPauliTensor({Pauli::Y, Pauli::Y}, 0.2)), {4, 5});
+    circ.add_box(
+        PauliExpBox(SymPauliTensor({Pauli::Z, Pauli::Z, Pauli::X}, 0.15)),
+        {2, 4, 5});
+    circ.add_box(
+        PauliExpBox(
+            SymPauliTensor({Pauli::X, Pauli::X, Pauli::X, Pauli::X}, 0.25)),
+        {2, 4, 5, 0});
+    circ.add_box(
+        PauliExpBox(
+            SymPauliTensor({Pauli::Y, Pauli::Z, Pauli::Z, Pauli::X}, 0.125)),
+        {1, 3, 5, 0});
+
+    circ.add_box(
+        PauliExpBox(SymPauliTensor(
+            {Pauli::X, Pauli::Z, Pauli::Y, Pauli::Y, Pauli::Z, Pauli::X},
+            0.125)),
+        {1, 3, 5, 0, 2, 4});
+
+    circ.add_box(
+        PauliExpBox(SymPauliTensor(
+            {Pauli::Z, Pauli::Y, Pauli::Y, Pauli::Z, Pauli::Z, Pauli::X},
+            0.125)),
+        {0, 1, 2, 3, 4, 5});
+
+    circ.add_box(
+        PauliExpBox(SymPauliTensor(
+            {Pauli::X, Pauli::Z, Pauli::Y, Pauli::Z, Pauli::Z, Pauli::X},
+            0.125)),
+        {5, 2, 4, 1, 3, 0});
+
+    circ.add_box(
+        PauliExpBox(SymPauliTensor(
+            {Pauli::X, Pauli::Z, Pauli::Y, Pauli::Y, Pauli::Z, Pauli::X},
+            0.125)),
+        {0, 5, 1, 4, 3, 2});
+
+    Circuit d(circ);
+    REQUIRE(!Transforms::greedy_pauli_optimisation(
+                 0.7, 0.3, 500, 500, 0, true, 0, 10)
+                 .apply(d));
+    REQUIRE(Transforms::greedy_pauli_optimisation(
+                0.7, 0.3, 500, 500, 0, true, 10, 10)
+                .apply(d));
+    REQUIRE(test_unitary_comparison(circ, d, true));
+  }
+}
 }  // namespace test_GreedyPauliSimp
 }  // namespace tket