Multistability of Driven-Dissipative Quantum Spins

TL;DR

This paper investigates the multistability phenomena in driven-dissipative quantum spin lattices, revealing how correlations influence steady states and proposing a scheme to go beyond meanfield approximations, with implications for experiments.

Contribution

It introduces an efficient method to include correlations beyond meanfield and demonstrates the impact of correlations on multistability in quantum spin models.

Findings

Correlations eliminate meanfield bistability in 1D, leading to a unique steady state.

Multistability persists in 2D and higher dimensions when considering the thermodynamic limit.

Parameter variations cause jumps between steady states with critical slowing down.

Abstract

We study the dynamics of lattice models of quantum spins one-half, driven by a coherent drive and subject to dissipation. Generically the meanfield limit of these models manifests multistable parameter regions of coexisting steady states with different magnetizations. We introduce an efficient scheme accounting for the corrections to meanfield by correlations at leading order, and benchmark this scheme using high-precision numerics based on matrix-product-operators in one- and two-dimensional lattices. Correlations are shown to wash the meanfield bistability in dimension one, leading to a unique steady state. In dimension two and higher, we find that multistability is again possible, provided the thermodynamic limit of an infinitely large lattice is taken first with respect to the long time limit. Variation of the system parameters results in jumps between the different steady states, each showing a critical slowing down in the convergence of perturbations towards the steady state. Experiments with trapped ions can realize the model and possibly answer open questions in the nonequilibrium many-body dynamics of these quantum systems, beyond the system sizes accessible to present numerics.

Equations

$${\partial }_t\rho =\hat{\mathcal{L}}[\rho ]\equiv -i[H,\rho ]+\mathcal{D}[\rho ],\hslash =1.$$

$$H=\sum _R\left[\frac{\Delta }{2}{\sigma }_R^z+\Omega {\sigma }_R^x\right]-\sum _{⟨R,R^{\prime }⟩}J\left({\sigma }_R^{+}{\sigma }_{R^{\prime }}^{-}+h.c.\right),$$

$$\mathcal{D}[\rho ]=\sum _R\left[{\sigma }_R^{-}\rho {\sigma }_R^{+}-\frac{1}{2}\left({\sigma }_R^{+}{\sigma }_R^{-}\rho +\rho {\sigma }_R^{+}{\sigma }_R^{-}\right)\right].$$

$${\partial }_t{\mu }_x$$

$${\partial }_t{\mu }_y$$

$${\partial }_t{\mu }_z$$

$$\begin{gathered} {\eta }_{ab}(R,R^{\prime },t)\equiv ⟨\left({\sigma }_R^a-{\mu }_a\right)\left({\sigma }_{R^{\prime }}^b-{\mu }_b\right)⟩ \\ =⟨{\sigma }_R^a{\sigma }_{R^{\prime }}^b⟩-{\mu }_a{\mu }_b,R\ne R^{\prime }, \end{gathered}$$

$$H=\sum _R\left[\frac{\Delta }{2}{\sigma }_R^z+\Omega {\sigma }_R^x\right]-\sum _{⟨R,R^{\prime }⟩}\frac{J_z}{2}{\sigma }_R^z{\sigma }_{R^{\prime }}^z,$$

$${\vartheta }_{ab}(R,R^{\prime },t)\equiv ⟨{\sigma }_R^a{\sigma }_{R^{\prime }}^b⟩,R\ne R^{\prime },$$

$${\eta }_{ab}(R,R^{\prime },t)={\vartheta }_{ab}(R,R^{\prime })-{\mu }_a{\mu }_b,$$

$${\zeta }_{abc}(R,R^{\prime },R^{\prime \prime },t)\equiv ⟨\left({\sigma }_R^a-{\mu }_a\right)\left({\sigma }_{R^{\prime }}^b-{\mu }_b\right)\left({\sigma }_{R^{\prime \prime }}^c-{\mu }_c\right)⟩,$$

$$\begin{gathered} {\partial }_t{\vartheta }_{ab}(R)= \\ \sum _d{\Pi }_{ad}{\vartheta }_{db}(R)+\sum _d{\Pi }_{bd}{\vartheta }_{ad}(R)+f_{ab}(\mu ,\vartheta )+g_{ab}(\mu ,\vartheta ), \end{gathered}$$

$$\Pi=\left(\begin{array}[]{ccc}0&-\Delta&0\\ \Delta&0&-2{\Omega}\\ 0&2{\Omega}&0\end{array}\right),$$

$${\partial }_t{\eta }_{ab}(R,t)={\partial }_t{\vartheta }_{ab}(R)-{\partial }_t\left[{\mu }_a{\mu }_b\right],$$

$$-\sum _{⟨R,R^{\prime }⟩}J\left({\sigma }_R^{+}{\sigma }_{R^{\prime }}^{-}+h.c.\right)=-\sum _{⟨R,R^{\prime }⟩}\frac{J}{2}\left({\sigma }_R^x{\sigma }_{R^{\prime }}^x+{\sigma }_R^y{\sigma }_{R^{\prime }}^y\right),$$

$$f_{xx}(R)=2J\left[2{\mu }_x{\mu }_y{\mu }_z-{\mu }_x{\vartheta }_{yz}(1)-{\mu }_y{\vartheta }_{xz}(R)\right]\left[\mathcal{Z}-{\delta }_{\Vert R\Vert ,1}\right]-2J\sum _{\begin{array}{c}{R}^{\prime }\ne 0\\ \Vert R^{\prime }-R\Vert =1\end{array}}{\mu }_z{\vartheta }_{xy}(R^{\prime }),$$

