qcr:2606.69981.1

Quantum Real-Time Evolution with Trotterization

This Qiskit tutorial implements quantum real-time evolution using Trotterization, the most direct approach to simulating how a quantum state evolves under a Hamiltonian on a quantum computer. The exact time-evolution operator is the exponential of the Hamiltonian, which is generally hard to implement directly; the Trotter-Suzuki decomposition approximates it by breaking the evolution into many small time steps and, within each step, applying the exponentials of the individual Hamiltonian terms in sequence, an approximation that becomes exact as the step size shrinks. The tutorial shows how to use Qiskit's TrotterQRTE (Trotter quantum real-time evolution) with the algorithms library: defining a Hamiltonian as a sum of Pauli terms, specifying an initial state and an evolution time, choosing the Trotter order and number of steps, and running the evolution through the Estimator primitive to track observables over time. It discusses the trade-off between Trotter step size (accuracy) and circuit depth (hardware feasibility), and how Trotterization compares to the variational time-evolution methods. It is a clear, foundational introduction to Hamiltonian simulation via Trotterization in Qiskit, the building block behind digital quantum simulation.

Quantum Simulation

Qubit

Circuit-based

Uploaded 1 day ago

Views

GitHub

183

Citing this entry? Use this QCR ID

Uploaded by

QCR Librarian

qiskit-community/qiskit-algorithms

183

In [ ]:

# --- Setup cell added by QCR (not part of the original tutorial) ---
# Source: qiskit-community/qiskit-algorithms @ 0.4.0, Apache License 2.0.
# Installs the example's dependencies. If a later cell still reports a missing
# package, restart the runtime/kernel and run again from the top.
%pip install -q qiskit-algorithms==0.4.0 qiskit-aer matplotlib pylatexenc qiskit-algorithms

Quantum Real Time Evolution using Trotterization

As a real time evolution technique, the Trotterization or Trotterized Real Time Evolution (RTE) consists in the successive application of a quantum gate, assumed to approximate the time evolution of a system for a time slice [1]. Following from the Schrödinger equation, the time evolution of a system initially in the state takes the form:

where is the time-independent Hamiltonian under which the system is submitted. We consider a Hamiltonian that can be written as a weighted sum of Pauli terms , with representing a tensor product of Pauli terms acting on qubits. Thus, the sequence of gates that approximate a single evolution step is found under the assumption that the product of matrix exponentials is a good approximation for the sum of matrix exponentials:

where is the time slice (evolution step) that we are choosing. As a result, a gate to be applied times is created. A smaller timestep leads to a more accurate approximation. However, this also leads to deeper circuits which, in practice, leads to more error accumulation, becoming a critical aspect especially on NISQ devices. For the purposes of this tutorial, all simulations are Statevector simulator, which doesn't include any shot or hardware noise.

In this tutorial, we will implement real time evolutions using the TrotterQRTE class. To illustrate this, we will study the time evolution of the Ising model on linear lattices of and sites. These lattices consist of an array of spins that interact only with their nearest neighbors. These spins can have two orientations: and , which correspond to a magnetization of and respectively.

where describes the interaction energy, and the magnitude of an external field. Let us write this expression using Pauli matrices, and considering that the external field has an angle with respect to the transversal direction,

In the computational basis, the system will be encoded as follows:

Quantum state	Spin representation

First, we will create a function that takes in the system parameters , , and , and returns our Hamiltonian as a SparsePauliOp. A SparsePauliOp is a sparse representation of an operator in terms of weighted Pauli terms.

In [1]:

from qiskit.quantum_info import SparsePauliOp
from math import sin, cos


def get_hamiltonian(L, J, h, alpha=0):

    # List of Hamiltonian terms as 3-tuples containing
    # (1) the Pauli string,
    # (2) the qubit indices corresponding to the Pauli string,
    # (3) the coefficient.
    ZZ_tuples = [("ZZ", [i, i + 1], -J) for i in range(0, L - 1)]
    Z_tuples = [("Z", [i], -h * sin(alpha)) for i in range(0, L)]
    X_tuples = [("X", [i], -h * cos(alpha)) for i in range(0, L)]

    # We create the Hamiltonian as a SparsePauliOp, via the method
    # `from_sparse_list`, and multiply by the interaction term.
    hamiltonian = SparsePauliOp.from_sparse_list([*ZZ_tuples, *Z_tuples, *X_tuples], num_qubits=L)
    return hamiltonian.simplify()

