Solvable Markovian dynamics of lattice quantum spin models

TL;DR

This paper presents an exact method for analyzing the real-time dynamics of lattice quantum spin models with Markovian dissipation, enabling the study of large systems and their steady states.

Contribution

The authors develop a closed hierarchy approach that allows exact solutions of quantum spin models in arbitrary dimensions with dissipation.

Findings

Exact solutions for driven-dissipative quantum spin models.

Analytical derivation of characteristic time scales and finite-size scaling.

Numerical confirmation of dissipative gap scaling exponents.

Abstract

We address the real-time dynamics of lattice quantum spin models coupled to single or multiple Markovian dissipative reservoirs using the method of closed hierarchies of correlation functions. This approach allows us to solve a number of quantum spin models exactly in arbitrary dimensions, which is illustrated explicitly with two examples of driven-dissipative systems. We investigate their respective nonequilibrium steady states as well as the full real-time evolution on unprecedented system sizes. Characteristic time scales are derived analytically, which allows us to understand the nontrivial finite-size scaling of the dissipative gap. The corresponding scaling exponents are confirmed by solving numerically for the full real-time evolution of two-point correlation functions.

Figures

Equations

$$\frac{d}{dt}A(t)=\frac{i}{\hslash }\left[H,A\right]+\sum _{\ell }{\gamma }_{\ell }{\mathcal{D}}_{\ell }(A).$$

$$\boldsymbol{\mathcal{D}}_{\ell}(\hskip 0.9pt\bullet\hskip 0.7pt)=\ell^{\dagger}\hskip 0.5pt\big{[}\hskip 0.9pt\bullet\hskip 0.7pt,\ell\big{]}+\big{[}\ell^{\dagger},\hskip 0.7pt\bullet\hskip 0.9pt\big{]}\hskip 0.9pt\ell,$$

$$A_m:{\mathtt{C}}_m\to \mathfrak{B}({\mathcal{H}}_N),$$

$$i\left[H,A_m(X)\right]+\sum _{\ell }{g}_{\ell }{\mathcal{D}}_{\ell }\left(A_m(X)\right)=c_0+\sum _{n=1}^{\infty }\sum _{Y\in {\mathtt{C}}_n}{c}_n(Y)A_n(Y),$$

$$c_n(Y)=0,\forall n>\mathcal{N}\left(A_m(X)\right).$$

$$\mathcal{N}\left(A_m(X)\right)=m,\forall X,$$

$$\mathcal{N}\left(A_m(X)\right)=m+k,k\ge 0,$$

$$\ell \in \left\{L_{xy}=\frac{1}{2}{\delta }_{⟨x,y⟩}(s_x^{+}+s_y^{+})(s_x^{-}-s_y^{-})\right\},$$

$$\frac{d}{dt}{A}_m(X)=\sum _{Y\in {\mathtt{C}}_m}{c}_m(Y)A_m(Y),$$

$$\frac{d}{dt}{C}_{xy}$$

$$\frac{d}{dt}{D}_{xy}$$

$$C_{xy}(t=0)=4D_{xy}(t=0)={\delta }_{xy}.$$

$$\frac{\partial }{\partial t}R=\frac{{\partial }^2}{\partial r^2}R+(d-1)\frac{1}{r}\frac{\partial }{\partial r}R,$$

