Matrix-product-operator approach to the nonequilibrium steady state of driven-dissipative quantum arrays

Abstract

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 methods. Accurate and numerically stable convergence was always achieved when applying the method to systems with a gapped Liouvillian and a nondegenerate steady state.

Linked cluster expansions for open quantum systems on a lattice

“…While in one dimension tensor-network approaches based on a straightforward generalization of matrix product states to operators can be effective [47][48][49] and alternative strategies have been proposed in order to improve their performances [50][51][52] , going to higher dimensions is much harder. Numerical strategies specifically suited for this purpose have been recently put forward, including cluster mean-field 53 , correlated variational Ansätze 54,55 , truncated correlation hierarchy schemes 56 , corner-space renormalization methods 57 , and even two-dimensional tensor-network structures 58 .…”

Spontaneous Beliaev-Landau scattering out of equilibrium

“…Because of the absence of a spectral gap and the convexity of the Bogoliubov dispersion of conservative continuum systems, criterium (20) in continuum models can only be satisfied in two or more spatial dimensions, while in 1D sytems phonons can only decay through higher-order scattering processes [44]. However, there exist specifically engineered 1D optical lattices with a nonconvex (but gapless) spectrum, such that energy and momentum conservation can be simultaneously satisfied [45,46].…”

“…The allowed wavevectors k and q that exactly satisfy the energy and momentum conservation condition (20) in one-dimension are indicated in Fig. 1b).…”

“…First of all the total photon number is not conserved, resulting in a much larger effective Hilbert space, and secondly the steady state is a mixed state, thus requiring the evaluation of a full density matrix rather than a single wavefunction. This led to the development of new numerical tools such as, among others, variational approaches based on matrix product states in 1D [17][18][19][20], resummation techniques [21], self-consistent projection operator theory [22], extensions of the variational principle [23] and the corner-space renormalization method for 2D lattices [24].…”

Gaussian Quantum Trajectories for the Variational Simulation of Open Quantum-Optical Systems

“…We are interested in a method that is variational and provides a description on the level of trajectories. Some current variational (Time-Dependent Variational Principle (TDVP) [19]) approaches to the description of open systems aim to describe the system at the level of the master equation [20][21][22][23][24], with ansatzes including the Gutzwiller density matrix- [25][26][27], matrix-product state [28,29] and matrix-product operator [30][31][32][33] methods. Variational descriptions on the trajectory level have so far focused on lattice-like systems: systems where the complexity arises mainly due to a large amount of modes, while the Hilbert space per mode remains small, with methods such as Gutzwiller Monte Carlo [34][35][36] or based on t-DMRG [37][38][39].…”