Fastest numerical method to solve Lindblad Master Equation?

Question

The Lindblad Master Equation is a generalization of the Schrodinger Equation for open quantum systems, given by

$$ \frac{\mathrm{d} \rho}{\mathrm{d}t} = -i \left[ H, \rho\right] + \sum_k \gamma_k \left( L_k \rho L_k^\dagger - \frac{1}{2} \left\{ L_k^\dagger L_k, \rho \right\}\right) = \mathcal{L}(\rho) $$ where $\rho$ is the density matrix describing the state of the system, $H$ is some Hamiltonian, and $L_i$ are the jump operators. All of these objects are of size $N \times N$ where $N$ is the Hilbert size of the system.

Question: Assuming that everything is time-independent (constant Hamiltonian and jump operators), what would be the fastest numerical method to solve this equation? More specifically, starting from a known initial density matrix $\rho_0$, how can I get the density matrix at time $t_f$, $\rho(t_f)$?

Starting point: I already know of two possible methods.

(1) Computing the complete Lindbladian propagator, i.e. $\Lambda(t_f) = \exp(t_f \mathcal{L})$ and then doing the simple superoperator-matrix product $\rho_f = \Lambda(t_f) \rho_0$. However, $\mathcal{L}$ is a superoperator of size $N \times N \times N\times N$, so computing its exponential does not seem numerically very fast.

(2) Making some Kraus operator propagation, i.e.

$$ \rho(t+\mathrm{d}t) = \sum_\nu M_\nu \rho(t) M_\nu^\dagger $$ where $M_\nu$ are the Kraus operators, given by $$ M_0 = I - \mathrm{d}t \left(i H + \frac{1}{2} \sum_k L_k^\dagger L_k\right) \quad \text{and} \quad M_\nu = \sqrt{\mathrm{d}t} L_\nu $$ which is a trace preserving method up to second order in $\mathrm{d}t$. The good thing about this is that it only requires matrix-matrix products. However, we do need to take sufficiently many time steps, and thus perform quite a lot of matrix-matrix products.

Is there any known fast method for such problems?

Answer 1

Are you familiar with the Stochastic simulation method (quantum trajectories)? It reduces the cost from evolving and $N\times N$ density matrix, to wavefunctions $|\psi\rangle$ of only $N$ elements. This must be repeated until convergence, but can be highly parallellized and overcomes memory problems.

Depending on the nature of the problem, additional approximations to the state can be made that overcome the still exponential scaling (both when working with the master equation directly or with trajectories) such as Gutzwiller/tensor network/neural network/Gaussian ansatzes.

In practice, people may study complicated systems instead with a different framework: e.g. using PDEs for phase-space distributions, Keldysh field theory or the renormalisation group flow.

EDIT: for the promised references, I'll do a shameless self-promotion to my thesis, which also contains lots of further references. Chapter 2 discusses open quantum systems (including Lindblad equation and trajectories), Sec 3.3 Lists a number of existing advanced methods to deal with large systems.

Answer 2

You can use a Python package called QuTiP which will allow you to solve the LME. It also has different numerical methods to solve your problem, you can take a look at them and compare. Here is some documentation (http://qutip.org/docs/3.1.0/guide/dynamics/dynamics-master.html) There are also plenty of examples to aid you in using the package.