diff --git a/docs/source/refs.bib b/docs/source/refs.bib
index 383049df0..bf7e889d9 100644
--- a/docs/source/refs.bib
+++ b/docs/source/refs.bib
@@ -669,6 +669,14 @@ @article{Sagastizabal_2019_PRA
   url       = {https://link.aps.org/doi/10.1103/PhysRevA.100.010302},
 }
 
+@misc{saki2023hypothesis,
+      title={Hypothesis Testing for Error Mitigation: How to Evaluate Error Mitigation}, 
+      author={Abdullah Ash Saki and Amara Katabarwa and Salonik Resch and George Umbrarescu},
+      year={2023},
+      eprint={2301.02690},
+      archivePrefix={arXiv},
+      primaryClass={quant-ph}
+}
 
 
 @article{Serbyn_2021_NatPhys,
diff --git a/mitiq/pt/__init__.py b/mitiq/pt/__init__.py
index 046ce9bca..e735c5ec8 100644
--- a/mitiq/pt/__init__.py
+++ b/mitiq/pt/__init__.py
@@ -4,7 +4,7 @@
 # LICENSE file in the root directory of this source tree.
 
 from mitiq.pt.pt import (
-    execute_with_pauli_twirling,
+    pauli_twirl_circuit,
     twirl_CNOT_gates,
     twirl_CZ_gates,
 )
diff --git a/mitiq/pt/pt.py b/mitiq/pt/pt.py
index 59a33fdd1..c7c07f390 100644
--- a/mitiq/pt/pt.py
+++ b/mitiq/pt/pt.py
@@ -4,12 +4,11 @@
 # LICENSE file in the root directory of this source tree.
 
 import random
-from typing import Callable, List, Optional, Union, cast
+from typing import List
 
 import cirq
-import numpy as np
 
-from mitiq import QPROGRAM, Executor, Observable, QuantumResult
+from mitiq import QPROGRAM
 from mitiq.interface import accept_qprogram_and_validate
 
 # P, Q, R, S from https://arxiv.org/pdf/2301.02690.pdf
@@ -51,38 +50,29 @@
 ]
 
 
-def execute_with_pauli_twirling(
+def pauli_twirl_circuit(
     circuit: QPROGRAM,
-    executor: Union[Executor, Callable[[QPROGRAM], QuantumResult]],
-    observable: Optional[Observable] = None,
-    *,
     num_circuits: int = 10,
-) -> float:
-    """Estimates the expectation value of the input circuit by averaging
-    expectation values obtained from Pauli twirled circuits.
+) -> List[QPROGRAM]:
+    r"""Return the Pauli twirled versions of the input circuit.
+
+    Only the $\mathrm{CZ}$ and $\mathrm{CNOT}$ gates in an
+    input circuit are Pauli twirled as specified in
+    :cite:`saki2023hypothesis`.
 
     Args:
         circuit: The input circuit to execute with twirling.
-        executor: A Mitiq executor that executes a circuit and returns the
-            unmitigated ``QuantumResult`` (e.g. an expectation value).
-        observable: Observable to compute the expectation value of. If
-            ``None``, the ``executor`` must return an expectation value
-            (float). Otherwise, the ``QuantumResult`` returned by ``executor``
-            is used to compute the expectation of the observable.
         num_circuits: Number of circuits to be twirled, and averaged.
 
     Returns:
         The expectation value estimated with Pauli twirling.
     """
-    executor = (
-        executor if isinstance(executor, Executor) else Executor(executor)
-    )
     CNOT_twirled_circuits = twirl_CNOT_gates(circuit, num_circuits)
     twirled_circuits = [
         twirl_CZ_gates(c, num_circuits=1)[0] for c in CNOT_twirled_circuits
     ]
-    expvals = executor.evaluate(twirled_circuits, observable)
-    return cast(float, np.average(expvals))
+
+    return twirled_circuits
 
 
 def twirl_CNOT_gates(circuit: QPROGRAM, num_circuits: int) -> List[QPROGRAM]:
diff --git a/mitiq/pt/tests/test_pt.py b/mitiq/pt/tests/test_pt.py
index 858f72c9a..e707d3b65 100644
--- a/mitiq/pt/tests/test_pt.py
+++ b/mitiq/pt/tests/test_pt.py
@@ -8,6 +8,7 @@
 import cirq
 import networkx as nx
 import numpy as np
+import pytest
 import qiskit
 
 from mitiq.benchmarks import generate_mirror_circuit
@@ -15,10 +16,11 @@
 from mitiq.pt.pt import (
     CNOT_twirling_gates,
     CZ_twirling_gates,
-    execute_with_pauli_twirling,
+    pauli_twirl_circuit,
     twirl_CNOT_gates,
     twirl_CZ_gates,
 )
+from mitiq.utils import _equal
 
 num_qubits = 2
 qubits = cirq.LineQubit.range(num_qubits)
@@ -116,7 +118,7 @@ def test_twirl_CNOT_increases_layer_count():
         assert num_gates_after == num_gates_before
 
 
-def test_execute_with_pauli_twirling():
+def test_pauli_twirl_circuit():
     num_qubits = 3
     num_layers = 20
     circuit, _ = generate_mirror_circuit(
@@ -124,7 +126,21 @@ def test_execute_with_pauli_twirling():
         two_qubit_gate_prob=1.0,
         connectivity_graph=nx.complete_graph(num_qubits),
     )
-    expval = execute_with_pauli_twirling(
-        circuit, amp_damp_executor, num_circuits=10
-    )
-    assert 0 <= expval < 0.5
+    num_circuits = 10
+    twirled_output = pauli_twirl_circuit(circuit, num_circuits)
+    assert len(twirled_output) == num_circuits
+
+
+@pytest.mark.parametrize(
+    "twirl_func", [pauli_twirl_circuit, twirl_CNOT_gates, twirl_CZ_gates]
+)
+def test_no_CNOT_CZ_circuit(twirl_func):
+    num_qubits = 2
+    qubits = cirq.LineQubit.range(num_qubits)
+    circuit = cirq.Circuit()
+    circuit.append(cirq.X.on_each(qubits))
+    twirled_output = twirl_func(circuit, 5)
+    assert len(twirled_output) == 5
+
+    for i in range(5):
+        assert _equal(circuit, twirled_output[i])