Let us get started, and create a Hamiltonian as an operator for , and a choice of , and :

In [2]:

from math import pi

H = get_hamiltonian(L=2, J=0.2, h=1.0, alpha=pi / 8)
H

SparsePauliOp(['ZZ', 'IZ', 'ZI', 'IX', 'XI'],
              coeffs=[-0.2       +0.j, -0.38268343+0.j, -0.38268343+0.j, -0.92387953+0.j,
 -0.92387953+0.j])

Let us create an instance of TimeEvolutionProblem. Conceptually, it contains all the information relevant on the physical problem. In our case, we will pass it the our Hamiltonian, an initial state, a final time. As initial state, we will take a spin up and a spin down.

In [3]:

from qiskit.quantum_info import Statevector
from qiskit_algorithms import TimeEvolutionProblem

final_time = 1.6

# First spin up, second spin down
# (remember that the labels are interpreted from right to left)
initial_state = Statevector.from_label("10")

problem = TimeEvolutionProblem(H, initial_state=initial_state, time=final_time)

We can now create an instance of TrotterQRTE. Conceptually, it contains the information relevant to solve a physical problem, which does by means of the method evolve.

In [4]:

from qiskit_algorithms import TrotterQRTE

trotter = TrotterQRTE()
result = trotter.evolve(problem)

Our evolved state is in the result's attribute evolved_state, which is a QuantumCircuit,

In [5]:

result.evolved_state.draw("mpl")

By sequentially decomposing the circuit, we can show it in terms of Qiskit's Circuit Library standard gates.

In [6]:

result.evolved_state.decompose(reps=5).decompose("disentangler_dg").decompose(
    "multiplex1_reverse_dg"
).draw("mpl")

The evolved state, like all QuantumCircuits, can be turned into a Statevector.

In [7]:

statevector = Statevector(result.evolved_state)
print(statevector)

Statevector([ 0.02895926+0.08738739j, -0.9411211 +0.31187756j,
              0.00811432-0.002689j  ,  0.02895926+0.08738739j],
            dims=(2, 2))

Let us find out the effect of the field direction after a certain a certain time . Note that the units of time are the inverse of the units of the Hamiltonian. For the purposes of this tutorial we leave all units dimensionless.

In [8]:

import numpy as np
import matplotlib.pyplot as plt

bar_width = 0.1
# We prepare an initial state ↑↓ (01).
# Note that Statevector and SparsePauliOp interpret the qubits from right to left
initial_state = Statevector.from_label("10")
trotter = TrotterQRTE()
final_time = 1.6
eps = 1e-5

# We create the list of angles in radians, with a small epsilon
# the exactly longitudinal field, which would present no dynamics at all
alphas = np.linspace(-np.pi / 2 + eps, np.pi / 2 - eps, 5)

for i, alpha in enumerate(alphas):
    H_alpha = get_hamiltonian(L=2, J=0.2, h=1.0, alpha=alpha)
    problem = TimeEvolutionProblem(H_alpha, initial_state=initial_state, time=1.6)
    result = trotter.evolve(problem)
    evolved_state = Statevector(result.evolved_state)
    # Dictionary of probabilities
    amplitudes_dict = evolved_state.probabilities_dict()
    labels = list(amplitudes_dict.keys())
    values = list(amplitudes_dict.values())
    # Convert angle to degrees
    alpha_str = f"$\\alpha={int(np.round(alpha * 180 / pi))}^\\circ$"
    plt.bar(np.arange(4) + i * bar_width, values, bar_width, label=alpha_str, alpha=0.7)

plt.xticks(np.arange(4) + 2 * bar_width, labels)
plt.xlabel("Measurement")
plt.ylabel("Probability")
plt.suptitle(
    f"Measurement probabilities at $t={final_time}$, for various field angles $\\alpha$\n"
    f"Initial state: 10, Linear lattice of size $L=2$"
)
plt.legend()
plt.show()

