Simulate Hamiltonian Time Evolution

A key task in quantum simulation is calculating the time evolution of a system under a given Hamiltonian, H. pyrauli is exceptionally well-suited for this task when the Hamiltonian can be expressed as a single, multi-qubit Pauli string (e.g., H = XZIX), which is common in physics models and algorithm building blocks.

This is because the evolution operator, \(U(t) = e^{-iHt}\), can be applied directly to the observable in the Heisenberg picture without needing to be decomposed into a complex sequence of one- and two-qubit gates.

The eiht Gate: Your Primary Tool

The primary method for Hamiltonian evolution is eiht(). It directly applies the evolution operator for a given time t.

circuit.eiht(pauli_axis, t)

pauli_axis: A list of strings defining the Pauli operators in the Hamiltonian term (e.g., ["X", "Z", "I", "X"]).
t: A float representing the evolution time.

The following example simulates an 8-qubit system evolving under a 4-local Wen Plaquette operator, \(H = X_0 Z_1 X_4 Z_5\). We then measure the expectation value of the \(Z_0\) observable.

n_qubits = 8
evolution_time = 0.5

# Define the Hamiltonian axis for H = X_0 Z_1 X_4 Z_5
hamiltonian_axis = ["X", "Z", "I", "I", "X", "Z", "I", "I"]

# 1. Build the circuit
circuit = pyrauli.Circuit(n_qubits)
circuit.eiht(hamiltonian_axis, evolution_time)

# 2. Define an observable to measure
observable = pyrauli.Observable("ZIIIIIII") # Z on qubit 0

# 3. Run the simulation and get the result
ev, err = circuit.expectation_value(observable)

print(f"Expectation value of Z_0 at t={evolution_time}: {ev:.4f}")

# Verify the result (for this case, EV = cos(2*t))
expected = math.cos(2 * evolution_time)
print(f"Expected theoretical value: {expected:.4f}")

The Rp Gate vs. eiht

pyrauli also provides the rp() gate, which applies a generalized rotation around a Pauli axis. It is important to understand the distinction between the two:

eiht() implements the standard physics definition of time evolution: \(U = e^{-i P t}\).
rp() implements the standard quantum gate rotation: \(U = e^{-i P \theta / 2}\).

This leads to a simple relationship between their parameters: a time evolution for time t is equivalent to a rotation by an angle theta = 2*t.

Tip

`circuit.eiht(axis, t)` is equivalent to `circuit.rp(axis, 2 * t)`.

The following code provides a concrete demonstration of this equivalence.

n_qubits = 4
time = 0.25
axis = ["X", "Y", "Z", "I"]

# Create one circuit using eiht
circuit_eiht = pyrauli.Circuit(n_qubits)
circuit_eiht.eiht(axis, time)

# Create another circuit using rp with theta = 2*t
circuit_rp = pyrauli.Circuit(n_qubits)
circuit_rp.rp(axis, 2 * time)

# They should produce the exact same result
observable = pyrauli.Observable("ZIII")
ev_eiht, _ = circuit_eiht.expectation_value(observable)
ev_rp, _ = circuit_rp.expectation_value(observable)

print(f"Result from eiht(t={time}): {ev_eiht:.4f}")
print(f"Result from rp(theta={2*time}):   {ev_rp:.4f}")
assert ev_eiht == pytest.approx(ev_rp)

Symbolic Time Evolution

A powerful feature of pyrauli is that the time parameter t can be a symbolic variable. This allows you to run the simulation just once to obtain a symbolic expression for an expectation value as a function of time. You can then evaluate this expression for any number of time points without re-running the simulation.

This is ideal for analyzing the dynamics of a system or for algorithms where time is a parameter to be optimized.

n_qubits = 4
axis = ["X", "Y", "Z", "I"]

# 1. Build a SymbolicCircuit with a variable 't' for time
symbolic_circuit = pyrauli.SymbolicCircuit(n_qubits)
symbolic_circuit.eiht(axis, pyrauli.SymbolicCoefficient("t"))

# 2. Run the simulation on a standard observable
observable = pyrauli.SymbolicObservable("ZIII")
final_observable = symbolic_circuit.run(observable)

# 3. The final expectation value is a symbolic expression of 't'
symbolic_ev = final_observable.expectation_value()
print(f"Symbolic expectation value: {symbolic_ev}")

# 4. Evaluate the expression for different concrete times
times_to_check = [0.0, 0.5, 1.0]
for t_val in times_to_check:
    numeric_ev = symbolic_ev.evaluate({"t": t_val})
    print(f"EV at t={t_val}: {numeric_ev:.4f}")