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 = 18
H = qu.ham_heis(n, sparse=True).real
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)

and \(\langle \psi | H | \psi \rangle\):

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

or \(\langle \psi | H^2 | \psi \rangle\):

# find the initial variance in energy
psi0.H @ H @ H @ psi0

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


en_low, en_high = qu.bound_spectrum(H)

print("Highest energy:", en_high)
print("Lowest energy:", en_low)
Highest energy: 4.500000000000008
Lowest energy: -8.022749087033766
CPU times: user 313 ms, sys: 249 ms, total: 562 ms
Wall time: 5.25 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.

100%|██████████| 100/100 [00:10<00:00, 10.18%/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)

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)