We have prepared a system initially with a sequence of spins , which corresponds to . After letting it evolve for under a transversal field (), we are almost guaranteed to measure , i.e. have a spin swap. (Note that the labels are interpreted from right to left). If the field is longitudinal (), we will have no evolution, therefore we will measure the system as it was initially prepared, . With intermediate angles, at , we will be able to measure all combinations will different probabilities, being a spin swap the most likely with a probability of 67%.

Auxiliary operators

Let us look into another feature of TrotterQRTE. We will now explore its ability to perform time evolutions of a system, while keeping track of some observables. The system that we now consider has a size of .

In [9]:

from math import pi

L = 6
H = get_hamiltonian(L=L, J=0.2, h=1.2, alpha=pi / 8)

The TrotterQRTE instance now will be created with a num_timesteps argument, and an Estimator primitive. The Qiskit's Estimator primitive estimates expectation values of observables, .

In [10]:

from qiskit_algorithms import TrotterQRTE
from qiskit.primitives import StatevectorEstimator

num_timesteps = 60
trotter = TrotterQRTE(num_timesteps=num_timesteps, estimator=StatevectorEstimator())

Let us define a magnetization operator , and a mean spin correlation operator .

In [11]:

magnetization_op = SparsePauliOp.from_sparse_list(
    [("Z", [i], 1.0) for i in range(0, L)], num_qubits=L
)
correlation_op = SparsePauliOp.from_sparse_list(
    [("ZZ", [i, i + 1], 1.0) for i in range(0, L - 1)], num_qubits=L
) / (L - 1)

Our new initial state will have the two middle spins facing down, and the TimeEvolutionProblem will incorporate some observables that will be kept track of:

the energy, i.e. the expectation value of the Hamiltonian,
the magnetization,
the mean spin correlation,

In [12]:

final_time = 30.0
initial_state = Statevector.from_label("001100")
problem = TimeEvolutionProblem(
    H,
    initial_state=initial_state,
    time=final_time,
    aux_operators=[H, magnetization_op, correlation_op],
)

Again, we let our TrotterQRTE evolve our newly created problem instance.

In [13]:

result = trotter.evolve(problem)

The new result also features an observables attribute. Let's extract the observables stored in it.

In [14]:

import numpy as np

observables = np.array(np.array(result.observables)[:, :, 0])
observables.shape

(61, 3)

In [15]:

import matplotlib.pyplot as plt

fig, axes = plt.subplots(3, sharex=True)
times = np.linspace(0, final_time, num_timesteps + 1)  # includes initial state
axes[0].plot(
    times, observables[:, 0], label="First order", marker="x", c="darkmagenta", ls="-", lw=0.8
)
axes[1].plot(
    times, observables[:, 1], label="First order", marker="x", c="darkmagenta", ls="-", lw=0.8
)
axes[2].plot(
    times, observables[:, 2], label="First order", marker="x", c="darkmagenta", ls="-", lw=0.8
)
axes[0].set_ylabel("Energy")
axes[1].set_ylabel("Magnetization")
axes[2].set_ylabel("Mean spin correlation")
axes[2].set_xlabel("Time")
fig.suptitle("Observable evolution")

Text(0.5, 0.98, 'Observable evolution')

Let us verify these results by comparing these expected values using Trotter evolution with the exact ones. To that end, we evaluate directly the expression discussed in the introduction,

on each one of the timesteps used by Trotter. We compute this exponential using SciPy's linalg.expm function, and then we let the initial system evolve using the Statevector's evolve method. This becomes too costly to be performed on larger systems very quickly.

In [16]:

import scipy as sc

H_array = H.to_matrix()

# We define a slightly denser time mesh
exact_times = np.linspace(0, final_time, 101)

# We compute the exact evolution using the exp
exact_evolution = [
    initial_state.evolve(sc.linalg.expm(-1j * time * H_array)) for time in exact_times
]

Having the exact state vectors, we compute the exact evolution of our operators' expectation values.

In [17]:

exact_energy = np.real([sv.expectation_value(H) for sv in exact_evolution])
exact_magnetization = np.real([sv.expectation_value(magnetization_op) for sv in exact_evolution])
exact_correlation = np.real([sv.expectation_value(correlation_op) for sv in exact_evolution])

We incorporate the exact evolution of the operators alongside the expected values resulting from the Trotterization.

In [18]:

