diff --git a/CHANGELOG.md b/CHANGELOG.md
index a0f40251f..23a096e40 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -98,7 +98,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - #468, #511: The `pyzx` module has been moved to a separate plugin: https://github.com/thierry-martinez/graphix-pyzx/
   Consequently, the `pyproject.toml` no longer defines an `extra` dependency group for the `pyzx` package.
 
-- #512: Method `Circuit.simulate_statevector` accepts a `backend: DenseStateBackend[_DenseStateT] | Literal["statevector", "densitymatrix"]` parameter.
+- #512: Method `Circuit.simulate` (formerly `Circuit.simulate_statevector`, see #567)  accepts a `backend: DenseStateBackend[_DenseStateT] | Literal["statevector", "densitymatrix"]` parameter.
+
+- #567:
+  - Method `Circuit.simulate_statevector` has been renamed `Circuit.simulate`.
+  - Method `Pattern.simulate_pattern` has been renamed `Pattern.simulate`.
 
 ## [0.3.5] - 2026-03-26
 
diff --git a/README.md b/README.md
index 86e6c5303..e1b00257f 100644
--- a/README.md
+++ b/README.md
@@ -74,7 +74,7 @@ pattern.draw_graph()
 ### simulating the pattern
 
 ```python
-state_out = pattern.simulate_pattern(backend="statevector")
+state_out = pattern.simulate(backend="statevector")
 ```
 
 ### and more..
diff --git a/benchmarks/statevec.py b/benchmarks/statevec.py
index 73e1c37c6..0944528b9 100644
--- a/benchmarks/statevec.py
+++ b/benchmarks/statevec.py
@@ -5,7 +5,7 @@
 Here we benchmark our statevector simulator for MBQC.
 
 The methods and modules we use are the followings:
-    1. :meth:`graphix.pattern.Pattern.simulate_pattern`
+    1. :meth:`graphix.pattern.Pattern.simulate`
         Pattern simulator with statevector backend.
     2. :mod:`paddle_quantum.mbqc`
         Pattern simulation using :mod:`paddle_quantum.mbqc`.
@@ -87,12 +87,12 @@ def simple_random_circuit(nqubit, depth):
     pattern.minimize_space()
     nqubit = len(pattern.extract_nodes())
     start = perf_counter()
-    pattern.simulate_pattern(max_qubit_num=30)
+    pattern.simulate(max_qubit_num=30)
     end = perf_counter()
     print(f"width: {width}, nqubit: {nqubit}, depth: {DEPTH}, time: {end - start}")
     pattern_time.append(end - start)
     start = perf_counter()
-    circuit.simulate_statevector()
+    circuit.simulate()
     end = perf_counter()
     circuit_time.append(end - start)
 
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 3f25f64cd..8ee360ba4 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -35,7 +35,9 @@
 autosectionlabel_prefix_document = True
 
 intersphinx_mapping = {
+    "python": ("https://docs.python.org/3", None),
     "networkx": ("https://networkx.github.io/documentation/stable/", None),
+    "sympy": ("https://docs.sympy.org/latest/", None),
 }
 
 sys.path.insert(0, os.path.abspath("../../"))
diff --git a/docs/source/generator.rst b/docs/source/generator.rst
index 97918ae39..da89a0a81 100644
--- a/docs/source/generator.rst
+++ b/docs/source/generator.rst
@@ -9,38 +9,7 @@ Pattern Generation
 .. currentmodule:: graphix.transpiler
 
 .. autoclass:: Circuit
-
-    .. automethod:: __init__
-
-    .. automethod:: transpile
-
-    .. automethod:: simulate_statevector
-
-    .. automethod:: cz
-
-    .. automethod:: cnot
-
-    .. automethod:: h
-
-    .. automethod:: s
-
-    .. automethod:: x
-
-    .. automethod:: y
-
-    .. automethod:: z
-
-    .. automethod:: rx
-
-    .. automethod:: ry
-
-    .. automethod:: rz
-
-    .. automethod:: j
-
-    .. automethod:: ccx
-
-    .. automethod:: m
+    :members:
 
 .. autoclass:: TranspiledPattern
 
diff --git a/docs/source/intro.rst b/docs/source/intro.rst
index 3fa3686c8..0d47c4367 100644
--- a/docs/source/intro.rst
+++ b/docs/source/intro.rst
@@ -132,7 +132,7 @@ Note that the input state has *teleported* to qubits 6 and 7 after the computati
         import networkx as nx
         og = OpenGraph(nx.Graph([0, 1]), output_nodes=[0, 1])
 
-    >>> print(og.to_pattern().simulate_pattern())
+    >>> print(og.to_pattern().simulate())
     Statevec object with statevector [[ 0.5+0.j  0.5+0.j]
     [ 0.5+0.j -0.5+0.j]] and length (2, 2).
 
diff --git a/docs/source/modifier.rst b/docs/source/modifier.rst
index 2d3b96073..2f3197419 100644
--- a/docs/source/modifier.rst
+++ b/docs/source/modifier.rst
@@ -7,84 +7,6 @@ Pattern Manipulation
 .. currentmodule:: graphix.pattern
 
 .. autoclass:: Pattern
-
-    .. automethod:: __init__
-
-    .. automethod:: add
-
-    .. automethod:: extend
-
-    .. automethod:: clear
-
-    .. automethod:: replace
-
-    .. automethod:: reorder_output_nodes
-
-    .. automethod:: reorder_input_nodes
-
-    .. automethod:: simulate_pattern
-
-    .. automethod:: compute_max_degree
-
-    .. automethod:: compose
-
-    .. automethod:: connected_nodes
-
-    .. automethod:: remove_input_nodes
-
-    .. automethod:: perform_pauli_pushing
-
-    .. automethod:: remove_pauli_measurements
-
-    .. automethod:: remove_local_clifford_commands
-
-    .. automethod:: to_ascii
-
-    .. automethod:: to_unicode
-
-    .. automethod:: to_latex
-
-    .. automethod:: standardize
-
-    .. automethod:: shift_signals
-
-    .. automethod:: is_standard
-
-    .. automethod:: extract_graph
-
-    .. automethod:: extract_nodes
-
-    .. automethod:: extract_causal_flow
-
-    .. automethod:: extract_gflow
-
-    .. automethod:: extract_opengraph
-
-    .. automethod:: extract_measurement_commands
-
-    .. automethod:: parallelize_pattern
-
-    .. automethod:: minimize_space
-
-    .. automethod:: draw
-
-    .. automethod:: max_space
-
-    .. automethod:: to_qasm3
-
-    .. automethod:: is_parameterized
-
-    .. automethod:: subs
-
-    .. automethod:: xreplace
-
-    .. automethod:: check_runnability
-
-    .. automethod:: map
-
-    .. automethod:: infer_pauli_measurements
-
-    .. automethod:: to_bloch
-
+    :members:
 
 .. autofunction:: shift_outcomes
diff --git a/docs/source/tutorial.rst b/docs/source/tutorial.rst
index a828719e7..431f345f5 100644
--- a/docs/source/tutorial.rst
+++ b/docs/source/tutorial.rst
@@ -51,9 +51,9 @@ This is the simplest of the `graph state
 <https://en.wikipedia.org/wiki/Graph_state>`_ with nodes = [0, 1] and edge = (0, 1). Any MBQC pattern has a corresponding resource graph state on which the computation occurs only with single-qubit measurements.
 
 We can use the :class:`~graphix.simulator.PatternSimulator` to classically simulate the pattern above and obtain the output state, for default input state of :math:`|+\rangle`.
-Alternatively, we can simply call :meth:`~graphix.pattern.Pattern.simulate_pattern` of :class:`~graphix.pattern.Pattern` object to do it in one line:
+Alternatively, we can simply call :meth:`~graphix.pattern.Pattern.simulate` of :class:`~graphix.pattern.Pattern` object to do it in one line:
 
->>> print(pattern.simulate_pattern(backend='statevector'))
+>>> print(pattern.simulate(backend='statevector'))
 Statevec object with statevector [1.+0.j 0.+0.j] and length (2,).
 
 Note again that we started with :math:`|+\rangle` state so the answer is correct.
@@ -240,7 +240,7 @@ We can simulate the MBQC pattern with various noise models to understand their e
 
 .. code-block:: python
 
-    out_state = pattern.simulate_pattern(backend="statevector")
+    out_state = pattern.simulate(backend="statevector")
 
 The statevector backend simulates the pattern ideally, without any noise.
 
diff --git a/examples/deutsch_jozsa.py b/examples/deutsch_jozsa.py
index 5f4d6dee4..dbae1e82c 100644
--- a/examples/deutsch_jozsa.py
+++ b/examples/deutsch_jozsa.py
@@ -90,8 +90,8 @@
 # So the preprocessing has done all the necessary computations, and all nodes are isolated with no further measurements required.
 # Let us make sure the result is correct:
 
-out_state = pattern.simulate_pattern()
-state = circuit.simulate_statevector().statevec
+out_state = pattern.simulate()
+state = circuit.simulate().statevec
 print("overlap of states: ", np.abs(np.dot(state.psi.flatten().conjugate(), out_state.psi.flatten())))
 
 # %%
diff --git a/examples/ghz_with_tn.py b/examples/ghz_with_tn.py
index 2834c4e0e..67d9e1ade 100644
--- a/examples/ghz_with_tn.py
+++ b/examples/ghz_with_tn.py
@@ -45,7 +45,7 @@
 # %%
 # Calculate the amplitudes of ``|00...0>`` and ``|11...1>`` states.
 
-tn = pattern.simulate_pattern(backend="tensornetwork")
+tn = pattern.simulate(backend="tensornetwork")
 print(f"The amplitude of |00...0>: {tn.basis_amplitude(0)}")
 print(f"The amplitude of |11...1>: {tn.basis_amplitude(2**n - 1)}")
 
diff --git a/examples/qaoa.py b/examples/qaoa.py
index da993fe37..da6212c28 100644
--- a/examples/qaoa.py
+++ b/examples/qaoa.py
@@ -47,8 +47,8 @@
 # %%
 # finally, simulate the QAOA circuit
 
-out_state = pattern.simulate_pattern()
-state = circuit.simulate_statevector().statevec
+out_state = pattern.simulate()
+state = circuit.simulate().statevec
 print("overlap of states: ", np.abs(np.dot(state.psi.flatten().conjugate(), out_state.psi.flatten())))
 # sphinx_gallery_thumbnail_number = 2
 
diff --git a/examples/qft_with_tn.py b/examples/qft_with_tn.py
index 5cb991c2a..40ad1cb5d 100644
--- a/examples/qft_with_tn.py
+++ b/examples/qft_with_tn.py
@@ -78,7 +78,7 @@ def qft(circuit: Circuit, n: int) -> None:
 import time
 
 t1 = time.time()
-tn = pattern.simulate_pattern(backend="tensornetwork")
+tn = pattern.simulate(backend="tensornetwork")
 value = tn.basis_amplitude(0)
 t2 = time.time()
 print("amplitude of |00...0> is ", value)
diff --git a/examples/rotation.py b/examples/rotation.py
index 0de287855..278af6751 100644
--- a/examples/rotation.py
+++ b/examples/rotation.py
@@ -68,7 +68,7 @@
 # Internally, we are executing the command sequence we inspected above on a statevector simulator.
 # We also have a tensornetwork simulation backend to handle larger MBQC patterns. see other examples for how to use it.
 
-out_state = pattern.simulate_pattern(backend="statevector")
+out_state = pattern.simulate(backend="statevector")
 print(out_state.flatten())
 
 # %%
