Unfolding a quantum master equation into a system of real-valued equations: Computationally effective expansion over the basis of $\mathrm{SU}\left(N\right)$ generators

Dynamics of an open $N$-state quantum system is often modeled with a Markovian master equation describing the evolution of the system density operator. By using generators of $\text{SU}\left(N\right)$ group as a basis, the density operator can be transformed into a real-valued “coherence-vector.” A generator of the dissipative evolution, so-called “Lindbladian,” can be expanded over the same basis and recast in the form of a real matrix. Together, these expansions result is a nonhomogeneous system of $N^2-1$ real-valued linear ordinary differential equations. Now one can, e.g., implement standard high-performance algorithms to integrate the system of equations forward in time while being sure in exact preservation of the trace (norm) and Hermiticity of the density operator. However, when performed in a straightforward way, the expansion turns to be an operation of the time complexity $O\left(N^{10}\right)$. The complexity can be reduced when the number of dissipative operators is independent of $N$, which is often the case for physically meaningful models. Here we present an algorithm to transform quantum master equation into a system of real-valued differential equations and propagate it forward in time. By using a specific scalable model, we evaluate computational efficiency of the algorithm and demonstrate that it is possible to handle the model system with $N={10}^3$ states on a single node of a computer cluster.

Figure 1

Sparsity patterns of tensors $\left\{d_{ijk}\right\}$ (a) and $\left\{f_{ijk}\right\}$ (b) for $N=3$. The red points indicate nonzero elements.

Figure 2

Error in the calculated density matrix by the quantum jump method and the proposed method as a function of computation time. For $N=200$, integration step $1E-4$ corresponds to a computation time of $18500$ trajectories using the quantum jump method and step $1.25E-5$ corresponds to $145300$ trajectories. For $N=400$, these values are equal to $\delta t=3E-4$ and 1600 trajectories, and $\delta t=5E-5$ and 9600 trajectories, respectively. The initial state is $\varrho \left(0\right)=|0\rangle \langle 0|$.

Figure 3

Computation time as a function of $N$. The total computation time (line; right $y$ axis) was measured for the propagation time $10T$. More detail analysis was performed for the propagation time $T$ (bars; left $y$ axis). To get estimates for longer propagation times, one has to scale linearly computation time of the ODE integration step.

Figure 4

Ratio between the actual computation time and the corresponding asymptotic predictions, Table 2 (time complexity for the dimer model) For a fixed $N$, the ratio is scaled to make the total time (to perform all steps) equal to one.

Figure 5

Scaling of memory used during the different steps of the algorithm as functions of $N$, for $N=200,400,600,800,1000$ (from bottom to top).

Figure 6

Ratio between the maximal memory use measured in computation experiments and the corresponding asymptotic predictions, Table 2 (space complexity for the dimer model). For a fixed $N$, the ratio is scaled to make the total time (to perform all steps) equal to one.