Lindblad Master Equation Solver¶
The dynamics of a closed (pure) quantum system is governed by the Schrödinger equation
where \(\Psi\) is the wave function, \(\hat H\) the Hamiltonian, and \(\hbar\) is Planck’s constant. In general, the Schrödinger equation is a partial differential equation (PDE) where both \(\Psi\) and \(\hat H\) are functions of space and time. For computational purposes it is useful to expand the PDE in a set of basis functions that span the Hilbert space of the Hamiltonian, and to write the equation in matrix and vector form
where \(\left|\psi\right>\) is the state vector and \(H\) is the matrix representation of the Hamiltonian. This matrix equation can, in principle, be solved by diagonalizing the Hamiltonian matrix \(H\). In practice, however, it is difficult to perform this diagonalization unless the size of the Hilbert space (dimension of the matrix \(H\)) is small. Analytically, it is a formidable task to calculate the dynamics for systems with more than two states. If, in addition, we consider dissipation due to the inevitable interaction with a surrounding environment, the computational complexity grows even larger, and we have to resort to numerical calculations in all realistic situations. This illustrates the importance of numerical calculations in describing the dynamics of open quantum systems, and the need for efficient and accessible tools for this task.
The Schrödinger equation, which governs the time-evolution of closed quantum systems, is defined by its Hamiltonian and state vector. In the previous section, Using Tensor Products and Partial Traces, we showed how Hamiltonians and state vectors are constructed in QuTiP. Given a Hamiltonian, we can calculate the unitary (non-dissipative) time-evolution of an arbitrary state vector \(\left|\psi_0\right>\) (
psi0) using the QuTiP function
qutip.mesolve. It evolves the state vector and evaluates the expectation values for a set of operators
expt_ops at the points in time in the list
times, using an ordinary differential equation solver. Alternatively, we can use the function
qutip.essolve, which uses the exponential-series technique to calculate the time evolution of a system. The
qutip.essolve functions take the same arguments and it is therefore easy switch between the two solvers.
For example, the time evolution of a quantum spin-1/2 system with tunneling rate 0.1 that initially is in the up state is calculated, and the expectation values of the \(\sigma_z\) operator evaluated, with the following code
In : H = 2 * np.pi * 0.1 * sigmax() In : psi0 = basis(2, 0) In : times = np.linspace(0.0, 10.0, 20) In : result = sesolve(H, psi0, times, [sigmaz()])
The brackets in the fourth argument is an empty list of collapse operators, since we consider unitary evolution in this example. See the next section for examples on how dissipation is included by defining a list of collapse operators.
The function returns an instance of
qutip.solver.Result, as described in the previous section Dynamics Simulation Results. The attribute
result is a list of expectation values for the operators that are included in the list in the fifth argument. Adding operators to this list results in a larger output list returned by the function (one array of numbers, corresponding to the times in times, for each operator)
In : result = sesolve(H, psi0, times, [sigmaz(), sigmay()]) In : result.expect Out: [array([ 1. , 0.78914057, 0.24548559, -0.40169513, -0.8794735 , -0.98636142, -0.67728219, -0.08258023, 0.54694721, 0.94581685, 0.94581769, 0.54694945, -0.08257765, -0.67728015, -0.98636097, -0.87947476, -0.40169736, 0.24548326, 0.78913896, 1. ]), array([ 0.00000000e+00, -6.14212640e-01, -9.69400240e-01, -9.15773457e-01, -4.75947849e-01, 1.64593874e-01, 7.35723339e-01, 9.96584419e-01, 8.37167094e-01, 3.24700624e-01, -3.24698160e-01, -8.37165632e-01, -9.96584633e-01, -7.35725221e-01, -1.64596567e-01, 4.75945525e-01, 9.15772479e-01, 9.69400830e-01, 6.14214701e-01, 2.77159958e-06])]
The resulting list of expectation values can easily be visualized using matplotlib’s plotting functions:
In : H = 2 * np.pi * 0.1 * sigmax() In : psi0 = basis(2, 0) In : times = np.linspace(0.0, 10.0, 100) In : result = sesolve(H, psi0, times, [sigmaz(), sigmay()]) In : fig, ax = subplots() In : ax.plot(result.times, result.expect); In : ax.plot(result.times, result.expect); In : ax.set_xlabel('Time'); In : ax.set_ylabel('Expectation values'); In : ax.legend(("Sigma-Z", "Sigma-Y")); In : show()
In : times = [0.0, 1.0] In : result = mesolve(H, psi0, times, , ) no collapse operator, using sesolve In : result.states