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:



In [1]:

    
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$:



In [2]:

    
# check normalization
qu.expectation(psi0, psi0)









    Out[2]:





0.9999999999999982

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



In [3]:

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









    Out[3]:





-0.452591067342989

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



In [4]:

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









    Out[4]:





qarray([[3.521528-7.077672e-16j]])

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



In [5]:

    
%%time

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.

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:



In [6]:

    
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:



In [7]:

    
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.



In [8]:

    
evo.update_to(5)









    



100%|██████████| 100/100 [00:10<00:00, 10.18%/s]

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



In [9]:

    
%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")

plt.show()

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:



In [10]:

    
qu.expec(H, evo.pt)









    Out[10]:





-0.45259091947115565