Product-formula (Trotter-Suzuki) simulation

What it fills

• Simulate Hamiltonian evolution

  Implement $e^{-iHt}$ to error $\varepsilon$ given some access model for $H$. It is an application in its own right and also the engine inside phase estimation and several linear-system solvers.

When it applies

Needs a decomposition of $H$ into efficiently exponentiable summands. The error is governed by commutators among those summands: Childs, Su, Tran, Wiebe and Zhu's analysis "directly exploits the commutativity of operator summands, producing tighter error bounds for both real- and imaginary-time evolutions", and they show local observables can be simulated with complexity independent of the system size for power-law interacting systems. Term ordering affects the error and is not fixed by the formula itself.

Requires

Every step this method names moves its route along, so there is nothing it needs alongside them.

Example

given  H = sum_j H_j  with each H_j efficiently exponentiable,
       evolution time t, step count r,  h = t/r

for step = 1 ... r:
    # first order (Lie-Trotter): one short evolution per summand, in order
    for j = 1 ... m:   apply exp(-i h H_j)
    # higher order: the same summands in the symmetrised sub-step sequence
    # of the chosen Suzuki formula

# no block-encoding is built, and there is no all-zeros flag to amplify

# the error is governed by commutators among the summands, and the term
# order affects it -- the formula itself does not fix that order

Cost, as the source states it

In the sparse-access oracle model, Berry, Ahokas, Cleve and Sanders: when $H$ acts on $n$ qubits, has at most a constant number of nonzero entries in each row/column, and $||H||$ is bounded by a constant, one may select any positive integer $k$ such that the simulation requires $O((\log* n)·t^{1+1/2k})$ accesses to matrix entries of $H$. They also show the temporal scaling cannot be significantly improved beyond this, because sublinear time scaling is not possible.

Implementations

