6. Time Evolution#

Time evolutions in quimb are handled by the class Evolution, which is initialized with a starting state and hamiltonian.

6.1. Basic Usage#

Set up the Evolution object with a initial state and hamiltonian.

from quimb import *

p0 = rand_ket(2**10)
h = ham_heis(10, sparse=True)
evo = Evolution(p0, h)

Update it in a single shot to a new time and get the state,

evo.update_to(1)
evo.pt
[[0.006074+0.011257j]
 [0.028495+0.016557j]
 [0.01617 -0.015705j]
 ...
 [0.018271+0.026071j]
 [0.018582-0.016482j]
 [0.027511+0.048377j]]

Lazily generate the state at multiple times:

for pt in evo.at_times([2, 3, 4]):
    print(expec(pt, p0))
0.003244693043705576
0.008360765055266119
0.009413166876677937

6.2. Methods of Updating#

There are three methods of updating the state:

  • Evolution(..., method='integrate'): use definite integration. Get system at each time step, only need action of Hamiltonian on state. Generally efficient. For pure and mixed states. The additional option int_small_step={False, True} determines whether a low or high order adaptive stepping scheme is used, giving naturally smaller or larger times steps. See scipy.integrate.ode for details, False corresponds to "dop853", True to "dopri5".

  • Evolution(..., method='solve'). Diagonalize the hamiltonian, which once done, allows quickly updating to arbitrary times. Supports pure and mixed states, recomended for small systems.

  • Evolution(..., method='expm'): compute the evolved state using the action of the matrix exponential in a ‘single shot’ style. Only needs action of Hamiltonian, for very large systems can use distributed MPI. Only for pure states.

6.3. Computing on the fly#

Sometimes, if integrating, it is best to just query the state at time-steps chosen dynamically by the adaptive scheme. This is achieved using the compute keyword supplied to Evolution. It can also just be a convenient way to set up calculations as well:

p0 = rand_product_state(10)
h = ham_heis(10, sparse=True)

dims = [2] * 10
sysa, sysb = (0, 1), (2, 3)

def calc_t_and_logneg(t, pt):
    ln = logneg_subsys(pt, dims, sysa, sysb)
    return t, ln

evo = Evolution(p0, h, compute=calc_t_and_logneg, progbar=True)
evo.update_to(1)

ts, lns = zip(*evo.results)
100%|##########| 100/100 [00:00<00:00, 20144.58%/s]
ts
(0.0,
 0.2494153162899183,
 0.4631400312076022,
 0.6768647461252861,
 0.9475557600771725,
 1.0)
lns
(0.0,
 0.2589846231770502,
 0.4423181326621942,
 0.592697104829058,
 0.7379099150663264,
 0.7601643146969251)

If a dict of callables is supplied to compute, (each should take two arguments, the time, and the state, as above), Evolution.results will itself be a dictionary containing the results of each function at each time step, under the respective key. This can be more convenient:

def calc_t(t, _):
    return t

def calc_logneg(_, pt):
    return logneg_subsys(pt, [2]*10, 0, 1)

evo = Evolution(p0, h, compute={'t': calc_t, 'ln': calc_logneg}, progbar=True)
evo.update_to(1)
100%|##########| 100/100 [00:00<00:00, 18631.41%/s]
evo.results
{'t': [0.0,
  0.2494153162899183,
  0.4631400312076022,
  0.6768647461252861,
  0.9475557600771725,
  1.0],
 'ln': [0.0,
  0.150158101589196,
  0.23732495484813082,
  0.2922696608578911,
  0.3242059171688551,
  0.3263498611123527]}

6.4. Time-Dependent Evolutions#

If you are using method='integrate' you can supply a callable to ham to evolve the state with a time dependent Hamiltonian. It should take a single argument t and return the Hamiltonian at the time. It probably makes sense to use a custom class here to avoid reconstructing as much of the Hamiltonian as possible at each step.

Here we’ll evolve the Neel state:

\[ | \psi(0) \rangle = | \uparrow \downarrow \uparrow \downarrow \uparrow \ldots \rangle \]

with the Hamiltonian:

\[ H(t) = \sum_{i = 0}^{L - 1} S^Z_{i} S^Z_{i + 1} + \cos(t) \sum_{i}^{L} S^X_i \]
class MyTimeDepIsingHam:
    
    def __init__(self, L):
        self.h_interaction = ham_ising(L, sparse=True, jz=1.0, bx=0.0, cyclic=False)
        self.h_field = ham_ising(L, sparse=True, jz=0.0, bx=1.0, cyclic=False)
    
    def __call__(self, t):
        return self.h_interaction + cos(t) * self.h_field
L = 16

# our initial state
psi0 = neel_state(L)

# instantiate the ham object, it's __call__ method will be used by Evolution
fn_ham_t = MyTimeDepIsingHam(L)

We still want to compute some properties during the evolution:

compute = {
    'time': lambda t, p: t,
    'entropy': lambda t, p: entropy_subsys(p, dims=[2] * L, sysa=range(L // 2))
}

Now we set up the evolution object again:

evo = Evolution(psi0, fn_ham_t, progbar=True, compute=compute)
evo.update_to(10)
100%|##########| 100/100 [00:07<00:00, 12.53%/s]

We can plot the half chain entropy that we computed on the fly:

%matplotlib inline

import matplotlib.pyplot as plt

with plt.style.context(NEUTRAL_STYLE):
    plt.plot(evo.results['time'], evo.results['entropy'])
_images/3f4a07fff8292fe0b3911ac5af80a01c3c36e2f1fe85d36a58cf43e898ce279d.svg

Or we can use the final state:

fidelity(psi0, evo.pt)
0.003302180752068547