$$f_{yy}(R)=-2J\left[2{\mu }_x{\mu }_y{\mu }_z-{\mu }_y{\vartheta }_{xz}(1)-{\mu }_x{\vartheta }_{yz}(R)\right]\left[\mathcal{Z}-{\delta }_{\Vert R\Vert ,1}\right]+2J\sum _{\begin{array}{c}{R}^{\prime }\ne 0\\ \Vert R^{\prime }-R\Vert =1\end{array}}{\mu }_z{\vartheta }_{xy}(R^{\prime }),$$

$$f_{zz}(R)=-2J\left[{\mu }_x{\vartheta }_{yz}(R)-{\mu }_y{\vartheta }_{xz}(R)\right]\left[\mathcal{Z}-{\delta }_{\Vert R\Vert ,1}\right]-2J\sum _{\begin{array}{c}{R}^{\prime }\ne 0\\ \Vert R^{\prime }-R\Vert =1\end{array}}\left[{\mu }_y{\vartheta }_{xz}(R^{\prime })-{\mu }_x{\vartheta }_{yz}(R^{\prime })\right].$$

$$\begin{gathered} f_{xy}(R)=J\left[2{\mu }_y^2{\mu }_z-2{\mu }_x^2{\mu }_z-{\mu }_y{\vartheta }_{yz}(1)-{\mu }_y{\vartheta }_{yz}(R)+{\mu }_x{\vartheta }_{xz}(1)+{\mu }_x{\vartheta }_{xz}(R)\right]\left[\mathcal{Z}-{\delta }_{\Vert R\Vert ,1}\right] \\ -J\sum _{\begin{array}{c}{R}^{\prime }\ne 0\\ \Vert R^{\prime }-R\Vert =1\end{array}}\left[{\mu }_z{\vartheta }_{yy}(R^{\prime })-{\mu }_z{\vartheta }_{xx}(R^{\prime })\right], \end{gathered}$$

$$\begin{gathered} f_{xz}(R)=-J{\mu }_y{\delta }_{\Vert R\Vert ,1}+J\left[2{\mu }_z^2{\mu }_y-{\mu }_z{\vartheta }_{yz}(1)-{\mu }_y{\vartheta }_{zz}(R)-{\mu }_x{\vartheta }_{xy}(R)+{\mu }_y{\vartheta }_{xx}(R)\right]\left[\mathcal{Z}-{\delta }_{\Vert R\Vert ,1}\right] \\ -J\sum _{\begin{array}{c}{R}^{\prime }\ne 0\\ \Vert R^{\prime }-R\Vert =1\end{array}}\left[{\mu }_z{\vartheta }_{yz}(R^{\prime })+{\mu }_y{\vartheta }_{xx}(R^{\prime })-{\mu }_x{\vartheta }_{xy}(R^{\prime })\right], \end{gathered}$$

$$\begin{gathered} f_{yz}(R)=J{\mu }_x{\delta }_{\Vert R\Vert ,1}+J\left[-2{\mu }_z^2{\mu }_x+{\mu }_z{\vartheta }_{xz}(1)+{\mu }_x{\vartheta }_{zz}(R)-{\mu }_x{\vartheta }_{yy}(R)+{\mu }_y{\vartheta }_{xy}(R)\right]\left[\mathcal{Z}-{\delta }_{\Vert R\Vert ,1}\right] \\ +J\sum _{\begin{array}{c}{R}^{\prime }\ne 0\\ \Vert R^{\prime }-R\Vert =1\end{array}}\left[{\mu }_z{\vartheta }_{xz}(R^{\prime })-{\mu }_y{\vartheta }_{xy}(R^{\prime })+{\mu }_x{\vartheta }_{yy}(R^{\prime })\right]. \end{gathered}$$

$$g_{aa}=-\Gamma {\vartheta }_{aa},g_{xy}=-\Gamma {\vartheta }_{xy},g_{xz}=-\Gamma \left[2{\vartheta }_{xz}+\frac{3}{2}{\mu }_x\right],g_{yz}=-\Gamma \left[2{\vartheta }_{yz}+\frac{3}{2}{\mu }_y\right].$$

$$-\sum _{⟨R,R^{\prime }⟩}\frac{J_z}{2}{\sigma }_R^z{\sigma }_{R^{\prime }}^z,$$

$$f_{xx}(R)=-2J_z\left[2{\mu }_x{\mu }_y{\mu }_z-{\mu }_x{\vartheta }_{yz}(1)-{\mu }_z{\vartheta }_{xy}(R)\right]\left[\mathcal{Z}-{\delta }_{\Vert R\Vert ,1}\right]+2J_z\sum _{\begin{array}{c}{R}^{\prime }\ne 0\\ \Vert R^{\prime }-R\Vert =1\end{array}}{\mu }_y{\vartheta }_{xz}(R^{\prime }),$$

$$f_{yy}(R)=2J_z\left[2{\mu }_x{\mu }_y{\mu }_z-{\mu }_y{\vartheta }_{xz}(1)-{\mu }_z{\vartheta }_{xy}(R)\right]\left[\mathcal{Z}-{\delta }_{\Vert R\Vert ,1}\right]-2J_z\sum _{\begin{array}{c}{R}^{\prime }\ne 0\\ \Vert R^{\prime }-R\Vert =1\end{array}}{\mu }_x{\vartheta }_{yz}(R^{\prime }),$$

$$f_{zz}(R)=0.$$

$$\begin{gathered} f_{xy}(R)=J_z\left[2{\mu }_x^2{\mu }_z-2{\mu }_y^2{\mu }_z-{\mu }_x{\vartheta }_{xz}(1)-{\mu }_z{\vartheta }_{xx}(R)+{\mu }_y{\vartheta }_{yz}(1)+{\mu }_z{\vartheta }_{yy}(R)\right]\left[\mathcal{Z}-{\delta }_{\Vert R\Vert ,1}\right] \\ +J_z\sum _{\begin{array}{c}{R}^{\prime }\ne 0\\ \Vert R^{\prime }-R\Vert =1\end{array}}\left[{\mu }_y{\vartheta }_{yz}(R^{\prime })-{\mu }_x{\vartheta }_{xz}(R^{\prime })\right], \end{gathered}$$