diff --git a/examples/tn_simulation.py b/examples/tn_simulation.py
index 84a5e8b19..1def877c6 100644
--- a/examples/tn_simulation.py
+++ b/examples/tn_simulation.py
@@ -93,7 +93,7 @@ def ansatz(
 # The graph_prep argument is optional,
 # but with 'parallel' the TensorNetworkBackend will prepeare the graph state faster.
 
-mbqc_tn = pattern.simulate_pattern(backend="tensornetwork", graph_prep="parallel")
+mbqc_tn = pattern.simulate(backend="tensornetwork", graph_prep="parallel")
 sv = mbqc_tn.to_statevector().flatten()
 print("Statevector after the simulation:", sv)
 
@@ -206,7 +206,7 @@ def cost(
     pattern.remove_input_nodes()
     pattern = pattern.infer_pauli_measurements()
     pattern.remove_pauli_measurements(standardize=True)
-    mbqc_tn = pattern.simulate_pattern(backend="tensornetwork", graph_prep="parallel")
+    mbqc_tn = pattern.simulate(backend="tensornetwork", graph_prep="parallel")
     exp_val: float = 0
     for op in ham:
         exp_val += np.real(mbqc_tn.expectation_value(op, range(n), optimize=opt))
@@ -248,7 +248,7 @@ def __call__(
 pattern = circuit.transpile().pattern
 pattern.standardize()
 pattern.shift_signals()
-mbqc_tn = pattern.simulate_pattern(backend="tensornetwork", graph_prep="parallel")
+mbqc_tn = pattern.simulate(backend="tensornetwork", graph_prep="parallel")
 
 # %%
 # Let's use the defined optimizer and find the most probable basis states.
diff --git a/graphix/pattern.py b/graphix/pattern.py
index ddfebcf3d..b04c36c86 100644
--- a/graphix/pattern.py
+++ b/graphix/pattern.py
@@ -389,7 +389,7 @@ def input_nodes(self) -> list[int]:
 
     @property
     def output_nodes(self) -> list[int]:
-        """List all nodes that are either `input_nodes` or prepared with `N` commands and that have not been measured with an `M` command."""
+        """List all nodes that are either ``input_nodes`` or prepared with ``N`` commands and that have not been measured with an ``M`` command."""
         return list(self.__output_nodes)  # copy for preventing modification
 
     def __len__(self) -> int:
@@ -1118,7 +1118,7 @@ def extract_pauli_flow(self) -> PauliFlow[Measurement]:
         return self.extract_xzcorrections().to_pauli_flow()
 
     def extract_xzcorrections(self) -> XZCorrections[Measurement]:
-        """Extract the XZ-corrections from the current measurement pattern.
+        r"""Extract the XZ-corrections from the current measurement pattern.
 
         Returns
         -------
@@ -1130,7 +1130,7 @@ def extract_xzcorrections(self) -> XZCorrections[Measurement]:
         XZCorrectionsError
             If the extracted correction dictionaries are not well formed.
         ValueError
-            If `N` commands in the pattern do not represent a |+⟩ state or if the pattern corrections form closed loops.
+            If ``N`` commands in the pattern do not represent a :math:`\ket{+}` state or if the pattern corrections form closed loops.
 
         Notes
         -----
@@ -1398,7 +1398,7 @@ def space_list(self) -> list[int]:
         return n_list
 
     @overload
-    def simulate_pattern(
+    def simulate(
         self,
         backend: StatevectorBackend | Literal["statevector"] = "statevector",
         input_state: State
@@ -1412,7 +1412,7 @@ def simulate_pattern(
     ) -> Statevec: ...
 
     @overload
-    def simulate_pattern(
+    def simulate(
         self,
         backend: DensityMatrixBackend | Literal["densitymatrix"],
         input_state: State
@@ -1426,7 +1426,7 @@ def simulate_pattern(
     ) -> DensityMatrix: ...
 
     @overload
-    def simulate_pattern(
+    def simulate(
         self,
         backend: TensorNetworkBackend | Literal["tensornetwork", "mps"],
         input_state: State
@@ -1439,7 +1439,7 @@ def simulate_pattern(
     ) -> MBQCTensorNet: ...
 
     @overload
-    def simulate_pattern(
+    def simulate(
         self,
         backend: Backend[_StateT_co],
         input_state: Data | None = ...,
@@ -1447,7 +1447,7 @@ def simulate_pattern(
         **kwargs: Any,
     ) -> _StateT_co: ...
 
-    def simulate_pattern(
+    def simulate(
         self,
         backend: Backend[_StateT_co] | _BackendLiteral = "statevector",
         input_state: Data | None = BasicStates.PLUS,
diff --git a/graphix/sim/statevec.py b/graphix/sim/statevec.py
index cb45e1eb5..a970976e0 100644
--- a/graphix/sim/statevec.py
+++ b/graphix/sim/statevec.py
@@ -574,7 +574,7 @@ def draw(
         >>> circuit = Circuit(2)
         >>> circuit.h(0)
         >>> circuit.cz(0, 1)
-        >>> print(circuit.simulate_statevector().statevec.draw())
+        >>> print(circuit.simulate().statevec.draw())
         √2/2(|00⟩ + |01⟩)
         """
         return statevec_to_str(
diff --git a/graphix/transpiler.py b/graphix/transpiler.py
index aa49d76bb..724e854f5 100644
--- a/graphix/transpiler.py
+++ b/graphix/transpiler.py
@@ -525,7 +525,7 @@ def transpile(self, *, transpile_swaps: bool = True) -> TranspiledPattern:
         return TranspiledPattern(result.pattern, classical_outputs)
 
     @overload
-    def simulate_statevector(
+    def simulate(
         self,
         backend: StatevectorBackend | Literal["statevector"] = ...,
         input_state: Data | None = None,
@@ -536,7 +536,7 @@ def simulate_statevector(
     ) -> SimulateResult[Statevec]: ...
 
     @overload
-    def simulate_statevector(
+    def simulate(
         self,
         backend: DensityMatrixBackend | Literal["densitymatrix"],
         input_state: Data | None = None,
@@ -547,7 +547,7 @@ def simulate_statevector(
     ) -> SimulateResult[DensityMatrix]: ...
 
     @overload
-    def simulate_statevector(
+    def simulate(
         self,
         backend: DenseStateBackend[_DenseStateT],
         input_state: Data | None = None,
@@ -557,7 +557,7 @@ def simulate_statevector(
         stacklevel: int = 1,
     ) -> SimulateResult[_DenseStateT]: ...
 
-    def simulate_statevector(
+    def simulate(
         self,
         backend: DenseStateBackend[_DenseStateT] | _DenseStateBackendLiteral = "statevector",
         input_state: Data | None = None,
@@ -649,24 +649,24 @@ def evolve(op: Matrix, qargs: Iterable[int]) -> None:
         return SimulateResult(_backend.state, tuple(classical_measures))
 
     def visit(self, visitor: InstructionVisitor) -> Circuit:
-        """Apply `visitor` to all instructions in the circuit."""
+        """Apply ``visitor`` to all instructions in the circuit."""
         result = Circuit(self.width)
         for instr in self.instruction:
             result.instruction.append(instr.visit(visitor))
         return result
 
     def map_angle(self, f: Callable[[ParameterizedAngle], ParameterizedAngle]) -> Circuit:
-        """Apply `f` to all angles that occur in the circuit."""
+        """Apply ``f`` to all angles that occur in the circuit."""
         return self.visit(_MapAngleVisitor(f))
 
     def is_parameterized(self) -> bool:
         """
-        Return `True` if there is at least one measurement angle that is not just an instance of `SupportsFloat`.
+        Return ``True`` if there is at least one measurement angle that is not just an instance of :class:`SupportsFloat`.
 
         A parameterized circuit is a circuit where at least one
         measurement angle is an expression that is not a number,
-        typically an instance of `sympy.Expr` (but we don't force to
-        choose `sympy` here).
+        typically an instance of :class:`sympy.Expr` (but we don't force to
+        choose ``sympy`` here).
 
         """
         for instr in self.instruction:
diff --git a/noxfile.py b/noxfile.py
index ab29a86c6..8523d0344 100644
--- a/noxfile.py
+++ b/noxfile.py
@@ -95,18 +95,21 @@ class ReverseDependency:
 @nox.parametrize(
     "package",
     [
-        ReverseDependency("https://github.com/TeamGraphix/graphix-stim-backend"),
-        ReverseDependency("https://github.com/TeamGraphix/graphix-symbolic"),
+        ReverseDependency("https://github.com/thierry-martinez/graphix-stim-backend", branch="rename-simulate"),
+        ReverseDependency("https://github.com/thierry-martinez/graphix-symbolic", branch="rename-simulate"),
         ReverseDependency("https://github.com/TeamGraphix/graphix-qasm-parser"),
-        ReverseDependency("https://github.com/TeamGraphix/graphix-ibmq", doctest_modules=False),
-        ReverseDependency("https://github.com/TeamGraphix/graphix-stim-compiler"),
-        ReverseDependency("https://github.com/TeamGraphix/graphix-pyzx"),
         ReverseDependency(
-            "https://github.com/qat-inria/veriphix",
+            "https://github.com/thierry-martinez/graphix-ibmq", doctest_modules=False, branch="rename-simulate"
+        ),
+        ReverseDependency("https://github.com/thierry-martinez/graphix-stim-compiler", branch="rename-simulate"),
+        ReverseDependency("https://github.com/thierry-martinez/graphix-pyzx", branch="rename-simulate"),
+        ReverseDependency(
+            "https://github.com/thierry-martinez/veriphix",
             doctest_modules=False,
             install_target=".[dev]",
+            branch="rename-simulate",
         ),
-        ReverseDependency("https://github.com/matulni/graphix-mqtbench", branch="minimal"),
+        ReverseDependency("https://github.com/thierry-martinez/graphix-mqtbench", branch="rename-simulate"),
     ],
 )
 def tests_reverse_dependencies(session: Session, package: ReverseDependency) -> None:
diff --git a/tests/test_branch_selector.py b/tests/test_branch_selector.py
index 1dfaa6b4a..bc4557cc8 100644
--- a/tests/test_branch_selector.py
+++ b/tests/test_branch_selector.py
@@ -62,7 +62,7 @@ def test_expectation_value(fx_rng: Generator, backend: _BackendLiteral) -> None:
     # Pattern that measures 0 on qubit 0 with probability 1.
     pattern = Pattern(cmds=[N(0), M(0)])
     branch_selector = CheckedBranchSelector(expected={0: 1.0})
-    pattern.simulate_pattern(backend, branch_selector=branch_selector, rng=fx_rng)
+    pattern.simulate(backend, branch_selector=branch_selector, rng=fx_rng)
 
 
 @pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.")
@@ -84,7 +84,7 @@ def test_random_branch_selector(fx_rng: Generator, backend: _BackendLiteral) ->
     pattern = Pattern(cmds=[N(0), M(0)])
     for _ in range(NB_ROUNDS):
         measure_method = DefaultMeasureMethod()
-        pattern.simulate_pattern(backend, branch_selector=branch_selector, measure_method=measure_method, rng=fx_rng)
+        pattern.simulate(backend, branch_selector=branch_selector, measure_method=measure_method, rng=fx_rng)
         assert measure_method.results[0] == 0
 
 
@@ -104,7 +104,7 @@ def test_random_branch_selector_without_pr_calc(fx_rng: Generator, backend: _Bac
     nb_outcome_1 = 0
     for _ in range(NB_ROUNDS):
         measure_method = DefaultMeasureMethod()
-        pattern.simulate_pattern(backend, branch_selector=branch_selector, measure_method=measure_method, rng=fx_rng)
+        pattern.simulate(backend, branch_selector=branch_selector, measure_method=measure_method, rng=fx_rng)
         if measure_method.results[0]:
             nb_outcome_1 += 1
     assert abs(nb_outcome_1 - NB_ROUNDS / 2) < NB_ROUNDS / 5
@@ -126,7 +126,7 @@ def test_fixed_branch_selector(backend: _BackendLiteral, outcome: list[Outcome])
     branch_selector = FixedBranchSelector(results1, default=FixedBranchSelector(results2))
     pattern = Pattern(cmds=[cmd for qubit in range(3) for cmd in (N(qubit), M(qubit, Measurement.XY(0.1)))])
     measure_method = DefaultMeasureMethod()
-    pattern.simulate_pattern(backend, branch_selector=branch_selector, measure_method=measure_method)
+    pattern.simulate(backend, branch_selector=branch_selector, measure_method=measure_method)
     for qubit, value in enumerate(outcome):
         assert measure_method.results[qubit] == value
 
@@ -146,7 +146,7 @@ def test_fixed_branch_selector_no_default(backend: _BackendLiteral) -> None:
     pattern = Pattern(cmds=[N(0), M(0, Measurement.XY(1e-5))])
     measure_method = DefaultMeasureMethod()
     with pytest.raises(ValueError):
-        pattern.simulate_pattern(backend, branch_selector=branch_selector, measure_method=measure_method)
+        pattern.simulate(backend, branch_selector=branch_selector, measure_method=measure_method)
 
 
 @pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.")
@@ -164,5 +164,5 @@ def test_const_branch_selector(backend: _BackendLiteral, outcome: Outcome) -> No
     pattern = Pattern(cmds=[N(0), M(0, Measurement.XY(1e-5))])
     for _ in range(NB_ROUNDS):
         measure_method = DefaultMeasureMethod()
-        pattern.simulate_pattern(backend, branch_selector=branch_selector, measure_method=measure_method)
+        pattern.simulate(backend, branch_selector=branch_selector, measure_method=measure_method)
         assert measure_method.results[0] == outcome
diff --git a/tests/test_circ_extraction.py b/tests/test_circ_extraction.py
index cdd6e1549..33005e34b 100644
--- a/tests/test_circ_extraction.py
+++ b/tests/test_circ_extraction.py
@@ -107,8 +107,8 @@ class TestPauliExponentialDAG:
     def test_to_circuit(self, test_case: PauliExpTestCase, fx_rng: Generator) -> None:
         qc = Circuit(len(test_case.pexp_dag.output_nodes))
         pexp_ladder_pass(test_case.pexp_dag, qc)
-        state = qc.simulate_statevector(rng=fx_rng).statevec
-        state_ref = test_case.qc.simulate_statevector(rng=fx_rng).statevec
+        state = qc.simulate(rng=fx_rng).statevec
+        state_ref = test_case.qc.simulate(rng=fx_rng).statevec
         assert state.isclose(state_ref)
 
     def test_to_circuit_outputs_order(self, fx_rng: Generator) -> None:
@@ -131,17 +131,17 @@ def test_to_circuit_outputs_order(self, fx_rng: Generator) -> None:
 
         qc_1 = Circuit(2)
         pexp_ladder_pass(pexp_dag_1, qc_1)
-        s_1 = qc_1.simulate_statevector(rng=fx_rng, input_state=[BasicStates.PLUS, BasicStates.MINUS]).statevec
+        s_1 = qc_1.simulate(rng=fx_rng, input_state=[BasicStates.PLUS, BasicStates.MINUS]).statevec
 
         qc_2 = Circuit(2)
         qc_2.swap(0, 1)  # We must swap before and after the Pauli exponential!
         pexp_ladder_pass(pexp_dag_2, qc_2)
 
-        s_2 = qc_2.simulate_statevector(rng=fx_rng, input_state=[BasicStates.PLUS, BasicStates.MINUS]).statevec
+        s_2 = qc_2.simulate(rng=fx_rng, input_state=[BasicStates.PLUS, BasicStates.MINUS]).statevec
         assert not s_1.isclose(s_2)
 
         qc_2.swap(0, 1)
-        s_2 = qc_2.simulate_statevector(rng=fx_rng, input_state=[BasicStates.PLUS, BasicStates.MINUS]).statevec
+        s_2 = qc_2.simulate(rng=fx_rng, input_state=[BasicStates.PLUS, BasicStates.MINUS]).statevec
 
         assert s_1.isclose(s_2)
 
@@ -431,8 +431,8 @@ def test_cm_berg_pass(self, qc_ref: Circuit, cm: CliffordMap, fx_rng: Generator)
         qc = Circuit(qc_ref.width)
         cm_berg_pass(cm, qc)
 
-        s_test = qc.simulate_statevector(rng=fx_rng).statevec
-        s_ref = qc_ref.simulate_statevector(rng=fx_rng).statevec
+        s_test = qc.simulate(rng=fx_rng).statevec
+        s_ref = qc_ref.simulate(rng=fx_rng).statevec
 
         assert s_test.isclose(s_ref)
 
@@ -448,8 +448,8 @@ def test_extract_rnd_circuit(self, fx_bg: PCG64, jumps: int) -> None:
 
         circuit = pattern.extract_opengraph().extract_circuit()
 
-        s_ref = circuit.simulate_statevector(rng=rng).statevec
-        s_test = circuit_ref.simulate_statevector(rng=rng).statevec
+        s_ref = circuit.simulate(rng=rng).statevec
+        s_test = circuit_ref.simulate(rng=rng).statevec
         assert s_ref.isclose(s_test)
 
     @pytest.mark.parametrize(
@@ -526,8 +526,8 @@ def test_extract_og(self, test_case: OpenGraph[Measurement], fx_rng: Generator)
         pattern = test_case.to_pattern()
         circuit = test_case.extract_circuit()
 
-        state = circuit.simulate_statevector(rng=fx_rng).statevec
-        state_ref = pattern.simulate_pattern(rng=fx_rng)
+        state = circuit.simulate(rng=fx_rng).statevec
+        state_ref = pattern.simulate(rng=fx_rng)
         assert state.isclose(state_ref)
 
     @pytest.mark.parametrize("infer_pauli", [True, False])
@@ -551,8 +551,8 @@ def test_extract_og_infer_pauli(self, infer_pauli: bool, fx_rng: Generator) -> N
 
         circuit = og.extract_circuit()
 
-        state = circuit.simulate_statevector(rng=fx_rng).statevec
-        state_ref = pattern.simulate_pattern(rng=fx_rng)
+        state = circuit.simulate(rng=fx_rng).statevec
+        state_ref = pattern.simulate(rng=fx_rng)
         assert state.isclose(state_ref)
 
     def test_extract_og_gflow(self, fx_rng: Generator) -> None:
@@ -570,8 +570,8 @@ def test_extract_og_gflow(self, fx_rng: Generator) -> None:
         pattern = og.to_pattern()
         circuit = og.extract_gflow().extract_circuit().to_circuit()
 
-        state = circuit.simulate_statevector(rng=fx_rng).statevec
-        state_ref = pattern.simulate_pattern(rng=fx_rng)
+        state = circuit.simulate(rng=fx_rng).statevec
+        state_ref = pattern.simulate(rng=fx_rng)
         assert state.isclose(state_ref)
 
     @pytest.mark.parametrize("test_case", [0.2, 0.5, 1.0])
@@ -592,13 +592,13 @@ def test_parametric_angles(self, test_case: float, fx_rng: Generator) -> None:
 
         # Substitute parameter at the level of the extracted circuit
         qc1 = og.extract_circuit()
-        s1 = qc1.subs(alpha, alpha_val).simulate_statevector(rng=fx_rng).statevec
+        s1 = qc1.subs(alpha, alpha_val).simulate(rng=fx_rng).statevec
 
         # Substitute parameter at the level of the open graph object
         # Calling `infer_pauli_measurements` is not necessary for the test to pass
         # (and it should not be), but it suppresses the warnings.
         qc2 = og.subs(alpha, alpha_val).infer_pauli_measurements().extract_circuit()
-        s2 = qc2.simulate_statevector(rng=fx_rng).statevec
+        s2 = qc2.simulate(rng=fx_rng).statevec
 
         assert s1.isclose(s2)
 
diff --git a/tests/test_clifford.py b/tests/test_clifford.py
index 06bd981bd..50ca15eb7 100644
--- a/tests/test_clifford.py
+++ b/tests/test_clifford.py
@@ -104,6 +104,6 @@ def test_to_pattern(self, fx_rng: Generator, c: Clifford) -> None:
         pattern = og.to_pattern()
         pattern_ref = Pattern(input_nodes=[0], cmds=[Command.C(0, c)])
         input_state = rand_state_vector(nqubits=1, rng=fx_rng)
-        state = pattern.simulate_pattern(input_state=input_state, rng=fx_rng)
-        state_ref = pattern_ref.simulate_pattern(input_state=input_state, rng=fx_rng)
+        state = pattern.simulate(input_state=input_state, rng=fx_rng)
+        state_ref = pattern_ref.simulate(input_state=input_state, rng=fx_rng)
         assert state.isclose(state_ref)
diff --git a/tests/test_db.py b/tests/test_db.py
index 90361cce8..5b983b351 100644
--- a/tests/test_db.py
+++ b/tests/test_db.py
@@ -85,7 +85,7 @@ def generate_clifford_pauli_decomposition(rng: Generator) -> tuple[tuple[PauliMe
         (
             clifford,
             tuple(
-                Pattern(input_nodes=[0], cmds=[Command.C(0, clifford)]).simulate_pattern(input_state=input_state)
+                Pattern(input_nodes=[0], cmds=[Command.C(0, clifford)]).simulate(input_state=input_state)
                 for input_state in input_states
             ),
         )
@@ -107,7 +107,7 @@ def explore(n: int) -> None:
                 if patterns[clifford.value] is not None:
                     continue
                 if all(
-                    pattern.simulate_pattern(input_state=input_state, rng=rng).isclose(output_state_ref)
+                    pattern.simulate(input_state=input_state, rng=rng).isclose(output_state_ref)
                     for input_state, output_state_ref in zip(input_states, output_states_ref, strict=True)
                 ):
                     patterns[clifford.value] = measurement_list
diff --git a/tests/test_flow_core.py b/tests/test_flow_core.py
index 4e52a9943..d881c9853 100644
--- a/tests/test_flow_core.py
+++ b/tests/test_flow_core.py
@@ -407,10 +407,10 @@ def test_corrections_to_pattern(self, test_case: XZCorrectionsTestCase, fx_rng:
 
             for plane in {Plane.XY, Plane.XZ, Plane.YZ}:
                 alpha = 2 * ANGLE_PI * fx_rng.random()
-                state_ref = test_case.pattern.simulate_pattern(input_state=PlanarState(plane, alpha), rng=fx_rng)
+                state_ref = test_case.pattern.simulate(input_state=PlanarState(plane, alpha), rng=fx_rng)
 
                 for _ in range(n_shots):
-                    state = pattern.simulate_pattern(input_state=PlanarState(plane, alpha), rng=fx_rng)
+                    state = pattern.simulate(input_state=PlanarState(plane, alpha), rng=fx_rng)
                     assert state.isclose(state_ref)
 
 
diff --git a/tests/test_noise_model.py b/tests/test_noise_model.py
index c3f3f9b58..a831e60b5 100644
--- a/tests/test_noise_model.py
+++ b/tests/test_noise_model.py
@@ -42,12 +42,12 @@ def test_noiseless_noise_model_simulation(fx_rng: Generator) -> None:
     nqubits = 5
     depth = 5
     circuit = rand_circuit(nqubits, depth, rng=fx_rng)
-    state = circuit.simulate_statevector().statevec
+    state = circuit.simulate().statevec
     pattern = circuit.transpile().pattern
     pattern.standardize()
     pattern.minimize_space()
     noise_model = NoiselessNoiseModel()
-    state_mbqc = pattern.simulate_pattern(backend="densitymatrix", noise_model=noise_model, rng=fx_rng)
+    state_mbqc = pattern.simulate(backend="densitymatrix", noise_model=noise_model, rng=fx_rng)
     assert np.abs(np.dot(state_mbqc.flatten().conjugate(), DensityMatrix(state).rho.flatten())) == pytest.approx(1)
 
 
@@ -94,13 +94,13 @@ def test_compose_noise_model_simulation(fx_rng: Generator) -> None:
     nqubits = 5
     depth = 5
     circuit = rand_circuit(nqubits, depth, rng=fx_rng)
-    state = circuit.simulate_statevector().statevec
+    state = circuit.simulate().statevec
     pattern = circuit.transpile().pattern
     pattern.standardize()
     pattern.minimize_space()
     # By default, `DepolarisingNoiseModel` is noiseless.
     noise_model = ComposeNoiseModel([NoiselessNoiseModel(), DepolarisingNoiseModel()])
-    state_mbqc = pattern.simulate_pattern(backend="densitymatrix", noise_model=noise_model, rng=fx_rng)
+    state_mbqc = pattern.simulate(backend="densitymatrix", noise_model=noise_model, rng=fx_rng)
     assert np.abs(np.dot(state_mbqc.flatten().conjugate(), DensityMatrix(state).rho.flatten())) == pytest.approx(1)
 
 
@@ -162,14 +162,12 @@ def test_confuse_result(fx_rng: Generator, noise_model: NoiseModel) -> None:
     # Pattern that measures 0 on qubit 0 with probability 1.
     pattern = Pattern(cmds=[N(0), M(0)])
     measure_method = DefaultMeasureMethod()
-    pattern.simulate_pattern(
+    pattern.simulate(
         backend="densitymatrix", noise_model=NoiselessNoiseModel(), rng=fx_rng, measure_method=measure_method
     )
     assert measure_method.results[0] == 0
     measure_method = DefaultMeasureMethod()
-    pattern.simulate_pattern(
-        backend="densitymatrix", noise_model=noise_model, rng=fx_rng, measure_method=measure_method
-    )
+    pattern.simulate(backend="densitymatrix", noise_model=noise_model, rng=fx_rng, measure_method=measure_method)
     assert measure_method.results[0] == 1
 
 
@@ -230,13 +228,13 @@ def test_compose_amplitude_damping_depolarising_simulation(fx_rng: Generator) ->
     nqubits = 5
     depth = 5
     circuit = rand_circuit(nqubits, depth, rng=fx_rng)
-    state = circuit.simulate_statevector().statevec
+    state = circuit.simulate().statevec
     pattern = circuit.transpile().pattern
     pattern.standardize()
     pattern.minimize_space()
     # both models default to noiseless (all probs 0)
     noise_model = ComposeNoiseModel([AmplitudeDampingNoiseModel(), DepolarisingNoiseModel()])
-    state_mbqc = pattern.simulate_pattern(backend="densitymatrix", noise_model=noise_model, rng=fx_rng)
+    state_mbqc = pattern.simulate(backend="densitymatrix", noise_model=noise_model, rng=fx_rng)
     assert np.abs(np.dot(state_mbqc.flatten().conjugate(), DensityMatrix(state).rho.flatten())) == pytest.approx(1)
 
 
diff --git a/tests/test_noisy_density_matrix.py b/tests/test_noisy_density_matrix.py
index e5521db9d..78af9397c 100644
--- a/tests/test_noisy_density_matrix.py
+++ b/tests/test_noisy_density_matrix.py
@@ -49,9 +49,9 @@ class TestNoisyDensityMatrixBackend:
     @pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.")
     def test_noiseless_noisy_hadamard(self, fx_rng: Generator) -> None:
         hadamardpattern = hpat()
-        noiselessres = hadamardpattern.simulate_pattern(backend="densitymatrix", rng=fx_rng)
+        noiselessres = hadamardpattern.simulate(backend="densitymatrix", rng=fx_rng)
         # noiseless noise model
-        noisynoiselessres = hadamardpattern.simulate_pattern(
+        noisynoiselessres = hadamardpattern.simulate(
             backend="densitymatrix",
             noise_model=NoiselessNoiseModel(),
             rng=fx_rng,
@@ -65,7 +65,7 @@ def test_noiseless_noisy_hadamard(self, fx_rng: Generator) -> None:
     # test measurement confuse outcome
     def test_noisy_measure_confuse_hadamard(self, fx_rng: Generator) -> None:
         hadamardpattern = hpat()
-        res = hadamardpattern.simulate_pattern(
+        res = hadamardpattern.simulate(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(measure_error_prob=1.0),
             rng=fx_rng,
@@ -80,7 +80,7 @@ def test_noisy_measure_confuse_hadamard_arbitrary(self, fx_rng: Generator, outco
         hadamardpattern = hpat()
         measure_error_pr = fx_rng.random()
         print(f"measure_error_pr = {measure_error_pr}, outcome = {outcome}")
-        res = hadamardpattern.simulate_pattern(
+        res = hadamardpattern.simulate(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(measure_error_prob=measure_error_pr),
             branch_selector=ConstBranchSelector(outcome),
@@ -100,7 +100,7 @@ def test_noisy_measure_channel_hadamard(self, fx_rng: Generator) -> None:
         measure_channel_pr = fx_rng.random()
         print(f"measure_channel_pr = {measure_channel_pr}")
         # measurement error only
-        res = hadamardpattern.simulate_pattern(
+        res = hadamardpattern.simulate(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(measure_channel_prob=measure_channel_pr),
             rng=fx_rng,
@@ -119,7 +119,7 @@ def test_noisy_x_hadamard(self, fx_rng: Generator, outcome: Outcome) -> None:
         # x error only
         x_error_pr = fx_rng.random()
         print(f"x_error_pr = {x_error_pr}, outcome = {outcome}")
-        res = hadamardpattern.simulate_pattern(
+        res = hadamardpattern.simulate(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(x_error_prob=x_error_pr),
             branch_selector=ConstBranchSelector(outcome),
@@ -141,7 +141,7 @@ def test_noisy_x_hadamard(self, fx_rng: Generator, outcome: Outcome) -> None:
     def test_noisy_entanglement_hadamard(self, fx_rng: Generator) -> None:
         hadamardpattern = hpat()
         entanglement_error_pr = fx_rng.uniform()
-        res = hadamardpattern.simulate_pattern(
+        res = hadamardpattern.simulate(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(entanglement_error_prob=entanglement_error_pr),
             rng=fx_rng,
@@ -174,7 +174,7 @@ def test_noisy_preparation_hadamard(self, fx_rng: Generator) -> None:
         hadamardpattern = hpat()
         prepare_error_pr = fx_rng.random()
         print(f"prepare_error_pr = {prepare_error_pr}")
-        res = hadamardpattern.simulate_pattern(
+        res = hadamardpattern.simulate(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(prepare_error_prob=prepare_error_pr),
             rng=fx_rng,
@@ -193,9 +193,9 @@ def test_noisy_preparation_hadamard(self, fx_rng: Generator) -> None:
     def test_noiseless_noisy_rz(self, fx_rng: Generator) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
-        noiselessres = rzpattern.simulate_pattern(backend="densitymatrix", rng=fx_rng)
+        noiselessres = rzpattern.simulate(backend="densitymatrix", rng=fx_rng)
         # noiseless noise model or DepolarisingNoiseModel() since all probas are 0
-        noisynoiselessres = rzpattern.simulate_pattern(
+        noisynoiselessres = rzpattern.simulate(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(),
             rng=fx_rng,
@@ -212,7 +212,7 @@ def test_noisy_preparation_rz(self, fx_rng: Generator) -> None:
         rzpattern = rzpat(alpha)
         prepare_error_pr = fx_rng.random()
         print(f"prepare_error_pr = {prepare_error_pr}")
-        res = rzpattern.simulate_pattern(
+        res = rzpattern.simulate(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(prepare_error_prob=prepare_error_pr),
             rng=fx_rng,
@@ -246,7 +246,7 @@ def test_noisy_entanglement_rz(self, fx_rng: Generator) -> None:
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
         entanglement_error_pr = fx_rng.uniform()
-        res = rzpattern.simulate_pattern(
+        res = rzpattern.simulate(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(entanglement_error_prob=entanglement_error_pr),
             rng=fx_rng,
@@ -299,7 +299,7 @@ def test_noisy_measure_channel_rz(self, fx_rng: Generator) -> None:
         measure_channel_pr = fx_rng.random()
         print(f"measure_channel_pr = {measure_channel_pr}")
         # measurement error only
-        res = rzpattern.simulate_pattern(
+        res = rzpattern.simulate(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(measure_channel_prob=measure_channel_pr),
             rng=fx_rng,
@@ -341,7 +341,7 @@ def test_noisy_x_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Outc
         m_nodes = (cmd.node for cmd in rzpattern if cmd.kind == CommandKind.M)
         results: dict[int, Outcome] = {next(m_nodes): z_outcome, next(m_nodes): x_outcome}
 
-        res = rzpattern.simulate_pattern(
+        res = rzpattern.simulate(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(x_error_prob=x_error_pr),
             branch_selector=FixedBranchSelector(results),
@@ -379,7 +379,7 @@ def test_noisy_z_rz(self, fx_rng: Generator, outcome_z: Outcome, outcome_x: Outc
         # M(0) determines Z, M(1) determines X
         results: dict[int, Outcome] = {0: outcome_z, 1: outcome_x}
 
-        res = rzpattern.simulate_pattern(
+        res = rzpattern.simulate(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(z_error_prob=z_error_pr),
             branch_selector=FixedBranchSelector(results),
@@ -420,7 +420,7 @@ def test_noisy_xz_rz(self, fx_rng: Generator, z_outcome: Outcome, x_outcome: Out
         # M(0) determines Z correction, M(1) determines X correction
         results: dict[int, Outcome] = {0: z_outcome, 1: x_outcome}
 
-        res = rzpattern.simulate_pattern(
+        res = rzpattern.simulate(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(x_error_prob=x_error_pr, z_error_prob=z_error_pr),
             branch_selector=FixedBranchSelector(results),
@@ -481,7 +481,7 @@ def test_noisy_measure_confuse_rz(self, fx_rng: Generator, z_outcome: Outcome, x
         results: dict[int, Outcome] = {0: z_outcome, 1: x_outcome}
 
         # Test with probability 1 to flip both outcomes
-        res = rzpattern.simulate_pattern(
+        res = rzpattern.simulate(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(measure_error_prob=1.0),
             branch_selector=FixedBranchSelector(results),
@@ -507,7 +507,7 @@ def test_noisy_measure_confuse_rz_arbitrary(
         # Test with arbitrary probability
         measure_error_pr = fx_rng.random()
         print(f"measure_error_pr = {measure_error_pr}, z_outcome = {z_outcome}, x_outcome = {x_outcome}")
-        res = rzpattern.simulate_pattern(
+        res = rzpattern.simulate(
             backend="densitymatrix",
             noise_model=DepolarisingNoiseModel(measure_error_prob=measure_error_pr),
             branch_selector=FixedBranchSelector(results),
@@ -535,8 +535,8 @@ class TestNoisyDensityMatrixAmplitudeDamping:
     def test_noiseless_amplitude_damping_hadamard(self, fx_rng: Generator) -> None:
         """A zero-parameter amplitude damping model reproduces the noiseless result."""
         hadamardpattern = hpat()
-        noiselessres = hadamardpattern.simulate_pattern(backend="densitymatrix", rng=fx_rng)
-        noisynoiselessres = hadamardpattern.simulate_pattern(
+        noiselessres = hadamardpattern.simulate(backend="densitymatrix", rng=fx_rng)
+        noisynoiselessres = hadamardpattern.simulate(
             backend="densitymatrix",
             noise_model=AmplitudeDampingNoiseModel(),
             rng=fx_rng,
@@ -551,8 +551,8 @@ def test_noiseless_amplitude_damping_rz(self, fx_rng: Generator) -> None:
         """A zero-parameter amplitude damping model reproduces the noiseless RZ result."""
         alpha = fx_rng.random()
         rzpattern = rzpat(alpha)
-        noiselessres = rzpattern.simulate_pattern(backend="densitymatrix", rng=fx_rng)
-        noisynoiselessres = rzpattern.simulate_pattern(
+        noiselessres = rzpattern.simulate(backend="densitymatrix", rng=fx_rng)
+        noisynoiselessres = rzpattern.simulate(
             backend="densitymatrix",
             noise_model=AmplitudeDampingNoiseModel(),
             rng=fx_rng,
@@ -569,8 +569,8 @@ def test_compose_amplitude_damping_with_depolarising_runs(self, fx_rng: Generato
         ad = AmplitudeDampingNoiseModel(prepare_error_prob=0.2)
         composed = ComposeNoiseModel([depol, ad])
 
-        r_depol = hadamardpattern.simulate_pattern(backend="densitymatrix", noise_model=depol, rng=fx_rng)
-        r_both = hadamardpattern.simulate_pattern(backend="densitymatrix", noise_model=composed, rng=fx_rng)
+        r_depol = hadamardpattern.simulate(backend="densitymatrix", noise_model=depol, rng=fx_rng)
+        r_both = hadamardpattern.simulate(backend="densitymatrix", noise_model=composed, rng=fx_rng)
         assert isinstance(r_depol, DensityMatrix)
         assert isinstance(r_both, DensityMatrix)
         assert np.isclose(r_both.rho.trace(), 1.0)
@@ -673,7 +673,7 @@ def test_amplitude_damping_step_matches_analytic_hadamard(
         self, param: str, outcome: Outcome, fx_rng: Generator
     ) -> None:
         gamma = fx_rng.random()
-        res = hpat().simulate_pattern(
+        res = hpat().simulate(
             backend="densitymatrix",
             noise_model=AmplitudeDampingNoiseModel(**{param: gamma}),
             branch_selector=ConstBranchSelector(outcome),
@@ -752,7 +752,7 @@ def test_amplitude_damping_step_matches_analytic_rz(
         # transpiler puts Z corrections before X corrections by
         # default.
         results: dict[int, Outcome] = {0: outcome_z, 1: outcome_x}
-        res = rzpattern.simulate_pattern(
+        res = rzpattern.simulate(
             backend="densitymatrix",
             noise_model=AmplitudeDampingNoiseModel(**{param: gamma}),
             branch_selector=FixedBranchSelector(results),
diff --git a/tests/test_opengraph.py b/tests/test_opengraph.py
index 3840c74d6..476e2cae0 100644
--- a/tests/test_opengraph.py
+++ b/tests/test_opengraph.py
@@ -1055,10 +1055,10 @@ def check_determinism(pattern: Pattern, fx_rng: Generator, n_shots: int = 3) ->
     """Verify if the input pattern is deterministic."""
     for plane in {Plane.XY, Plane.XZ, Plane.YZ}:
         alpha = 2 * ANGLE_PI * fx_rng.random()
-        state_ref = pattern.simulate_pattern(input_state=PlanarState(plane, alpha), rng=fx_rng)
+        state_ref = pattern.simulate(input_state=PlanarState(plane, alpha), rng=fx_rng)
 
         for _ in range(n_shots):
-            state = pattern.simulate_pattern(input_state=PlanarState(plane, alpha), rng=fx_rng)
+            state = pattern.simulate(input_state=PlanarState(plane, alpha), rng=fx_rng)
             if not state.isclose(state_ref):
                 return False
 
@@ -1171,8 +1171,8 @@ def test_pflow_focused(self, test_case: OpenGraphFlowTestCase) -> None:
     def test_double_entanglement(self) -> None:
         pattern = Pattern(input_nodes=[0, 1], cmds=[E((0, 1)), E((0, 1))])
         pattern2 = pattern.extract_opengraph().to_pattern()
-        state = pattern.simulate_pattern()
-        state2 = pattern2.simulate_pattern()
+        state = pattern.simulate()
+        state2 = pattern2.simulate()
         assert state.isclose(state2)
 
     def test_from_to_pattern(self, fx_rng: Generator) -> None:
@@ -1186,8 +1186,8 @@ def test_from_to_pattern(self, fx_rng: Generator) -> None:
 
         for plane in {Plane.XY, Plane.XZ, Plane.YZ}:
             alpha = 2 * ANGLE_PI * fx_rng.random()
-            state_ref = pattern_ref.simulate_pattern(input_state=PlanarState(plane, alpha), rng=fx_rng)
-            state = pattern.simulate_pattern(input_state=PlanarState(plane, alpha), rng=fx_rng)
+            state_ref = pattern_ref.simulate(input_state=PlanarState(plane, alpha), rng=fx_rng)
+            state = pattern.simulate(input_state=PlanarState(plane, alpha), rng=fx_rng)
             assert state.isclose(state_ref)
 
     def test_isclose_measurement(self) -> None:
@@ -1337,8 +1337,8 @@ def test_compose_clifford(self, fx_rng: Generator) -> None:
         og_c, _ = og1.compose(og2, mapping)
         pc_test = og_c.to_pattern()
 
-        sv = pc.simulate_pattern(rng=fx_rng)
-        sv_test = pc_test.simulate_pattern(rng=fx_rng)
+        sv = pc.simulate(rng=fx_rng)
+        sv_test = pc_test.simulate(rng=fx_rng)
 
         assert sv.isclose(sv_test)
 
diff --git a/tests/test_optimization.py b/tests/test_optimization.py
index 5c386d8d3..39093e73e 100644
--- a/tests/test_optimization.py
+++ b/tests/test_optimization.py
@@ -52,8 +52,8 @@ def test_standardize_clifford_entanglement(fx_rng: Generator) -> None:
                 assert p[2].kind == CommandKind.C
                 assert p[3].kind == CommandKind.C
 
-                state_ref = p_ref.simulate_pattern(input_state=PlanarState(Plane.XY, alpha))
-                state_p = p.simulate_pattern(input_state=PlanarState(Plane.XY, alpha))
+                state_ref = p_ref.simulate(input_state=PlanarState(Plane.XY, alpha))
+                state_p = p.simulate(input_state=PlanarState(Plane.XY, alpha))
                 assert state_p.isclose(state_ref)
 
 
@@ -87,8 +87,8 @@ def test_remove_useless_domains(fx_bg: PCG64, jumps: int) -> None:
     pattern.remove_pauli_measurements()
     pattern2 = remove_useless_domains(pattern)
     pattern2 = StandardizedPattern.from_pattern(pattern2).to_space_optimal_pattern()
-    state = pattern.simulate_pattern(rng=rng)
-    state2 = pattern2.simulate_pattern(rng=rng)
+    state = pattern.simulate(rng=rng)
+    state2 = pattern2.simulate(rng=rng)
     assert state.isclose(state2)
 
 
@@ -109,8 +109,8 @@ def test_to_space_optimal_pattern(fx_rng: Generator) -> None:
         output_nodes=[17, 18],
     )
     pattern2 = StandardizedPattern.from_pattern(pattern).to_space_optimal_pattern()
-    state = pattern.simulate_pattern(rng=fx_rng)
-    state2 = pattern2.simulate_pattern(rng=fx_rng)
+    state = pattern.simulate(rng=fx_rng)
+    state2 = pattern2.simulate(rng=fx_rng)
     assert state.isclose(state2)
 
 
@@ -149,8 +149,8 @@ def test_remove_local_clifford_commands(fx_bg: PCG64, jumps: int) -> None:
     new_pattern = pattern.remove_local_clifford_commands(copy=True)
     assert not any(cmd.kind == CommandKind.C for cmd in new_pattern)
     input_state = rand_state_vector(nqubits, rng=rng)
-    state_ref = pattern.simulate_pattern(input_state=input_state, rng=rng)
-    state = new_pattern.simulate_pattern(input_state=input_state, rng=rng)
+    state_ref = pattern.simulate(input_state=input_state, rng=rng)
+    state = new_pattern.simulate(input_state=input_state, rng=rng)
     assert state.isclose(state_ref)
 
 
diff --git a/tests/test_parameter.py b/tests/test_parameter.py
index 85d0e7ecf..910bd869c 100644
--- a/tests/test_parameter.py
+++ b/tests/test_parameter.py
@@ -68,7 +68,7 @@ def test_instantiated_pattern_simulation(fx_rng: Generator) -> None:
     pattern.add(graphix.command.M(1, Measurement.XY(alpha)))
     pattern0 = pattern.subs(alpha, 0)
     # Instantied patterns can be simulated.
-    pattern0.simulate_pattern(rng=fx_rng)
+    pattern0.simulate(rng=fx_rng)
     pattern1 = pattern.subs(alpha, 1)
     assert not pattern1.is_parameterized()
     assert list(pattern1) == [graphix.command.M(node=0), graphix.command.M(1, Measurement.XY(1))]
@@ -76,7 +76,7 @@ def test_instantiated_pattern_simulation(fx_rng: Generator) -> None:
         graphix.command.M(node=0),
         graphix.command.M(1, -Measurement.X),
     ]
-    pattern1.simulate_pattern(rng=fx_rng)
+    pattern1.simulate(rng=fx_rng)
 
 
 def test_multiple_parameters(fx_rng: Generator) -> None:
@@ -104,7 +104,7 @@ def test_multiple_parameters(fx_rng: Generator) -> None:
         graphix.command.N(node=2),
         graphix.command.M(2, -Measurement.X),
     ]
-    pattern23.simulate_pattern(rng=fx_rng)
+    pattern23.simulate(rng=fx_rng)
 
 
 def test_parallel_substitution() -> None:
@@ -175,11 +175,11 @@ def test_random_circuit_with_parameters(fx_bg: PCG64, jumps: int, use_xreplace:
     pattern.minimize_space()
     assignment: dict[Parameter, float] = {alpha: rng.uniform(high=2), beta: rng.uniform(high=2)}
     if use_xreplace:
-        state = circuit.xreplace(assignment).simulate_statevector().statevec
-        state_mbqc = pattern.xreplace(assignment).simulate_pattern(rng=rng)
+        state = circuit.xreplace(assignment).simulate().statevec
+        state_mbqc = pattern.xreplace(assignment).simulate(rng=rng)
     else:
-        state = circuit.subs(alpha, assignment[alpha]).subs(beta, assignment[beta]).simulate_statevector().statevec
-        state_mbqc = pattern.subs(alpha, assignment[alpha]).subs(beta, assignment[beta]).simulate_pattern(rng=rng)
+        state = circuit.subs(alpha, assignment[alpha]).subs(beta, assignment[beta]).simulate().statevec
+        state_mbqc = pattern.subs(alpha, assignment[alpha]).subs(beta, assignment[beta]).simulate(rng=rng)
     assert state_mbqc.isclose(state)
 
 
@@ -199,4 +199,4 @@ def test_simulation_exception(fx_rng: Generator) -> None:
     alpha = Placeholder("alpha")
     pattern.add(graphix.command.M(1, Measurement.XY(alpha)))
     with pytest.raises(PlaceholderOperationError):
-        pattern.simulate_pattern(rng=fx_rng)
+        pattern.simulate(rng=fx_rng)
diff --git a/tests/test_pattern.py b/tests/test_pattern.py
index 6e16278d1..e230be685 100644
--- a/tests/test_pattern.py
+++ b/tests/test_pattern.py
@@ -74,9 +74,9 @@ def test_standardize(self, fx_rng: Generator) -> None:
         pattern = pattern.infer_pauli_measurements()
         pattern.standardize()
         assert pattern.is_standard()
-        state = circuit.simulate_statevector().statevec
+        state = circuit.simulate().statevec
         pattern.remove_pauli_measurements()
-        state_mbqc = pattern.simulate_pattern(rng=fx_rng)
+        state_mbqc = pattern.simulate(rng=fx_rng)
         assert state_mbqc.isclose(state)
 
     def test_minimize_space(self, fx_rng: Generator) -> None:
@@ -86,8 +86,8 @@ def test_minimize_space(self, fx_rng: Generator) -> None:
         pattern = circuit.transpile().pattern
         pattern.standardize()
         pattern.minimize_space()
-        state = circuit.simulate_statevector().statevec
-        state_mbqc = pattern.simulate_pattern(rng=fx_rng)
+        state = circuit.simulate().statevec
+        state_mbqc = pattern.simulate(rng=fx_rng)
         assert state_mbqc.isclose(state)
 
     # https://github.com/TeamGraphix/graphix/issues/157
@@ -158,8 +158,8 @@ def test_minimize_space_with_gflow(self, fx_bg: PCG64, jumps: int) -> None:
         pattern = pattern.infer_pauli_measurements()
         pattern.remove_pauli_measurements()
         pattern.minimize_space()
-        state = circuit.simulate_statevector().statevec
-        state_mbqc = pattern.simulate_pattern(rng=rng)
+        state = circuit.simulate().statevec
+        state_mbqc = pattern.simulate(rng=rng)
         assert state_mbqc.isclose(state)
 
     @pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.")
@@ -202,8 +202,8 @@ def test_parallelize_pattern(self, fx_rng: Generator) -> None:
         pattern = circuit.transpile().pattern
         pattern.standardize()
         pattern.parallelize_pattern()
-        state = circuit.simulate_statevector().statevec
-        state_mbqc = pattern.simulate_pattern(rng=fx_rng)
+        state = circuit.simulate().statevec
+        state_mbqc = pattern.simulate(rng=fx_rng)
         assert state_mbqc.isclose(state)
 
     @pytest.mark.parametrize("jumps", range(1, 11))
@@ -217,8 +217,8 @@ def test_shift_signals(self, fx_bg: PCG64, jumps: int) -> None:
         pattern.shift_signals(method="mc")
         assert pattern.is_standard()
         pattern = StandardizedPattern.from_pattern(pattern).to_space_optimal_pattern()
-        state = circuit.simulate_statevector().statevec
-        state_mbqc = pattern.simulate_pattern(rng=rng)
+        state = circuit.simulate().statevec
+        state_mbqc = pattern.simulate(rng=rng)
         assert state_mbqc.isclose(state)
 
     @pytest.mark.parametrize("jumps", range(1, 11))
@@ -238,8 +238,8 @@ def test_pauli_measurement_random_circuit(
         pattern = pattern.infer_pauli_measurements()
         pattern.remove_pauli_measurements()
         pattern.minimize_space()
-        state = circuit.simulate_statevector().statevec
-        state_mbqc: Statevec | DensityMatrix = pattern.simulate_pattern(backend, rng=rng)
+        state = circuit.simulate().statevec
+        state_mbqc: Statevec | DensityMatrix = pattern.simulate(backend, rng=rng)
         assert compare_backend_result_with_statevec(state_mbqc, state) == pytest.approx(1)
 
     @pytest.mark.parametrize("jumps", range(1, 11))
@@ -268,8 +268,8 @@ def test_pauli_measurement_single(self, pm: PauliMeasurement) -> None:
         pattern_ref = pattern.copy()
         pattern.remove_pauli_measurements()
         branch_selector = ConstBranchSelector(0)
-        state = pattern.simulate_pattern(branch_selector=branch_selector)
-        state_ref = pattern_ref.simulate_pattern(branch_selector=branch_selector)
+        state = pattern.simulate(branch_selector=branch_selector)
+        state_ref = pattern_ref.simulate(branch_selector=branch_selector)
         assert state.isclose(state_ref)
 
     def test_pauli_measurement(self) -> None:
@@ -296,8 +296,8 @@ def test_pauli_measurement(self) -> None:
         assert isolated_nodes == set()
         pattern.minimize_space()
         pattern_opt.minimize_space()
-        state = pattern.simulate_pattern()
-        state_opt = pattern.simulate_pattern()
+        state = pattern.simulate()
+        state_opt = pattern.simulate()
         assert state.isclose(state_opt)
 
     @pytest.mark.parametrize("jumps", range(1, 6))
@@ -311,8 +311,8 @@ def test_pauli_measured_against_nonmeasured(self, fx_bg: PCG64, jumps: int) -> N
         pattern1 = copy.deepcopy(pattern)
         pattern1 = pattern1.infer_pauli_measurements()
         pattern1.remove_pauli_measurements()
-        state = pattern.simulate_pattern(rng=rng)
-        state1 = pattern1.simulate_pattern(rng=rng)
+        state = pattern.simulate(rng=rng)
+        state1 = pattern1.simulate(rng=rng)
         assert state.isclose(state1)
 
     def test_extract_measurement_commands(self) -> None:
@@ -374,10 +374,10 @@ def test_shift_signals_plane(self, plane: Plane, method: str) -> None:
         for outcomes_ref_list in itertools.product(*([zero_one] * 3)):
             outcomes_ref = dict(enumerate(outcomes_ref_list))
             branch_selector = FixedBranchSelector(results=outcomes_ref)
-            state_ref = pattern_ref.simulate_pattern(branch_selector=branch_selector)
+            state_ref = pattern_ref.simulate(branch_selector=branch_selector)
             outcomes_p = shift_outcomes(outcomes_ref, signal_dict)
             branch_selector = FixedBranchSelector(results=outcomes_p)
-            state_p = pattern.simulate_pattern(branch_selector=branch_selector)
+            state_p = pattern.simulate(branch_selector=branch_selector)
             assert state_p.isclose(state_ref)
 
     @pytest.mark.parametrize("jumps", range(1, 11))
@@ -390,8 +390,8 @@ def test_standardize_direct(self, fx_bg: PCG64, jumps: int) -> None:
         pattern.standardize()
         assert pattern.is_standard()
         pattern.minimize_space()
-        state_p = pattern.simulate_pattern(rng=rng)
-        state_ref = circuit.simulate_statevector().statevec
+        state_p = pattern.simulate(rng=rng)
+        state_ref = circuit.simulate().statevec
         assert state_p.isclose(state_ref)
 
     @pytest.mark.parametrize("jumps", range(1, 11))
@@ -404,8 +404,8 @@ def test_shift_signals_direct(self, fx_bg: PCG64, jumps: int) -> None:
         pattern.standardize()
         pattern.shift_signals(method="direct")
         pattern.minimize_space()
-        state_p = pattern.simulate_pattern(rng=rng)
-        state_ref = circuit.simulate_statevector().statevec
+        state_p = pattern.simulate(rng=rng)
+        state_ref = circuit.simulate().statevec
         assert state_p.isclose(state_ref)
 
     @pytest.mark.parametrize("jumps", range(1, 11))
@@ -419,8 +419,8 @@ def test_pauli_measurement_then_standardize(self, fx_bg: PCG64, jumps: int) -> N
         pattern.remove_pauli_measurements()
         pattern.standardize()
         pattern.minimize_space()
-        state = circuit.simulate_statevector().statevec
-        state_mbqc = pattern.simulate_pattern(rng=rng)
+        state = circuit.simulate().statevec
+        state_mbqc = pattern.simulate(rng=rng)
         assert compare_backend_result_with_statevec(state_mbqc, state) == pytest.approx(1)
 
     @pytest.mark.parametrize("jumps", range(1, 11))
@@ -432,8 +432,8 @@ def test_standardize_two_cliffords(self, fx_bg: PCG64, jumps: int) -> None:
         pattern.add(C(node=0, clifford=Clifford(c1)))
         pattern_ref = pattern.copy()
         pattern.standardize()
-        state_ref = pattern_ref.simulate_pattern()
-        state_p = pattern.simulate_pattern()
+        state_ref = pattern_ref.simulate()
+        state_p = pattern.simulate()
         assert state_p.isclose(state_ref)
 
     # Simple pattern composition
@@ -670,8 +670,8 @@ def test_compose_6(self, fx_bg: PCG64, jumps: int) -> None:
         )
         p.minimize_space()
         p_compose.minimize_space()
-        s = p.simulate_pattern(rng=rng)
-        s_compose = p_compose.simulate_pattern(rng=rng)
+        s = p.simulate(rng=rng)
+        s_compose = p_compose.simulate(rng=rng)
         assert s.isclose(s_compose)
 
     # Test warning composition after standardization
@@ -760,7 +760,7 @@ def test_check_runnability_failures(self) -> None:
 
         pattern = Pattern(cmds=[N(0), M(0, s_domain={0})])
         with pytest.raises(RunnabilityError) as exc_info:
-            pattern.simulate_pattern()
+            pattern.simulate()
         assert exc_info.value.node == 0
         assert exc_info.value.reason == RunnabilityErrorReason.DomainSelfLoop
 
@@ -943,8 +943,8 @@ def test_extract_causal_flow_rnd_circuit(self, fx_bg: PCG64, jumps: int) -> None
         p_ref.remove_pauli_measurements()
         p_test.remove_pauli_measurements()
 
-        s_ref = p_ref.simulate_pattern(rng=rng)
-        s_test = p_test.simulate_pattern(rng=rng)
+        s_ref = p_ref.simulate(rng=rng)
+        s_test = p_test.simulate(rng=rng)
         assert s_ref.isclose(s_test)
 
     # Extract gflow from random circuits
@@ -963,8 +963,8 @@ def test_extract_gflow_rnd_circuit(self, fx_bg: PCG64, jumps: int) -> None:
         p_ref.remove_pauli_measurements()
         p_test.remove_pauli_measurements()
 
-        s_ref = p_ref.simulate_pattern(rng=rng)
-        s_test = p_test.simulate_pattern(rng=rng)
+        s_ref = p_ref.simulate(rng=rng)
+        s_test = p_test.simulate(rng=rng)
         assert s_ref.isclose(s_test)
 
     # Extract Pauli flow from random circuits
@@ -981,18 +981,18 @@ def test_extract_pauli_flow_rnd_circuit(self, fx_bg: PCG64, jumps: int) -> None:
         p_ref.remove_pauli_measurements()
         p_test.remove_pauli_measurements()
 
-        s_ref = p_ref.simulate_pattern(rng=rng)
-        s_test = p_test.simulate_pattern(rng=rng)
+        s_ref = p_ref.simulate(rng=rng)
+        s_test = p_test.simulate(rng=rng)
         assert s_ref.isclose(s_test)
 
     @pytest.mark.parametrize("test_case", PATTERN_FLOW_TEST_CASES)
     def test_extract_causal_flow(self, fx_rng: Generator, test_case: PatternFlowTestCase) -> None:
         if test_case.has_cflow:
             alpha = 2 * np.pi * fx_rng.random()
-            s_ref = test_case.pattern.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
+            s_ref = test_case.pattern.simulate(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
 
             p_test = test_case.pattern.to_bloch().extract_causal_flow().to_corrections().to_pattern()
-            s_test = p_test.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
+            s_test = p_test.simulate(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
 
             assert s_ref.isclose(s_test)
         else:
@@ -1004,10 +1004,10 @@ def test_extract_gflow(self, fx_rng: Generator, test_case: PatternFlowTestCase)
         """Tests the round trip Pattern -> XZCorrections -> GFlow -> XZCorrections -> Pattern."""
         if test_case.has_gflow:
             alpha = 2 * np.pi * fx_rng.random()
-            s_ref = test_case.pattern.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
+            s_ref = test_case.pattern.simulate(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
 
             p_test = test_case.pattern.to_bloch().extract_gflow().to_corrections().to_pattern()
-            s_test = p_test.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
+            s_test = p_test.simulate(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
 
             assert s_ref.isclose(s_test)
         else:
@@ -1023,10 +1023,10 @@ def test_extract_pauli_flow(self, fx_rng: Generator, test_case: PatternFlowTestC
         if not test_case.has_gflow:
             pytest.skip("no gflow; Pauli-flow existence is not determined by has_gflow")
         alpha = 2 * np.pi * fx_rng.random()
-        s_ref = test_case.pattern.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
+        s_ref = test_case.pattern.simulate(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
 
         p_test = test_case.pattern.to_bloch().extract_pauli_flow().to_corrections().to_pattern()
-        s_test = p_test.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
+        s_test = p_test.simulate(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
 
         assert s_ref.isclose(s_test)
 
@@ -1046,10 +1046,10 @@ def test_extract_cflow_og(self, fx_rng: Generator) -> None:
             },
         )
         p_ref = og.extract_causal_flow().to_corrections().to_pattern()
-        s_ref = p_ref.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
+        s_ref = p_ref.simulate(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
 
         p_test = p_ref.extract_causal_flow().to_corrections().to_pattern()
-        s_test = p_test.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
+        s_test = p_test.simulate(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
 
         assert s_ref.isclose(s_test)
 
@@ -1070,10 +1070,10 @@ def test_extract_gflow_og(self, fx_rng: Generator) -> None:
         )
 
         p_ref = og.extract_gflow().to_corrections().to_pattern()
-        s_ref = p_ref.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
+        s_ref = p_ref.simulate(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
 
         p_test = p_ref.extract_gflow().to_corrections().to_pattern()
-        s_test = p_test.simulate_pattern(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
+        s_test = p_test.simulate(input_state=PlanarState(Plane.XZ, alpha), rng=fx_rng)
 
         assert s_ref.isclose(s_test)
 
@@ -1094,8 +1094,8 @@ def test_extract_xzc_rnd_circuit(self, fx_bg: PCG64, jumps: int) -> None:
         p_test = p_test.infer_pauli_measurements()
         p_test.remove_pauli_measurements()
 
-        s_ref = p_ref.simulate_pattern(rng=rng)
-        s_test = p_test.simulate_pattern(rng=rng)
+        s_ref = p_ref.simulate(rng=rng)
+        s_test = p_test.simulate(rng=rng)
         assert s_ref.isclose(s_test)
 
     def test_extract_xzc_empty_domains(self) -> None:
@@ -1202,8 +1202,8 @@ def test_extract_opengraph_standardization(self) -> None:
     def test_extract_opengraph_roundtrip(self, pattern: Pattern, fx_rng: Generator) -> None:
         pattern_test = pattern.extract_opengraph().to_pattern()
 
-        sv = pattern.simulate_pattern(rng=fx_rng)
-        sv_test = pattern_test.simulate_pattern(rng=fx_rng)
+        sv = pattern.simulate(rng=fx_rng)
+        sv_test = pattern_test.simulate(rng=fx_rng)
 
         assert sv.isclose(sv_test)
 
@@ -1310,8 +1310,8 @@ def test_standardize(self, fx_bg: PCG64, jumps: int) -> None:
         assert pattern.is_standard()
         pattern.minimize_space()
         pattern_mc.minimize_space()
-        state_d = pattern.simulate_pattern(rng=rng)
-        state_ref = pattern_mc.simulate_pattern(rng=rng)
+        state_d = pattern.simulate(rng=rng)
+        state_ref = pattern_mc.simulate(rng=rng)
         assert state_d.isclose(state_ref)
 
     @pytest.mark.parametrize("jumps", range(1, 11))
@@ -1329,8 +1329,8 @@ def test_shift_signals(self, fx_bg: PCG64, jumps: int) -> None:
         assert pattern.is_standard()
         pattern.minimize_space()
         pattern_mc.minimize_space()
-        state_d = pattern.simulate_pattern(rng=rng)
-        state_ref = pattern_mc.simulate_pattern(rng=rng)
+        state_d = pattern.simulate(rng=rng)
+        state_ref = pattern_mc.simulate(rng=rng)
         assert state_d.isclose(state_ref)
 
     @pytest.mark.parametrize("jumps", range(1, 11))
@@ -1344,8 +1344,8 @@ def test_standardize_and_shift_signals(self, fx_bg: PCG64, jumps: int) -> None:
         pattern.shift_signals()
         assert pattern.is_standard()
         pattern.minimize_space()
-        state_p = pattern.simulate_pattern(rng=rng)
-        state_ref = circuit.simulate_statevector().statevec
+        state_p = pattern.simulate(rng=rng)
+        state_ref = circuit.simulate().statevec
         assert state_p.isclose(state_ref)
 
     @pytest.mark.parametrize("jumps", range(1, 4))
@@ -1364,7 +1364,7 @@ def test_mixed_pattern_operations(self, fx_bg: PCG64, jumps: int) -> None:
         nqubits = 3
         depth = 2
         circuit = rand_circuit(nqubits, depth, rng)
-        state_ref = circuit.simulate_statevector().statevec
+        state_ref = circuit.simulate().statevec
         for process in processes:
             pattern = circuit.transpile().pattern
             for operation in process:
@@ -1374,7 +1374,7 @@ def test_mixed_pattern_operations(self, fx_bg: PCG64, jumps: int) -> None:
                     pattern.shift_signals(method=operation[1])
             assert pattern.is_standard()
             pattern.minimize_space()
-            state_p = pattern.simulate_pattern(rng=rng)
+            state_p = pattern.simulate(rng=rng)
             assert state_p.isclose(state_ref)
 
     def test_pauli_measurement_end_with_measure(self) -> None:
@@ -1393,8 +1393,8 @@ def test_arbitrary_inputs(
         rand_planes = fx_rng.choice(np.array(Plane), nqb)
         states = [PlanarState(plane=i, angle=j) for i, j in zip(rand_planes, rand_angles, strict=True)]
         randpattern = rand_circ.transpile().pattern
-        out: Statevec | DensityMatrix = randpattern.simulate_pattern(backend=backend, input_state=states, rng=fx_rng)
-        out_circ = rand_circ.simulate_statevector(input_state=states).statevec
+        out: Statevec | DensityMatrix = randpattern.simulate(backend=backend, input_state=states, rng=fx_rng)
+        out_circ = rand_circ.simulate(input_state=states).statevec
         assert compare_backend_result_with_statevec(out, out_circ) == pytest.approx(1)
 
     def test_arbitrary_inputs_tn(self, fx_rng: Generator, nqb: int, rand_circ: Circuit) -> None:
@@ -1403,9 +1403,7 @@ def test_arbitrary_inputs_tn(self, fx_rng: Generator, nqb: int, rand_circ: Circu
         states = [PlanarState(plane=i, angle=j) for i, j in zip(rand_planes, rand_angles, strict=True)]
         randpattern = rand_circ.transpile().pattern
         with pytest.raises(NotImplementedError):
-            randpattern.simulate_pattern(
-                backend="tensornetwork", graph_prep="sequential", input_state=states, rng=fx_rng
-            )
+            randpattern.simulate(backend="tensornetwork", graph_prep="sequential", input_state=states, rng=fx_rng)
 
 
 def assert_equal_edge(edge: Sequence[int], ref: Sequence[int]) -> bool:
diff --git a/tests/test_qasm3_exporter_to_qiskit.py b/tests/test_qasm3_exporter_to_qiskit.py
index 77c5bab76..31fdfb5f5 100644
--- a/tests/test_qasm3_exporter_to_qiskit.py
+++ b/tests/test_qasm3_exporter_to_qiskit.py
@@ -76,7 +76,7 @@ def check_qasm3(pattern: Pattern) -> None:
     # Reorder qubits to match the pattern's expected output ordering.
     backend.finalize(pattern.output_nodes)
     state_qiskit = backend.state
-    state_mbqc = pattern.simulate_pattern(branch_selector=branch_selector)
+    state_mbqc = pattern.simulate(branch_selector=branch_selector)
     assert state_mbqc.isclose(state_qiskit)
 
 
diff --git a/tests/test_remove_pauli_measurements.py b/tests/test_remove_pauli_measurements.py
index 10870ab36..aad694748 100644
--- a/tests/test_remove_pauli_measurements.py
+++ b/tests/test_remove_pauli_measurements.py
@@ -199,8 +199,8 @@ def check_pattern_equivalence(pattern: Pattern, pattern2: Pattern, rng: Generato
     pattern2.minimize_space()
     for _ in range(4):
         input_state = rand_state_vector(len(pattern.input_nodes), rng=rng)
-        state = pattern.simulate_pattern(input_state=input_state, rng=rng)
-        state2 = pattern2.simulate_pattern(input_state=input_state, rng=rng)
+        state = pattern.simulate(input_state=input_state, rng=rng)
+        state2 = pattern2.simulate(input_state=input_state, rng=rng)
         assert state.isclose(state2)
 
 
@@ -288,7 +288,7 @@ def test_pattern_remove_pauli_measurements_output_nodes() -> None:
     )
     pattern = og.to_pattern()
     pattern.remove_pauli_measurements()
-    pattern.simulate_pattern()
+    pattern.simulate()
 
 
 def test_try_pivot_x_with_output_node_after_pivot() -> None:
diff --git a/tests/test_simulator.py b/tests/test_simulator.py
index 6c9ecbb81..e3979d186 100644
--- a/tests/test_simulator.py
+++ b/tests/test_simulator.py
@@ -13,14 +13,14 @@
 def test_no_explicit_input_state(hadamardpattern: Pattern, fx_rng: Generator) -> None:
     # No explicit input state: the default initial state is |+⟩.
     # H|+⟩ = |0⟩, so we expect the final state to be |0⟩.
-    state = hadamardpattern.simulate_pattern(rng=fx_rng)
+    state = hadamardpattern.simulate(rng=fx_rng)
     assert state.isclose(Statevec(BasicStates.ZERO))
 
 
 def test_explicit_input_state_zero(hadamardpattern: Pattern, fx_rng: Generator) -> None:
     # Provide an explicit input state |0⟩.
     # H|0⟩ = |+⟩, so the final state should be |+⟩.
-    state = hadamardpattern.simulate_pattern(input_state=BasicStates.ZERO, rng=fx_rng)
+    state = hadamardpattern.simulate(input_state=BasicStates.ZERO, rng=fx_rng)
     assert state.isclose(Statevec(BasicStates.PLUS))
 
 
@@ -30,7 +30,7 @@ def test_backend_prepared_zero(hadamardpattern: Pattern, fx_rng: Generator) -> N
     # therefore H|0⟩ = |+⟩.
     backend = StatevectorBackend()
     backend.add_nodes(hadamardpattern.input_nodes, BasicStates.ZERO)
-    state = hadamardpattern.simulate_pattern(backend=backend, input_state=None, rng=fx_rng)
+    state = hadamardpattern.simulate(backend=backend, input_state=None, rng=fx_rng)
     assert state.isclose(Statevec(BasicStates.PLUS))
 
 
@@ -39,7 +39,7 @@ def test_no_prepared_qubits_and_input_state_none(hadamardpattern: Pattern, fx_rn
     # This is ambiguous, so a ValueError must be raised.
     backend = StatevectorBackend()
     with pytest.raises(ValueError, match="the backend is expected to have 1 input nodes already prepared"):
-        hadamardpattern.simulate_pattern(backend=backend, input_state=None, rng=fx_rng)
+        hadamardpattern.simulate(backend=backend, input_state=None, rng=fx_rng)
 
 
 def test_prepared_qubits_and_input_state(hadamardpattern: Pattern, fx_rng: Generator) -> None:
@@ -50,11 +50,11 @@ def test_prepared_qubits_and_input_state(hadamardpattern: Pattern, fx_rng: Gener
     backend = StatevectorBackend()
     backend.add_nodes(hadamardpattern.input_nodes, BasicStates.ZERO)
     with pytest.raises(ValueError, match="the backend is expected to have no pre-allocated qubits"):
-        hadamardpattern.simulate_pattern(backend=backend, rng=fx_rng)
+        hadamardpattern.simulate(backend=backend, rng=fx_rng)
 
 
 def test_node_index_after_finalize() -> None:
     pattern = Pattern(input_nodes=[0, 1], output_nodes=[1, 0])
     backend = StatevectorBackend()
-    pattern.simulate_pattern(backend=backend)
+    pattern.simulate(backend=backend)
     assert list(backend.node_index) == [1, 0]
diff --git a/tests/test_statevec.py b/tests/test_statevec.py
index a91e56f7c..b1c657ee0 100644
--- a/tests/test_statevec.py
+++ b/tests/test_statevec.py
@@ -163,7 +163,7 @@ def test_nqubits(self) -> None:
 
     def test_nqubits_pattern(self) -> None:
         p = Pattern(input_nodes=[0, 1, 2])
-        sv = p.simulate_pattern(backend="statevector")
+        sv = p.simulate(backend="statevector")
         assert sv.nqubit == 3
 
 
diff --git a/tests/test_tnsim.py b/tests/test_tnsim.py
index 2bb76a21b..8c3cab71d 100644
--- a/tests/test_tnsim.py
+++ b/tests/test_tnsim.py
@@ -53,7 +53,7 @@ def test_entangle_nodes(self, fx_rng: Generator) -> None:
         circuit = Circuit(2)
         pattern = circuit.transpile().pattern
         pattern.add(E(nodes=(0, 1)))
-        tn = pattern.simulate_pattern(backend="tensornetwork", graph_prep="sequential", rng=fx_rng)
+        tn = pattern.simulate(backend="tensornetwork", graph_prep="sequential", rng=fx_rng)
         dummy_index = [gen_str() for _ in range(2)]
         for qubit_index, n in enumerate(tn._dangling):
             ind = tn._dangling[n]
@@ -82,7 +82,7 @@ def test_entangle_nodes(self, fx_rng: Generator) -> None:
     #     pattern.results[15] = 1  # X&Z operator will be applied.
     #     for cmd in cmds:
     #         pattern.add(cmd)
-    #     tn = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+    #     tn = pattern.simulate(backend="tensornetwork", rng=fx_rng)
     #     dummy_index = gen_str()
     #     ind = tn._dangling.pop("0")
     #     tensor = tn.tensor_map[tn._get_tids_from_inds(ind).popleft()]
@@ -102,9 +102,9 @@ def test_entangle_nodes(self, fx_rng: Generator) -> None:
 
     def test_expectation_value1(self, fx_rng: Generator) -> None:
         circuit = Circuit(1)
-        state = circuit.simulate_statevector().statevec
+        state = circuit.simulate().statevec
         pattern = circuit.transpile().pattern
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op1 = random_op(1, fx_rng)
         value1 = state.expectation_value(random_op1, [0])
         value2 = tn_mbqc.expectation_value(random_op1, [0])
@@ -112,9 +112,9 @@ def test_expectation_value1(self, fx_rng: Generator) -> None:
 
     def test_expectation_value2(self, fx_rng: Generator) -> None:
         circuit = Circuit(2)
-        state = circuit.simulate_statevector().statevec
+        state = circuit.simulate().statevec
         pattern = circuit.transpile().pattern
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op2 = random_op(2, fx_rng)
         input_list = [0, 1]
         for qargs in itertools.permutations(input_list):
@@ -124,9 +124,9 @@ def test_expectation_value2(self, fx_rng: Generator) -> None:
 
     def test_expectation_value3(self, fx_rng: Generator) -> None:
         circuit = Circuit(3)
-        state = circuit.simulate_statevector().statevec
+        state = circuit.simulate().statevec
         pattern = circuit.transpile().pattern
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op3 = random_op(3, fx_rng)
         input_list = [0, 1, 2]
         for qargs in itertools.permutations(input_list):
@@ -136,9 +136,9 @@ def test_expectation_value3(self, fx_rng: Generator) -> None:
 
     def test_expectation_value3_sequential(self, fx_rng: Generator) -> None:
         circuit = Circuit(3)
-        state = circuit.simulate_statevector().statevec
+        state = circuit.simulate().statevec
         pattern = circuit.transpile().pattern
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", graph_prep="sequential", rng=fx_rng)
+        tn_mbqc = pattern.simulate(backend="tensornetwork", graph_prep="sequential", rng=fx_rng)
         random_op3 = random_op(3, fx_rng)
         input_list = [0, 1, 2]
         for qargs in itertools.permutations(input_list):
@@ -148,9 +148,9 @@ def test_expectation_value3_sequential(self, fx_rng: Generator) -> None:
 
     def test_expectation_value3_subspace1(self, fx_rng: Generator) -> None:
         circuit = Circuit(3)
-        state = circuit.simulate_statevector().statevec
+        state = circuit.simulate().statevec
         pattern = circuit.transpile().pattern
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op1 = random_op(1, fx_rng)
         input_list = [0, 1, 2]
         for qargs in itertools.permutations(input_list, 1):
@@ -160,9 +160,9 @@ def test_expectation_value3_subspace1(self, fx_rng: Generator) -> None:
 
     def test_expectation_value3_subspace2(self, fx_rng: Generator) -> None:
         circuit = Circuit(3)
-        state = circuit.simulate_statevector().statevec
+        state = circuit.simulate().statevec
         pattern = circuit.transpile().pattern
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op2 = random_op(2, fx_rng)
         input_list = [0, 1, 2]
         for qargs in itertools.permutations(input_list, 2):
@@ -172,9 +172,9 @@ def test_expectation_value3_subspace2(self, fx_rng: Generator) -> None:
 
     def test_expectation_value3_subspace2_sequential(self, fx_rng: Generator) -> None:
         circuit = Circuit(3)
-        state = circuit.simulate_statevector().statevec
+        state = circuit.simulate().statevec
         pattern = circuit.transpile().pattern
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", graph_prep="sequential", rng=fx_rng)
+        tn_mbqc = pattern.simulate(backend="tensornetwork", graph_prep="sequential", rng=fx_rng)
         random_op2 = random_op(2, fx_rng)
         input_list = [0, 1, 2]
         for qargs in itertools.permutations(input_list, 2):
@@ -187,8 +187,8 @@ def test_hadamard(self, fx_rng: Generator) -> None:
         circuit.h(0)
         pattern = circuit.transpile().pattern
         pattern.standardize()
-        state = circuit.simulate_statevector().statevec
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        state = circuit.simulate().statevec
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op1 = random_op(1, fx_rng)
         value1 = state.expectation_value(random_op1, [0])
         value2 = tn_mbqc.expectation_value(random_op1, [0])
@@ -198,8 +198,8 @@ def test_s(self, fx_rng: Generator) -> None:
         circuit = Circuit(1)
         circuit.s(0)
         pattern = circuit.transpile().pattern
-        state = circuit.simulate_statevector().statevec
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        state = circuit.simulate().statevec
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op1 = random_op(1, fx_rng)
         value1 = state.expectation_value(random_op1, [0])
         value2 = tn_mbqc.expectation_value(random_op1, [0])
@@ -209,8 +209,8 @@ def test_x(self, fx_rng: Generator) -> None:
         circuit = Circuit(1)
         circuit.x(0)
         pattern = circuit.transpile().pattern
-        state = circuit.simulate_statevector().statevec
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        state = circuit.simulate().statevec
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op1 = random_op(1, fx_rng)
         value1 = state.expectation_value(random_op1, [0])
         value2 = tn_mbqc.expectation_value(random_op1, [0])
@@ -220,8 +220,8 @@ def test_y(self, fx_rng: Generator) -> None:
         circuit = Circuit(1)
         circuit.y(0)
         pattern = circuit.transpile().pattern
-        state = circuit.simulate_statevector().statevec
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        state = circuit.simulate().statevec
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op1 = random_op(1, fx_rng)
         value1 = state.expectation_value(random_op1, [0])
         value2 = tn_mbqc.expectation_value(random_op1, [0])
@@ -231,8 +231,8 @@ def test_z(self, fx_rng: Generator) -> None:
         circuit = Circuit(1)
         circuit.z(0)
         pattern = circuit.transpile().pattern
-        state = circuit.simulate_statevector().statevec
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        state = circuit.simulate().statevec
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op1 = random_op(1, fx_rng)
         value1 = state.expectation_value(random_op1, [0])
         value2 = tn_mbqc.expectation_value(random_op1, [0])
@@ -243,8 +243,8 @@ def test_rx(self, fx_rng: Generator) -> None:
         circuit = Circuit(1)
         circuit.rx(0, theta)
         pattern = circuit.transpile().pattern
-        state = circuit.simulate_statevector().statevec
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        state = circuit.simulate().statevec
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op1 = random_op(1, fx_rng)
         value1 = state.expectation_value(random_op1, [0])
         value2 = tn_mbqc.expectation_value(random_op1, [0])
@@ -255,8 +255,8 @@ def test_ry(self, fx_rng: Generator) -> None:
         circuit = Circuit(1)
         circuit.ry(0, theta)
         pattern = circuit.transpile().pattern
-        state = circuit.simulate_statevector().statevec
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        state = circuit.simulate().statevec
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op1 = random_op(1, fx_rng)
         value1 = state.expectation_value(random_op1, [0])
         value2 = tn_mbqc.expectation_value(random_op1, [0])
@@ -267,8 +267,8 @@ def test_rz(self, fx_rng: Generator) -> None:
         circuit = Circuit(1)
         circuit.rz(0, theta)
         pattern = circuit.transpile().pattern
-        state = circuit.simulate_statevector().statevec
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        state = circuit.simulate().statevec
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op1 = random_op(1, fx_rng)
         value1 = state.expectation_value(random_op1, [0])
         value2 = tn_mbqc.expectation_value(random_op1, [0])
@@ -278,8 +278,8 @@ def test_i(self, fx_rng: Generator) -> None:
         circuit = Circuit(1)
         circuit.i(0)
         pattern = circuit.transpile().pattern
-        state = circuit.simulate_statevector().statevec
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        state = circuit.simulate().statevec
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op1 = random_op(1, fx_rng)
         value1 = state.expectation_value(random_op1, [0])
         value2 = tn_mbqc.expectation_value(random_op1, [0])
@@ -290,8 +290,8 @@ def test_cz(self, fx_rng: Generator) -> None:
         circuit.cz(0, 1)
         pattern = circuit.transpile().pattern
         pattern.standardize()
-        state = circuit.simulate_statevector().statevec
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        state = circuit.simulate().statevec
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op2 = random_op(2, fx_rng)
         value1 = state.expectation_value(random_op2, [0, 1])
         value2 = tn_mbqc.expectation_value(random_op2, [0, 1])
@@ -302,8 +302,8 @@ def test_cnot(self, fx_rng: Generator) -> None:
         circuit.cnot(0, 1)
         pattern = circuit.transpile().pattern
         pattern.standardize()
-        state = circuit.simulate_statevector().statevec
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        state = circuit.simulate().statevec
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op2 = random_op(2, fx_rng)
         value1 = state.expectation_value(random_op2, [0, 1])
         value2 = tn_mbqc.expectation_value(random_op2, [0, 1])
@@ -316,8 +316,8 @@ def test_ccx(self, fx_bg: PCG64, jumps: int, fx_rng: Generator) -> None:
         circuit.ccx(0, 1, 2)
         pattern = circuit.transpile().pattern
         pattern.minimize_space()
-        state = circuit.simulate_statevector().statevec
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        state = circuit.simulate().statevec
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op3 = random_op(3, rng)
         value1 = state.expectation_value(random_op3, [0, 1, 2])
         value2 = tn_mbqc.expectation_value(random_op3, [0, 1, 2])
@@ -332,8 +332,8 @@ def test_with_graphtrans(self, fx_bg: PCG64, jumps: int, fx_rng: Generator) -> N
         pattern.shift_signals()
         pattern = pattern.infer_pauli_measurements()
         pattern.remove_pauli_measurements()
-        state = circuit.simulate_statevector().statevec
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        state = circuit.simulate().statevec
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op3 = random_op(3, rng)
         input_list = [0, 1, 2]
         for qargs in itertools.permutations(input_list):
@@ -351,8 +351,8 @@ def test_with_graphtrans_sequential(self, fx_bg: PCG64, jumps: int, fx_rng: Gene
         pattern.shift_signals()
         pattern = pattern.infer_pauli_measurements()
         pattern.remove_pauli_measurements()
-        state = circuit.simulate_statevector().statevec
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", graph_prep="sequential", rng=fx_rng)
+        state = circuit.simulate().statevec
+        tn_mbqc = pattern.simulate(backend="tensornetwork", graph_prep="sequential", rng=fx_rng)
         random_op3 = random_op(3, rng)
         input_list = [0, 1, 2]
         for qargs in itertools.permutations(input_list):
@@ -368,9 +368,9 @@ def test_coef_state(self, fx_bg: PCG64, jumps: int, fx_rng: Generator) -> None:
         pattern = circuit.transpile().pattern
         pattern.standardize()
         pattern.shift_signals()
-        statevec_ref = circuit.simulate_statevector().statevec
+        statevec_ref = circuit.simulate().statevec
 
-        tn = pattern.simulate_pattern("tensornetwork", rng=fx_rng)
+        tn = pattern.simulate("tensornetwork", rng=fx_rng)
         for number in range(len(statevec_ref.flatten())):
             coef_tn = tn.basis_coefficient(number)
             coef_sv = statevec_ref.flatten()[number]
@@ -384,9 +384,9 @@ def test_to_statevector(self, fx_bg: PCG64, nqubits: int, jumps: int, fx_rng: Ge
         pattern = circuit.transpile().pattern
         pattern.standardize()
         pattern.shift_signals()
-        statevec_ref = circuit.simulate_statevector().statevec
+        statevec_ref = circuit.simulate().statevec
 
-        tn = pattern.simulate_pattern("tensornetwork", rng=fx_rng)
+        tn = pattern.simulate("tensornetwork", rng=fx_rng)
         statevec_tn = tn.to_statevector()
 
         assert Statevec(data=statevec_tn).isclose(statevec_ref)
@@ -400,8 +400,8 @@ def test_evolve(self, fx_bg: PCG64, jumps: int, fx_rng: Generator) -> None:
         pattern.shift_signals()
         pattern = pattern.infer_pauli_measurements()
         pattern.remove_pauli_measurements()
-        state = circuit.simulate_statevector().statevec
-        tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
+        state = circuit.simulate().statevec
+        tn_mbqc = pattern.simulate(backend="tensornetwork", rng=fx_rng)
         random_op3 = random_op(3, rng)
         random_op3_exp = random_op(3, rng)
 
diff --git a/tests/test_transpiler.py b/tests/test_transpiler.py
index f3ae5f181..76487f383 100644
--- a/tests/test_transpiler.py
+++ b/tests/test_transpiler.py
@@ -64,8 +64,8 @@ def test_instructions(self, fx_bg: PCG64, jumps: int, instruction: InstructionTe
         circuit = Circuit(3, instr=[instruction(rng)])
         pattern = circuit.transpile().pattern
         input_state = rand_state_vector(3, rng=rng)
-        state = circuit.simulate_statevector(input_state=input_state).statevec
-        state_mbqc = pattern.simulate_pattern(input_state=input_state, rng=rng)
+        state = circuit.simulate(input_state=input_state).statevec
+        state_mbqc = pattern.simulate(input_state=input_state, rng=rng)
         assert state_mbqc.isclose(state)
 
     def test_transpiled(self, fx_rng: Generator) -> None:
@@ -74,8 +74,8 @@ def test_transpiled(self, fx_rng: Generator) -> None:
         pairs = [(i, np.mod(i + 1, nqubits)) for i in range(nqubits)]
         circuit = rand_gate(nqubits, depth, pairs, fx_rng, use_rzz=True)
         pattern = circuit.transpile().pattern
-        state = circuit.simulate_statevector(rng=fx_rng).statevec
-        state_mbqc = pattern.simulate_pattern(rng=fx_rng)
+        state = circuit.simulate(rng=fx_rng).statevec
+        state_mbqc = pattern.simulate(rng=fx_rng)
         assert state_mbqc.isclose(state)
 
     @pytest.mark.parametrize("backend", ["statevector", "densitymatrix"])
@@ -91,11 +91,11 @@ def test_measure(
         circuit.m(0, axis)
         input_state = rand_state_vector(2, rng=rng)
         branch_selector = ConstBranchSelector(outcome)
-        state = circuit.simulate_statevector(
+        state = circuit.simulate(
             rng=rng, input_state=input_state, branch_selector=branch_selector, backend=backend
         ).statevec
         pattern = circuit.transpile().pattern
-        state_mbqc = pattern.simulate_pattern(
+        state_mbqc = pattern.simulate(
             rng=rng, input_state=input_state, branch_selector=branch_selector, backend=backend
         )
         if isinstance(state_mbqc, Statevec) and isinstance(state, Statevec):
@@ -113,9 +113,9 @@ def test_measure_early(self, fx_bg: PCG64, jumps: int, axis: Axis, outcome: Outc
         circuit.cnot(1, 2)
         input_state = rand_state_vector(3, rng=rng)
         branch_selector = ConstBranchSelector(outcome)
-        state = circuit.simulate_statevector(rng=rng, input_state=input_state, branch_selector=branch_selector).statevec
+        state = circuit.simulate(rng=rng, input_state=input_state, branch_selector=branch_selector).statevec
         pattern = circuit.transpile().pattern
-        state_mbqc = pattern.simulate_pattern(rng=rng, input_state=input_state, branch_selector=branch_selector)
+        state_mbqc = pattern.simulate(rng=rng, input_state=input_state, branch_selector=branch_selector)
         assert state_mbqc.isclose(state)
 
     @pytest.mark.parametrize("input_axis", [Axis.X, Axis.Y, Axis.Z])
@@ -141,7 +141,7 @@ def test_measurement_expectation_value(
         circuit.m(0, measurement_axis)
         expectation_value0 = 0.5 if input_axis != measurement_axis else 1 if input_sign == Sign.PLUS else 0
         branch_selector = CheckedBranchSelector(expected={0: expectation_value0}, abs_tol=1e-15)
-        circuit.simulate_statevector(input_state=input_state, branch_selector=branch_selector, rng=fx_rng)
+        circuit.simulate(input_state=input_state, branch_selector=branch_selector, rng=fx_rng)
 
     @pytest.mark.parametrize("jumps", range(1, 11))
     @pytest.mark.parametrize("axis", [Axis.X, Axis.Y, Axis.Z])
@@ -152,12 +152,10 @@ def test_transpile_measurements_to_z_axis(self, fx_bg: PCG64, jumps: int, axis:
         circuit.m(0, axis)
         input_state = rand_state_vector(2, rng=rng)
         branch_selector = ConstBranchSelector(outcome)
-        state = circuit.simulate_statevector(rng=rng, input_state=input_state, branch_selector=branch_selector).statevec
+        state = circuit.simulate(rng=rng, input_state=input_state, branch_selector=branch_selector).statevec
         circuit_z = circuit.transpile_measurements_to_z_axis()
         assert all(instr.axis == Axis.Z for instr in circuit_z.instruction if instr.kind == InstructionKind.M)
-        state_z = circuit.simulate_statevector(
-            rng=rng, input_state=input_state, branch_selector=branch_selector
-        ).statevec
+        state_z = circuit.simulate(rng=rng, input_state=input_state, branch_selector=branch_selector).statevec
         assert state_z.isclose(state)
 
     @pytest.mark.parametrize("jumps", range(1, 11))
@@ -170,8 +168,8 @@ def test_transpile_j_to_rzh(self, fx_bg: PCG64, jumps: int) -> None:
         assert any(instr.kind == InstructionKind.J for instr in circuit.instruction)
         circuit2 = circuit.transpile_j_to_rzh()
         assert not any(instr.kind == InstructionKind.J for instr in circuit2.instruction)
-        state = circuit.simulate_statevector(rng=rng).statevec
-        state2 = circuit2.simulate_statevector(rng=rng).statevec
+        state = circuit.simulate(rng=rng).statevec
+        state2 = circuit2.simulate(rng=rng).statevec
         assert state.fidelity(state2) == pytest.approx(1)
 
     @pytest.mark.parametrize("jumps", range(1, 11))
@@ -191,10 +189,8 @@ def test_transpile_swaps_with_measurements(self, fx_bg: PCG64, jumps: int, axis:
         assert I(2) in circuit2.instruction
         input_state = rand_state_vector(3, rng=rng)
         branch_selector = ConstBranchSelector(outcome)
-        state = circuit.simulate_statevector(rng=rng, input_state=input_state, branch_selector=branch_selector).statevec
-        state2 = circuit2.simulate_statevector(
-            rng=rng, input_state=input_state, branch_selector=branch_selector
-        ).statevec
+        state = circuit.simulate(rng=rng, input_state=input_state, branch_selector=branch_selector).statevec
+        state2 = circuit2.simulate(rng=rng, input_state=input_state, branch_selector=branch_selector).statevec
         assert transpiled_swaps.outputs == (
             OutputIndex(OutputKind.Qubit, 2),
             OutputIndex(OutputKind.Bit, 0),
@@ -211,9 +207,9 @@ def test_cz_ccx(self, fx_rng: Generator) -> None:
         circuit = Circuit(width=3)
         circuit.cz(2, 0)
         circuit.ccx(0, 1, 2)
-        ref_state = circuit.simulate_statevector(rng=fx_rng).statevec
+        ref_state = circuit.simulate(rng=fx_rng).statevec
         pattern = circuit.transpile().pattern
-        state = pattern.simulate_pattern(rng=fx_rng)
+        state = pattern.simulate(rng=fx_rng)
         assert state.isclose(ref_state)
 
     def test_ccx_decomposition(self) -> None:
@@ -223,16 +219,16 @@ def test_ccx_decomposition(self) -> None:
         circuit2 = Circuit(width=3)
         circuit2.cz(2, 0)
         circuit2.extend(decompose_ccx(instruction.CCX(controls=(0, 1), target=2)))
-        state = circuit.simulate_statevector().statevec
-        state2 = circuit2.simulate_statevector().statevec
+        state = circuit.simulate().statevec
+        state2 = circuit2.simulate().statevec
         assert state.isclose(state2)
 
     def test_cnot_cz(self, fx_rng: Generator) -> None:
         """Test regression about output node reordering."""
         circuit = Circuit(width=3, instr=[instruction.CNOT(0, 1), instruction.CZ((0, 1))])
-        state = circuit.simulate_statevector(rng=fx_rng).statevec
+        state = circuit.simulate(rng=fx_rng).statevec
         pattern = circuit.transpile().pattern
-        state_mbqc = pattern.simulate_pattern(rng=fx_rng)
+        state_mbqc = pattern.simulate(rng=fx_rng)
         assert state.isclose(state_mbqc)
 
     @pytest.mark.parametrize("jumps", range(1, 6))
@@ -257,12 +253,12 @@ def test_classical_outputs_consistency(self, fx_bg: PCG64, jumps: int, axes: lis
         results_pattern: dict[int, Outcome] = {node: m_outcomes.get(node, 0) for node in non_output_nodes}
         input_state = rand_state_vector(width, rng=rng)
         measure_method = DefaultMeasureMethod()
-        circuit_result = circuit.simulate_statevector(
+        circuit_result = circuit.simulate(
             rng=rng,
             input_state=input_state,
             branch_selector=FixedBranchSelector(results=results_circuit),
         )
-        pattern.simulate_pattern(
+        pattern.simulate(
             rng=rng,
             input_state=input_state,
             branch_selector=FixedBranchSelector(results=results_pattern),
@@ -280,7 +276,7 @@ def test_classical_outputs_empty(self) -> None:
         circuit.h(0)
         result = circuit.transpile()
         assert len(result.classical_outputs) == 0
-        assert len(circuit.simulate_statevector().classical_measures) == 0
+        assert len(circuit.simulate().classical_measures) == 0
 
 
 class TestCircuits:
@@ -318,8 +314,8 @@ def test_instructions(self, fx_bg: PCG64, jumps: int, instruction: InstructionTe
         circuit = Circuit(3, instr=[instruction(rng)])
         pattern = circuit.transpile().pattern
         input_state = rand_state_vector(3, rng=rng)
-        state = circuit.simulate_statevector(input_state=input_state).statevec
-        state_mbqc = pattern.simulate_pattern(input_state=input_state, rng=rng)
+        state = circuit.simulate(input_state=input_state).statevec
+        state_mbqc = pattern.simulate(input_state=input_state, rng=rng)
         assert state_mbqc.isclose(state)
 
     def test_simple(self) -> None:
@@ -328,8 +324,8 @@ def test_simple(self) -> None:
         pattern = circuit.transpile().pattern
         pattern.minimize_space()
         input_state = rand_state_vector(3, rng=rng)
-        state = circuit.simulate_statevector(input_state=input_state).statevec
-        state_mbqc = pattern.simulate_pattern(input_state=input_state, rng=rng)
+        state = circuit.simulate(input_state=input_state).statevec
+        state_mbqc = pattern.simulate(input_state=input_state, rng=rng)
         assert state_mbqc.isclose(state)
 
     @pytest.mark.parametrize("jumps", range(1, 3))
@@ -340,8 +336,8 @@ def test_dm_backend(self, fx_bg: PCG64, jumps: int) -> None:
         pattern = circuit.transpile().pattern
         pattern.minimize_space()
         input_state = rand_state_vector(nqubits, rng=rng)
-        state = circuit.simulate_statevector(input_state=input_state, backend="densitymatrix").statevec
-        state_mbqc = pattern.simulate_pattern(input_state=input_state, backend="densitymatrix", rng=rng)
+        state = circuit.simulate(input_state=input_state, backend="densitymatrix").statevec
+        state_mbqc = pattern.simulate(input_state=input_state, backend="densitymatrix", rng=rng)
         assert np.allclose(state_mbqc.rho, state.rho)
 
 
@@ -355,8 +351,8 @@ def test_transpile_swaps(fx_bg: PCG64, jumps: int) -> None:
     transpiled_swaps = transpile_swaps(circuit)
     circuit2 = transpiled_swaps.circuit
     assert not any(instr.kind == InstructionKind.SWAP for instr in circuit2.instruction)
-    state = circuit.simulate_statevector(rng=rng).statevec
-    state2 = circuit2.simulate_statevector(rng=rng).statevec
+    state = circuit.simulate(rng=rng).statevec
+    state2 = circuit2.simulate(rng=rng).statevec
     state2.permute(transpiled_swaps.extract_output_node_indices())
     assert state.isclose(state2)
 
@@ -378,8 +374,8 @@ def test_transpile_swaps_with_measurements(fx_bg: PCG64, jumps: int, axis: Axis,
     assert I(2) in circuit2.instruction
     input_state = rand_state_vector(3, rng=rng)
     branch_selector = ConstBranchSelector(outcome)
-    state = circuit.simulate_statevector(rng=rng, input_state=input_state, branch_selector=branch_selector).statevec
-    state2 = circuit2.simulate_statevector(rng=rng, input_state=input_state, branch_selector=branch_selector).statevec
+    state = circuit.simulate(rng=rng, input_state=input_state, branch_selector=branch_selector).statevec
+    state2 = circuit2.simulate(rng=rng, input_state=input_state, branch_selector=branch_selector).statevec
     assert transpiled_swaps.outputs == (
         OutputIndex(OutputKind.Qubit, 2),
         OutputIndex(OutputKind.Bit, 0),
@@ -405,12 +401,12 @@ def test_transpile_swaps_vs_no_transpile_swaps() -> None:
     circuit.swap(0, 1)
     pattern_without_swap = circuit.transpile().pattern
     pattern_with_swap = circuit.transpile(transpile_swaps=False).pattern
-    state_without_swap = pattern_without_swap.simulate_pattern()
-    state_with_swap = pattern_with_swap.simulate_pattern()
+    state_without_swap = pattern_without_swap.simulate()
+    state_with_swap = pattern_with_swap.simulate()
     assert state_without_swap.isclose(state_with_swap)
 
 
 def test_backend_branch_selector() -> None:
     circ = Circuit(1)
     with pytest.raises(ValueError, match="already instantiated"):
-        circ.simulate_statevector(backend=StatevectorBackend(), branch_selector=ConstBranchSelector(0))
+        circ.simulate(backend=StatevectorBackend(), branch_selector=ConstBranchSelector(0))