Computing propagators

Sometime the evolution of a single state is not sufficient and the full propagator is desired. QuTiP has the propagator function to compute them:

>>> H = sigmaz() + np.pi *sigmax()
>>> psi_t = sesolve(H, basis(2, 1), [0, 0.5, 1]).states
>>> prop = propagator(H, [0, 0.5, 1])

>>> print((psi_t[1] - prop[1] @ basis(2, 1)).norm())
2.455965272327082e-06

>>> print((psi_t[2] - prop[2] @ basis(2, 1)).norm())
2.0071900004562142e-06

The first argument is the Hamiltonian, any time dependent system format is accepted. The function also accepts an optional c_ops argument for collapse operators. When used, a propagator for density matrices is computed: \(\rho(t) = U(t)(\rho(0))\):

>>> rho_t = mesolve(H, fock_dm(2, 1), [0, 0.5, 1], c_ops=[sigmam()]).states
>>> prop = propagator(H, [0, 0.5, 1], c_ops=[sigmam()])

>>> print((rho_t[1] - prop[1](fock_dm(2, 1))).norm())
7.23009476734681e-07

>>> print((rho_t[2] - prop[2](fock_dm(2, 1))).norm())
1.2666967766644768e-06

The propagator function is also available as a class:

>>> U = Propagator(H, c_ops=[sigmam()])

>>> state_0_5 = U(0.5)(fock_dm(2, 1))
>>> state_1 = U(1., t_start=0.5)(state_0_5)

>>> print((rho_t[1] - state_0_5).norm())
7.23009476734681e-07

>>> print((rho_t[2] - state_1).norm())
8.355518501351504e-07

The Propagator can take options and args as a solver instance.

Using a solver to compute a propagator

Many solvers accept an operator as the initial state. When an identity matrix is passed as the initial state, the propagator is computed. This can be used to compute a propagator for Bloch-Redfield or Floquet equations:

>>> delta = 0.2 * 2*np.pi
>>> eps0 = 1.0 * 2*np.pi
>>> gamma1 = 0.5

>>> H = - delta/2.0 * sigmax() - eps0/2.0 * sigmaz()

>>> def ohmic_spectrum(w):
>>>   if w == 0.0: # dephasing inducing noise
>>>     return gamma1
>>>   else: # relaxation inducing noise
>>>     return gamma1 / 2 * (w / (2 * np.pi)) * (w > 0.0)

>>> prop = brmesolve(H, qeye(2), [0, 1], a_ops=[[sigmax(), ohmic_spectrum]]).final_state