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

Figure 1

Diagrammatic representations of matrix products. (a) Diagram of the vectorized density matrix as an MPS in the mixed canonical form (5). (b) Diagram of the Liouvillian operator in the MPO representation (6). (c) Diagram representing the on-site Liouvillian operator ${\mathcal{L}}_l$. In all diagrams, triangles pointing right represent the left-normalized matrices $A$, while triangles pointing left denote the right-normalized matrices $B$, both entering the mixed canonical form of the MPS, Eq. (5). The local representation of the density matrix $\varrho$ is depicted as a circle at site $l$. Thin lines represent physical indices, while thick lines denote bond indices. The MPO matrices $W$ in Eq. (6) are represented as squares. (d) Diagrammatic scheme illustrating the computational complexity associated with index contractions at one site in the case where the simple Liouvillian $\mathcal{L}$ is used in the variational approach. (e) Same as (d) in the case where the squared Liouvillian ${\mathcal{L}}^{†}\mathcal{L}$ is instead used. In this second case, contraction of the MPO bond indices bears an additional $O\left(D_W\right)$ computational cost.

Figure 2

Comparison between the spatial correlations $\langle {\widehat{X}}_m{\widehat{X}}_{m+l}\rangle$ for $l=1,...,4$, computed by Trotter dynamics and through the direct variational MPS determination of the NESS. Parameters were set as $V=0$ and $\gamma =J$, and the array length was $N=15$. For the Trotter dynamics, the time step was set as $dt=0.1J$, and time integration was carried out until $T=10/\gamma$. For both the Trotter dynamics and the variational NESS methods the bond dimension was at most $D=20$. These results also displayed perfect agreement with MCWF calculations (not shown).

Figure 3

Comparison between the spatial correlations $\langle {\widehat{X}}_m{\widehat{X}}_{m+l}\rangle$ for $l=1,...,4$, computed by MCWF and through the direct variational MPS determination of the NESS. Parameters were set as $V=0.5J$ and $\gamma =J$, and the array length was $N=15$. The MCWF simulation was performed with 1000 trajectories, and time integration was carried out until $T=10/\gamma$. For the variational NESS calculations, the bond dimension was at most $D=30$.

Figure 4

Comparison between the spatial correlations $\langle {\widehat{X}}_m{\widehat{X}}_{m+l}\rangle$ for $l=1,...,4$, computed by MCWF and through the direct variational MPS determination of the NESS, in the case of weak dissipation. Parameters were set as $\gamma =0.1J$, and the array length was $N=15$. The MCWF simulation was performed with 16 trajectories, and time averages were taken from $t=10/\gamma$ to $t=100/\gamma$. For the variational NESS the bond dimension was at most $D=50$. The inset shows the correlations $\langle {\widehat{X}}_m{\widehat{X}}_{m+l}\rangle$ as a function of $l$ with $h=J$ for a system of $N=50$ sites and bond dimension $D=60$. This system size lies beyond the computational reach of the MCWF method.