A solvable family of driven-dissipative many-body systems

TL;DR

This paper introduces a family of non-equilibrium, dissipative many-body models that are exactly solvable across dimensions, enabling new insights into steady states and decoherence effects relevant for quantum systems.

Contribution

It presents a class of non-equilibrium models analogous to the transverse-field Ising model that are exactly solvable, and uses these solutions to analyze phase transitions and decoherence in quantum systems.

Findings

Broad class of models solvable in any dimension

Proved a no-go theorem on steady-state phase transitions

Analyzed decoherence effects in quantum computing architectures

Abstract

Exactly solvable models have played an important role in establishing the sophisticated modern understanding of equilibrium many-body physics. And conversely, the relative scarcity of solutions for non-equilibrium models greatly limits our understanding of systems away from thermal equilibrium. We study a family of non-equilibrium models, some of which can be viewed as dissipative analogues of the transverse-field Ising model, in that an effectively classical Hamiltonian is frustrated by dissipative processes that drive the system toward states that do not commute with the Hamiltonian. Surprisingly, a broad and experimentally relevant subset of these models can be solved efficiently in any number of spatial dimensions. We leverage these solutions to prove a no-go theorem on steady-state phase transitions in a many-body model that can be realized naturally with Rydberg atoms or trapped ions, and to compute the effects of decoherence on a canonical trapped-ion-based quantum computation architecture.

Equations

$$\hat{H}=\sum _{j\ne k}{h}_{j,k}{\hat{\sigma }}_j^z{\hat{\sigma }}_k^z.$$

$$\dot{\hat{\rho }}$$

$$\mathcal{D}(\star )$$

$$\displaystyle{\rm Tr}\big{(}\hat{\sigma}^{z}_{j}\mathcal{D}(\hat{\sigma}^{\pm}_{j})\big{)}=0~{}~{}~{}({\rm for~{}all}~{}j).$$

$$\displaystyle{\rm Tr}\big{(}\hat{\sigma}^{z}_{j}\mathcal{L}(\hat{\sigma}^{\pm}_{j})\big{)}=0~{}~{}~{}({\rm for~{}all}~{}j).$$