$$\begin{gathered} f_{xz}(R)=J_z{\mu }_y{\delta }_{\Vert R\Vert ,1}-J_z\left[2{\mu }_z^2{\mu }_y-{\mu }_z{\vartheta }_{yz}(1)-{\mu }_z{\vartheta }_{yz}(R)\right]\left[\mathcal{Z}-{\delta }_{\Vert R\Vert ,1}\right] \\ +J_z\sum _{\begin{array}{c}{R}^{\prime }\ne 0\\ \Vert R^{\prime }-R\Vert =1\end{array}}{\mu }_y{\vartheta }_{zz}(R^{\prime }), \end{gathered}$$

$$\begin{gathered} f_{yz}(R)=-J_z{\mu }_x{\delta }_{\Vert R\Vert ,1}+J_z\left[2{\mu }_z^2{\mu }_x-{\mu }_z{\vartheta }_{xz}(1)-{\mu }_z{\vartheta }_{xz}(R)\right]\left[\mathcal{Z}-{\delta }_{\Vert R\Vert ,1}\right] \\ -J_z\sum _{\begin{array}{c}{R}^{\prime }\ne 0\\ \Vert R^{\prime }-R\Vert =1\end{array}}{\mu }_x{\vartheta }_{zz}(R^{\prime }). \end{gathered}$$

Full text

Coherent control over quantum single- and few-body dynamics is continuously improving, spanning atomic, optical, and solid-state systems Haroche and Raimond (2006); Yamamoto et al. (2000); Schoelkopf and Girvin (2008). An ongoing effort is focused on assembling many individually tunable systems and studying the ensuing many-body dynamics. A significant challenge lies in realizing unitary dynamics, however, the inevitable presence of dissipative processes can be utilized in different scenarios, such as by reservoir engineering Bardyn et al. (2013). Coherent time-periodic driving is a useful tool Oka and Kitamura (2019), and rich dynamics are observed with systems of strong light-matter interactions at the interface between quantum optics and condensed matter Greentree et al. (2006); Hartmann et al. (2006); Chang et al. (2008); Koch et al. (2010); Petrescu et al. (2012); Cho et al. (2008); Umucalılar and Carusotto (2011); Hafezi et al. (2013); Rechtsman et al. (2013); Otterbach et al. (2013); Carusotto and Ciuti (2013); Le Hur et al. (2016); Noh and Angelakis (2017); Hartmann (2016); Schiró et al. (2016); Fitzpatrick et al. (2017); Fink et al. (2017). Systems with a competition between interactions, nonlinearity, coherent external driving and dissipative dynamics include arrays of coupled circuit quantum electrodynamic units Houck et al. (2012); Schmidt and Koch (2013), cold atoms Bloch et al. (2012), and ions Bohnet et al. (2016). Critical phenomena and dissipative phase transitions in these open systems often come with new properties and novel dynamic universality classes Diehl et al. (2010); Sieberer et al. (2013); Lee et al. (2013); Jin et al. (2013a); Marino and Diehl (2016); Scarlatella et al. (2019); Munoz et al. (2019).

The state of an open quantum system is defined by a density matrix $\rho$, with the dynamics often treated using a Lindblad master equation, describing a memory-less bath, and the time evolution generated by the Liouvillian superoperator acting on $\rho$ Breuer and Petruccione (2002). The theoretical tools available for open quantum many-body systems are relatively limited. For driven-dissipative lattice models the meanfield (MF) approach is often employed, with $\rho$ approximated as a product of single-site density matrices. The dynamics of local observables are described by nonlinear equations, studied, e.g., for lattice Rydberg atoms Lee et al. (2011); Qian et al. (2012); Marcuzzi et al. (2014); Carr et al. (2013); Letscher et al. (2017); Parmee and Cooper (2018), coupled quantum-electrodynamics cavities and circuits Jin et al. (2013b, 2014); Schiró et al. (2016), nonlinear photonic models Foss-Feig et al. (2017); Biondi et al. (2017a), and spin lattices Chan et al. (2015); Wilson et al. (2016). A key feature of the MF phase diagrams are multistable parameter regions where two or more steady states coexist.

However, the Lindblad equation converges in general to a unique steady state in finite systems Spohn (1977); Minganti et al. (2018), making the status of the MF approximation unclear. Indeed, significant deviations from MF have been found using approximation schemes accounting for quantum correlations Weimer (2015); Biondi et al. (2017b); Jin et al. (2016), and also using exact numerical methods (quantum trajectories Daley (2014) and Matrix Product Operators (MPO) Prosen and Žnidarič (2009)). In one-dimensional (1D) lattices with nearest-neighbour (NN) interactions, the MF bistability is found to be replaced by a crossover driven by large quantum fluctuations Weimer (2015); Mendoza-Arenas et al. (2016); Foss-Feig et al. (2017); Vicentini et al. (2018). In contrast, in certain 2D NN models, MF bistability has been found by approximate methods to be replaced by a first-order phase transition between two states, for nonlinear bosons using a truncated Wigner approximation Foss-Feig et al. (2017); Vicentini et al. (2018), and for Ising spins using a variational ansatz acounting for short-range correlations Weimer (2015), a cluster MF approach Jin et al. (2018) and two-dimensional tensor network states Kshetrimayum et al. (2017). In a parameter region around the jump, the convergence towards the steady state slows down Weimer (2015); Vicentini et al. (2018), a phenomenon related to a gap closing in the spectrum of the Liouvillian Cai and Barthel (2013); Macieszczak et al. (2016); Minganti et al. (2018).