• Qiskit's product-formula evolution synthesis (LieTrotter and SuzukiTrotter)

  IBM's Qiskit ships product formulas not as circuits but as synthesis strategies plugged into one gate. `PauliEvolutionGate` is documented as the "Time-evolution of an operator consisting of Paulis": for a Hermitian $H$ consisting of Pauli terms and real evolution time $t$, "this gate represents the unitary" $U(t) = e^{-itH}$. The release note that introduced it states the separation outright — "The synthesis of this gate is performed by `EvolutionSynthesis` and is decoupled from the gate itself" — and the gate's own docstring gives the reason for it: "This gate represents the exact evolution $U(t)$. Implementing this operation exactly, however, generally requires an exponential number of gates. The compiler therefore typically implements an *approximation* of the unitary $U(t)$, e.g. using a product formula such as defined by `LieTrotter`." The gate and three strategies arrived together in Qiskit 0.19, whose release note lists them as `LieTrotter`, "first order Trotterization"; `SuzukiTrotter`, "higher order Trotterization"; and `MatrixExponential`, "exact, matrix-based evolution". Only the first two are product formulas.

  `LieTrotter` is not a separate construction: it subclasses `SuzukiTrotter` and calls `super().__init__(1, ...)`, so first order is the order-1 case of the same recursion. That recursion is `SuzukiTrotter._recurse(order, grouped_paulis)`. At order 1 it returns the groups untouched. At order 2 it halves the coefficients of every group but the last and returns `halves + full + list(reversed(halves))` — the last summand at full time between two mirrored halves. Above that it sets `reduction = 1 / (4 - 4 ** (1 / (order - 1)))`, which is the coefficient this record's hop writes $u_{k} := 1/(4 - 4^{1/(2k-1)})$ for a formula of order $2k$, and returns `outer + inner + outer`, where `outer` is the expansion of order $2k-2$ at `reduction` taken twice and `inner` is the same expansion at `1 - 4 * reduction`. Odd orders above 1 are refused in the constructor: "Suzuki product formulae are symmetric and therefore only defined for when the order is 1 or even". `expand` returns the sequence as `(pauli, indices, angle)` triples with the angle `coeff * time * 2 / self.reps`, the docstring warning that "the rotation angle contains a factor of 2". By default each Pauli term becomes its own non-commuting factor, and the source says so against itself — "Assume no commutativity here. If we were to group commuting Paulis, here would be the location to do so." — a caller wanting grouping passes a list of operators whose terms already commute. The term ordering this record's `conditions` says the formula does not fix is exposed as `preserve_order`: set to `False`, `reorder_paulis` builds a graph joining any two terms that "act on the same qubit", colours it with rustworkx's `graph_greedy_color` under `ColoringStrategy.Saturation`, and emits terms colour by colour so that same-colour subcircuits "can be run in parallel". It is documented as "deterministic and invariant under permutation of the Pauli term in ``paulis``". Two further knobs are gate-level rather than formula-level: `cx_structure`, `"chain"` for "next neighbor connections" or `"fountain"` where "all qubits are connected to one", and `insert_barriers`.

  Nothing in the repository runs the module on a physical system; what exercises it is `test/python/circuit/library/test_evolution_gate.py`. `test_lie_trotter` and `test_suzuki_trotter` both evolve `(X ^ X ^ X) + (Y ^ Y ^ Y) + (Z ^ Z ^ Z)` for `time = 0.123` with `reps = 4` — the first at first order, the second at orders 2, 4 and 6 — and assert CX counts derived from the recursion depth, the test's own comment reading "recurse (order - 2) / 2 times, base case has 5 blocks with 4 CX each". `test_suzuki_trotter_manual_no_reorder` evolves `X + Y` on one qubit at order 4 against a circuit built by hand from `p_4 = 1 / (4 - 4 ** (1 / 3))`, commented "coefficient for reduced time from Suzuki paper". `test_suzuki_trotter_manual` evolves `(X ^ X ^ I ^ I) + (I ^ Y ^ Y ^ I) + (I ^ I ^ Z ^ Z)` at order 2 with `preserve_order = False` and asserts an exact five-gate circuit. The accuracy helper is narrower than it looks: `assertSuzukiTrotterIsCorrect` compares the synthesised circuit against a `SuzukiTrotter` of the same `order` and `reps` built with an `exact_atomic_evolution`, so it checks that the decomposition reproduces the product formula, not that the product formula is close to $e^{-itH}$. `scipy.linalg.expm` does appear in six of that file's tests, but wherever a product-formula circuit is the thing compared against it, the instance is one on which the formula carries no error: `test_global_phase` binds the default-synthesis `evo.definition` and asserts equality with `expm` under its own stated condition, "since all terms in the Pauli operator commute, we can compare to an exact matrix exponential". `test_matrix_decomposition` runs two subtests and neither exercises the recursion — the first synthesises with `MatrixExponential`, and the second calls `PauliEvolutionGate.to_matrix()`, which is itself `sc.sparse.linalg.expm(-1j * time * spmatrix)` on the operator and never consults the synthesis. The one test that measures an approximation against `expm` is `test_qdrift_evolution`, and its subject is `QDrift`. The largest instance named anywhere is not in the tests but in the Qiskit 1.3 release note: a 100-qubit Heisenberg Hamiltonian with 10 timesteps at fourth order. No molecule, lattice or device is named.

  The module is `qiskit/synthesis/evolution/` in https://github.com/Qiskit/qiskit, Apache License 2.0. `product_formula.py` holds the `ProductFormula` base and `reorder_paulis`, `suzuki_trotter.py` the recursion, `lie_trotter.py` the order-1 subclass; beside them sit `qdrift.py` (Campbell's randomised `QDrift`, a different construction), `matrix_synthesis.py` (`MatrixExponential`), `evolution_synthesis.py` and `pauli_network.py`. Seven of the directory's eight Python files are headed "(C) Copyright IBM 2021"; the eighth, `pauli_network.py`, is headed "(C) Copyright IBM 2024". The gate is `qiskit/circuit/library/pauli_evolution.py`, and with no `synthesis` argument it constructs `LieTrotter()` — first order, one repetition. `ProductFormula.synthesize` has two paths. With `wrap=True` or a caller-supplied `atomic_evolution` it goes through `_custom_evolution`, where those two cases are themselves disjoint: given a custom `atomic_evolution` the loop builds a `SparseObservable` for each Pauli rotation and calls `self.atomic_evolution(circuit, operator, time=1)`, appending whatever the caller's function queues and creating no gate object of its own; only the other branch, commented "this means self._wrap is True", builds each Pauli's circuit Rust-side and wraps it with `evo.to_gate(label=f"exp(it {pauli_string})")`. Otherwise the whole network is built Rust-side in one call through `from qiskit._accelerate.circuit_library import pauli_evolution`. Since Qiskit 1.3 the transpiler reaches all of it through high-level-synthesis plugins in `qiskit/transpiler/passes/synthesis/hls_plugins.py`. The second of those, `PauliEvolutionSynthesisRustiq`, is not an alternative product formula: it calls `algo.expand(evo)` on the same `ProductFormula`, warns "Cannot apply Rustiq if the evolution synthesis does not implement ``expand``", and differs only in how the resulting Pauli network becomes gates. Read at tag 2.5.2, released 2026-08-13. The same recursion also exists in Rust at `crates/synthesis/src/evolution/suzuki_trotter.rs`, where the coefficient is written `1.0 / (4.0 - 4_f64.powf(1.0 / (order as f64 - 1.0)))`, reached from `crates/circuit_library/src/suzuki_trotter.rs`.

  One number is reported for this module and it is a construction time, not an accuracy and not a gate count. The Qiskit 1.3 release note says the default plugin "constructs circuit as before, but faster as it internally uses Rust", that "the larger the circuit (e.g. by the Hamiltonian size, the number of timesteps, or the Suzuki-Trotter order), the higher the speedup", and gives a single instance: "a 100-qubit Heisenberg Hamiltonian with 10 timesteps and a 4th-order Trotter formula is now constructed ~9.4x faster". The baseline is named only as constructing "as before", no machine is stated, and the quantity measured is the time to build the circuit rather than to run it. The other quantities the repository pins are exact circuit contents, not measurements: the order-2 expansion of `H = IX + ZZ` at one repetition is documented as `("X", [0], t), ("ZZ", [0, 1], 2t), ("X", [0], t)`, and the four-qubit reordering test fixes the output to `rzz(time, 0, 1)`, `rxx(time, 2, 3)`, `ryy(2 * time, 1, 2)`, `rxx(time, 2, 3)`, `rzz(time, 0, 1)`.

