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'])

Or we can use the final state:

fidelity(psi0, evo.pt)

0.003302180752068547