In this Letter we study a driven-dissipative model of spins one-half with XY (flip-flop) interactions in presence of coherent drive and dissipation, using a combination of MPO simulations and an approximation scheme which accounts for quantum fluctuations beyond meanfield (MFQF). For one dimensional lattices we confirm the existence of a unique steady state in the thermodynamic limit. As our main result, we find that in dimension two and higher multistability (in particular, bistability) is again possible, with jumps between the different steady states accompanied by a critical slowing down, provided that the thermodynamic limit of an infinitely large lattice is taken first with respect to the long time limit. We argue that this order of limits is physically plausible, and we link the bistability to the fact that for finite size and time the probability distribution of relevant observables develops a strong bimodal structure. Depending on the order of limits, bimodality leads either to a first-order dissipative phase transition (as usually discussed when the long time limit is taken first), or to a bistable regime. We thus provide a theoretical scenario reconciling our results with the literature cited above, and finding similar dynamics in a model with Ising interactions (see below and App. B) indicates the generality of our results. We suggest that in an experimental platform based on trapped-ion quantum simulators, such a question can be addressed.

Model. We consider a quantum system with $N$ sites $R\in\mathbb{Z}^{D}$ on a hypercubic lattice in $D$ spatial dimensions, for which the connectivity is $\mathcal{Z}=2D$. The master equation for $\rho$ is defined using the Liouvillian $\hat{\mathcal{L}}$,

[TABLE]

The Hamiltonian describing Rabi oscillations of two-level systems with a drive detuned by $\Delta$ from the resonant transition frequency and a Rabi frequency $\Omega$, is given in a frame rotating with the drive by

[TABLE]

where the second sum extends over all pairs of NN sites, describing hopping with amplitude $J$, with spin-$\frac{1}{2}$ operators (Pauli matrices) $\sigma_{R}^{a}$, $a=\{x,y,z\}$, and $\sigma^{\pm}_{R}=({\sigma^{x}_{R}\pm i\sigma^{y}_{R}})/{2}$. For spin losses occurring independently at each site with rate $\Gamma=1$ (which fixes the frequency and time units),

[TABLE]

Aside from translation invariance, this model has no manifest microscopic symmetries. Its MF phase diagram displays bistability Mendoza-Arenas et al. (2016); Wilson et al. (2016) in a region of parameters that terminates at a second order point where an emerging $Z_{2}$ symmetry spontaneously breaks Marcuzzi et al. (2014). In 1D the steady state is unique as obtained by MPO simulations Mendoza-Arenas et al. (2016). Here we focus on higher dimensions, which we find to manifest bistability in the thermodynamic limit.

Dynamics of Observables. From the master equation one can derive a hierarchy of equations of motion for $n$-points expectation values of the form $\langle\sigma_{R_{1}}^{a}\sigma_{R_{2}}^{b}\cdots\sigma_{R_{n}}^{c}\rangle\equiv{\rm tr}\{\rho\sigma_{R_{1}}^{a}\sigma_{R_{2}}^{b}\cdots\sigma_{R_{n}}^{c}\}$, which depend on the value of correlators at the next order, $n+1$. Assuming a translationally-invariant density matrix, we define the uniform vector mean magnetization, $\mu_{a}(t)=\left\langle\frac{1}{N}\sum_{R}\sigma_{R}^{a}\right\rangle=\left\langle\sigma_{R}^{a}\right\rangle$, and its equations of motion

[TABLE]

with the connected two-point correlation functions,

[TABLE]

and setting $R^{\prime}=0$ using the translation invariance, ${\eta}_{ab}(1)$ is the correlator at a NN of the origin.

Equations (4)-(6) are exact. The limit ${\eta}\to 0$ reduces $\rho$ to a product of identical on-site states, leading to the MF equations, whose steady state and dynamics are studied in detail in Landa et al. . We present an approximate scheme going beyond MF, formally based on an expansion in $1/\mathcal{Z}$ (with a related approach in Biondi et al. (2017b)). Neglecting the connected three-point correlators $\left\langle\left(\sigma_{R}^{a}-\mu_{a}\right)\left(\sigma_{R^{\prime}}^{b}-\mu_{b}\right)\left(\sigma_{R^{\prime\prime}}^{c}-\mu_{c}\right)\right\rangle\approx 0$, allows us to derive (see App. E) coupled equations for $\eta_{ab}(R,t)$, which we solve numerically together with their feedback into Eqs. (4)-(6). Since the short-range correlators $\eta_{ab}(1)$ appearing in Eqs. (4)-(6) are dynamically coupled to all distances in the lattice, the MFQF method accounts for the spatial structure of correlation functions. The simulations have been verified to converge as a function of $N$, and hence we can approximate the system dynamics as a function of time with the limit $N\to\infty$ taken first.

Correlations wash away Bistability in 1D. We start our analysis with numerically exact MPO calculations of large lattices in 1D. The density matrix $\rho$ can be considered as a pure state in an enlarged Hilbert space with four states per site Mascarenhas et al. (2015), allowing us to solve the Lindblad evolution using a method formally similar to pure state unitary evolution encoded using well-established matrix product states (see Schollwöck (2011); Zaletel et al. (2015) and references therein). We evolve $\rho(t)$ in a 1D chain with open BC (translation invariance is not enforced), using an MPO algorithm Prosen and Žnidarič (2009); Benenti et al. (2009); Mascarenhas et al. (2015), with an implementation based on the iTensor library ite (n 21), a Trotter decomposition of order four Zaletel et al. (2015); Bidzhiev and Misguich (2017), and bond dimension $\chi=300$. With up to 200 spins we checked that observables measured in the central region of the chain had negligible finite-size effects and truncation errors at the scale of the plots, allowing us to obtain their steady-state bulk values corresponding to the thermodynamic limit.