• PennyLane's TrotterProduct operation

  Xanadu's PennyLane ships the same family as a single operation rather than as a compiler strategy. `TrotterProduct` arrived in v0.33.0, added by pull request 4661, whose release note introduces it as "Higher-order Trotter-Suzuki methods are now easily accessible through a new operation called `TrotterProduct`" and tells the caller: "Simply specify the `order` of the approximation and the evolution `time`." What makes it worth recording beside the Qiskit module is not the recursion, which is the same one, but that in the release read here — v0.45.1 — it ships the error analysis with it: there the class is declared `TrotterProduct(ErrorOperation, ResourcesOperation)` and carries an `error()` method returning "an *upper-bound* on the spectral norm error", implemented from the paper this record's `conditions` credits for the commutator bounds. That method is a property of the version read rather than of the operation for all time, and every claim about it below is scoped to v0.45.1.

  The class docstring states the recursion in the same symbols this record's hop uses: symmetrised products $S_{m}(t) = S_{m-2}(p_{m}t)^{2} \cdot S_{m-2}((1-4p_{m})t) \cdot S_{m-2}(p_{m}t)^{2}$ "where the coefficient is" $p_{m} = 1 / (4 - \sqrt[m - 1]{4})$, and the $m$-th order, $n$-step approximation $e^{iHt} \approx [S_{m}(t / n)]^{n}$. The sign convention is documented opposite to Qiskit's: `time` is "the parameter $t$ in $e^{iHt}$" where `PauliEvolutionGate` documents $e^{-itH}$, and the docstring's own way of recovering the older `ApproxTimeEvolution` is `qp.adjoint(qp.TrotterProduct(hamiltonian, time, order=1, n=n))`. In code, `_scalar(order)` returns `(4 - 4**root) ** -1` with `root = 1 / (order - 1)`, and `_recursive_expression` returns `(2 * ops_lst_1) + ops_lst_2 + (2 * ops_lst_1)`. Its order-2 base case is `[Evolution(op, -x * 0.5) for op in ops + ops[::-1]]` — every term at half time forwards then backwards, so the last term appears as two adjacent half-time exponentials of the same operator where Qiskit's `_recurse` emits one at full time; the sequences differ by that merge, not by the formula. Two constraints on the input have no counterpart in Qiskit. The terms need not be Paulis: the argument is documented as "the Hamiltonian written as a linear combination of operators with known matrix exponentials", and the constructor requires a `Sum` (converting a `LinearCombination` or `SProd` first) with at least two terms — "There should be at least 2 terms in the Hamiltonian. Otherwise use `qp.exp`". And `check_hermitian=True` by default verifies every operand. On term order the class does not choose for you and says so in a standing warning: "The Trotter-Suzuki decomposition depends on the order of the summed observables. Two mathematically identical `LinearCombination` objects may undergo different time evolutions due to the order in which those observables are stored" — `compute_decomposition` reads `kwargs["base"].operands` in stored order, where Qiskit optionally recolours it. `error()` offers two bounds, both attributed by section: "one-norm-bound" to "Section 2.3 (lemma 6, equation 22 and 23)" and the default "commutator-bound" to "Appendix C (equation 189)" of Childs et al. (2021), https://arxiv.org/abs/1912.08854. `_one_norm_error` computes `(h_one_norm * t) ** (p + 1) / (math.factorial(p + 1) * n**p)` times `(upsilon ** (p + 1) + 1)`; `_commutator_error` sums nested commutators over every combination of exponents summing to `p`, with prefactor `(2 * upsilon * t ** (p + 1)) / ((p + 1) * n**p)`. Both use $\Upsilon$ = `_compute_repetitions`, documented against "the definition of upsilon from section 2.3 (equation 15)" of the same paper and returning 1 at order 1 and `(5 ** (k - 1)) * 2` otherwise. One qualifier the default hides: `fast=True` replaces each spectral norm with an upper bound — the sum of the absolute Pauli coefficients where the operator has a Pauli representation, otherwise the Frobenius norm — so the number returned by default is a bound computed from bounds.

  The documented runs are all small and all in software. The v0.33.0 release note builds `H = qml.dot([0.25, 0.75], [qml.PauliX(0), qml.PauliZ(0)])`, applies a Hadamard and then `qml.TrotterProduct(H, time=2.4, order=2)` on `qml.device("default.qubit", wires=2)`, and returns the state. The `error()` examples use `qp.dot([1.0, 0.5, -0.25], [qp.X(0), qp.Y(0), qp.Z(0)])` at `time=0.01, order=2` and compute a bound rather than running a circuit. The usage-details example is a four-term two-qubit Hamiltonian, `[qp.X(0), qp.Y(1), qp.Y(0) @ qp.Z(1), qp.X(0) @ qp.Y(1)]` with coefficients `[0.5, 0.2, 0.1, -0.6]`, used only to show that grouping the operands changes the decomposition. `default.qubit` is the only device named in the `TrotterProduct` docstring; no dataset, molecule or lattice appears, and no hardware run is reported.

  `pennylane/templates/subroutines/time_evolution/trotter.py` in https://github.com/PennyLaneAI/pennylane — Python, Apache License 2.0, headed "Copyright 2018-2023 Xanadu Quantum Technologies Inc." and described in its module docstring as "templates for Suzuki-Trotter approximation based subroutines". The path moved: at v0.33.0, where the class was added, the same file was `pennylane/templates/subroutines/trotter.py`. The two error bounds are not defined in it — the module imports `_one_norm_error` and `_commutator_error` `from pennylane.resource.error` for `error()` to call, and both live in `pennylane/resource/error/trotter_error.py`, which holds nine functions and no classes: those two plus `_compute_repetitions`, `_spectral_norm`, `_generate_combinations`, `_recursive_nested_commutator` and the flattening helpers `_recursive_flatten`, `_simplify` and `_flatten_trotter`. Back in `trotter.py` itself, beneath `TrotterProduct`, sit `TrotterizedQfunc` and the `trotterize` function, which apply the identical recursion to an arbitrary quantum function instead of to a Hamiltonian — the docstring's framing is that a caller supplies "the first order expansion $S_{1}(t)$" and the machinery "expands it to any higher order". Read at tag v0.45.1, released 2026-06-26.

  What the sources report are documented outputs, not benchmarks. The v0.33.0 release note prints the evolved state for its two-term Hamiltonian as `[-0.13259524+0.59790098j 0. +0.j -0.13259524-0.77932754j 0. +0.j ]`. The `error()` docstring prints, for `qp.dot([1.0, 0.5, -0.25], [qp.X(0), qp.Y(0), qp.Z(0)])` at `time=0.01, order=2`, `SpectralNormError(8.039062500000003e-06)` under `method="one-norm-bound"` and `SpectralNormError(6.166666666666668e-06)` under `method="commutator-bound"` — on that example the commutator bound is the smaller of the two, and it is the method's default. Both are upper bounds computed classically from the operator, not measured errors, and both come from the default `fast=True` path. No gate count, runtime or hardware result is reported for this operation.

