Symbolic Heisenberg Evolution

In the Heisenberg picture of quantum mechanics, operators evolve in time while states remain fixed. Let's see how to simulate this evolution in Aleph for a simple spin chain.

The Heisenberg equation of motion

For an operator $A(t)$ evolving under a Hamiltonian $H$, the equation of motion is:

$$\frac{dA}{dt} = i[H, A(t)]$$

where $[H, A] = HA - AH$ is the commutator. We'll solve this numerically using Euler's method and see how to manage the growth of operator expressions.

Building the Hamiltonian

We consider the following Hamiltonian for the transverse-field Ising model. For $N$ qubits, we have

$$H = J \sum_{j=0}^{N-1} Z_j Z_{j+1} + h \sum_{i=0}^{N-1} X_i$$

where $J$ is the coupling strength and $h$ is the transverse field strength. We can build this in Aleph!

var N = 2;              // number of qubits
var J = 1.0;            // coupling strength
var h = 0.5;            // transverse field

var Hmt = operator_sum(); // empty ComplexOperatorSum

// ZZ interaction terms
for(j : [0..N-1]) {
    var k_site = (j + 1) % N;
    Hmt += J * (Z(j) * Z((j+1)%N));
}

// X field terms
for (var m = 0; m < N; ++m) {
    Hmt += h * X(m);
}

// Simplify the Hamiltonian
simplify_paulis(Hmt);

The simplify_paulis function is a specialized transformation for Pauli operators that combines Pauli operators according to their algebra rules, merges terms, and prunes the sum according to some tolerance. This will be the main driver for efficiently computing the Pauli expressions.

Setting up the initial operator

We start with a simple observable, the $Z$ operator on the first qubit:

var A = operator_sum(Z(0)); // ComplexOperatorSum

Time evolution with Euler's method

To solve the Heisenberg equations of motion for the observable $A$, we will use a Euler's method which approximates the time derivative with a finite difference:

$$A(t + \Delta t) \approx A(t) + i\Delta t \cdot [H, A(t)]$$

Let's implement this step by step.

var Nt = 5;            // number of time steps
var dt = 0.1;           // time step size
var epsilon = 1e-3;     // pruning threshold
var rhs = operator_sum() 

for (step : [0..Nt]) {
    // Compute the commutator [H, A(t)]
    rhs = operator_sum(commutator(Hmt,A));
    
    // Expand the products into sums of products. 
    rhs = flattened(rhs);
    
    // Euler step: A(t+dt) = A(t) + i*dt*[H,A(t)]
    A = A + (dt * 1i )* rhs;
    
    // Clean up the new operator, and prune terms with coefficients < epsilon
    simplify_paulis(A,epsilon);
    
    print(A);
}

For each step, we can see the result

((1+0i) * Z(0) + (0.1+0i) * Y(0))
((0.99+0i) * Z(0) + (0.2+0i) * Y(0) + (0.02+0i) * X(0) * Z(1))
((0.97+0i) * Z(0) + (0.295+0i) * Y(0) + (0.06+0i) * X(0) * Z(1) + (0.002+0i) * X(0) * Y(1))
((0.9405+0i) * Z(0) + (0.38+0i) * Y(0) + (0.1188+0i) * X(0) * Z(1) + (0.008+0i) * X(0) * Y(1))
((0.9025+0i) * Z(0) + (0.45029+0i) * Y(0) + (0.194+0i) * X(0) * Z(1) + (0.01988+0i) * X(0) * Y(1))

Notice how the operator evolves:

• Step 0: A Y(0) component appears due to the transverse field h*X(0)
• Step 1: The coupling J*Z(0)*Z(1) creates the two-qubit term X(0)*Z(1)
• Later steps: More complex terms emerge, showing operator spreading across the chain
• Coefficients: The original Z(0) weight decreases as the operator delocalizes

Without transformations, the operator expression grows exponentially with each time step. Each commutator doubles the number of terms, and nested structures accumulate. The transformation pipeline manages this growth:

