# 2. Quenching a Random Product State¶

First we’ll set up a 20 qubit random product state, and a sparse heisenberg hamiltonian to evolve it with:

import quimb as qu

n = 20
H = qu.ham_heis(n, sparse=True)
psi0 = qu.rand_product_state(n)


We can do a few checks on the system like $$\langle \psi | \psi \rangle$$:

# check normalization
qu.expectation(psi0, psi0)

0.9999999999999987


and $$\langle \psi | H | \psi \rangle$$:

# find the initial energy
qu.expec(H, psi0)

0.3002436137794859


or $$\langle \psi | H^2 | \psi \rangle$$:

# find the initial variance in energy
Hpsi = H @ psi0
qu.expec(Hpsi, Hpsi)

11.30918741947443


Let’s compare this to the total energy spectrum of the Hamiltonian H:

%%time

en_low, en_high = qu.bound_spectrum(H)

print("Highest energy:", en_high)
print("Lowest energy:", en_low)

Highest energy: 4.749999999999985
Lowest energy: -8.682473334398967
CPU times: user 2min 48s, sys: 2min 16s, total: 5min 5s
Wall time: 45.9 s


From which we can infer that our initial state has roughly has overlap between many states in the centre of the spectrum.

## 2.1. Evolution¶

Now let’s set up some things we want to compute while we evolve with the hamiltonian. Namely the logarithmic negativity and mutual information between neighbouring qubits. Since we’ll use the default adaptive integrating scheme, we’ll return the current time as well:

def compute(t, pt):
"""Perform computation at time t with state pt.
"""
dims = [2] * n
lns = [qu.logneg_subsys(pt, dims, i, i + 1) for i in range(n - 1)]
mis = [qu.mutinf_subsys(pt, dims, i, i + 1) for i in range(n - 1)]
return t, lns, mis


Set up the evolution with the initial state, hamiltonian and the compute dict:

evo = qu.Evolution(psi0, H, compute=compute, progbar=True)


Update the evolution to t=5. The functions in compute will be called at each step the integrator uses. If we had set method='solve' or method='expm', we should use the generator evo.at_times(ts) to specify the time steps when computation takes place.

evo.update_to(5)

100%|##########| 100/100 [00:37<00:00,  2.66%/s]


We can extract the results of the computation from evo.results and plot them:

%matplotlib inline
import matplotlib.pyplot as plt

ts, lns, mis = zip(*evo.results)

with plt.style.context(qu.NEUTRAL_STYLE):
fig, axs = plt.subplots(2, 1, sharex=True)
axs[0].plot(ts, lns, '-');
axs[0].set_title("Logarithmic Negativity")
axs[1].plot(ts, mis, '-');
axs[1].set_title("Mutual Information")


We can see that the classical correlations outlast the quantum correlations.

Finally, let’s check that energy has been conserved in the current state at t=5:

qu.expec(H, evo.pt)

0.3002436962453416