Figure 1(a) shows the $x$ component of the steady-state magnetization, $\vec{\mu}^{S}\equiv\lim_{t\to\infty}\vec{\mu}(t)$, in 1D, for $J\mathcal{Z}=4$ as a function of $\Delta$. In MF, $\vec{\mu}^{S}$ is unique except for $1.3\lesssim\Delta\lesssim 1.9$, where there are two co-existing stable solutions in addition to an unstable solution. At the presence of quantum correlations (in MPO and MFQF), the magnetization departs significantly from the MF prediction, with a crossover between the two limiting regimes, in the range $1.5\lesssim\Delta\lesssim 5$. We define the inverse correlation lengths $\lambda_{ab}$ by fitting the six correlation functions to $\eta_{ab}(R)\sim\exp\{-\lambda_{ab}R\}$. For simplicity, we present in Fig. 1(b) only one correlation length, and Fig. 1(c) shows the corresponding total correlation measured by $\Sigma_{ab}=\sum_{R}\eta_{ab}(R)$. The spatial structure of the two-point correlation undergoes a sharp change within the crossover region, from relatively small but widely extended correlations for low $\Delta$, to much larger but very short-ranged correlations, for high $\Delta$. A separate analysis of the correlations shows that at the same time, the correlations change nature from periodic modulations (spin density-wave character), to being overdamped in space.

As Fig. 1 shows, the MFQF approximation correctly captures the uniqueness of the steady state and the disappearance of bistability in 1D. On both sides of the crossover region, the results are quantitatively accurate. In its center, the approximation reaches too large values for $\eta_{ab}(R)$ and the correlation length. More generally, we find that as $J$ is increased in 1D, the MFQF approach loses its accuracy (for parameters of strong correlations), plausibly because of the role of higher-order correlation functions that are neglected, which can lead at much larger $J$ to the breakdown of the approximation. However, the MFQF approach is easy to generalize to higher dimensions, and quantitatively accurate in regions with moderate correlations.

Bistability in Higher dimensions. Figure 2(a) shows the results of simulations with large 2D and 3D lattices, for $J\mathcal{Z}=4$. The MFQF theory, that allows simulating large lattices (with $N\to\infty$), is compared with 2D-MPO calculations, limited to a finite-size system, for which, as in 1D, $\rho$ is encoded as a product of matrices. The matrix product runs over a snake-like path visiting all the sites of a cylinder of length $L_{x}=12$ and perimeter $L_{y}=4$ (see App. C). Such an approach has been applied in ground-state calculations of 2D models Stoudenmire and White (2012), but we are not aware of previous 2D-MPO Lindblad calculations. For $\Delta\lesssim 1.5$ and $\Delta\gtrsim 2.5$ the agreement between MPO and MFQF is almost perfect, giving a nontrivial check of the ability of MFQF to capture significant correlation effects (that result in $\vec{\mu}$ strongly departing from MF). The computational cost of guaranteeing a high accuracy in 2D-MPO calculation is exponential in $L_{y}$ (see App. C), limiting the present MPO calculations to relatively small systems, which cannot show bistability (and a possible discontinuity would also be smeared out).

As our main result, using MFQF we find in 2D two stable $\vec{\mu}^{S}$ branches, that in 3D extend over larger ranges of $\Delta$, converging towards the MF bistability region and magnetization values. Figure 2(b) shows that $\Sigma_{xx}$ increases by two orders of magnitude for one of the bistable states, and Fig. 2(c) shows the asymptotic relaxation rate associated to the convergence to $\vec{\mu}^{S}$. It is obtained by fitting $\partial_{t}\vec{\mu}^{2}\sim e^{-\kappa t}$ at large times ($t\sim 100$). The fact that $\kappa\to 0$ at the branch edges in 2D indicates a critical slowing down when approaching the end of the bistability region in the phase that is about to disappear, leading to a discontinuous jump. The MFQF approach does however not always predict bistability in 2D. Replacing each hopping term in Eq. (2) by the Ising coupling $J_{z}\sigma_{R}^{z}\sigma_{R^{\prime}}^{z}$, we find a smooth crossover for moderate $J_{z}$ (as obtained using cluster meanfield Jin et al. (2018)), and a small bistability region for stronger couplings (again as in Jin et al. (2018)); See App. B.

Bistability, Liouvillian spectrum and bimodality. We henceforth return to the question raised in the introduction: how to reconcile the uniqueness of the steady state in finite systems, with bistability seen when taking first the thermodynamic limit of infinite size, and then the long-time limit?