• flatten expands nested products into a flat sum
• simplify_paulis applies Pauli algebra rules ($X^2 = I$, $XY = iZ$, etc.), merges like terms, and truncates the sum based on the coefficient magnitudes.

Scaling to larger systems: While our 2-qubit, 5-step example remains manageable, the computational cost grows rapidly. For $N=10$ qubits and 50 time steps, the expression can grow to millions of terms without pruning. The pruning threshold $\varepsilon$ becomes critical for controlling this growth - smaller values preserve accuracy but increase memory usage, while larger values keep expressions compact at the cost of approximation error.

Understanding the results

After evolution, A contains the time-evolved operator $A(t)$ as a sum of Pauli strings with complex coefficients. You can:

• Analyze support: See which qubits the operator acts on (it spreads from qubit 0)
• Compute observables: Use the evolved operator to calculate expectation values

Benchmarking performance

Now that we've established all the pieces, we can combine them into a single function. This lets us easily test how performance scales with the number of qubits $N$.

By creating a single function that takes N as input, we can benchmark different system sizes:

var my_Heisenberg = fun(integer N)
{
    // Hamiltonian parameters
    var J = 1.0; 
    var h = 0.5;
    // Time evolution parameters
    var Nt = 30;
    var dt = 0.1;
    var epsilon = 1e-3;
    // Build Hamiltonian
    var Hmt = operator_sum();
    // ZZ interaction terms
    for (j : [0..N-1]) {
        Hmt += J * (Z(j) * Z((j + 1) % N));
    }
    // X field terms
    for (m : [0..N]) {
        Hmt += h * X(m);
    }
    // Setup initial operator
    var A = operator_sum(Z(0));
    var rhs = operator_sum();
    // Time evolution loop
    for (step : [0..Nt]) {
        rhs = operator_sum(commutator(Hmt,A));
        rhs = flattened(rhs);
        A = A + (dt * 1i )* rhs;
        simplify_paulis(A,epsilon);
    }
}

To run benchmark for different system sizes, we can build the following for-loop

// Run benchmark for various N
for (N : [16, 18, 20, 22, 24]) {
    var results = benchmark(bind(my_Heisenberg, N), 10);
    print(results.string());
}

This runs the full evolution 10 times for each system size (16, 18, 20, 22 and 24 qubits) and reports timing statistics. You'll see how the computational cost increases as you add more qubits and lower $\epsilon$.

BenchmarkResult(total=675.670 ms, mean=67.567 ms, best=62.649 ms, worst=74.341 ms, stddev=2.912 ms) 
BenchmarkResult(total=759.219 ms, mean=75.922 ms, best=71.885 ms, worst=91.282 ms, stddev=5.631 ms) 
BenchmarkResult(total=832.948 ms, mean=83.295 ms, best=79.803 ms, worst=90.024 ms, stddev=3.657 ms) 
BenchmarkResult(total=945.349 ms, mean=94.535 ms, best=88.209 ms, worst=102.571 ms, stddev=4.826 ms) 
BenchmarkResult(total=1.011 s, mean=101.086 ms, best=96.198 ms, worst=119.028 ms, stddev=6.316 ms) 

Computing expectation values

Now that we have the time-evolved operator $A(t)$, we now may want to compute expectation values $\langle \psi | A(t) | \psi \rangle$ for a given quantum state $|\psi\rangle$.

There are two equivalent ways to compute expectation values:

1. Schrödinger picture: The state evolves in time while operators remain fixed $\langle A(t) \rangle = \langle \psi(t) | A | \psi(t) \rangle$

2. Heisenberg picture: Operators evolve in time while the state remains fixed $\langle A(t) \rangle = \langle \psi(0) | A(t) | \psi(0) \rangle$

We'll implement both methods and compare their results to verify our Heisenberg evolution, tracking how observables evolve in many-body systems.

info

Aleph provides a built-in expval(A, state) function that computes expectation values $\langle \psi | A | \psi \rangle$ efficiently. See the expval for details.

Schrödinger picture: Evolving the state