axes[0].plot(exact_times, exact_energy, c="k", ls=":", label="Exact")
axes[1].plot(exact_times, exact_magnetization, c="k", ls=":", label="Exact")
axes[2].plot(exact_times, exact_correlation, c="k", ls=":", label="Exact")
# Select the labels of only the first axis
legend = fig.legend(
    *axes[0].get_legend_handles_labels(),
    bbox_to_anchor=(1.0, 0.5),
    loc="center left",
    framealpha=0.5,
)
fig.tight_layout()
fig

We see that, as an approximation, a Pauli-Trotter evolution isn't too far from the exact solution, but the accuracy is limited. Let's see how to find higher order formulas to address this.

Product formula overview

If it isn't specified, the default product formula that TrotterQRTE uses is the Lie product formula [2], which is at first order. In Qiskit this is implemented in the LieTrotter class. A first order formula consists of the approximation stated in the introduction, where the matrix exponential of a sum is approximated by a product of matrix exponentials:

Knowing this, we can have a look at the circuit that performs a single Trotter step.

In [19]:

from qiskit.circuit import QuantumCircuit
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import LieTrotter

dt = final_time / num_timesteps
trotter_step_first_order = PauliEvolutionGate(H, dt, synthesis=LieTrotter())

# We create an empty circuit
circuit = QuantumCircuit(H.num_qubits)
circuit.append(trotter_step_first_order, range(H.num_qubits))
circuit = circuit.decompose(reps=2)

# Let us print some stats
print(
    f"""
Trotter step with Lie-Trotter
-----------------------------

                  Depth: {circuit.depth()}
             Gate count: {len(circuit)}
    Nonlocal gate count: {circuit.num_nonlocal_gates()}
         Gate breakdown: {", ".join([f"{k.upper()}: {v}" for k, v in circuit.count_ops().items()])}
"""
)

# And finally draw the circuit
circuit.draw("mpl")

Trotter step with Lie-Trotter
-----------------------------

                  Depth: 17
             Gate count: 27
    Nonlocal gate count: 10
         Gate breakdown: CX: 10, U1: 6, R: 6, RZ: 5

There exists a second-order formula, called the Suzuki-Trotter decomposition [3], and can be used in Qiskit by means of the SuzukiTrotter class. Using this formula, a second order decomposition is:

By means of recursions, higher-order approximations can be found [1].

In [20]:

from qiskit.synthesis import SuzukiTrotter

second_order_formula = SuzukiTrotter()  # if not passed, order defaults to 2
trotter_step_second_order = PauliEvolutionGate(H, dt, synthesis=second_order_formula)
circuit = QuantumCircuit(H.num_qubits)
circuit.append(trotter_step_second_order, range(H.num_qubits))
circuit = circuit.decompose(reps=2)

# Let us print some stats
print(
    f"""
Trotter step with Suzuki Trotter (2nd order)
--------------------------------------------

                  Depth: {circuit.depth()}
             Gate count: {len(circuit)}
    Nonlocal gate count: {circuit.num_nonlocal_gates()}
         Gate breakdown: {", ".join([f"{k.upper()}: {v}" for k, v in circuit.count_ops().items()])}

"""
)

# And finally
circuit.draw("mpl")

Trotter step with Suzuki Trotter (2nd order)
--------------------------------------------

                  Depth: 34
             Gate count: 53
    Nonlocal gate count: 20
         Gate breakdown: CX: 20, U1: 12, R: 11, RZ: 10

In [21]:

fourth_order_formula = SuzukiTrotter(order=4)
trotter_step_fourth_order = PauliEvolutionGate(H, dt, synthesis=fourth_order_formula)
circuit = QuantumCircuit(H.num_qubits)
circuit.append(trotter_step_fourth_order, range(H.num_qubits))
circuit = circuit.decompose(reps=2)

# Let us print some stats
print(
    f"""
Trotter step with Suzuki Trotter (4th order)
--------------------------------------------

                  Depth: {circuit.depth()}
             Gate count: {len(circuit)}
    Nonlocal gate count: {circuit.num_nonlocal_gates()}
         Gate breakdown: {", ".join([f"{k.upper()}: {v}" for k, v in circuit.count_ops().items()])}

"""
)