Considering the Liouvillian $\hat{\mathcal{L}}$ of Eq. (1), the unique thermodynamic steady state corresponds to $\rho_{ss}=\lim_{N\to\infty}\lim_{t\to\infty}\rho(t)$, which is independent of the initial conditions, and is an eigenstate of $\hat{\mathcal{L}}$ at any $N$, $\hat{\mathcal{L}}\rho_{ss}=0$. Assuming bistability, we define $\rho_{1}$ and $\rho_{2}$ as the two distinct density matrices obtained from $\lim_{t\to\infty}\lim_{N\to\infty}\rho(t)$, which depend on the initial conditions. For $N$ large but finite the bistability should be replaced by long-lived metastable states, in which case $\rho_{1}$ and $\rho_{2}$ are defined at times $t\gg 1/\Gamma$, but small compared to the lifetime of these metastable states. As in the model studied here, these two states have different local properties: a local observable (i.e. sum of local terms, e.g. $M_{x}=\frac{1}{N}\sum_{R}\sigma^{x}_{R}$) has a probability distribution $P_{1}(m)$ with a single peak centred around $m_{1}$ in $\rho_{1}$, and a distribution $P_{2}(m)$ peaked around $m_{2}$($\neq m_{1}$) in $\rho_{2}$. At the same time, metastability implies some relaxation time diverging with $N$, and the spectrum of $\hat{\mathcal{L}}$ must have at least one nonzero eigenvalue $\lambda_{E}$ with a vanishingly small real part, $\lim_{N\to\infty}{\rm Re}(\lambda_{E})=0$, other eigenvalues being separated by a gap $\mathcal{O}(\Gamma)$ Macieszczak et al. (2016); Rose et al. (2016). We assume for simplicity that $\hat{\mathcal{L}}$ has a unique such small eigenvalue (therefore real), and denote by $\rho_{E}$ the associated eigenstate (or eigenmatrix). $\rho_{ss}$ and $\rho_{E}$ are the eigenstates from which all long-lived states can be constructed, since for $t$ much larger than $1/\Gamma$ we can ignore higher “excited” eigenstates. So, for $\rho_{1}$ and $\rho_{2}$ to be long-lived, they must be linear combinations of $\rho_{ss}$ and $\rho_{E}$. As physical states have a trace equal to 1, and since ${\rm Tr}\rho_{ss}=1$ and ${\rm Tr}\rho_{E}=0$, there must exist two distinct scalars $a_{1}$ and $a_{2}$ such that $\rho_{i}=\rho_{ss}+a_{i}\rho_{E}$ with $i=1,2$ Macieszczak et al. (2016); Rose et al. (2016); Minganti et al. (2018); Letscher et al. (2017). Inverting these relations we get $\rho_{ss}=(a_{2}\rho_{1}-a_{1}\rho_{2})/(a_{2}-a_{1})$. So, if $a_{1}$ and $a_{2}$ are both nonzero (which may not be always the case) $\rho_{ss}$ is a “cat state” (with correlation functions extending over the system size), being a linear combination of two uniform physical states with different local properties. In $\rho_{ss}$, the probability distribution $P_{\rm ss}=(a_{2}P_{1}-a_{1}P_{2})/(a_{2}-a_{1})$ of a local observable is bimodal, peaked around the two mean values $m_{1,2}$ realized in the states $\rho_{i=1,2}$.

Using exact diagonalization on small systems we have computed such distributions for the fully-connected (FC) version of the present XY model, which is bistable in the thermodynamic limit (where MF becomes exact), and for the 1D and 2D cases. We find (App. A), that the magnetization becomes bimodal in parts of the MF bistability region for the FC and 2D cases, whereas it stays mono-modal in 1D. The scaling with $N$ of the bimodal peaks is beyond the scope of the current work, however, the mean value of an observable computed with $\rho_{ss}$ may become discontinuous as $N\to\infty$ at some value of the parameter Foss-Feig et al. (2017); Casteels et al. (2017); Vicentini et al. (2018); Minganti et al. (2018). This could correspond, in the discussion above, to smoothly varying $\rho_{i=1,2}$ but discontinuous jumps of $a_{1}$ and $a_{2}$. Hence, a unique steady state with a discontinuous jump is a-priori compatible with bistability and hysteresis and, in the present scenario, finding one or the other in a theory calculation is a matter of order of limits.

Moreover, since the support of $P_{1}$ has essentially no overlap with that of $P_{2}$ (for a large enough system), any density matrix which is not a convex combination of $\rho_{1}$ and $\rho_{2}$ would give some (unphysical) negative probability density. This means that all physical long-lived states are convex combinations of the mono-modal states $\rho_{i}$, and the latter thus coincide with the extreme states of Macieszczak et al. (2016). The above discussion therefore connects our results both with the theories of first-order phase transitions, and the theory based on the extreme metastable states. The lifetimes of the many-body metastable states would diverge with $N$, plausibly $\propto e^{N}$, and for large enough $N$, exceed the time accessible in numerical or experimental realizations. We conjecture that an initial state with a finite correlation length will lead, in the time window $1/\Gamma\ll t\ll 1/|\lambda_{E}|$, to one of the two mono-modal states $\rho_{i}$, and not to an arbitrary combination of the two. A heuristic argument is given in App. D. A product state is a natural reproducible initial state in an experiment, allowing to explore the metastability. As a parameter is swept back and forth across the bistability region in an experimental setup, observables will show hysteretis loops – unless the sweep is unrealistically slow ($\propto e^{-N}$).

Experimental feasibility. In addition to possible realizations with circuit-QED arrays Mendoza-Arenas et al. (2016), driven-dissipative spin models can be realized in current experiments with a few tens to a few hundreds of trapped ions. Ising and XY interactions can be implemented by laser beams inducing spin-motion coupling along one or two orthogonal directions Porras and Cirac (2004); Friedenauer et al. (2008); Schneider et al. (2012), with an additional laser for the on-site Hamiltonian. As recently demonstrated experimentally, the interaction can be varied from being almost independent of distance to a dipolar power-law, and therefore short-range in 1D Islam et al. (2013); Smith et al. (2016) and 2D lattices Bohnet et al. (2016). The interaction strength in these works is of order $J/\hbar\sim 10^{4}s^{-1}$, one to two orders of magnitude larger than the qubit dephasing rates, and the rate of spin-flip processes in Eq. (3) can be potentially controlled as well.

To conclude, studying lattices of driven-dissipative interacting spins using state-of-the-art 1D MPO simulations, for the parameters presented here and in further parameter regimes Landa et al. , we have found no phase transition but a crossover between two regimes with different characteristics. On the other hand, using a new approach that accounts for the leading-order lattice correlations and their feedback onto the mean magnetization, bistability appears to be possible in driven-dissipative quantum systems already in 2D. Thus, the present exact and approximate calculations suggest that $D=2$ is a lower critical dimension for bistability in this problem. This conclusion is consistent with works done in the context of Rydberg atoms on related models Marcuzzi et al. (2014), pointing toward a model-A dynamic universality class (whose lower critical dimension is known to be two) for the second order phase transition at the ending point of the bistability regime. This implies that in one dimension fluctuations destroy the critical point and with it the entire bistability region, in line with our results.

