Lindblad equation implementation using QuTiP.
QuTiP's mesolve() integrates the full Lindblad master equation when the c_ops list is non-empty:
$$\frac{\mathrm{d}\rho}{\mathrm{d}t} = -\frac{i}{\hbar}[H,\,\rho] + \sum_k\!\left(L_k\rho L_k^\dagger - \frac{1}{2}\{L_k^\dagger L_k,\,\rho\}\right)$$
Each element of c_ops is a collapse operator $C_k$; QuTiP constructs the full dissipator term internally as $C_k\rho C_k^\dagger - \frac{1}{2}\{C_k^\dagger C_k, \rho\}$ for each $k$. The collapse operator encodes both the jump operator and its rate: passing $C_k = \sqrt{\gamma}\,L_k$ is equivalent to a jump operator $L_k$ with rate $\gamma$.
The example below uses amplitude damping ($C = \sqrt{\gamma}\,\sigma^-$), which models spontaneous emission: the qubit loses energy to the environment at rate $\gamma$. Starting from $\lvert +\rangle\langle +\rvert$, the coherences ($\langle X\rangle$, $\langle Y\rangle$) decay exponentially while the population settles into $\lvert 0\rangle$.
# run: python3 lindblad_qutip.py # desc: evolve rho=|+><+| under precession H=(omega/2)*Z with amplitude damping at rate gamma import numpy as np import qutip as qt # Initial density matrix: rho = |+><+| psi0 = (qt.basis(2, 0) + qt.basis(2, 1)).unit() rho0 = qt.ket2dm(psi0) omega = 2 * np.pi H = 0.5 * omega * qt.sigmaz() # Amplitude damping: L = sqrt(gamma) * sigma_minus = sqrt(gamma) * destroy(2) # qt.basis(2, 0) is |0> (ground state); destroy(2) lowers |1> -> |0> gamma = 1.5 c_ops = [np.sqrt(gamma) * qt.destroy(2)] tlist = np.linspace(0, 3.0, 400) # Integrate the full Lindblad equation result = qt.mesolve(H, rho0, tlist, c_ops=c_ops, e_ops=[qt.sigmax(), qt.sigmaz()]) # <X> -> 0 as coherences are destroyed (T2 decay) # <Z> -> +1 as population settles in |0> (T1 decay; sigmaz|0> = +|0>) print(result.expect[0][-1]) # ~0.0 (decoherence) print(result.expect[1][-1]) # ~+1.0 (ground state)