$$R(\epsilon <r<\xi (t))\sim \{\begin{array}{cc}1-r/\xi (t),& d=1,\\ 1-\ln \left(r/\epsilon \right)/\ln \left(\xi (t)/\epsilon \right),& d=2,\end{array}$$

$$R(r>\epsilon ,t)\sim (\epsilon /r{)}^{d-2},d>2.$$

$$t_N\sim \{\begin{array}{cc}{N}^2,& d=1,\\ N\ln N,& d=2,\\ N,& d>2.\end{array}$$

$$\big{(}L_{1}^{\alpha}\big{)}_{xy}=\frac{1}{2}\delta_{\langle x,y\rangle}\begin{pmatrix}\sqrt{2}s_{x}^{+}s_{y}^{+}(s_{x}^{-}-s_{y}^{-})\\[5.0pt] ~{}-\sqrt{2}s_{x}^{-}s_{y}^{-}(s_{x}^{+}-s_{y}^{+})~{}\\[5.0pt] (s_{x}^{+}+s_{y}^{+})(s_{x}^{-}-s_{y}^{-})\end{pmatrix},$$

$$\big{(}L_{2}^{\beta}\big{)}_{xy}=\frac{1}{2}\delta_{\langle x,y\rangle}\begin{pmatrix}(\mathbbm{1}+2s_{x}^{3})s_{y}^{-}\\[5.0pt] s_{x}^{-}(\mathbbm{1}+2s_{y}^{3})\\[5.0pt] ~{}s_{x}^{+}(\mathbbm{1}-2s_{y}^{3})~{}\\[5.0pt] (\mathbbm{1}-2s_{x}^{3})s_{y}^{+}\end{pmatrix}.$$

$$\frac{d}{dt}M$$

$$\frac{d}{dt}{M}_{\mathrm{s}}$$

$$\frac{d}{dt}{C}_{xy}$$

$$-{\delta }_{⟨x,y⟩}\left[(1-\lambda )(C_{xx}+C_{xy}+8D_{xy})-4\lambda C_{xy}\right],$$

$$\frac{d}{dt}{D}_{xy}$$

$$-(1-\lambda ){\delta }_{⟨x,y⟩}\left[D_{xx}-D_{xy}+C_{xy}\right],$$

$$\frac{d}{dt}{\tilde{D}}_{xy}$$

$$+(1-\lambda ){\delta }_{⟨x,y⟩}\left[{\tilde{D}}_{xx}+{\tilde{D}}_{xy}+C_{xy}\right],$$

Full text

Introduction. The dynamics of quantum spin systems is of fundamental importance to atomic, molecular, and optical physics, condensed matter physics, quantum chemistry, and beyond. Yet, computing their real-time dynamics poses a tremendous challenge for classical computers, since the required resources scale exponentially with the system size. In general, standard methods for quantum spin systems in equilibrium as, e.g., Monte Carlo importance sampling Suzuki et al. (1977); Suzuki (1993), cannot be applied out of equilibrium due to the absence of a positive-definite probability measure. The density matrix renormalization group White (1992), on the other hand, which can be understood as a variational method on the class of matrix product states Verstraete et al. (2008); Schollwöck (2011), manages to reduce the exponential complexity by following the time-dependence of quantum states with low entanglement Cazalilla and Marston (2002); White and Feiguin (2004); Vidal (2004); Daley et al. (2004). Unfortunately, this restriction largely limits its applicability to one-dimensional systems and typically does not allow to follow the time evolution of the system to arbitrarily large times. Other time-dependent variational methods Haegeman et al. (2011); Transchel et al. rely on similar ansatz states and hence also face similar constraints. In view of these limitations, exactly solvable models are immensely useful. Not only do they provide important insights into intriguing nonequilibrium phenomena, but they also serve as benchmarks for systematic approximations and numerical techniques.

In this Letter, we investigate quantum spin models that are coupled to a Markovian reservoir, which can be solved exactly in the framework of closed hierarchies of correlation functions Caspar et al. (2016a, b) (see also earlier work Ref. Žunkovič (2014); Klich in the context of many-body boson and fermion systems). Their solution relies on a set of closure conditions that need to be satisfied in order for correlation functions of quantum spins (multi-point functions) to close among themselves. One is bound to ask whether these models are still sufficiently interesting to warrant further discussion – are they maybe “too simple to be real”? Here, we argue that the latter is not the case. Indeed, several protocols for the dissipative preparation of long-range entangled states Kraus et al. (2008); Diehl et al. (2008); Verstraete et al. (2009); Weimer et al. (2010) lie within this class of exactly solvable models, some of which have also been realized recently using ion-trap experiments Barreiro et al. (2011); Schindler et al. (2013). Other engineered dissipative systems with competing dissipative couplings Weimer (2017) are also directly accessible within our approach. This gives us the unique opportunity to investigate these systems on sufficiently large system sizes in order to determine those points in parameter space where the dissipative gap closes. There, the relevant length- and timescales diverge, which is characteristic of a dissipative phase transition Verstraete et al. (2009); Diehl et al. (2010); Eisert and Prosen ; Lesanovsky et al. (2013). In particular, we calculate the finite-size scaling of the dissipative gap, which characterizes the nature of the long-time relaxation towards the nonequilibrium steady state and essentially defines the nonequilibrium dynamical universality class. In combination with exact diagonalization of the quantum master equation on small system sizes, closed hierarchies allow us to construct a detailed picture of the full nonequilibrium dynamics. These insights might also help to elucidate the dynamics of systems that are outside of scope of these methods. In this spirit, we conclude this work with some general arguments regarding the properties of driven-dissipative quantum systems.

Closed hierarchies of correlations. Before we turn to specific examples, we need to make the nature and properties of closed hierarchies precise. We will be concerned with open systems of quantum spin-$1/2$ particles (qubits) on a regular $d$-dimensional hypercubic lattice with periodic boundary conditions. The latter assumption is neither necessary for the following arguments, nor essential for our conclusions, but it allows for exploiting the translation symmetry of the problem. We assume that the spins are weakly coupled to an external reservoir and that there is a separation of timescales between the dynamics of the quantum spin system and the reservoir. This second assumption allows us to integrate out the “fast” degrees of freedom of the bath and to consider only the reduced density matrix $\rho_{N}$ (i.e., the statistical ensemble) of the $N$-spin system. Proceeding along these lines, we arrive at an effective description of the dynamics of the reduced system $\rho_{N}$, which is described by a Markovian quantum master equation Breuer and Petruccione (2002).

To address the time evolution of observables $A$, which are elements of the set $\mathfrak{B}(\mathcal{H}_{N})$ of bounded linear operators on the Hilbert space $\mathcal{H}_{N}\simeq\mathbbm{C}^{2^{N}}$, it is convenient to work in the Heisenberg picture, where their time evolution is governed by the (adjoint) Markovian quantum master equation

[TABLE]

$H$ is the system Hamiltonian, and $\gamma_{\ell}=\gamma g_{\ell}$ defines the rate of dissipation, where $\gamma$ sets the scale and $g_{\ell}$ denotes the dimensionless dissipative coupling. In the following, we will always set $\hbar=\gamma=1$. The “dissipator” $\boldsymbol{\mathcal{D}}_{\ell}$ is conveniently represented in Lindblad form Gorini et al. (1976); Lindblad (1976)

[TABLE]

where $\ell$ are (Lindblad) jump operators that describe the effective interaction with the environment. If one considers individual realizations of the statistical ensemble $\rho_{N}$, these jump operators yield a stochastic time evolution of observables Gardiner (1991); Breuer and Petruccione (2002).

Specifically, we will be interested in the time evolution of multi-point functions

[TABLE]

which can be viewed as maps from the configuration space of $m$ spin-$1/2$ particles, $\mathtt{C}_{m}$, to the space of bounded operators on the $N$-particle Hilbert space, $\mathcal{H}_{N}$. We define the map (3) for $m>0$ such that $A_{m}(X)=\prod_{i=1}^{m}s_{x_{i}}^{a_{i}}$ for each point $X=\{(x_{i},a_{i})\}_{1\leq i\leq m}\in\mathtt{C}_{m}$ (i.e., $x_{i}\neq x_{j}$, $\forall i\neq j$) and $A_{0}$ is identified with the identity operator. Inserting the operator $A_{m}(X)$ into Eq. (1), we observe that the commutators on the right hand side of this equation typically induce operators of higher order, i.e.,

[TABLE]

with coefficients $c_{n}\in\mathbb{R}$. For most systems of interest, Eq. (4) yields an infinite hierarchy of equations of motion that cannot be solved without truncation. One may ask, however, under what conditions the hierarchy closes exactly, such that

[TABLE]

$\mathcal{N}\!\left(A_{m}(X)\right)$ is an integer that, in general, might depend on the particular observable under consideration. Stated differently, we are asking the following question: Are there quantum spin systems that yield a closed hierarchy, and if yes, what are their properties? The strongest form of the closure conditions

[TABLE]

was derived by us in previous work Caspar et al. (2016b). This condition constrains the form of the Hamiltonian $H$ and the jump operators $\ell$ and defines a class of efficiently solvable models (in the sense that the problem of computing the time evolution of $m$-point functions scales only polynomially in the system size $N$). The general idea behind condition (6) can be easily summarized as follows: Inserting $m$ quantum spin operators $s_{x_{i}}^{a_{i}}$ on lattice sites $x_{i}$, $i=1,2,\ldots,m$ (which defines the $m$-point operator, $A_{m}(X)$, of interest), the infinitesimal generator of the dynamical semigroup $\boldsymbol{\mathcal{G}}(\hskip 0.9pt\bullet\hskip 0.7pt)=i\left[H,\hskip 0.9pt\bullet\hskip 0.7pt\right]+\sum_{\ell}g_{\ell}\boldsymbol{\mathcal{D}}_{\ell}\left(\hskip 0.9pt\bullet\hskip 0.7pt\right)$ is allowed only (i) to move around single quantum spin operators from one site to another, (ii) to convert them into different “species”, i.e., $s_{x}^{a}\mapsto s_{y}^{b}$, or (iii) to completely remove them from the lattice. Thus, systems that satisfy the closure condition (6) are characterized by hopping of quantum spins from one site to another, with some probability for single spins being eliminated from the lattice. In a diagrammatic language, these processes can be viewed either as $m\rightarrow m-k$ particle scattering ($A_{m}\rightarrow A_{m-k}$, $0\leq k\leq m-1$), or decay $m\rightarrow 0$ ($A_{m}\rightarrow\varnothing$). The opposite process, however, namely particle creation ($A_{m}\rightarrow A_{m+k}$, $k>0$) is not allowed. Therefore, if we define the set of operators $A|_{\mathtt{C}_{m}}\equiv A_{m}(\mathtt{C}_{m})$, we may write Eq. (6) in short as $\boldsymbol{\mathcal{G}}(A|_{\mathtt{C}_{m}})=\bigoplus_{n=0}^{m}c_{n}A|_{\mathtt{C}_{n}}$.

One might imagine less restrictive conditions on the dynamics, e.g.,

[TABLE]

with $k$ finite. However, it is easy to see that the hierarchy is unbounded if transitions of the type $A_{m}\rightarrow A_{m+k}$, $k>0$, are allowed in general. The key to construct closed hierarchies in this case is the following observation: The action of the infinitesimal generator $\boldsymbol{\mathcal{G}}$ should map the set of operators $A|_{\mathtt{C}_{m}}$ to a subset $A|_{\mathtt{S}}\subseteq\bigoplus_{n=0}^{m+k}c_{n}A|_{\mathtt{C}_{n}}$, which itself defines an invariant subspace for the dynamics, i.e., $\boldsymbol{\mathcal{G}}(A|_{\mathtt{S}})\subseteq A|_{\mathtt{S}}$. In this case, only the knowledge of a subset of operators, namely $A|_{\mathtt{S}}$, is needed to determine the time evolution of $m$-point functions $A_{m}$. Finding efficiently solvable models that satisfy Eqs. (6) or (7) hence corresponds to the problem of identifying invariant subspaces $A|_{\mathtt{S}}$ which are small (i.e., whose size should scale only polynomially with the system size). The dissipative transverse-field Ising model Foss-Feig et al. is one example of a system that satisfies Eq. (7), which is of interest in the context of “engineered” quantum systems, e.g., using ion trap experiments Gärttner et al. or Rydberg atoms Weimer et al. (2010); Glaetzle et al. (2015); Roghani and Weimer . Here, we will limit ourselves to the condition Eq. (6) and discuss two examples within this class of problems that have been at the focus of interest recently.

Dissipative quantum state preparation. We consider the set of non-Hermitian jump operators Diehl et al. (2008)

[TABLE]

that act uniformly on adjacent sites with $g_{\ell}=1$; $\delta_{\langle x,y\rangle}$ is nonzero and equal to one if and only if $x$ and $y$ are nearest neighbor sites, i.e., $||x-y||=1$ (here and in the following the lattice spacing is set to one). Note that $L_{xy}$ acts on pairs of nearest-neighbor sites $(x,y)$ and maps the two-particle spin singlet state to the spin triplet, while it conserves $S^{3}=\sum_{x}s_{x}^{3}$. At the same time it annihilates the spin triplet state, i.e., $L_{xy}^{2}=0$.

In the absence of a Hamiltonian, $H=0$, we obtain a linear system of equations of motion

[TABLE]

which can be solved exactly for arbitrary $m$-point functions.111Strictly speaking, the assumption $H=0$ is not necessary since Eq. (6) also allows for the possibility of including a coupling to an external field, i.e., $H=\sum_{x,a}h_{x}^{a}s_{x}^{a}$, with $h_{x}^{a}\in\mathbbm{R}$. The dynamics of the (ensemble-averaged) two-point correlation functions C_{xy}=\operatorname{tr}\rho_{N}\big{(}s_{x}^{+}s_{y}^{-}+s_{x}^{-}s_{y}^{+}\big{)} and $D_{xy}=\operatorname{tr}\rho_{N}s_{x}^{3}s_{y}^{3}$ is governed by the evolution equations ($x\neq y$)

[TABLE]

$\Delta_{xy}\equiv\Delta_{x}+\Delta_{y}$ is the discrete Laplacian on the hypercubic lattice and Eqs. (10) and (11) are subject to the boundary condition $C_{xx}=4D_{xx}=1$. We choose to initialize the system in a fully disordered state with the reduced density matrix $\rho_{N}=2^{-N}\mathbbm{1}$, i.e.,

[TABLE]

This set up allows us to study the build-up of correlations by putting the spin system in contact with the “engineered” reservoir. In fact, these initial conditions are rather convenient, since Eq. (11) can be solved trivially in this case – note that the two-point correlation function is invariant in time $D_{xy}(t)=(1/4)\delta_{xy}$ so that the growth of correlations is fully described by the function $C_{xy}$ alone.

In order to study its long-time behavior, it is useful to consider the limit of large separations $||x-y||\gg 1$, in which case the contact term on the right hand side of Eq. (10) vanishes. If we furthermore assume that the initial state $\rho_{N}$ is translationally invariant and isotropic (which, of course, is true for the fully disordered state), then Eq. (10) can be mapped onto the continuum diffusion equation, which in spherical coordinates reads

[TABLE]

i.e., $R=R(r,t)$, where $r$ is the radial coordinate and the boundary condition, $R(r=\varepsilon,t)=1$, $\varepsilon>0$, is imposed. Note that Eq. (13) describes the behavior of the correlation function $C_{xy}$ exactly in the limit of large separations in arbitrary spatial dimensions $d$. Solving this equation is a standard exercise in statistical physics Krapivsky et al. (2010). At late times, one obtains

[TABLE]

where $\xi(t)\sim t^{1/2}$, while

[TABLE]

Thus, asymptotically the system develops long-range correlations. For $d\leq 2$ we find $R(r,t\rightarrow\infty)=1$, while correlations decay algebraically for $d>2$. We may estimate the timescale $t_{N}$ on which the system reaches the final, long-range correlated state by requiring $(1/N)\sum_{x,y}C_{xy}\sim\lim_{\varepsilon\rightarrow 0}\int_{\varepsilon}^{\xi(t)}dr\,r^{d-1}R(r,t)\sim N$, which yields

[TABLE]

These results are confirmed explicitly by solving numerically for the ensemble-averaged two-point function $C_{xy}$ on the hypercubic lattice starting from the fully disordered state Caspar et al. (2016a). Similarly, the spectral gap $\delta$ of the infinitesimal generator of the quantum dynamical semigroup $\boldsymbol{\mathcal{G}}$ scales as $\delta\sim t_{N}^{-1}$. The vanishing of $\delta$ in the limit $N\rightarrow\infty$ signals the presence of a diverging time- and lengthscale, $t_{N}\sim\xi^{2}$. This means that the system becomes critical in the infinite-volume limit.

Competing dissipative couplings. As a second example we consider a model that was introduced recently Weimer (2017) to investigate the possibility of purely dissipation-driven phase transitions in lattice quantum spin systems with Hamiltonian $H=0$. The non-Hermitian jump operators $\ell$ are now taken from two distinct sets of bilocal operators, i.e., \mathcal{S}_{1}=\big{\{}\big{(}L_{1}^{\alpha}\big{)}_{xy}\big{\}} with

[TABLE]

and \mathcal{S}_{2}=\big{\{}(L_{2}^{\beta})_{xy}\big{\}}, with

[TABLE]

The associated dissipative couplings are given by $g_{\ell\in\mathcal{S}_{1}}\equiv g_{1}=1-\lambda$ and $g_{\ell\in\mathcal{S}_{2}}\equiv g_{2}=\lambda$, with $\lambda\in[0,1]$, such that the infinitesimal generator can be decomposed into two distinct contributions $\boldsymbol{\mathcal{G}}=\boldsymbol{\mathcal{G}}_{1}+\boldsymbol{\mathcal{G}}_{2}$, with $\boldsymbol{\mathcal{G}}_{i}\equiv g_{i}\sum_{\ell\in\mathcal{S}_{i}}\boldsymbol{\mathcal{D}}_{\ell}$. Note that the infinitesimal generator $\boldsymbol{\mathcal{G}}_{1}$ conserves $S^{3}$, i.e., \boldsymbol{\mathcal{G}}_{1}\big{(}S^{3}\big{)}=0, while $\boldsymbol{\mathcal{G}}_{2}$ conserves $\widetilde{S}^{3}=\sum_{x}(-1)^{||x||}s_{x}^{3}$. Nevertheless, the total generator for the dynamics $\boldsymbol{\mathcal{G}}$ conserves neither $S^{3}$ or $\widetilde{S}^{3}$, which makes it clear that the points $\lambda=0$ and $\lambda=1$ in parameter space are special.

Again, we find that the hierarchy of correlation functions closes, arriving at an equation of the form (9). In particular, the evolution equations for the average magnetization $M=(1/N)\operatorname{tr}\rho_{N}S^{3}$ and the average staggered magnetization $M_{\rm s}=(1/N)\operatorname{tr}\rho_{N}\widetilde{S}^{3}$ are given by

[TABLE]

From the equations of motion, we confirm that the average (staggered) magnetization is conserved for $\lambda=0$ ($\lambda=1$). On the other hand, in the absence of conserved charges $0<\lambda<1$ both $M$ and $M_{\rm s}$ decay to zero, which characterizes a unique fixed point for the dynamics that is independent of the initial conditions.

The equations of motion for the ensemble-averaged two-point correlation functions $C_{xy}$ and $D_{xy}$ (defined above Eq. (10) and (11)), are given by ($x\neq y$)

[TABLE]

with $C_{xx}=4D_{xx}=1$. It is useful to employ a different choice of basis and define the two-point function $\widetilde{D}_{xy}=(-1)^{||x-y||}\operatorname{tr}\rho_{N}s_{x}^{3}s_{y}^{3}$. This function signals the presence of antiferromagnetic correlations and its equation of motion ($x\neq y$) reads

[TABLE]

with $\widetilde{D}_{xx}\equiv D_{xx}=1/4$. Initializing the system in the fully disordered initial state $\rho_{N}=2^{-N}\mathbbm{1}$, we determine the presence of long-range correlations in the nonequilibrium steady state (Fig. 1), and calculate the characteristic timescales associated with the growth of correlations induced by the “engineered” dissipative reservoir (Fig. 2). This problem is discussed separately for two distinct cases, i.e., either $\lambda\in\{0,1\}$ or $0<\lambda<1$. In the former case, the infinitesimal generator $\boldsymbol{\mathcal{G}}=\sum_{\ell}g_{\ell}\boldsymbol{\mathcal{D}}_{\ell}$ conserves either $M$ ($\lambda=0$) or $M_{\rm s}$ ($\lambda=1$) and the finite-size scaling of $t_{N}$ is identical to that derived in (16). This is clear from the diffusive nature of the dynamics, which is ultimately a consequence of the presence of a (global) conserved charge. Long-range correlations in $C_{xy}$ ($\widetilde{D}_{xy}$) are hence only observed infinitesimally close to the special points $\lambda=0$ and $\lambda=1$ (cf. Fig. 1). On the other hand, for $0<\lambda<1$, the long-time behavior is fully determined by the presence of the local relaxation. Consequently, the characteristic time scale $t_{N}$ and the dissipative gap $\delta$ are finite for all $N$ (cf. Fig. 2). This explicitly rules out the possibility of a driven-dissipative phase transition between the two distinct long-range correlated phases at $\lambda=0$ and $\lambda=1$, which has been suggested previously Weimer (2017). However, we want to point out that the dissipative gap $\delta$ nevertheless displays a nonanalytic behavior at $\lambda=0.6$, which indicates a dynamic transition between a regime where the late-time dynamics is characterized either by the set of jump operators $\mathcal{S}_{1}$ or $\mathcal{S}_{2}$.

Conclusions. We have shown that closed hierarchies of correlation functions provide a unique perspective on the properties of open quantum spin systems. They allow us to reach definite conclusions regarding the vanishing of the dissipative gap, which has been suggested as a signature for steady-state phase transitions. In this way, specific proposals for quantum-state preparation via “engineered” dissipation can be solved exactly on large system sizes. Furthermore, this allows us to directly contrast our results to those of other numerical methods as, e.g., functional integral methods Hebenstreit (2016) or variational methods Haegeman et al. (2011); Transchel et al. ; Weimer (2015a, b), which can be tested on large system sizes in arbitrary dimensions. In particular, having discussed a recently introduced purely dissipative model for quantum magnetism Weimer (2017), we find no indications of a nonequilibrium steady-state phase transition, which was suggested by variational calculations. We want to point out however that this does not exclude the possibility of singular (or nonanalytic) behavior in systems that exhibit closed hierarchies (see, e.g., Linzner et al. (2016)).

In closing, we want to mention related work on “solvable” models in which local jump operators compete with the Hamiltonian dynamics Foss-Feig et al. . While these models feature a closed hierarchy in the sense described in this work, they are different from the purely dissipative quantum spin models considered here, where jump operators act on pairs of neighboring lattice sites. Nevertheless, these two examples demonstrate that closed hierarchies provide a much broader framework for classifying and addressing “solvable” models with increasing computational complexity in a systematic way. It will certainly be interesting to study more elaborate closure conditions and their physical realizations in future work.

Acknowledgments. We want to thank S. Diehl, M. Fleischhauer, I. Lesanovsky, and U.-J. Wiese for illuminating discussions. This research is funded by the Swiss National Science Foundation and the European Research Council under the European Union’s Seventh Framework Programme, FP7/2007-2013, 339220.