A Matrix-Product-Operator Approach to the Nonequilibrium Steady State of Driven-Dissipative Quantum Arrays

We develop a numerical procedure to efficiently model the nonequilibrium steady state of one-dimensional arrays of open quantum systems, based on a matrix-product operator ansatz for the density matrix. The procedure searches for the null eigenvalue of the Liouvillian superoperator by sweeping along the system while carrying out a partial diagonalization of the single-site stationary problem. It bears full analogy to the density-matrix renormalization group approach to the ground state of isolated systems, and its numerical complexity scales as a power law with the bond dimension. The method brings considerable advantage when compared to the integration of the time-dependent problem via Trotter decomposition, as it can address arbitrarily long-ranged couplings. Additionally, it ensures numerical stability in the case of weakly dissipative systems thanks to a slow tuning of the dissipation rates along the sweeps. We have tested the method on a driven-dissipative spin chain, under various assumptions for the Hamiltonian, drive, and dissipation parameters, and compared the results to those obtained both by Trotter dynamics and Monte-Carlo wave function. Accurate and numerically stable convergence was always achieved when applying the method to systems with a gapped Liouvillian and a non-degenerate steady-state.