In the Schrödinger picture, we evolve the initial state $|\psi(0)\rangle$ forward in time using the unitary operator $U(t) = e^{-iHt}$, then compute expectation values with the fixed operator $A(0)$:

$$|\psi(t)\rangle = e^{-iHt}|\psi(0)\rangle$$ $$\langle A(t) \rangle = \langle \psi(t) | A(0) | \psi(t) \rangle$$

First we choose an initial state to evolve, a simple (and common) choice is the all zero computational basis state $|\psi\rangle = |00...0\rangle$. We can then make our $2^N$ dimensional many body state from this basis state:

// Create empty state vector container
var state_0 = zeros([2**N], as_matrix | as_complex); 
var q = Qbit(N); // |00...0⟩ → index 0 
state_0[q] = 1.0; // Set amplitude

We'll use exact diagonalization to compute the time evolution by obtaining obtaining computing the matrix exponential directly:

def schrodingerPicture(N, J, h, dt, Nt, A_0, Hmt, state_0)
{
    // Diagonalize the Hamiltonian H = V D V†
    var Hmat = matrix(Hmt);
    var es = eigensolver(Hmat, ["is_hermitian" : true]);
    var evals = es.eigenvalues();      // Eigenvalues 
    var evecs = es.eigenvectors();     // Eigenvectors

    var expecvals = [];
    
    // Time evolution loop
    for (var step = 0; step <= Nt; ++step) {
        var t = dt * step;

        // Compute e^(-iHt) 
        var expdt_arr = exp(-1i * t * array(evals));
        var expdt_mat = matrix(expdt_arr);
        
        // Evolve the state |ψ(t)⟩ = e^(-iHt)|ψ(0)⟩
        var state_t = evecs * as_diagonal(expdt_mat) * evecs.adjointed() * state_0;
        
        // Compute expectation value ⟨ψ(t)|A_0|ψ(t)⟩
        var expecval_t = expval(A_0, state_t);
        
        expecvals.push_back(expecval_t.real());
    }

    return expecvals;
}

Using the Schrödinger picture function:

We can call this function with our parameters to compute the time evolution of expectation values:

// Use the parameters and objects we defined earlier
var expecvals_S = schrodingerPicture(N, J, h, dt, Nt, A_0, Hmt, state_0);

Heisenberg picture: Evolving the operator

In the Heisenberg picture, the state $|\psi(0)\rangle$ remains fixed while the operator evolves according to the Heisenberg equation of motion. We compute the expectation value using the evolved operator and the initial state:

$$\langle A(t) \rangle = \langle \psi(0) | A(t) | \psi(0) \rangle$$

In our above Euler's method evolution, we already have access to the evolved operator, all we need to do is specify our initial state. Since Qbits are in terms of orthonormal basis state, we need only check if our operator expressions leave our state unchanged after being acted upon:

// Symbolic expectation value
def symbolic_expval(ComplexOperatorSum S, Qbit state) {
    var sum = complex(0.0)
    for (var k = 0; k < S.get_coefficients().size(); k += 1) {
        var term = S [k]
        var out = term.second() * state
        if (integer(out) == integer(state)) {
            sum += (term.first)
        }
    }
    return sum
}

note

The symbolic_expval function assumes that the input state has coefficient 1.0 and is complex. If instead the coefficient can be anything, the addition to the sum should be changed to (term.first * state.coefficient()).

After the Heisenberg evolution has been performed and the operator $A(t)$ is obtained, the expectation value can calculated with

var basis_state = Qbit(N)
var exp_val = symbolic_expval(A,basis_state)

Higher Order Solvers: The RK4 Method

While we used Euler's method earlier for simplicity, the Runge-Kutta 4th order (RK4) method offers significantly better accuracy and stability, particularly for longer time evolutions [1, 2]. The RK4 method requires the derivative $i[H, A]$ at four intermediate points:

$$k_1 = i\Delta t[H, A(t)]$$ $$k_2 = i\Delta t[H, A(t) + \frac{k_1}{2}]$$ $$k_3 = i\Delta t[H, A(t) + \frac{k_2}{2}]$$ $$k_4 = i\Delta t[H, A(t) + k_3]$$