The question of the existence of a lower critical dimension for bistability, bimodality and hysteresis and the accompanied dissipative phase transitions in this model can be directly addressed experimentally. If a definite answer is found, it would constitute the first demonstration of deciding a question currently intractable classically, by a controlled quantum simulation. It could ascertain the status of the meanfield approximation in these systems, and shed light on the differences between equilibrium and nonequilibrium phase transitions.

Appendix A Bimodality

We start studying the model of Eqs.(1)-(3) of the main text on a finite FC lattice with $N$ sites, i.e. a graph where all sites are linked (with connectivity $\mathcal{Z}=N-1$). The FC version of the model is interesting because: (i) it has a unique steady state for any finite $N$, (ii) its meanfield (MF) solution gives a set of non-linear differential equations which can support bistability, and (iii) the MF approximation becomes asymptotically exact in the limit $N\to\infty$. This is due to the fact that each site is coupled to the $x$ and $y$ components of the total magnetization of the $N-1$ other sites, and the fluctuations of this magnetization generically become small compared to its mean when $N\to\infty$. As the thermodynamic limit $N\to\infty$ is approached, the bistability is thus expected to appear in some way and the FC model can be viewed as a playground to investigate how a unique steady state at finite $N$ can be reconciled with bistability in the thermodynamic limit.

We thus compute exactly the unique stationary state density matrix $\rho_{ss}=\lim_{t\to\infty}\rho(t)$ on small systems. We focus on the probability distribution of the $x$ magnetization per site, i.e. $P_{x}=P(M_{x}=m)$, where $M_{x}=\frac{1}{N}\sum_{R}\sigma_{R}^{x}$, that we plot in Fig. 3(a) at fixed $\Omega$ and $J\mathcal{Z}$ (the rescaling by $\mathcal{Z}$ allows to compare lattices of a different connectivity), and for different values of $\Delta$. Quite interestingly, such probability distribution evolves from a single peak centered around $m>0$ to a single peak centered around $m<0$, with an intermediate $\Delta$ range where $P_{x}$ has two separated peaks. To quantify the bimodality of the distribution we define an index $b_{x}=2(P_{\max,2}-P_{\min})/(P_{\max,1}+P_{\max,2})$, where $P_{\max,2}\leq P_{\max,1}$ are the two maxima of the distribution, and $P_{\min}$ is the minimum between the two. As can be seen in Fig. 3(b), the extent of the bimodality region in $\Delta$ and the maximum of $b_{x}$ increase with $J\mathcal{Z}$ at a fixed $\Omega$, following the MF bistability region.

While in the finite size limit we consider here the stationary state is still unique and the average value of the observable will either evolve smoothly with $\Delta$ or give rise to a step, depending on the scaling with $N$ of the peaks, the emergence of a bimodal probability distribution suggests a possible scenario for bistability to survive. Indeed for a FC lattice taking the thermodynamic limit first leads to a set of non-linear differential equations which can support multiple stable solutions. This should translate, in the probability distribution of the magnetization at finite size and finite time, and depending on the initial condition, in the emergence of a second peak on time scales exponentially in the system size such that, in the thermodynamic limit the switching from one solution to the other becomes exponentially suppressed.

Repeating this analysis for a 2D square lattice with periodic boundary conditions (BC) we find again a bimodal parameter region, shifted towards higher $\Delta$, with and the maximum $b_{x}$ smaller than in the FC model [Fig. 3(c)]. With an increasing interaction strength $J$, the magnitude and extent of correlated fluctuations in the lattice grow, explaining the gradual decrease of the bimodal parameter region for a small lattice. As further consistency checks, we verified that when increasing $N$ (up to 10) at specific parameter values, the bimodality range and maximum increase, and for a 1D chain we find (not shown) that $b_{x}$ remains strictly zero.

Figure 4 shows the bimodality index $b_{x}$ calculated exactly for three increasing lattice sizes (in the 2D case shown in panel (b), different parallelograms are constructed with periodic boundary conditions to incorporate $N$ sites), for the driven-dissipative XY model discussed in the main main text. Although sensitive to the finite system size, the width of the $\Delta$ ranges of bimodality increase with $N$ and the maximal $b_{x}$ value increases accordingly, which is consistent with bistability in the $N\to\infty$ limit.

Appendix B Bistability in a driven-dissipative Ising model in 2D

We have employed the MFQF method to the study of the driven dissipative Ising model, obtained after replacing the XY (flip-flop) interaction term in the Hamiltonian,

[TABLE]

while keeping the dissipator identical. The equations of motion for $\vec{\mu}$ that are obtained from this model are identical to those in Eqs. (4)-(6) of the main text with just the replacement $J\to-J_{z}$. The correlator equations (given in the following), are different.

Setting $\Delta=0$ we have studied the resulting MFQF steady state on a 2D lattice. A necessary condition for bistability in the meanfield limit Landa et al. , is $|J_{z}\mathcal{Z}|>2$. For $J_{z}\mathcal{Z}=-4$, Fig. 5 shows a smooth crossover of the steady-state magnetization $\vec{\mu}^{S}$ as a function of $\Omega$ across the parameter region of meanfield bistability. The $z$ magnetization curve of Fig. 5(b) can be compared with the results plotted in Fig. 1(a) of Jin et al. (2018), which were obtained using a cluster MF approach for different cluster sizes. In the notation of Ref. Jin et al. (2018), $J_{z}\mathcal{Z}=-V/\gamma$ and $\Omega=(h_{x}/\gamma)/2$, making the parameters in the plots identical (with only a somewhat larger range taken here for the abscissa). The MFQF curve that we obtain contains a further noticeable feature around $\Omega\approx 1.5$, with a relatively sharper change in magnetization.

