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: ρ(t) = U(t)(ρ(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.
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