which are used to construct the full step, combining them with a weighted average:

$$A(t + \Delta t) = A(t) + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)$$

A single RK4 step can be performed with:

def evolve_RK4_step(ComplexOperatorSum Hmt, ComplexOperatorSum A, integer N, real dt, real epsilon)
{
    // k1 = i*dt*[H, A]
    var k1 = 1.0i * dt * flattened(operator_sum(commutator(H, A)));
    simplify_paulis(k1,eps)

    // k2 = i*dt*[H, A + k1/2]
    var k2 = 1.0i * dt * flattened(operator_sum(commutator(H, A + 0.5 * k1)));
    simplify_paulis(k2,eps)

    // k3 = i*dt*[H, A + k2/2]
    var k3 = 1.0i * dt * flattened(operator_sum(commutator(H, A + 0.5 * k2)));
    simplify_paulis(k3,eps)

    // k4 = i*dt*[H, A + k3]
    var k4 = 1.0i * dt * flattened(operator_sum(commutator(H, A + k3)));
    simplify_paulis(k4,eps)

    // A_new = A + (k1 + 2*k2 + 2*k3 + k4)/6.0
    A = A + (1.0/6.0)*k1 + (2.0/6.0)*k2 + (2.0/6.0)*k3 + (1.0/6.0)*k4
    simplify_paulis(A,epsilon)

    return A;
}

The complete time evolution is then just looping over this step and saving the expectation value at each desired time.

def heisenbergPicture(N, J, h, dt, Nt, A, Hmt, state_0, epsilon)
{
    var expecvals = [];
    
    for (var step = 0; step <= Nt; ++step) {

        // Compute expectation value ⟨ψ(0)|A(t)|ψ(0)⟩
        var expecval_t = symbolic_expval(A, state_0);
        expecvals.push_back(expecval_t.real());
        
        // Evolve operator to next time step
        if (step < Nt) {
            A_t = evolve_RK4_step(Hmt, A, N, dt, epsilon);
        }
    }

    return expecvals;
}

Now we can call this function to compute the time evolution of expectation values using operator evolution:

var expecvals_H = heisenbergPicture(N, J, h, dt, Nt, A_0, Hmt, state_0, epsilon);

Key takeaways

1. Two equivalent pictures: Schrödinger (evolving states) and Heisenberg (evolving operators) give identical physical predictions for expectation values
2. Built-in expval function: Aleph provides expval(A, state) to compute $\langle \psi | A | \psi \rangle$. This works for dense states and product states (Qbit)
3. Schrödinger picture uses exact diagonalization: Compute $H = VDV^\dagger$ once, then evolve states with $e^{-iHt}$
   • exact but limited to small systems
4. Heisenberg picture uses RK4 integration: Evolve operators via $\frac{dA}{dt} = i[H, A(t)]$
   • approximate but offers flexibility for larger systems through expression simplification and pruning
5. RK4 over Euler time integration: Fourth-order method achieves $O(\Delta t^4)$ global error compared to Euler's $O(\Delta t)$, while maintaining numerical stability over longer time evolutions
6. Pruning is essential: The epsilon parameter controls operator expression growth - smaller values preserve accuracy but increase computational cost
7. Validation through comparison: Comparing both methods verifies implementation correctness and reveals numerical accuracy limits

Recommended next steps

• Try different Hamiltonians (Ising model, $XY$ model)
• Explore adaptive pruning strategies
• Export expectation values and visualize the comparison between both methods (see Data Files)

References

[1] Butcher, J.C. (2016). Numerical Methods for Ordinary Differential Equations (3rd ed.) . John Wiley & Sons. DOI:10.1002/97811191215344

[2] Griffiths, D.F., & Higham, D.J. (2010). *Numerical Methods for Ordinary Differential Equations: Initial Value Problems. Springer. DOI:10.1007/978-0-85729-148-6

Contributors

Eunji Yoo

Sebastien J. Avakian

Jonathon Riddell