$$\displaystyle\begin{matrix}{\text{\small Schr\"{o}dinger picture:}}&{\text{\small Heisenberg picture:}}\\ O(t)={\rm Tr}\big{(}\hat{O}\exp(t\mathcal{L})\,\hat{\rho}_{0}\big{)},~{}&O(t)={\rm Tr}\big{(}\hat{\rho}_{0}\exp(t\mathcal{L}^{\ddagger})\,\hat{O}\big{)}.\end{matrix}$$

$$\displaystyle O(t)={\rm Tr}_{\mathscr{A}\bigcup\mathscr{B}}\big{(}{\rm Tr}_{\mathscr{C}}\big{(}\hat{\rho}_{0}\sum_{n=0}^{\infty}\frac{t^{n}}{n!}\hat{O}_{n}\big{)}\big{)},$$

$$\displaystyle O(t)={\rm Tr}_{\mathscr{A}\bigcup\mathscr{B}}\big{(}\hat{\rho}_{\scriptscriptstyle\!\mathscr{A}\!\mathscr{B}}\sum_{n=0}^{\infty}\frac{t^{n}}{n!}(\mathcal{L}^{\ddagger}_{\mathscr{A\!B}})^{n}\hat{O}\big{)},$$

$$\displaystyle O(t)={\rm Tr}_{\mathscr{A}\bigcup\mathscr{B}}\big{(}\hat{O}\exp(t\mathcal{L}_{\mathscr{A\!B}})\hat{\rho}_{\scriptscriptstyle\!\mathscr{A}\!\mathscr{B}}\big{)}.$$

$${\hat{H}}_{\mathrm{TFIM}}=J\sum _{⟨j,k⟩}{\hat{\sigma }}_j^x{\hat{\sigma }}_k^x+\Delta \sum _j{\hat{\sigma }}_j^z,{\hat{J}}_j={\hat{\sigma }}_j^{-},$$

$${\hat{H}}_{\mathrm{fTFIM}}(t)=J\sum _{⟨j,k⟩}{\hat{\sigma }}_j^x{\hat{\sigma }}_k^x+\sum _j{\Delta }_j(t){\hat{\sigma }}_j^z.$$

$$\dot{\hat{\rho }}=-i[J\sum _{⟨j,k⟩}{\hat{\sigma }}_j^x{\hat{\sigma }}_k^x,\hat{\rho }]+\gamma \sum _j({\hat{\sigma }}_j^z\hat{\rho }{\hat{\sigma }}_j^z-\hat{\rho }).$$

$$\hat{\rho }=\sum _{\{{\mu }_1,\dots ,{\mu }_N\}}{\rho }_{{\mu }_1,\dots {\mu }_N}{\hat{\sigma }}_1^{{\mu }_1}\otimes \cdots \otimes {\hat{\sigma }}_N^{{\mu }_N};{\mu }_j\in \{1,x,y,z\},$$

$$\mathcal{L}=⟮\begin{array}{cccccc}{\mathcal{D}}^{00}& 0& 0& \cdots & 0& 0\\ {\mathcal{D}}^{10}& {\mathcal{D}}^{11}+i{\mathcal{H}}^{11}& 0& \cdots & 0& 0\\ 0& {\mathcal{D}}^{21}& {\mathcal{D}}^{22}+i{\mathcal{H}}^{22}& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & {\mathcal{D}}^{N,N-1}& {\mathcal{D}}^{NN}+i{\mathcal{H}}^{NN}\end{array}⟯.$$

$$\mathrm{eigs}(\mathcal{L})=\mathrm{eigs}({\mathcal{D}}^{00})\cup \mathrm{eigs}({\mathcal{D}}^{11}+i{\mathcal{H}}^{11})\cup \cdots \cup \mathrm{eigs}({\mathcal{D}}^{NN}+i{\mathcal{H}}^{NN}).$$

$$\mathrm{Tr}[{\widehat{\text{mathds}1}}_j{\mathcal{D}}_j({\hat{\sigma }}_j^x)]=\mathrm{Tr}[{\widehat{\text{mathds}1}}_j{\mathcal{D}}_j({\hat{\sigma }}_j^y)]=\mathrm{Tr}[{\hat{\sigma }}_j^z{\mathcal{D}}_j({\hat{\sigma }}_j^x)]=\mathrm{Tr}[{\hat{\sigma }}_j^z{\mathcal{D}}_j({\hat{\sigma }}_j^y)]=0.$$

$${\mathcal{D}}_{\mu \nu }=⟮\begin{array}{cccc}[cc|cc]0& 0& 0& 0\\ {\mathcal{D}}_{z1}& {\mathcal{D}}_{zz}& 0& 0\\ & & & \\ {\mathcal{D}}_{x1}& {\mathcal{D}}_{xz}& {\mathcal{D}}_{xx}& {\mathcal{D}}_{xy}\\ {\mathcal{D}}_{y1}& {\mathcal{D}}_{yz}& {\mathcal{D}}_{yx}& {\mathcal{D}}_{yy}\end{array}⟯.$$

$$\mathrm{Tr}[{\hat{\sigma }}_j^y{\mathcal{D}}_j({\hat{\sigma }}_j^x)]=\mathrm{Tr}[{\hat{\sigma }}_j^x{\mathcal{D}}_j({\hat{\sigma }}_j^y)],\mathrm{and}\mathrm{Tr}[{\hat{\sigma }}_j^z{\mathcal{D}}_j({\widehat{\text{mathds}1}}_j)]=0.$$

$$\displaystyle{\rm eigs}(\mathcal{D}^{dd})=\bigcup_{\mathscr{J}_{d}}\bigg{(}\sum_{j\in\mathscr{J}_{d}}\mathscr{E}_{j}^{1}+\sum_{j\notin\mathscr{J}_{d}}\mathscr{E}_{j}^{0}\bigg{)}.$$

$$⟨⟨d,v|{\mathcal{L}}^{dd}|d,v⟩⟩=⟨⟨d,v|{\mathcal{D}}^{dd}|d,v⟩⟩+i⟨⟨d,v|{\mathcal{H}}^{dd}|d,v⟩⟩=x+iy.$$

$$x=⟨⟨d,v|{\mathcal{D}}^{dd}|d,v⟩⟩\le -d\Gamma ,$$

Full text

The understanding of equilibrium many-body physics in both classical and quantum systems has relied heavily on exact solutions of highly simplified models. For example, solutions of the classical and quantum (transverse-field) Ising models have elucidated the structure of classical and quantum phase transitions, respectively, and the deep connections between them Sondhi et al. (1997); Sachdev (2011). For quantum systems that are not in thermal equilibrium such solutions are comparatively scarce, though important progress has been made in numerous specialized models that either: (a) are one dimensional and isolated Calabrese and Cardy (2005); Fagotti and Calabrese (2008), (b) impose nontrivial dissipation on otherwise free bosons or fermions Prosen (2008); Horstmann et al. (2013); Huffman et al. (2016); Žnidarič (2010, 2011); Eisler (2011); Žunkovič (2014); Guo and Poletti (2016a, b), (c) impose highly fine-tuned dissipation Žunkovič (2014), or (d) are coupled to an environment only at a boundary Prosen and Pižorn (2008); Prosen (2011a, b, 2015). On the other hand, interacting quantum spin systems that are dissipative in the bulk and driven so that they do not thermalize are expected to exhibit a variety of rich behaviors not found in their equilibrium counterparts, including unusual multi-critical points Lee et al. (2013), critical exponents not found in or near equilibrium Sieberer et al. (2016), and the existence of zero-entropy entangled steady states Kraus et al. (2008); Diehl et al. (2008); Verstraete et al. (2009); Dengis et al. (2014). And, naturally, such systems play an important role in the theory of quantum computation in the presence of decoherence. While recent experimental advances have enabled the controlled study of such physics in systems ranging from exciton-polariton fluids Deng et al. (2002); Kasprzak et al. (2006); Szymańska et al. (2006); Byrnes et al. (2014); Rodriguez et al. (2016), to trapped ions Bohnet et al. (2016); Schindler et al. (2013), to Rydberg gases Peyronel et al. (2012); Carr et al. (2013); Malossi et al. (2014), the minimal microscopic models expected to capture the essential qualitative physics—many-body quantum master equations Sieberer et al. (2016)—continue to pose severe challenges to existing theoretical techniques.

In this manuscript we investigate a broad and experimentally relevant class of driven-dissipative many-body spin models, some of which can be viewed as dissipative analogues of the transverse-field Ising model (TFIM), and show that they can be solved efficiently. The TFIM captures a characteristic feature of low-temperature quantum systems more generally: even at zero temperature, the ordering associated with a classical Hamiltonian (the Ising model) can be frustrated by the persistence of quantum fluctuations (induced by the transverse field, left panel of Fig. 1). In particular, the transverse field favors a zero-temperature density matrix that does not commute with the Ising Hamiltonian. As a result, energy minimization forces the system to develop quantum correlations even between arbitrarily distant spins, generically prohibiting exact calculations of its ground-state properties (except in certain very well-known circumstances, e.g. nearest-neighbor interactions in 1D). The models we consider realize a natural non-equilibrium analogue of this scenario, in which a classical Hamiltonian is frustrated by the presence of dissipative fluctuations (right panel, Fig. 1). Strikingly, we find that the inclusion of a broad class of dissipative processes favoring a steady-state density matrix that does not commute with the Hamiltonian is relatively benign, even in more than one spatial dimension. In particular, correlations remain dynamically localized whenever the Hamiltonian is of finite range, enabling the time-dependence of finite-order correlation functions to be efficiently obtained (even for an infinite system) by solving a finite-dimensional system of equations. This structure exists, for example, even when the dissipation alone drives the system towards the ground state of a transverse field. It also exists when the dissipation is explicitly derived from a temporally fluctuating transverse field; evidently, even though the Ising model is not (in general) expected to be solvable in the presence of a particular (static or time-varying) transverse field, the dynamics of the transverse-field Ising model averaged over all such fields can be computed exactly.

To illustrate the utility and generality of these solutions, we use them to prove that an experimentally relevant subset of the models considered possess a finite dissipative gap, and therefore cannot undergo dissipative phase transitions. We then apply this result to a model of interacting Rydberg atoms studied in Refs. Maghrebi and Gorshkov (2016); Overbeck et al. (2016), thereby confirming the structure of the phase diagram inferred from recently developed approximate techniques Weimer (2015). From a more applied perspective, we emphasize that Ising-like Hamiltonians, despite being “classical”, induce truly quantum dynamics and play a central role in the production of states sought for quantum information processing and quantum metrology Raussendorf and Briegel (2001); Kitagawa and Ueda (1993); Mølmer and Sørensen (1999). Their solvability in the presence of dissipation thus affords numerous exciting opportunities to investigate the effect of decoherence on the generation of useful entangled states. As an illustrative example, we compute the effects of a fluctuating transverse field (dephasing) on Ising dynamics generated by Mølmer-Sørensen gates Mølmer and Sørensen (1999); Lee et al. (2005); Kim et al. (2010), which can impose important technical limitations in trapped-ion based approaches to quantum computation Gaebler et al. (2016).

*Dissipative models.—*The models we consider can be constructed in close analogy to the transverse-field Ising model. There, one starts with a many-body Hamiltonian that is “classical” in the sense that it can be diagonalized by a choice of local basis. Thus its eigenstates have well defined local properties, and in this sense there are no quantum fluctuations FN (1). To be concrete we restrict our focus to the spin-$\frac{1}{2}$ Ising model,

[TABLE]

However, the scope of our results is significantly broader, and includes generic Hamiltonians that can be diagonalized by a choice of local basis or, with some additional constraints, commuting Hamiltonians that are not diagonalized by a choice of local basis (e.g. the toric code Kitaev (2003)). Extensions along these lines will be discussed at the end of the manuscript. In equilibrium, a natural question to ask is how one can modify $\hat{H}$ such that, in the low temperature limit where thermal fluctuations vanish, quantum fluctuations remain. One simple strategy is to add single-body terms to the Hamiltonian that are not diagonalized by the eigenstates of $\hat{H}$; the usual culprit is a transverse field, resulting in an additional term $B\sum_{j}\hat{\sigma}_{j}^{x}$. Even at zero temperature, i.e. in the quantum ground state, such a term induces fluctuations between the classical eigenstates of the interaction part of the Hamiltonian.

In a driven-dissipative setting, we are no longer interested in properties of the ground state, which control the low-temperature equilibrium physics, but rather the properties of the steady state, which control the long-time, non-equilibrium physics. A natural dissipative generalization of the procedure used above to frustrate the zero-temperature ordering associated with $\hat{H}$ is to introduce single-body dissipative processes that drive the system toward a steady state that does not commute with $\hat{H}$. Assuming that dissipation can be treated in the Born-Markov approximation, the dynamics of an open quantum system with Hamiltonian $\hat{H}$ is governed by a Markovian master equation of the form

[TABLE]

The dissipation is induced by jump operators $\hat{J}_{j\alpha}$, each of which we assume to be supported on a single site $j$. The index $\alpha$ may take on multiple values in order to describe multiple dissipative channels on a given site, though in what follows we will generally consider only one jump operator on each site, dropping the index $\alpha$. Consider the dissipative dynamics in the absence of the Hamiltonian, described by $\dot{\hat{\rho}}=\mathcal{D}(\hat{\rho})$. The steady-state solution of the purely dissipative dynamics is determined implicitly by solving $\mathcal{D}(\hat{\rho}_{\rm dis})=0$. If $\hat{\rho}_{\rm dis}$ commutes with the Hamiltonian, then it is automatically also a solution of $\mathcal{L}(\hat{\rho})=0$ and thus is a proper steady state of the complete dynamics including both coherent evolution and dissipation. If it does not, then in close analogy to the ground state of the transverse-field Ising model, we expect the steady state to possess fluctuations in the sense of admixing of classical states.

Surprisingly, a large class of dissipators that frustrate the classical Hamiltonian, i.e. for which $[\hat{\rho}_{\rm dis},\hat{H}]\neq 0$, nevertheless admit exact solutions for the dynamics of observables. In particular, suppose that FN (2)

[TABLE]

For $\hat{H}$ of the form in Eq. (1), we also have ${\rm Tr}(\hat{\sigma}^{z}_{j}[\hat{H},\hat{\sigma}^{\pm}_{j}])=0$, and therefore whenever the jump operators are chosen to satisfy Eq. (3), the complete Liouvillian also obeys

[TABLE]

Taken together with the equality ${\rm Tr}(\hat{\mathds{1}}_{j}\mathcal{L}(\hat{\sigma}^{\pm}_{j}))=0$, which follows trivially from the definition of $\mathcal{L}$ (and ensures the conservation of probability) regardless of the form of the jump operators, Eq. (4) can be understood colloquially as a statement that $\mathcal{L}$ does not map coherences onto populations (Fig. 2a). In what follows, we will show that finite-range Hamiltonians of the form in Eq. (1), subjected to dissipation obeying Eq. (3), produce correlations that are localized in space and can therefore be solved efficiently in the thermodynamic limit.

There are many natural jump operators satisfying Eq. (3). For example, dephasing, spontaneous emission, and incoherent pumping along the $z$ axis all do, and the models studied here therefore subsume finite-range versions of the models studied in Foss-Feig et al. (2013a, b) (and realized experimentally in Bohnet et al. (2016)) as special cases. However, all the examples of dissipation just mentioned lead to steady states $\hat{\rho}_{\rm dis}$ that commute with $\hat{H}$, and thus the steady states in the presence of $\hat{H}$ remain trivial. Examples of jump operators satisfying Eq. (3) but producing steady states that do not commute with $\hat{H}$ include dephasing along any direction in the $xy$ plane; $\hat{J}=\cos(\theta)\hat{\sigma}^{x}+\sin(\theta)\hat{\sigma}^{y}$, or spontaneous emission along any axis in the $xy$ plane; $\hat{J}=\hat{\sigma}^{z}+i(\cos(\theta)\hat{\sigma}^{x}+\sin(\theta)\hat{\sigma}^{y})$. Note that in the latter example, $\hat{\rho}_{\rm dis}$ is the ground state of a transverse field.

Localization of correlations.—The time-dependent expectation value of an arbitrary operator $\hat{O}$, initially supported on a set of sites $\mathscr{A}$, can be written as

[TABLE]

The Heisenberg picture expression utilizes the adjoint Liouvillian $\mathcal{L}^{\ddagger}(\hat{O})=i[\hat{H},\hat{O}]+\mathcal{D}^{\ddagger}(\hat{O})$, with the adjoint dissipator given by \mathcal{D}^{\ddagger}(\star)=\sum_{j,\alpha}\frac{\gamma_{j\alpha}}{2}\big{(}2\hat{J}^{\dagger}_{j\alpha}\star\hat{J}^{\phantom{\dagger}}_{j\alpha}-\{\hat{J}^{\dagger}_{j\alpha}\hat{J}^{\phantom{\dagger}}_{j\alpha},\star\}\big{)}. The operator $\hat{O}$ could be, for example, a product of two spin operators on different sites, in which case $O(t)$ is an equal-time correlation function (generalizations to unequal time correlation functions are also possible using the quantum regression formula). Equation (4) can be recast in terms of the Heisenberg-picture Liouvillian as ${\rm Tr}[\hat{\sigma}^{\pm}_{j}\mathcal{L}^{\ddagger}(\hat{\sigma}^{z}_{j})]=0$; thus, $\mathcal{L}^{\ddagger}$ does not map populations onto coherences (Fig. 2a). This condition on $\mathcal{L}^{\ddagger}$ can be restated in the following useful way: Referring to an operator as “diagonal on the set of sites $\mathscr{S}$” if it commutes with all operators $\hat{\sigma}_{j}^{z}$ for $j\in\mathscr{S}$, one can show that the set of operators that are diagonal on $\mathscr{S}$ is closed under the action of $\mathcal{L}^{\ddagger}$. With this structure in mind, we refer to $\mathcal{L}^{\ddagger}$ as diagonality preserving (in the $z$ basis), a property which will play a key role in the solution of Eq. (2).

To compute $O(t)$, it is useful to define the set of sites that are nearest neighbors of $\mathscr{A}$, denoted $\mathscr{B}$, as all sites outside of $\mathscr{A}$ for which a term in the Hamiltonian has simultaneous support on $\mathscr{A}$ and $\mathscr{B}$. The complement of $\mathscr{A}\bigcup\mathscr{B}$, denoted $\mathscr{C}$, is therefore not directly connected to $\mathscr{A}$ by terms in the Hamiltonian (see Fig. 2b). We then decompose the Hamiltonian as $\hat{H}=\hat{H}_{\mathscr{A}\!\mathscr{B}}+\hat{H}_{V}$, where $\hat{H}_{\mathscr{A\!B}}$ contains all terms in $\hat{H}$ that have support on $\mathscr{A}$, and $\hat{H}_{V}$ contains all terms that do not (note that, by the definition of $\mathscr{B}$, $\hat{H}_{\mathscr{A\!B}}$ is supported on $\mathscr{A}\bigcup\mathscr{B}$, while $\hat{H}_{V}$ is supported on $\mathscr{B}\bigcup\mathscr{C}$). Similarly, we decompose the Heisenberg-picture dissipator as $\mathcal{D}^{\ddagger}=\mathcal{D}^{\ddagger}_{\mathscr{A\!B}}+\mathcal{D}^{\ddagger}_{V}$, where $\mathcal{D}^{\ddagger}_{\mathscr{A\!B}}$ only contains jump operators supported on $\mathscr{A}\bigcup\mathscr{B}$ and $\mathcal{D}_{V}^{\ddagger}$ only contains jump operators supported on $\mathscr{C}$. Writing $\mathcal{L}^{\ddagger}=\mathcal{L}^{\ddagger}_{\mathscr{A\!B}}+\mathcal{L}^{\ddagger}_{V}$, with $\mathcal{L}^{\ddagger}_{\mathscr{A\!B}}(\star)=i[\hat{H}_{\mathscr{A\!B}},\star]+\mathcal{D}_{\mathscr{A\!B}}^{\ddagger}(\star)$ and $\mathcal{L}^{\ddagger}_{V}(\star)=i[\hat{H}_{V},\star]+\mathcal{D}^{\ddagger}_{V}(\star)$, we can write $O(t)$ as

[TABLE]

where $\hat{O}_{n}=(\mathcal{L}^{\ddagger}_{\mathscr{A\!B}}+\mathcal{L}^{\ddagger}_{V})^{n}\hat{O}$. The structure of Eq. (6) can be greatly simplified by induction. In particular, suppose that $\hat{O}_{n}$ satisfies the following two conditions:

(A) $\hat{O}_{n}$ is supported on $\mathscr{A}\bigcup\mathscr{B}$.

(B) $\hat{O}_{n}$ is diagonal on $\mathscr{B}$.

Note that both conditions are satisfied trivially for $\hat{O}_{0}=\hat{O}$, since we have assumed that $\hat{O}$ is fully supported on $\mathscr{A}$. We will show that the operator $\hat{O}_{n+1}=(\mathcal{L}^{\ddagger}_{\mathscr{A\!B}}+\mathcal{L}^{\ddagger}_{V})\hat{O}_{n}$ also satisfies these conditions, thereby proving that all $\hat{O}_{n}$ do.

Toward this end, we first show that $\mathcal{L}^{\ddagger}_{V}(\hat{O}_{n})=i[\hat{H}_{V},\hat{O}_{n}]+\mathcal{D}^{\ddagger}_{V}(\hat{O}_{n})$ vanishes by the assumptions. The first term vanishes because $\hat{H}_{V}$ is diagonal and supported on $\mathscr{B}\bigcup\mathscr{C}$, while $\hat{O}_{n}$ is diagonal on $\mathscr{B}\bigcup\mathscr{C}$ by conditions (A,B), giving $[\hat{H}_{V},\hat{O}_{n}]=0$. The second term vanishes because $\mathcal{D}^{\ddagger}_{V}$ only contains jump operators supported on $\mathscr{C}$, implying (together with condition A) that $\mathcal{D}^{\ddagger}_{V}(\hat{O}_{n})=0$. Thus $\mathcal{L}^{\ddagger}_{V}(\hat{O}_{n})=0$ as claimed, giving $\hat{O}_{n+1}=\mathcal{L}^{\ddagger}_{\mathscr{A\!B}}(\hat{O}_{n})$. Now we note that, because $\hat{H}_{\mathscr{A\!B}}$ and the jump operators in $\mathcal{D}^{\ddagger}_{\mathscr{A\!B}}$ are supported on $\mathscr{A}\bigcup\mathscr{B}$, $\mathcal{L}^{\ddagger}_{\mathscr{A\!B}}(\hat{O}_{n})$ must also be supported on $\mathscr{A}\bigcup\mathscr{B}$, and $\hat{O}_{n+1}$ therefore satisfies condition (A). Furthermore, because $\mathcal{L}^{\ddagger}_{\mathscr{A}\!\mathscr{B}}$ is diagonality preserving, $\hat{O}_{n+1}$ continues to satisfy condition (B). Therefore, by induction, we deduce that $\hat{O}_{n}=(\mathcal{L}^{\ddagger}_{\mathscr{A\!B}})^{n}\hat{O}$, and is supported on $\mathscr{A}\bigcup\mathscr{B}$. Plugging this result into Eq. (6) and using the isolated support of $\hat{O}_{n}$ to carry out the trace over $\mathscr{C}$, we obtain

[TABLE]

where $\hat{\rho}_{\scriptscriptstyle\!\mathscr{A}\!\mathscr{B}}={\rm Tr}_{\mathscr{C}}(\hat{\rho}_{0})$ is the reduced density matrix of the initial state upon tracing out all sites in $\mathscr{C}$. Returning to the Schrödinger picture we find our primary result,

[TABLE]

Equation (8) can be read as follows: The time-dependent expectation value of $\hat{O}$ can be computed by first calculating the reduced density matrix on $\mathscr{A}\bigcup\mathscr{B}$, and then evolving it using a Hamiltonian and jump operators supported only on $\mathscr{A}\bigcup\mathscr{B}$. It follows immediately that correlations are localized for finite-ranged Hamiltonians. For example, consider a correlation function $\langle\hat{O}\rangle=\langle\hat{\sigma}^{\mu}_{j}\hat{\sigma}^{\nu}_{k}\rangle$ ($\mu,\nu\in\{x,y,z\}$), and suppose the system starts in a product state. If the Hamiltonian is of finite range $r$, and $j$ and $k$ are separated by a distance $d_{jk}>2r$ (so that sites $j$ and $k$ do not share any neighbors, and $\mathscr{A}\bigcup\mathscr{B}$ decomposes into two disjoint sets), then it follows from Eq. (8) that $\langle\hat{\sigma}^{\mu}_{j}\hat{\sigma}^{\nu}_{k}\rangle=\langle\hat{\sigma}^{\mu}_{j}\rangle\langle\hat{\sigma}^{\nu}_{k}\rangle$ at all times FN (3).

Applications.—Just as many-body ground states can suddenly change in character as a parameter in the Hamiltonian is continuously adjusted, signaling a quantum phase transition, the steady state of a master equation can exhibit a sudden change when a parameter in the Liouvillian is adjusted continuously, signaling a dissipative phase transition. While a quantum phase transition is associated with the closing of an energy gap in the spectrum of the Hamiltonian, a dissipative phase transition is associated with the closing of a dissipative gap in the spectrum of the Liouvillian. Consider the dissipative transverse-field Ising model studied in Refs. Maghrebi and Gorshkov (2016); Overbeck et al. (2016),

[TABLE]

in which fluctuations due to both dissipation (strength $\gamma$) and a transverse field (strength $\Delta$) are considered; note that the state which minimizes the energy with respect to the transverse field (for $\Delta>0$) is identical to the steady state of the dissipation alone. In Eq. (9) we have chosen to make the Hamiltonian diagonal along the $x$ direction for consistency with Ref. Overbeck et al. (2016). However, for $\Delta=0$ and upon making a rotation of the coordinate system that takes the Ising part of $\hat{H}_{\rm TFIM}$ into the form of Eq. (1), it can easily be verified that this model satisfies Eq. (3). In Ref. Overbeck et al. (2016), a careful analysis of the variational techniques developed in Weimer (2015) suggests that the system is disordered for any value of $\gamma$ at $\Delta=0$ and in high dimensions. Here, we can make this conclusion rigorous in any dimension: In the supplement, we prove that this model at $\Delta=0$ [and also a more general subclass of the models obeying Eq. (3)] has a finite dissipative gap $\geq\gamma$. Since the model is disordered at $\gamma=\infty$ and a gap exists at $\Delta=0$ for any $\gamma>0$, it must be disordered along the entire $\Delta=0$ axis. In light of the Liouvillian stability results proved in Ref. Cubitt et al. (2015), the existence of a dissipative gap at $\Delta=0$ may actually imply something much stronger: the gap persists for small $\Delta>0$, ruling out a phase transition for $\Delta\lesssim\gamma$ and implying that the system is disordered for sufficiently small $\Delta$ at any finite $\gamma$.

The solutions developed here can also be used to calculate dynamics of the fluctuating-transverse-field Ising model,

[TABLE]

Here, the transverse fields $\Delta_{j}(t)$ are Gaussian random variables with white-noise spectrum $\overline{\Delta_{j}(t_{1})\Delta_{k}(t_{2})}=\gamma\delta_{j,k}\delta(t_{1}-t_{2})$. This model—interpreted such that observables are obtained by averaging dynamics induced by $\hat{H}_{\rm fTFIM}(t)$ over realizations of the fluctuating transverse field—can be viewed as a continuous unraveling of the master equation describing dephasing dynamics in the Ising model Gisin and Percival (1993); Barchielli and Gregoratti (2009),

[TABLE]

[See also Refs. Budini (2001); Chenu et al. (2016) for a constructive approach to Eq. (11) starting from Eq. (10)]. A master equation of this form arises naturally in ion trap experiments, where it captures the effects of qubit dephasing during Ising dynamics induced by Mølmer-Sørensen gates Mølmer and Sørensen (1999); Lee et al. (2005); Kim et al. (2010); Gaebler et al. (2016), and hence it plays a prominent role in the description of decoherence effects on trapped-ion based approaches to quantum computation. Naively one would not expect there to be an efficient solution to Eq. (11), given that its unraveling [either continuous or discrete, with the continuous case leading back to Eq. (10)] requires one to solve the Ising model in the presence of a time-dependent transverse field. Nevertheless, after a rotation by $\pi/2$ about the $y$ axis Eq. (11) falls into the class of solvable models studied here, and observables can be computed efficiently. In Fig. 3(a), we show simulations of the stochastic Schrödinger equation associated with $\hat{H}_{\rm fTFIM}$; while these simulations are computationally very expensive, and not feasible for more than $\sim 10$ spins, the average over sufficiently many trajectories (lighter to darker red lines indicating more averaging) can be seen to converge to the exact solution (black-dashed line).

Outlook.— We have identified an exactly solvable class of driven-dissipative quantum spin models in which dissipation frustrates the ordering of a classical Hamiltonian. The underlying algebraic structure that enables an exact solution, ensured by Eq. (3), can also be identified in a much more general class of Hamiltonians that are solvable in the absence of dissipation. For example, the results can easily be generalized to systems with arbitrary finite-dimensional Hilbert spaces on each site and time-dependent Hamiltonians. Long-range interactions do spoil some aspects of the solution, but the absence of dissipative phase transitions can still be proven. Hamiltonians that are not diagonalized by a local choice of basis but can still be written as a sum of local commuting terms, such as the toric code Kitaev (2003); Dengis et al. (2014), are also amenable to similar solution techniques, and will be the topic of future research.

Supplemental Material

This supplemental material contains a proof that, under conditions specified in the manuscript and reiterated below, the Liouvillians we consider possesses a finite dissipative gap. In particular, the dissipative gap of the Liouvillian $\mathcal{L}$ is bounded below by that of the dissipator $\mathcal{D}$ in the absence of a Hamiltonian. The first section briefly introduces relevant formalism, the second section discusses the structure of the Liovillians considered in the manuscript in the language of this formalism, and the final section contains the proof.

S0.1 Parametrization of density matrices and superoperators

An arbitrary density matrix of $N$ spin-$\frac{1}{2}$’s can be expanded as

[TABLE]

where the coefficients $\rho_{\mu_{1},\dots,\mu_{N}}$ are real numbers and $\hat{\sigma}^{1}_{j}=\hat{\mathds{1}}_{j}$ is the identity operator on site $j$. For convenience we also introduce a more compact notation in which the set of indices $\{\mu_{1}\dots\mu_{N}\}$ is replaced by a single collective index $\mu$ (which runs over all possible assignments of values $1,x,y,z$ to all $N$ sites), the basis element $\hat{\sigma}_{1}^{\,\mu_{1}}\otimes\dots\otimes\hat{\sigma}_{N}^{\,\mu_{N}}$ is denoted simply by $\hat{\mu}$, and the density operator $\hat{\rho}$ (or any other operator) is denoted with a modified ket notation as $|\rho\rangle\!\rangle$. In this notation, Eq. (S1) becomes $|\rho\rangle\!\rangle=\sum_{\mu}\rho_{\mu}|\mu\rangle\!\rangle$. The $\rho_{\mu}$ must be real because the density matrix is Hermitian, and the basis elements all are Hermitian, but it is nevertheless useful to consider complex coefficients. In this notation a general (not-necessarily Hermitian) operator $\hat{A}$ and its Hermitian conjugate $\hat{A}^{\dagger}$ become, respectively, $|A\rangle\!\rangle=\sum_{\mu}A_{\mu}|\mu\rangle\!\rangle$ and $\langle\!\langle A|=\sum_{\mu}A^{*}_{\mu}\langle\!\langle\mu|$. The basis elements $|\mu\rangle\!\rangle$ form a vector space (Liouville space) endowed with the inner product $\langle\!\langle A|B\rangle\!\rangle\equiv 2^{-N}{\rm Tr}(\hat{A}^{\dagger}\hat{B})=2^{-N}\sum_{\mu\nu}A^{*}_{\mu}B_{\nu}{\rm Tr}(\hat{\mu}\hat{\nu})=\sum_{\mu}A^{*}_{\mu}B_{\mu}$, where we have used the fact that $\hat{\mu}^{\dagger}=\hat{\mu}$ and ${\rm Tr}(\hat{\mu}\hat{\nu})=2^{N}\delta_{\mu,\nu}$. We denote the action of a superoperator on an operator, $\hat{B}=\mathcal{S}(\hat{A})$, with the notation $|B\rangle\!\rangle=\mathcal{S}|A\rangle\!\rangle$. The superoperator $\mathcal{S}$ inherets a matrix representation from the equation $|B\rangle\!\rangle=\mathcal{S}|A\rangle\!\rangle$ by expanding $|A\rangle\!\rangle=\sum_{\nu}A_{\nu}|\nu\rangle\!\rangle$ and left multiplying both sides of the equation by $\langle\!\langle\mu|$, giving $B_{\mu}=\sum_{\nu}\mathcal{S}_{\mu\nu}A_{\nu}$, where $\mathcal{S}_{\mu\nu}\equiv\langle\!\langle\mu|\mathcal{S}|\nu\rangle\!\rangle$. Note that a physical Liouvillian must be hermiticity preserving, which implies that its matrix elements in a basis of Hermitian operators (such as the operators $\hat{\mu}$ defined above) are always real.

S0.2 Structure of the Liouvillian

For our purposes below, an important property of each basis element $|\mu\rangle\!\rangle$ is the number of operators $\hat{\sigma}^{x}_{j}$ and $\hat{\sigma}^{y}_{j}$ that appear [see Eq. (S1)], which we denote by $d$. In literature on nuclear-magnetic resonance Liouville space can be decomposed into a direct sum of subspaces with fixed $d$, and we denote the basis elements of these subspaces by $|d,\mu_{d}\rangle\!\rangle$, with the understanding that the index $\mu_{d}$ enumerates only the basis elements in this subspace. An arbitrary vector within this subspace is denoted by $|d,v\rangle\!\rangle=\sum_{\mu_{d}}v_{\mu_{d}}|d,\mu_{d}\rangle\!\rangle$. If we define $\mathcal{P}_{d}$ as a superoperator that projects onto the subspace spanned by $|d,\mu_{d}\rangle\!\rangle$, then we can expand $\mathcal{S}=\sum_{d,d^{\prime}}\mathcal{P}_{d}\mathcal{S}\mathcal{P}_{d^{\prime}}\equiv\sum_{d,d^{\prime}}\mathcal{S}^{dd^{\prime}}$, which has a matrix representation $\mathcal{S}^{dd^{\prime}}_{\mu_{d}\nu_{d^{\prime}}}=\langle\!\langle d,\mu_{d}|\mathcal{S}^{dd^{\prime}}|d^{\prime},\nu_{d^{\prime}}\rangle\!\rangle$.

The Liouvillian superoperator can be decomposed as $\mathcal{L}=\mathcal{D}+i\mathcal{H}$, where $\mathcal{H}(\star)=-[\hat{H},\star]$. Because $\mathcal{D}$ and $i\mathcal{H}$ are each (independently) hermiticity preserving, they both have real-valued matrix representations in the chosen basis. The dissipators we consider in the manuscript obey the condition {\rm Tr}\big{(}\hat{\sigma}^{z}_{j}\mathcal{D}(\hat{\sigma}^{\pm}_{j})\big{)}={\rm Tr}\big{(}\hat{\mathds{1}}_{j}\mathcal{D}(\hat{\sigma}^{\pm}_{j})\big{)}=0 (Eq. (3) of the manuscript), which in the language used here implies that $\mathcal{D}^{dd^{\prime}}=0$ whenever $d<d^{\prime}$. In fact, because the dissipator is constructed of local terms, it is straightforward to see that $\mathcal{D}^{dd^{\prime}}$ vanishes unless either $d=d^{\prime}$ or $d=d^{\prime}+1$. It is also readily verified that $\mathcal{H}$ is block diagonal in the $d$-subspace decomposition, and that $\mathcal{H}^{00}=0$ (the latter condition following because $\hat{H}$ commutes with all basis elements in the $d=0$ sector). Hence the complete Liouvillian has the block-lower-triangular structure (dropping subscripts)

[TABLE]

Note that the block-lower-triangular structure of $\mathcal{L}$ of ensures that the time evolution within the $d=0$ subspace, which completely determines the populations in the $z$-basis, is governed by a closed set of equations. The dynamics within this sector is therefore equivalent to the dynamics of a classical master equation, a situation which sometimes arises in more general quantum systems whenever it is possible to eliminate the contribution of coherences from the equations of motion for populations Fischer et al. (2016); Hoening et al. (2014). The eigenvalues of a block-lower triangular matrix are given by the eigenvalues of the diagonal blocks, and therefore

[TABLE]

Here it is to be understood that the eigenvalues of a projected operator are computed only with respect to vectors supported on that subspace, i.e. the eigenvalues of $\mathcal{D}^{00}$ are given by the eigenvalues of the matrix $\mathcal{D}^{00}_{\mu_{0},\nu_{0}}$.

S0.3 Dissipative gap

The dissipators we consider are constructed as sums over single-site dissipators, $\mathcal{D}=\sum_{j}\mathcal{D}_{j}$, each of which independently obeys the constraints (equivalent to Eq. (3) in the manuscript)

[TABLE]

Here it is understood that $\mathcal{D}_{j}$ is an operator in the full Liouville space, but because it only acts nontrivially on site $j$ we will (in a slight abuse of notation) temporarily use the symbol $\mathcal{D}_{j}$ to represent an operator acting on just the Liouville space associated with site $j$, spanned by vectors $|\mu_{j}\rangle\!\rangle$; regardless of the interpretation of which space $\mathcal{D}_{j}$ acts on, its eigenvalues are the same up to degeneracies. With this interpretation in mind, we can represent $\mathcal{D}_{j}$ as a $4\times 4$ real matrix matrix $\langle\!\langle\mu_{j}|\mathcal{D}_{j}|\nu_{j}\rangle\!\rangle$. Using the shorthand $\mathcal{D}_{\mu\nu}\equiv\langle\!\langle\mu_{j}|\mathcal{D}_{j}|\nu_{j}\rangle\!\rangle$ (note the temporary suppression of the index $j$), ordering our basis as $\mu=1,z,x,y$, and taking into account the constraints on $\mathcal{D}_{j}$ imposed in Eq. (S4), we have

[TABLE]

In order to prove the existence of a dissipative gap for $\mathcal{L}$, we also make the following two assumptions:

(1) Each local dissipator $\mathcal{D}_{j}$ has a unique steady state and a dissipative gap $\Gamma_{j}$, and ${\rm min}_{j}\Gamma_{j}=\Gamma>0$. Since all eigenvalues of the dissipator must have non-positive real parts, this means that for any nonzero eigenvalue $\lambda$ of $\mathcal{D}_{j}$, we must have ${\rm Re}(\lambda)\leq-\Gamma_{j}$.

(2) In addition to Eq. (S4), the local dissipators obey the further constraints:

[TABLE]

In the notation just introduced, we can restate these conditions as: (1) The matrix $\mathcal{D}_{\mu\nu}$ has a unique zero eigenvalue and a dissipative gap $\Gamma_{j}$, and (2) The diagonal blocks of $\mathcal{D}_{\mu\nu}$ are symmetric ($\mathcal{D}_{z1}=0$ and $\mathcal{D}_{yx}=\mathcal{D}_{xy}$). Note that the Hamiltonian in Eq. (9) of the manuscript, when rotated about the $y$ axis and at $\Delta=0$, becomes of the form in Eq. (1) of the manuscript, while the jump operators become $\hat{J}_{j}=\hat{\sigma}_{j}^{y}-i\hat{\sigma}_{j}^{z}$. It is straightforward to verify that conditions (1) and (2) are satisfied for the corresponding dissipator, with $\Gamma=\gamma$.

As mentioned above, a block-lower-triangular matrix of this form shares its eigenvalues with the diagonal blocks. Let’s denote the set of eigenvalues of the upper-left block by $\mathscr{E}_{j}^{0}=\{0,\varepsilon^{0}\}$, and the set of eigenvalues of the lower-right block by $\mathscr{E}_{j}^{1}=\{\varepsilon^{1}_{a},\varepsilon^{1}_{b}\}$ (it is straightforward to see that the upper-left block must have a zero eigenvalue). By conditions (1) and (2), we know that the three non-zero eigenvalues are real, and that $\varepsilon^{0},\varepsilon^{1}_{a},\varepsilon^{1}_{b}\leq-\Gamma_{j}$. Because each $\mathcal{D}_{j}$ is supported on a single site, the eigenvalues of the full dissipator $\mathcal{D}=\sum_{j}\mathcal{D}_{j}$ are simply sums of eigenvalues from each $\mathcal{D}_{j}$. Eigenvalues of the projected dissipator $\mathcal{D}^{dd}$ can be constructed as follows: Defining $\mathscr{J}_{d}$ to be a set containing $d$ sites, and also defining the notation $\mathscr{A}+\mathscr{B}=\{a+b:~{}a\in\mathscr{A},~{}b\in\mathscr{B}\}$, we have

[TABLE]

We can draw two immediate conclusions:

(A) $\mathcal{D}^{00}$ has precisely one zero eigenvalue, obtained when the zero eigenvalue in $\mathscr{E}_{j}^{0}$ is chosen for all $j$, and any other eigenvalue $\varepsilon$ satisfies $\varepsilon\leq-\Gamma$.

(B) Any eigenvalue $\varepsilon$ of $\mathcal{D}^{dd}$ (with $d>0$) is real and satisfies $\varepsilon<-d\Gamma$.

Now consider an arbitrary normalized eigenvector of $\mathcal{L}^{dd}$ in the $d\neq 0$ subspace, $|d,v\rangle\!\rangle$, with eigenvalue $\varepsilon=x+iy$, and write

[TABLE]

It is straightforward to show that the superoperator $\mathcal{H}^{dd}$ is Hermitian (for any Hamiltonian), and thus $x={\rm Re}(\langle\!\langle d,v|\mathcal{D}^{dd}|d,v\rangle\!\rangle)$. In general, there is not a simple relation between the real part of the expectation value of an operator and the real parts of its eigenvalues. However, assumption (2) guarantees that $\mathcal{D}^{dd}$ is a real symmetric matrix. Therefore, we have

[TABLE]

where the inequality follows because the expectation value of a real symmetric matrix (more generally a Hermitian matrix) is bounded between its smallest and largest eigenvalues. Returning to Eq. (S3), keeping in mind conclusion (A), and using Eq. (S9), we find that $\mathcal{L}$ has a single zero eigenvalue, and a dissipative gap of at least $\Gamma$.