Trotter step with Suzuki Trotter (4th order)
--------------------------------------------

                  Depth: 170
             Gate count: 265
    Nonlocal gate count: 100
         Gate breakdown: CX: 100, U1: 60, R: 55, RZ: 50

Finally, let us perform a simulation at fourth order.

In [22]:

from qiskit.synthesis import SuzukiTrotter

trotter = TrotterQRTE(
    SuzukiTrotter(order=4), num_timesteps=num_timesteps, estimator=StatevectorEstimator()
)
problem = TimeEvolutionProblem(
    H,
    initial_state=initial_state,
    time=final_time,
    aux_operators=[H, magnetization_op, correlation_op],
)
result = trotter.evolve(problem)
observables_order4 = np.array(np.array(result.observables)[:, :, 0], dtype=np.float64)

and plot the results.

In [23]:

axes[0].plot(
    times, observables_order4[:, 0], label="Fourth Order", marker="x", c="limegreen", ls="-", lw=0.8
)
axes[1].plot(
    times, observables_order4[:, 1], label="Fourth Order", marker="x", c="limegreen", ls="-", lw=0.8
)
axes[2].plot(
    times, observables_order4[:, 2], label="Fourth Order", marker="x", c="limegreen", ls="-", lw=0.8
)

# Replace the legend
legend.remove()
legend = fig.legend(
    *axes[0].get_legend_handles_labels(),
    bbox_to_anchor=(1.0, 0.5),
    loc="center left",
    framealpha=0.5,
)
fig

As it is to expect, we can directly see that a higher-order product formula leads to more accurate expectation values.

Magnetization evolution plot

In this last step, we aim at visualizing the time evolution of each one of the sites of the lattice individually. Let us present the expected value of the magnetization of each one of the sites as a function of time in a color plot. As the initial state was (spin sequence ), at time we expect to see magnetization in the two middle signs and elsewhere. During the evolution, as the plot above suggests, we can expect to see oscillations which will eventually damp.

In [24]:

from matplotlib import cm

# An inner list comprehension loops over the terms of the SparsePauliOp magnetization_op,
# which corresponds to the magnetization of each one of the sites
magnetizations = np.real(
    [[sv.expectation_value(term) for term in magnetization_op] for sv in exact_evolution]
)
# The shape of magnetizations is (101, 6), containing <Z>(t) for each site 0, 1, ..., 5
plt.figure(figsize=(14, 2))
# Create the 2-dim xx and yy arrays tiling the grid with the x and y values
xx, yy = np.meshgrid(exact_times, np.arange(L))
plt.pcolor(xx, yy, magnetizations.T, vmin=-1, vmax=+1, cmap="RdBu")
# Force the figure to have all y ticks from 0 to 5
plt.yticks(np.arange(L))
plt.ylabel("Site $i$")
plt.xlabel("Time")
plt.colorbar(label="$\\langle Z_i \\rangle$", aspect=1.8)
plt.show()

References

[1] Hatano, Naomichi, and Masuo Suzuki. "Finding exponential product formulas of higher orders." Quantum annealing and other optimization methods. Berlin, Heidelberg: Springer Berlin Heidelberg, 2005. 37-68.

[2] Varadarajan, Veeravalli S. Lie groups, Lie algebras, and their representations. Vol. 102. Springer Science & Business Media, 2013.

[3] Magnus, Wilhelm. "On the exponential solution of differential equations for a linear operator." Communications on pure and applied mathematics 7.4 (1954): 649-673.

In [25]:

import tutorial_magics

%qiskit_version_table
%qiskit_copyright

Version Information

Software	Version
`qiskit`	2.0.2
`qiskit_algorithms`	0.4.0
System information
Python version	3.13.3
OS	Linux
Tue Jun 17 10:45:23 2025 CEST

This code is a part of a Qiskit project

This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.

Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.

Join the Discussion

Comments (0)

No comments yet. Be the first to share your thoughts!

Indexed by QCR Librarian

This entry was created automatically from publicly available records. QCR links to public sources and only stores repository content where the license permits redistribution.

Claim this entry →

Versions

v1 Latest

Jun 17, 2026

qcr:2606.69981.1

Cite all versions? Use the base QCR ID to always reference the latest version of this entry.

Tools used

Qiskit

Keywords

trotterization

time-evolution

hamiltonian-simulation

qiskit

dynamics