What it needs

Nothing below this — it bottoms out here.

Slots it makes unnecessary

These do not fill the slot. They replace the span it belongs to, so this layer is not on their path at all.

• Block-encode a matrix

  Wrap an operator $A$ inside a larger unitary $U$ so that $A/α$ sits in $U$'s top-left block, giving every routine above it one uniform way to touch the matrix. The subnormalization $α$ and the ancilla count are outputs of this layer, not free parameters.

Other ways to fill the same slot

Different approaches

• Truncated-Taylor LCU simulation

  Truncate the Taylor series of $e^{-iHt}$ over short segments and implement the truncated sum as a linear combination of unitaries — PREPARE loads the coefficients, SELECT applies the terms, PREPARE$†$ unprepares — with the all-zeros ancilla flag amplified.

• Qubitization walk simulation

  From a block-encoding pair $(U, |G⟩)$ with $H = (⟨G|⊗I)U(|G⟩⊗I)$, build a walk operator $W$ that splits the Hilbert space into invariant two-dimensional $SU(2)$ subspaces, one per eigenvalue of $H$, with eigenvalues $e^{±i·arccos(H/α)}$. Quantum signal processing phases applied to $W$ then produce $e^{-iHt}$.

In the Atlas

• Trotter–Suzuki Hamiltonian simulation

  Approximates the time evolution of a Hamiltonian with non-commuting terms by alternating their individual exponentials.

Sources

• Efficient quantum algorithms for simulating sparse Hamiltonians↗Dominic W. Berry, Graeme Ahokas, Richard Cleve, Barry C. Sanders · 2005theory yes · simulation ? · hardware no
• A Theory of Trotter Error↗Andrew M. Childs, Yuan Su, Minh C. Tran, Nathan Wiebe, Shuchen Zhu · 2019theory yes · simulation ? · hardware no
• Simulating Hamiltonian dynamics with a truncated Taylor series↗Dominic W. Berry, Andrew M. Childs, Richard Cleve, Robin Kothari, Rolando D. Somma · 2014theory yes · simulation ? · hardware no