# Fastest numerical method to solve Lindblad Master Equation?

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?

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.