In Fig. 6, ${\mu}_{x}^{S}$ and ${\mu}_{z}^{S}$ are shown for $J_{z}\mathcal{Z}=-6$, in a larger region of $\Omega$ across the parameter region of meanfield bistability. The $z$ magnetization curve of Fig. 6(b) can be compared with the curves plotted in Fig. 1(b) of Jin et al. (2018). In the region where the cluster MF of Ref. Jin et al. (2018) predicts a shrinking of the bistability (but not its complete disappearance up to the largest available cluster size), MFQF predicts a rather sharp but smooth crossover. However in MFQF approach bistability remains in a small region ($2.06\lesssim\Omega\lesssim 2.13$), at the right edge of MF bistability region. The structure of these the two co-existing stable phases is distinctly different in terms of the correlation functions $\eta_{ab}(R)$ defined in Eq. (7) of the main text. Figure 7 presents the characteristics of $\eta_{zz}(R)$ as an example, showing the large differences in the correlation length, and – as a result – the total correlation, between the two bistable states.

Appendix C 2D-MPO

The 2D MPO calculations presented in Fig. 3 of the main text and in Fig. 8, have been carried out on cylinders of length $L_{x}$ from 8 to 16, and fixed $L_{y}=4$. At fixed precision the required bond dimension $\chi$ is expected to be constant with $L_{x}$, but exponentially large in $L_{y}$. The data is obtained with a bond dimension $\chi=400$ (pushed to $\chi=600$ at $\Delta=1.75$). As in the 1D-MPO result, the steady state is obtained by evolving $\rho$ in time from an initial state where all the spins are pointing down. In practice a total time between $15$ and $25$ was used. Since the path of interacting spins associated to the MPO artificially breaks the translation invariance of the lattice in the $y$ direction, it is important to check that the bond dimension is large enough to restore the translation symmetry in the observables. In our case the magnetization was found to be translation invariant up to relative errors of the order of $\mathcal{O}(10^{-4})$.

Appendix D Finite correlation length and mono-modal states

Consider a system where the Liouvilian spectrum has a steady state $\rho_{ss}$, associated to the eigenvalue 0, and an eigenstate $\rho_{E}$ associated to an eigenvalue $\lambda_{E}$ which goes to zero when increasing the system size. We assume that all the other eigenstates are gapped and can be ignored for times $t\gg 1/\Gamma$. As discussed in the main text, there exits two special linear combinations $\rho_{i=1,2}$ of $\rho_{ss}$ and $\rho_{E}$, which are the extreme states [56], and in which local observables have probability distributions with a single peak. On the other hand, nontrivial linear combinations of $\rho_{1}$ and $\rho_{2}$ host bimodal distributions. A natural question is then: what initial conditions lead to $\rho_{1}$ and $\rho_{2}$, and what initial conditions instead lead to some convex combination of the two? We give below a simple heuristic argument suggesting that initial conditions where connected correlation functions have a finite correlation length, should, for a large enough system, lead to one of the two extreme states and not to arbitrary convex combinations of the two. Take an initial state $\rho(t=0)$ where all connected correlation functions are short-ranged. Thanks to the assumption of a gap in the Liouvillian spectrum (above $|\lambda_{E}|$), the state $\rho(t)$ will reach the metastable manifold (spanned by $\rho_{ss}$ and $\rho_{E}$) rather quickly, after a time scale of order $\mathcal{O}(N^{0})$, *i.e. *a time that is not growing with the system size. Because lattice models with short-range interactions typically have a maximum (Lieb-Robinson) velocity, beyond which the information cannot propagate, it should not be possible for the system to develop true long-range correlations (spanning the whole system) in a finite amount of time. And if correlations remain negligible at distances comparable with the system size, then this implies that the probability distributions of local observable cannot be bi-modal. So, the Lindblad dynamics of a large system initialized from a state with short-range correlations cannot lead to a bimodal state in the time range $1/\Gamma\ll t\ll N^{1/D}$ ($N^{1/D}$ is proportional to the linear size of the system). We thus conclude that in this time window the state that is reached is close to one of the two extreme states, since any other convex combination would be bimodal.

Appendix E Meanfield with Quantum Fluctuations

To derive the equations of the MFQF approach, we define a two-point correlation function (correlator),

[TABLE]

which is a function of the difference $R-R^{\prime}$ alone, symmetric in $a,b$ (because $\sigma_{R^{\prime}}^{a}$ and $\sigma_{R}^{b}$ commute). Using Eq. (9), the connected two-point correlator is defined as in Eq. (7) of the main text, for $R\neq R^{\prime}$, by

[TABLE]

The connected three-point correlator is defined for $R\neq R^{\prime}\neq R^{\prime\prime}$ by

[TABLE]

which is again a function of the differences only.

The approximation of the following treatment is based on assuming that $\zeta$ (and higher order connected correlators) can be neglected in comparison to $\eta$. The e.o.m of $\vartheta_{ab}$, setting $R^{\prime}=0$, is

[TABLE]

where the local Hamiltonian terms are described using the matrix

[TABLE]

while $f_{ab}(\mu,\vartheta)\propto J$ comes from the kinetic terms, and $g_{ab}(\mu,\vartheta)\propto\Gamma$ comes from the Lindbladian part, and both are given below. By using Eq. (10) we get the e.o.m system for $\eta(R,t)$,

[TABLE]

which we solve numerically together with the coupled system for $\vec{\mu}(t)$.

For the Hamiltonian in Eq. (2) of the main text where the kinetic term is

[TABLE]

we find that $f_{ab}(\mu,\vartheta)$ of Eq. (12) is given by,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The components of $g(\mu,\vartheta)$ in Eq. (12) are given by

[TABLE]

For the Ising model of Eq. (8) where the kinetic term is

[TABLE]

we get instead of the above $f_{ab}(\mu,\vartheta)$, the following expressions;

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]