Cluster mean-field approach to the steady-state phase diagram of dissipative spin systems

We show that short-range correlations have a dramatic impact on the steady-state phase diagram of quantum driven-dissipative systems. This effect, never observed in equilibrium, follows from the fact that ordering in the steady state is of dynamical origin, and is established only at very long times, whereas in thermodynamic equilibrium it arises from the properties of the (free) energy. To this end, by combining the cluster methods extensively used in equilibrium phase transitions to quantum trajectories and tensor-network techniques, we extend them to nonequilibrium phase transitions in dissipative many-body systems. We analyze in detail a model of spin-1=2 on a lattice interacting through an XYZ Hamiltonian, each of them coupled to an independent environment that induces incoherent spin flips. In the steady-state phase diagram derived from our cluster approach, the location of the phase boundaries and even its topology radically change, introducing reentrance of the paramagnetic phase as compared to the single-site mean field where correlations are neglected. Furthermore, a stability analysis of the cluster mean field indicates a susceptibility towards a possible incommensurate ordering, not present if short-range correlations are ignored.

We show that short-range correlations have a dramatic impact on the steady-state phase diagram of quantum driven-dissipative systems. This effect, never observed in equilibrium, follows from the fact that ordering in the steady state is of dynamical origin, and is established only at very long time, whereas in thermodynamic equilibrium it arises from the properties of the (free-)energy. To this end, by combining the cluster methods extensively used in equilibrium phase transitions to quantum trajectories and tensor-network techniques, we extend them to non-equilibrium phase transitions in dissipative many-body systems. We analyze in detail a model of spins-1/2 on a lattice interacting through an XYZ Hamiltonian, each of them coupled to an independent environment which induces incoherent spin flips. In the steady-state phase diagram derived from our cluster approach, the location of the phase boundaries and even its topology radically change, introducing re-entrance of the paramagnetic phase as compared to the single-site mean field where correlations are neglected. Furthermore a stability analysis of the cluster mean-field indicates a susceptibility towards a possible incommensurate ordering, not present if short-range correlations are ignored.

I. INTRODUCTION
In thermodynamic equilibrium a transition to a state with a spontaneous broken symmetry can be induced by a change in the external conditions (such as temperature or pressure) or in the control parameters (such as an external applied field). The most widely studied examples are for systems at non-zero temperature, in the framework of classical phase transitions [1]. Here, equilibrium thermal fluctuations are responsible for the critical behavior associated with the discontinuous change of the thermodynamic properties of the system. Transitions may also occur at zero temperature, as a function of some coupling constant [2]; in that case, since there are no thermal fluctuations, quantum fluctuations play a prominent role. For many decades, the study of phase transitions and critical phenomena has attracted the attention of a multitude of scientists from the most diverse fields of investigations: Phase transitions are present at all energy scales, in cosmology and high-energy physics, as well as in condensed matter.
Moving away from the thermodynamic equilibrium, collective phenomena and ordering also appear in open systems, upon tuning the rate of transitions caused by the environment [3]. For example, they emerge in most diverse situations [4] ranging from the synchronous flashing of fireflies [5] to the evolution of financial markets [6]. The classical statistical mechanics of such driven systems (including traffic models, active matter, and flocking) has attracted an increasing attention over the years, see e.g.
Refs. [7,8]. Such interest is in part due to the remarkable possibility of achieving ordered states that are not possible in equilibrium systems, displaying for example long-range order in two-dimensional flocking [9], something forbidden by the Mermin-Wagner theorem [10] in equilibrium.
Thanks to the recent impressing experimental progresses, see e.g. Refs. [11][12][13], the investigation of nonequilibrium properties of driven-dissipative systems has entered the quantum world. Rydberg atoms in optical lattices [14], systems of trapped ions [14], excitonpolariton condensates [15], cold atoms in cavities [16], arrays of coupled QED cavities [17,18], are probably the most intensively investigated experimental platforms in relation to this aim. The predicted steady-state phase diagram of these driven dissipative systems becomes incredibly rich, displaying a variety of phenomena. Just as for classical statistical mechanics, phases, which are not possible in an equilibrium phase diagram, may appear [19]. The steady state itself needs not be timeindependent and the system may end up in a limit cycle [20][21][22][23][24]. Renormalisation group (RG) calculations using the Keldysh formalism have been performed [25]; in some cases the universality class of the transitions may be modified both by the presence of the external environment and by non-equilibrium effects [26,27]. A judicious engineering of the system-bath couplings can lead to non-trivial many-body states in the stationary regime [28,29]. The field of dissipative many-body open systems embraces a much wider class of problems, rang-ing from transport to relaxation dynamics to quantum information processing (just to mention few examples). A more comprehensive panorama of the recent literature can be also found in Refs. [30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49] and citations therein.
In condensed matter systems, most notably in Josephson junction arrays, the impact of an external bath on the phase diagram, and the relative critical properties was thoroughly studied over the last twenty years, see e.g. [50][51][52][53]. In all those studies the system and the bath were in an overall equilibrium situation at a given (possibly zero) temperature. In quantum driven-dissipative systems, such as that considered here, non-equilibrium conditions and the flow of energy through the system play a major role.
Our work focuses on an important aspect of the physics of many-body open systems: the determination of the steady-state phase diagram. We are going to consider systems in which the coupling to the environment leads to a Markovian dynamics. In these cases the evolution of the corresponding density matrix ρ(t) obeys the Lindblad equation The first term in the r.h.s. describes the coherent unitary time evolution (ruled by the system HamiltonianĤ). The second term, corresponding to a sum of Lindbladian superoperators L j [ρ], takes into account the coupling to the external bath(s). The steady-state phase diagram is obtained by looking at the long-time limit (t → ∞) of the solution to Eq. (1) and computing appropriate averages Ô = Tr Ô ρ t→∞ ≡ Tr Ô ρ SS of local observablesÔ, in order to determine the (possible) existence of phases with broken symmetries (space, time, spin, . . . ) [54].
Nearly all the results obtained so far on the phase diagram (with the notable exception of the works based on Keldysh RG mentioned above) rely on the (single-site) mean-field approximation, where all the correlations are ignored. Very little is known beyond that limit about the interplay between many-body correlations and dissipation, although there are some contributions in this direction [55][56][57]. While quasi-exact numerical methods exist for open one-dimensional systems, unfortunately no true phase transitions are expected to occur in that context. Beyond one dimension, such methods are much harder to apply. However it is well known that the mean-field decoupling, while important to grasp the salient features of the system, is not at all accurate in locating the phase boundaries.
An improvement in the determination of the phase diagram can be obtained by a systematic inclusion of shortrange correlations (up to a given cluster size). In equilibrium, this has been achieved within the cluster mean-field approximation [58][59][60], and using linked cluster expansions [61]. In the cluster mean-field approach, the accuracy of the diagram is obviously related to the size of the considered cluster. Even though it is still mean-field in nature, a suitable scheme that combines it with finitesize scaling may in principle allow to extract non-classical critical exponents [62]. In higher dimensions (above the lower critical one) where one expects spontaneous symmetry breaking, cluster methods lead only to quantitative corrections (a mere shift) to the mean-field predictions. These corrections become smaller on increasing the dimensionality.
For equilibrium phase transitions, the topology of the phase diagram is well captured at the mean-field level, and the short-range fluctuations considered by cluster methods only lead to shifts in the location of the transition lines/points. Normally, they do not cut an ordered phase into two separate parts, divided by a disordered region. The possibility to have a radical change of topology is however permitted out of equilibrium, where the spontaneous breaking of symmetry is of pure dynamical nature: terms which are formally irrelevant in the RG sense can nonetheless modify the flow of RG-relevant terms, so as to move a point in parameter space from one side to the other side of a phase boundary. Such a scenario is rarely, if ever, seen in equilibrium.
We demonstrate that the above picture is indeed verified in the open many-body context, and ordering with a non-trivial spatial pattern may emerge (see Fig. 1). The most natural way to show this is to include correlations through a cluster mean-field analysis which, to the best of our knowledge, has never been systematically applied in the open many-body context. Although the general strategy is the same as for equilibrium systems, there are several peculiarities emerging in this scenario, which need to be carefully addressed. The steady-state solution typically needs to be obtained dynamically via Eq. (1) (and not through a solution of a self-consistent equation [63]). To increase the cluster size, we introduce a new approach that combines the cluster mean field with quantum trajectories [64] and with matrix-product-operators [65,66].
We apply our technique to a spin-1/2 XYZ-model with relaxation (as previously studied by Lee et al. [19]), and show that the short-range correlations captured by the cluster approach can have a dramatic effect on the phase diagram. This last point is exemplified in Fig. 1 (that summarises one of our main results). A mean-field analysis predicts a transition from a paramagnet to a ferromagnet (upper panel) in the whole region of large couplings J y > J c y . The lower panel sketches the outcome of the cluster analysis. The ferromagnetic regime has shrunk to a finite region disappearing in the limit of large couplings. For an equilibrium system, such behavior would be very strange: large coupling strengths increase the tendency toward ferromagnetic order, yet here we find that the ordered state is destroyed by strong couplings. Furthermore, indications from a stability analysis hint at a different type of ordering at large values of J y .
The paper is organized as follows. In the next section we define the spin-1/2 model with nearest-neighbor XYZ interactions coupled to a local bath, which will be considered in the following. We then introduce the -0. 8 (3), for Jz = 1. The single-site mean-field as worked out in Ref. [19] would predict the emergence of different phases: paramagnetic (PM), ferromagnetic (FM), antiferromagnetic (AFM), and spin-density-wave (SDW) (inset to the top panel). Here we focus on the region highlighted by the red box, which displays a transition from PM to FM states (magnified in the top panel). A proper inclusion of shortrange correlations (through the cluster mean-field) shrinks the ferromagnetic region to a small "island", thus suppressing the order at large couplings and hinting to a possible incommensurate (inc) ordering (bottom panel).
cluster mean-field approach to driven dissipative systems and show how to combine it with quantum trajectories (Sec. III B) and with the matrix-product-operator (Sec. III C) formalism. We will see this method at work by looking at the steady-state phase diagram and comparing its rich features with those pointed out in Ref. [19] at the single-site mean-field level. Specifically, in Sec. IV we discuss how the location of the transition lines is qualitatively changed in the cluster approach. Our aim is to highlight the key role of short-range correlations in driven-dissipative systems. For this purpose, we will concentrate on a specific region of the diagram where a paramagnetic to ferromagnetic transition takes place. In one dimension (Sec. IV A) the cluster approach with appro-priate scaling restores the absence of symmetry breaking. While the one-dimensional results presented here are as could be expected, we believe they are however useful as a benchmark of the numerical methods employed in the rest of this paper. Surprises appear in the two-dimensional case (Sec. IV B), where a ferromagnetic phase is possible. Including cluster correlations gives rise to a phase diagram radically different from what derived within single-site mean-field. The extent of the ferromagnetic region becomes finite. The nature of the such transitions is discussed in Sec. IV C, where a stability analysis around the mean-field solution is performed. The finite extent of the ordered phase appears to persist in higher-dimensional systems (Sec. IV D), even though the mean field progressively becomes, as expected, more accurate. The underlying dynamical mechanism responsible for such dramatic modifications in the phase diagram will be discussed in Sec. IV E where we will provide a more physical intuition of the results obtained in this work. Finally in Sec. V we conclude with a brief summary of our results.

II. THE MODEL
We consider a spin-1/2 lattice system whose coherent internal dynamics is governed by an anisotropic XYZ-Heisenberg Hamiltonian, σ α j (α = x, y, z) denoting the Pauli matrices on the j-th site of the system. The Lindbladian for this model reads where γ is the rate of the dissipative processes that tend to flip all the spins down independently [σ ± j = 1 2 σ x j ±σ y j stand for the corresponding raising and lowering operators along the z axis]. In the rest of the paper we set = 1 and work in units of γ. The (single-site) mean-field phase diagram of the model defined in Eqs. (2) and (3) has been worked out in Ref. [19]; for orientation we summarise the main results of this analysis here.
It is important to remark that an in-plane XY anisotropy (J x = J y ) is fundamental to counteract the dissipative spin flips along the orthogonal direction [19]. In the case in which J x = J y , Eq. (2) reduces to an XXZ Heisenberg model. Since this latter conserves the global magnetisation along the z axis, the steady-state solution ρ SS of Eq. (1) would trivially coincide with the pure product state having all the spins aligned and pointing down along the z direction. This corresponds to a paramagnetic state where the dissipation is dominant, and such that σ x j SS = σ y j SS = 0 and σ z j SS = −1, where Ô SS = Tr Ô ρ SS denotes the expectation value of a given observableÔ on the steady state.
The steady-state phase diagram presented in Ref. [19] is particularly rich and includes, for strongly anisotropic spin-spin interactions, ferromagnetic, antiferromagnetic, spin-density-wave, and staggered-XY states. Hereafter we concentrate in the regime of parameters J x , J y ≥ 1 and J z = 1, where the single-site mean field predicts a single ferromagnetic (FM) to paramagnetic (PM) phase transition. Indeed by changing the various coupling constants, the PM phase may become unstable and the system can acquire a finite magnetisation along the xy plane ( σ x j SS , σ y j SS = 0), thus entering a FM phase. This fact is associated to the spontaneous breaking of the Z 2 symmetry which is present in the model, and corresponds to a π rotation along the z axis (σ x → −σ x ,σ y → −σ y ). The picture changes dramatically when local correlations are included.
As already mentioned in the introduction, in an open system the stationary state may also break timetranslational invariance (the steady state is time periodic) [20][21][22][23][24]. Our numerics suggests that a timeindependent solution exists for all parameters we study, and so we will not consider this last case and concentrate on stationary time-independent solutions. This corresponds to the stationary point of Eq. (1), ∂ t ρ SS = 0, irrespective of the initial condition. In the remainder of the paper, we will always implicitly refer to this occurrence.

III. METHODS
Solving Eq. (1) for a many-body system is a formidable task, even from a numerical point of view. The exponential increase of the Hilbert space makes a direct integration of the master equation unfeasible already for relatively small system sizes. Indeed one needs to manipulate a density matrix of dimensions 2 L × 2 L , which becomes a computationally intractable task already for quite small number of sites (L 10). In order to access systems as large as possible and to perform finite-size scaling up to reasonable sizes, we employ a combination of strategies.
In this section we discuss how to use cluster mean-field methods for driven-dissipative systems; these will be employed to determine the phase diagram of the model defined by Eqs. (2) and (3). In order to keep the notation as simple as possible, we will describe the cluster approach in the spin-1/2 language for nearest neighbor Hamiltonians. A straightforward extension of our formalism allows to consider generic short-range Hamiltonians of the formĤ = iĥ . . (with the various terms respectively including on-site, nearestneighbor, next nearest-neighbor, . . . , couplings) and a generic dissipator containing more than one Lindblad operator on each site.
FIG. 2: Sketch of the cluster mean-field approach in a dissipative system of interacting spin-1/2 particles. The figure refers to 2 × 2 cluster on a two-dimensional square lattice.

A. Cluster mean-field
Let us isolate a given subset C of contiguous lattice sites, hereafter called cluster, from the rest of the lattice forming the system (which is supposed to be at the thermodynamic limit). This is pictorially shown in Fig. 2. The decoupled cluster mean-field (CMF) Hamiltonian with respect to the cluster can be written aŝ faithfully describes the interactions inside the cluster, whileĤ effectively represents the mean-field interactions of the cluster C with its neighbors [σ j = (σ x j ,σ y j ,σ z j )]. The sum is restricted to the sites on the boundary B(C) of the cluster. The parameter B eff j = (B x j , B y j , B z j ) in Eq. (6) is related to the average magnetisation of the neighboring spins of i belonging to the cluster C adjacent to C. The effective field needs to be computed self-consistently in time.
This reduced description arises from a factorized Ansatz for the global density matrix where ρ C is the density matrix of the C-th cluster. Inserting such Ansatz into Eq. (1) and exploiting the translational invariance with respect to the cluster periodicity (ρ C = ρ C , ∀ C, C ) we get an effective master equation of the form: We recall that the standard mean-field treatment derives from assuming that the cluster is formed by a single site.
The mean-field approach represents a crude approximation for a many-body interacting system, since all the correlations are effectively neglected. The decoupling on a larger structure described above partially overcomes this problem: the idea is that interactions among the sites inside a cluster are treated exactly [see Eq. (5)], while those among neighboring clusters are treated at the mean-field level [see Eq. (6)]. As a consequence, shortrange correlations inside the cluster are safely taken into account. The full problem is eventually simplified into the evolution of the density matrix ρ C of the cluster in the presence of a time-dependent effective field B eff j (t). So far what we discussed equally applies to any cluster mean-field approximation, either classical or quantum. The only non trivial modification in the present case is that one has to study the evolution of Eq. (8) in the presence of a time-dependent field that has to be determined self-consistently. In order to improve its accuracy and to have a reliable scaling of the correlations, clusters of sufficiently large dimensions need to be considered. For small clusters a direct integration of the cluster master equation is feasible, while larger clusters can be faithfully treated by combining the above explained approach with specific techniques designed to deal with open systems. Specifically, we are going to integrate the cluster mean-field approximation together with quantum trajectories and with tensor-network approaches. The idea and procedure is straightforward, but some practical details require stating explicitly. We present such details in the next sections.

B. Quantum trajectories
There is a simple procedure that allows to avoid simulating the mixed time evolution of the full master equation (1) [which would need to store and evolve a 2 L × 2 L matrix ρ(t)]. Indeed it can be shown that one can equivalently perform a stochastic evolution protocol of a pure state vector of size 2 L , according to the quantumtrajectory (QT) approach [64] [which requires to manipulate N × 2 L elements, N being the number of trajectories (typically N 2 L is sufficient to get reliable results)]. The unitary time evolution part of Eq. (1), together with the anti-commutator term in Eq. (3), can be regarded as if the evolution were performed by means of an effective non-Hermitian HamiltonianĤ eff =Ĥ + iK, witĥ The remaining term in Eq. (3) gives rise to the so-called quantum jumps. It can be shown that, if the density matrix at some reference time t 0 is given by the pure state ρ(t 0 ) = |ψ 0 ψ 0 |, after an infinitesimal amount of time δt it will evolve into the statistical mixture of the pure states { ψ 0 , ψ j } j=1,...,L (the : Quantum trajectories combined with mean-field / cluster mean-field method. Coloured boxes along a given line stand for the time-evolved state of the k-th trajectory, which is stochastically chosen among the set of pure states {|ψ0(t) k , |ψj(t) k } according to Eq. (9). For each of those states, one finds the corresponding mean fields on each site j inside the considered cluster, σ j (t) k . The mean field B eff j (t) parametrizing the effective Hamiltonian in Eq. (13), to be used in the master equation for the cluster density matrix, is obtained by averaging over all the N trajectories.
tilde indicates states at time t 0 + δt): where dp j = γ ψ 0 |σ + jσ − j |ψ 0 and Therefore, with probability dp j a jump to the state ψ j occurs, while with probability 1 − j dp j there are no jumps and the system evolves according toĤ eff . Assuming that there exists a single steady state ρ SS for Eq. (1), one has [64]: for any observableÔ and reference time T 0 . The state |ψ(t) is stochastically chosen among those in Eq. (10), according to the statistical mixture (9), after iterating the above algorithm for (t − t 0 )/δt times, where the time interval δt has to be much smaller than the relevant dynamical time scales. It is possible to combine the QT method with the above described CMF approach at the cost of some moderate modifications. In order to do that, it is necessary to perform a simulation of a sufficiently large number N of trajectories in parallel. For each trajectory k, the meanfield expectation value σ j (t) k ≡ k ψ(t)|σ j |ψ(t) k on each site j of the considered cluster C is computed iteratively in time. The average over all the trajectories gives the correct mean field at time t, which has to be self-consistently used to describe effective interactions between adjacent clusters [see Eq. (6)]. Note that this approach corresponds to performing the stochastic unravelling of the cluster mean field theory. Such an approach is different from performing a cluster mean-field decoupling of a stochastic unravelling of the original equation (i.e., each trajectory would evolve according to its own mean field). Eventually one gets an effective non-Hermitian cluster mean-field Hamiltonian which, together with the possibility of having quantum jumps, governs the time evolution of each trajectory for the next time step, as in Eqs. (9)- (10). The idea of this combined approach is schematically depicted in Fig. 3, and turns out to be effective to deal with clusters containing L 10 sites.

C. Matrix product operators
Quantum trajectories are not the only method which can be fruitfully combined to cluster mean-field techniques. Tensor networks are also ideally suited to this aim. Below we will consider Matrix Product Operators (MPO) that work very well for one-dimensional (1D) systems. It would be highly desirable to have tensor-network approaches also in higher dimensions. We believe that in combination with cluster mean-field, this will represent a significant step forward in an accurate analysis of this class of non-equilibrium critical points.
For 1D systems, the long-time limit of Eq. (1) can be faithfully addressed using a MPO Ansatz for the density matrix [65,66]. The solution ρ SS is reached dynamically by following the time evolution according to Eq. (1), using an algorithm based on the time-evolving block decimation (TEBD) scheme [67] adapted to open systems.
The starting point is based on the fact that a generic many-body mixed state on a L-site lattice, ρ = can be written as a matrix product state in the enlarged Hilbert space of dimension d L ⊗ d L , where d is the dimension of the onsite Hilbert space. By means of repeated singular value decompositions of the tensor C i1···i L , j1···j L , it is possible to obtain FIG. 4: One-dimensional TEBD scheme for 1D systems with open boundaries, combined with the cluster mean-field method. Circles denote the sites of the lattice. The manybody state corresponding to the OBC cluster made up of L black circles inside the orange box is written in an MPO representation, and evolved in time with the TEBD scheme. The cluster is coupled to the rest of the system (gray circles) through the mean field at the edges. At regular small time intervals, the mean fields B eff 1 (t) = Tr σ L ρ(t) on the leftmost site and B eff L (t) = Tr σ 1 ρ(t) on the rightmost site, are self consistently evaluated and used to construct the Hamiltonian for the next TEBD iteration.
where the super-ket ||i 1 · · · i L , j 1 · · · j L = L a=1 |i a j a | is used in order to deal with the super-operator formalism, i.e. with linear operators acting on vector spaces of linear operators. The bond-link dimension χ of the MPO (14) can be kept under a given threshold by cutting the smallest singular values, and is proportional to the amount of quantum correlations between the system sites that can be encoded in ρ MPO . Starting from χ = 1 (separable state) and increasing χ, quantum correlations can be taken into account at increasing distance.
The TEBD scheme can be naturally embedded in the Ansatz given in Eq. (14), by performing a Trotter decomposition of the Liouvillian super-operator [67] which describes the master equation (1) [this can be easily handled for Hamiltonian and Lindbladian written as sums of local terms, as in Eqs. (2)-(3)]. In the case of translationally invariant systems, it is even possible to adopt an infinite version of the TEBD (the i-TEBD), using the same approach that has been successfully applied to pure states [68]. Indeed this can be generalized to encompass arbitrary 1D evolution operators that can be expressed as a (translationally invariant) tensor network [69]. The TEBD method has been proven to be very effective in many different open 1D quantum systems, as for example coupled cavity arrays [44,70], Bose-Hubbard chains with bond dissipation [71] and driven/dissipative spin systems [72]. Alternative approaches based on the variational search of the Liouvillian super-operator [73,74] or on the local purification of the density matrix [75] have been recently proposed.
The description of 1D dissipative many-body systems in terms of MPO and the search for the steady state by time-evolving the Liouvillian super-operator can be combined with the CMF approach in a natural way (see Fig. 4 for a sketch of the idea). We consider a linear cluster of L sites with open boundary conditions (OBC); its master-equation dynamics can be simulated by means of the TEBD scheme. The only novel ingredient is provided by the mean fields which have to be applied only at the two edge sites of the chain (the leftmost and the rightmost site). These can be easily evaluated in a self consistent way in time, by computing the average expectation values: respectively on site 1 and site L of the chain, at regular time intervals, as outlined above for the other methods. Such fields are inserted in the effective Hamiltonian (6), which is used to build up the Liouvillian operator for the time evolution up to the next iterative step.
As mentioned at the beginning of this subsection, the extension of all these ideas to two dimensional systems would be very intriguing. For example, one could think to combine the cluster mean-field approach with MPOs using a mapping of the lattice to a one-dimensional structure with long-range interactions, through an appropriate wiring-up strategy. This has been already successfully employed in the context of equilibrium systems, where impressive results on wide strips have been obtained (see e.g. Ref. [76]). In higher dimensions, these methods suffer of problems related to the computational cost of the tensor network contraction [77], that is common to all planar structures. The presence of dissipation could help in reducing the amount of correlations in the steady state, so that it might be possible that relatively good accuracies will be reached even with small bond links.

IV. RESULTS
Let us now put into practice the methods outlined above, and study the PM-FM dissipative phase transition of the interacting spin model described by Eqs. (1), (2)-(3). As detailed in Sec. II, this is associated to a Z 2 -symmetry breaking mechanism, whose location in the phase diagram we would like to accurately unveil.
The full phase diagram at the single-site MF level has been already obtained in Ref. [19]. By writing the meanfield equations of motion for the magnetisation along the different axes, it is possible to analytically evaluate the critical point separating the PM from the FM phase. For fixed values of J x , J z , it is located at where z is the coordination number of the lattice, i.e. the number of nearest neighbors of each lattice site. As in any single-site mean field, the only effect of the system dimensionality enters through the integer z. From the theory of critical phenomena, we know that the role of dimensionality is crucial, particularly in low dimensions. Below we show that, under a more careful treatment of the short-range correlations, cluster mean-field produces important qualitative and quantitative changes to the phase diagram. In the different subsections, we will address the cases of increasing dimensionality. In onedimensional systems, where we do not expect any phase transition, the cluster mean-field together with quantum trajectories and MPO allows to recover this result.

A. One dimension
The 1D case represents the most suitable situation to benchmark our methods. Here, due to the reduced dimensionality (the system has a coordination number z = 2), the MF predictions are known to fail and no symmetry-breaking mechanism should occur (as already stated in Ref. [19]). Using a combination of strategies as described in Sec. III, we numerically verify the absence of symmetry breaking, thus gaining confidence on how accurate our methods can be for driven-dissipative systems.
We are able to perform a direct integration of the master equation (1) for systems with up to L = 9 spins, by employing a standard fourth-order Runge-Kutta (RK) method, without applying consistent MF terms at the boundaries. For larger systems, with 10 ≤ L ≤ 16, we use the quantum trajectory approach (the time evolution of each trajectory is computed by means of a fourthorder RK method) obtaining reliable results already with a number of trajectories not exceeding N = 500, for all the values of the parameters we have probed. For even larger clusters (L 40) we resort to an MPO approach combined with the cluster mean field.
In order to check for the (possible) existence of an ordered FM phase, we calculate the steady-state ferromagnetic spin-structure factor S xx A non-zero value of S xx SS (0) indicates the stabilization of a FM ordering in the thermodynamic limit. We do not look directly at the order parameter σ x j SS , since we are studying finite-size systems and the Z 2 symmetry may not spontaneously break [78].
In 3), for which we provide a finite-size scaling (Fig. 6), and an analysis of the two-point correlation functions (Fig. 7). Here we have set Jx = 0.9 and Jz = 1, and work in units of γ. Note that for Jy = 0.9 and Jy = 1 the spin structure factor is rigorously zero.
J y > J c y ). In striking contrast with this, our numerics displays a decrease of the xx correlations with the system size. We also observe a non-monotonic behavior with J y , and the fact that S xx SS (0) vanishes for J y = 0.9 and J y = 1. This can be explained as follows. For J x = J y , the Hamiltonian (2) conserves the magnetization along z. Since the dissipative spin-flip processes occur along the same direction, it is clear that they cannot be counteracted by the unitary dynamics and so the steady state is a pure product state, having all the spins aligned and pointing down along the z direction, making the xx and yy correlations vanishing at any distance. On the contrary for J y = J z , the total magnetization along the xaxis is conserved by the Hamiltonian. In this case, due to the different privileged axis with respect to the dissipation process, the steady state is not a product state. The correlators are generally different from zero, however the spin-structure factor of Eq. (17) at k = 0 sums to zero. It is worth noticing that, on the contrary, S yy SS (k = 0) is not affected by the Hamiltonian symmetry and is different from zero (not shown).
Coming back to the data in Fig. 5 on the spin-structure factor S xx SS (k = 0), we can pinpoint the emergence of two peaks at J y ≈ 0.4 and 1.3. Before commenting on the behavior of the spin-structure factor in proximity of such peaks, let us analyze more in detail their dependence with L by performing a finite-size scaling of our data. This is provided in Fig. 6. Black data sets correspond to those in Fig. 5. We observe a systematically drop of the correlations with L for both peaks, that can be nicely fitted with a power law where the exponent α depends on the value of J y as indicated in the various panels. We were able to reach longer sizes by employing a MPO approach for considerably larger chain lengths (L ≤ 40). We applied a cluster mean field at the edges of the chain, in order to better mimic the thermodynamic limit. The results obtained with this method are displayed in Fig. 6 by the blue sets of data, and qualitative agree with the previous results without mean field (black data). In particular, an analogous power-law behavior (18) emerges. Notice that, in correspondence to the peak that is remnant of the ferromagnetic phase (J y = 1.3), a nonmonotonic behavior in the combined MPO-CMF approach emerges. This has to be ascribed to the meanfield corrections that become very effective for very short clusters.
Further evidence of the remnants of the Z 2 -symmetry breaking predicted at the mean-field level is provided by analyzing the two-point correlation functions σ x j σ x j+r SS as a function of the distance r. Figure 7 shows results for parameters corresponding to the two distinct phases predicted by the mean-field theory. In particular we observe that, in the cases where the symmetry is not broken in MF, correlations of the order parameter exhibit a clear exponential decay with the distance, as one can recognize in the upper panel (J y = 0.4). This is evident already at very small sizes L ∼ 12. A more intriguing situation occurs in the case where the MF would predict a symmetry-broken phase (see the lower panel for J y = 1.3). In such case an instability of the PM phase at short lengths emerges, in the sense that a bump in the correlators clearly emerges at r 10 and the exponential suppression of correlations is not immediately visible. Longer sizes are needed to observe the absence of quasi-long-range correlations.
To corroborate our analysis, we also performed simulations by directly addressing the thermodynamic limit. We employed a TEBD numerical approach based on a translationally invariant Ansatz for the MPO [69]. Here the mean field need not to be used. The results are in perfect agreement with those obtained with the cluster mean field, thus validating our approach. In all the cases that we considered, we clearly see an emergence of exponential decay at large distances, thus signalling the absence of ferromagnetic order in any parameter range. Remarkably, the data obtained with MPO simulations (both in the finite and the infinite case), converge with a relatively small bond-link dimension (χ ≤ 120).

B. Two dimensions
We now proceed with the discussion of the model in a two-dimensional square lattice (z = 4). Here there is no chance to solve Eq. (1) exactly for any thermodynamically relevant system size, therefore we resort to approximate techniques combined with a CMF approach. In this framework we are able to highlight a number of significant modifications to the steady-state phase diagram predicted by the single-site MF. Clearly such differences must arise from taking into account the effect of shortrange correlations inside the cluster. The shape of the considered clusters always respects the square-lattice geometry (i.e. they have a number of sites L = × ). With the numerical capabilities at our disposal, we are able to deal with clusters up to size = 4. The ≤ 3 data have been computed by numerically integrating the time evolution of the cluster master equation with a standard RK method. In order to address the case = 4, we employed the quantum trajectories approach explained in Sec. III B.
Our main result is reported in Fig. 8, which displays the phase diagram in a region of the parameter space where the MF analysis would predict the occurrence of a Z 2 -symmetry breaking mechanism. It is immediately visible that, under a CMF treatment of the system, the extent of the FM phase is drastically reduced. Specifically we shall contrast the single-site MF predictions (black line) with the results obtained using a 3 × 3 cluster size (blue circles). On the one side, the single-site MF analysis predicts a symmetry-broken phase in a large and extended portion of the phase space [fixing J z = 1, for −1 J x 1 the ferromagnet extends for any J y 1 according to Eq. (16), and disappears only in the asymptotic limit J y → ∞]. On the other side, the latter analysis indicates a tendency to confine the FM phase into a finite-size region in the parameter phase, which is surrounded by the PM phase, thus modifying the topology of the diagram.
Our CMF numerics shows that the disappearance of the ordered phase at large J y is accompanied by the progressive shrinkage of the Bloch vector for the single-site density matrix, with increasing the coupling strength. This effect can be already seen from the Bloch equations of the single-site MF [19], which predict a saturation of the spins in the limit of infinite coupling-see Eqs. (20), (21) and analogous for [M y SS ] MF . It is however important to remark that, even though the left and upper boundaries of the FM phase shrink with the cluster size while the right and bottom ones are almost unaffected, our results support the existence of a finite region for the symmetry-broken phase even in the thermodynamic limit → ∞, as we will detail below. Since the calculations with large clusters are very demanding, we considered few (representative) couplings. Our analysis performed with clusters of size up to 4 × 4 indicates that the ferromagnet will survive in the limit → ∞, for fixed J x = 0.9 and for 1.04 J y 1.4 (see Fig. 10). We do expect that for other values of J x the behavior will be similar.
Before commenting on the scaling with the cluster size, let us point out the fact that the CMF data for = 2 (red squares in Fig. 8) evidence an intermediate situation. Indeed taking into account only nearest-neighbor interactions, the extent of the FM phase is slightly reduced as compared to the single-site MF, yet it is not sufficient to confine the symmetry-broken phase into a finite-size region surrounded by the PM phase. Nonetheless, after a more careful analysis of the magnitude of the order parameter, we are able to detect a clear tendency toward a topological modification of the diagram. Specifically we fixed several values of the coupling J x , while varying J y , and investigated the FM-PM phase transition by looking at the steady-state on-site magnetisation along the x-axis: so to explore the phase diagram of Fig. 8 along certain vertical cuts. Notice that we do not need to calculate the correlators S xx SS (0) of Eq. (17) as we did in the 1D geometry, since the self-adaptive mean-field method automatically breaks the symmetry in the FM phase.
The different panels of Fig. 9 refer to four values of J x , as indicated by the first four green arrows on the left in Fig. 8, and display M x SS as a function of J y , for different cluster sizes . The 1 × 1 MF data (black lines) can be found by working out the steady-state limit of the MF Bloch equations for the magnetization [19], giving the following result: with These curves exhibit a finite magnetization for all J y 1, with a maximum at a given value of J y (dependent of J x ) and they eventually go to zero in the limit J y → +∞. This vanishing order at strong coupling is similar to the absence of ordering on resonance in the Dicke model [12], and the suppression of ordering in the degenerate limit of the Rabi model [24]. The non-monotonicity of M x SS as a function of J y emerges also in the CMF analysis: the 2×2 data signal a strong suppression of the order parameter for J y ∼ 2, which however remains finite. Going further with a 3×3 cluster, we see the sharp disappearance of the FM in an intermediate extended region where M x SS = 0 (for 1.5 J y 3, depending on the value of J x , the system is not ferromagnetically ordered along x or y). The revival of the FM phase at large values of J y (J y 3) is outside the parameter range of Fig. 8. We will analyze this feature later in Sec. IV C.
Let us now have a closer look at the vertical cut of Fig. 8 for J x = 0.9; the magnetization is shown in Fig. 10 for clusters up to = 4. Both for the 3 × 3 and the 4 × 4 CMF analysis, we do not see any reappearance of the FM ordering at large J y (we numerically checked this statement up to J y = 10). The symmetry-broken phase is confined to a finite-size region which shrinks with increasing . While the left boundary is basically unaffected by the role of correlations (J c (left) y ≈ 1.04 ± 0.01), the right boundary is strongly sensitive to . Our simulations indicate a transition point J c(right) y ≈ 2.04 ± 0.005, 1.67 ± 0.01, 1.57 ± 0.03, for clusters respectively with = 2, 3, 4. A scaling with of these data for the right boundary indicates a behavior that is compatible with J c (right) y ≈ 1.40 + 2.54 −2 , thus which supports the existence of the FM phase in the limit of large cluster size → ∞, for 1.04 J y 1.40. In the data for = 2, a discontinuity of M x SS seem to appear immediately before the right transition point (at J y ≈ 2), which requires a further analysis (a similar behavior is observed in the lower right panel of Fig. 9, for J y = 0.5). We will return to this point in Sec. IV C.
We also checked that, close to the transition, our numerics predicts a growth of the order parameter that is well approximated by . We repeated a similar analysis for other vertical (fixed J y ) and horizontal (fixed J x ) cuts, and obtained qualitatively analogous results. This evidences the fact that the CMF remains a mean-field analysis, and leads to the same critical exponents as those of its singlesite version. In order to get the correct exponents, one would need a more careful finite-size analysis [62] which requires slightly larger values of and is unfortunately out of reach for the present computational capabilities.
The stability of the symmetry-broken phase for J c (left) y < J y < J c (right) y up to 4 × 4 clusters is corroborated by the behavior of the correlation functions σ x j σ x j+r SS and σ y j σ y j+r SS with the distance r, as reported in Fig. 11 for three different values of J y . As discussed in Section IV A for the 1D case, in the parameter region where we predict a PM, the correlators decay exponentially with r (black data set at J y = 1). On the opposite side, the point at J y = 1.2 (red data) displays a marked distance-independence of correlations with the distance, thus signalling the presence of a FM phase (notice that this point lies well inside the closed region in Fig. 8). The case J y = 1.7 (blue data) shows a subtler behavior and corresponds to a point for which the single-site and the 2 × 2 mean-field analysis would predict a symmetry-broken phase, contrary to our ≥ 3 CMF calculations which display no evidence of this type. The reminiscence of a kind of quasi-ordering at short distances is indeed forecast by a slow decay of correlations. While we are not able to see a clear exponential decay with r, due to our limited numerical capabilities, we expect that this would be visible for clusters appreciably longer than = 4. Nonetheless we stress that xx corre-lations here are one order of magnitude smaller than in the FM point.
A sketch of the phase diagram summarising all our results is provided in Fig. 1.

C. Two dimensions -stability analysis
As anticipated in the previous subsection, the 2 × 2 analysis reveals a discontinuity of the order parameter inside the first FM phase, very close to the transition point J c (right) y to the disordered phase. Such a jump, between two symmetry-broken states, is known as a metamagnetic transition. The jump is visible for certain values of J x , and seems to vanish quickly with increasing the cluster size (already for = 3 it is barely recognizable from our data). On the one hand, the latter observation suggests that this jump could be an artifact of the CMF analysis. On the other hand, a deeper investigation is required to understand its origins.
To highlight the existence of this feature, in Fig. 12 we show a magnification of the relevant parameter region of Fig. 10. We only consider the 2 × 2 case, since this is the situation where it is mostly relevant. We observe the presence of a first-order phase transition within the first ordered phase, where the order parameter exhibits a discontinuity. This is corroborated by a bistability effect: specifically, we calculated the magnetization M x SS starting from different initial states and we observed a slight difference in proximity to the jump, as is visible from the figure [79].
At this point we perform a linear stability analysis, in order to check whether and how the system becomes unstable in correspondence of the jump. We start from the CMF factorization Ansatz given in Eq. (7), where each cluster density matrix ρ C obeys the mean-field master equation (8). The stability analysis is performed directly on the factorized density matrix, as detailed in Ref. [33]. Let us first rewrite the equation of motion for a single cluster, say the n-th one, in the superoperator formalism as: where we omitted the index C . Here ||ρ n denotes a super ket, i.e.. a vectorised form of the density matrix, and M i denote superoperators. In this equation M 0 represents all the on-cluster terms, while M j is the on-cluster part of an off-cluster term, and E j the corresponding offcluster expectation. For example, in the term J x σ x N σ x 1 , we have M j = −iJ x ||σ x 1 , and E j = E[σ x N ] is the superoperator form of the expectation. Moreover e j is the direction to the neighboring cluster involved.
When performing linear stability analysis, we expand the fluctuations in terms of plane waves ||ρ n = ||ρ 0 + k e ik·rn ||δρ k (24) so that the resulting equation of motion for ||δρ k is The last term is a sum of rank-one matrices (since M j ||ρ 0 is a vector, like E j ). Thus we obtain where in the second part, we have grouped terms with the same vector e together, as these all get the same k dependent factor. We then numerically compute the eigenvalues of the effective superoperator in Eq. (26) for each value of k = (k x , k y ). The most unstable eigenvalue is the one with the largest positive real part. Since for a × cluster the vectors e j must be times the elementary lattice vectors, the range of lattice momenta coming from the × cluster stability analysis is restricted to the first Brillouin zone of the superlattice, |k j | < π/ .
In Fig. 13 we plot the real part of the most unstable eigenvalue as a function of the momentum k x and the coupling J y (for fixed J x = 0.9). We notice that the jump inside the FM phase occurs when there is an instability at finite k, around |k| = π/4. This suggests that the finite cluster size is responsible for the particular metamagnetic transition seen, and explains why the extent of the ordered phase reduces as larger clusters (capable of describing such short-range fluctuations) are used. The transition to the normal state also occurs from a finite momentum instability, at small |k|. We also see an incommensurate finite-momentum instability at the rebirth of the FM phase, for large J y , thus signaling that probably the reappearance of the ordered phase is an artifact of the translationally invariant CMF Ansatz. Finally we checked that the dispersion is almost isotropic in k x , k y .

D. Three-and four-dimensional systems
For completeness we consider also the case of higher dimensions. Although not relevant for direct experimental realisations it helps in completing the picture achieved so far; it may also be possible to study (finite sized) four dimensional systems by using synthetic dimensions as proposed recently by Ozawa et al. [80]. Mean-field results are expected to improve their validity with increasing the coordination number z in the lattice. It is therefore tempting to investigate systems in higher dimensionality by means of CMF techniques. Obviously, on increasing d, our ability in considering larger clusters goes drastically down. We checked the dependence on d of the PM-FM transition by means of a mean-field analysis with clusters of size L = 2 d . In these cases we looked again at the average on-site magnetisation along the x axis. The results are displayed in Fig. 14.
Naturally, the extent of the symmetry-broken phase region is increasing with the dimensionality, as it is apparent from Fig. 14 (despite the value of the order parameter does not necessarily become larger). This supports the common wisdom of the validity of single-site mean field in high dimensions. What is surprising from Fig. 14 is that even the four-dimensional system shows a critical value of J y beyond which the phase is paramagnetic. This is in sharp contrast with the mean-field result that does not capture this second critical point. Our limited analysis up to four-dimensions and for very small clusters does not allow to draw conclusions in determining how/if the second critical point moves in higher dimensions. It is however an interesting point to be understood. E. Short-range correlations and Lindblad dynamics: Origin of the re-entrant paramagnetic phase As discussed in detail in Sec IV B, our calculations show that, on improving the Ansatz for the steady-state density matrix by including short-range correlations, the critical points may shift from J y = +∞ to a finite value of J y . This situation may appear as counterintuitive. It is indeed unlikely to occur at equilibrium, where the inclusion of short-range fluctuations may only lead to a shift of the boundary position of the order of the energy fluctuations inside the cluster [O(zJ y ) in this case]. In this section we want to explain the mechanism responsible for this behavior. This will also help us to elucidate the nature of the PM phase observed at large J y within the CMF approach and, consequently, the re-entrance to a disordered phase. To this aim it is sufficient to compare the single-site with the 2 × 1 (two-site) cluster cases. We consider this minimal cluster dimension for simplicity, since taking larger clusters would not add new ingredients to the understanding of the mechanism.
First, it is important to stress that, already at the single-site MF, a steady state with vanishing spontaneous magnetization in all the directions is predicted in the limit J y → +∞. As shown in the top panel of Fig. 15 (for fixed J x = 0.9 and J z = 1), two phases emerge: a PM for J y < J c y and a FM for J y > J c y , with magnetization along y (or equivalently along x) initially increasing, but then decreasing asymptotically toward zero as J y is increased. This phenomenon is related to the progressive deterioration of the purity of the steady-state density matrix, P = Tr[ρ 2 SS ], for J y > J c y . This comes as a consequence of the out-of-equilibrium nature of the steady-state resulting from the interplay of driving and dissipation and cannot occur at equilibrium, where an increasing coupling typically stabilizes the ordering. A similar kind of behavior can be seen in a driven two-level system [81,82] where increasing driving enhances the population, but suppresses the purity of the system, leading to a suppression of the homodyne ampltitude | σ − | and of the purity when driven on resonance. We see that P = 1 in the PM phase and then it decreases toward its minimal value (P = 1/2 in the case of a single-site cluster) as J y is increased beyond the critical value J c y . This suggests the fact that the disordered phase detected for J y < J c y is different in nature compared to the one reached in the large-J y limit for the cluster mean-field simulation. The former is due to the stabilization of a fully polarized along the z direction, which coincides with the single-site MF solution for any J y < J c y (while it is the exact solution to the problem only for J x = J y ). The latter PM phase is a consequence of the fact that the steady-state for J y → +∞ is fully mixed.
What is the effect of including short range correlations? In order to understand this point, let us consider more in detail the smallest cluster where this feature can be observed, namely a 2 × 1 plaquette. As shown in the , that has to be contrasted with a nearly fully mixed state in the right PM region, J y > J c(right) y (the minimal value for a two-site cluster is P = 1/4). Thus, the exact inclusion of the nearest neighbor correlations allows for the reentrance of a PM phase for J y > J c(right) y . At singlesite MF level, such PM phase appears only in the limiting case J y → +∞ and then is never detectable for any finite value of the couplings. We remark that the decreasing of the purity, and the consequent reduction of the magnetizations, as J y is increased is a common nonequilibrium feature of the two Ansätze (1 × 1 and 2 × 1). While in the 1 × 1 case the purity reduction is not enough to kill the FM order, in the 2×1 plaquette this reduction of purity is more prominent and the latter phenomenon (suppression of magnetization) occurs.
The equations of motion in the Heisenberg picture for magnetization σ β j (β = x, y, z) are where αβγ is the Levi-Civita symbol and δ αβ is the Kro-necker delta. The steady-state density matrix in the 2×1 plaquette for J y > J c(right) y can be analytically computed, and is almost fully mixed. Therefore it can can be written as ρ [2×1] SS ≈ ρ [1] ⊗ ρ [2] . The two-point spin correlator appearing in Eq. (27) can be thus decomposed as σ γ j σ α j+1 = σ γ j σ α j+1 + Σ γ,α j,j+1 , where | Σ γ,α j,j+1 | 1. Inserting this expression into Eq. (27) and exploiting translational invariance we get where L β [1×1] are the terms one would get from the singlesite MF. Equation (29) shows that spin-spin correlations correct the single-site MF equations of motion only through the small term Σ γ,α j,j+1 . On the other hand we know that, for J y > J c(right) y , the steady-state properties can change dramatically when considering a single site or a plaquette as a cluster: in the former case one gets a ferromagnet, while in the latter case one gets a paramagnet. Spin-spin correlations, even if very weak, cannot be neglected and drastically modify the structure of the density matrix at long times. These conspire with the dynamically induced reduction of purity at large J y , already visible for the single siste mean-field, to suppress the ordering altogether. This is the key to understand the dramatic changes in the phase boundaries we presented in the previous sections.
We believe that the mechanism is generic and should be relevant for other driven-dissipative models as well.

V. CONCLUSIONS
In this work we introduced a cluster mean-field approach combined with quantum trajectories and tensornetwork techniques to the study of the steady-state phase diagram in driven-dissipative systems. This approach allowed us to analyze the effect of short-range correlations. The result is somewhat unexpected. The whole structure of the phase diagram is radically modified in clear opposition to what typically happens in equilibrium phase transitions. In particular, we observed that the location of critical points may shift from infinite to finite values of the system parameters. The reason underlying this behavior is related to the fact that, differently from equilibrium, spontaneous symmetry breaking is of pure dynamical nature and is not determined through a free-energy analysis. It is already known that in dissipative systems, energy-minimizing ferromagnetic phases may be destabilized, and replaced by incommensurate or antiferromagnetic order. Such behaviour has been extensively studied in classical pattern forming systems [4], including examples such as active matter and flocking [7][8][9]. As such, short range correlations can be expected to play a much greater role in dissipative than in equilibrium systems. Accordingly, the topology of the phase diagram can significantly change. This appears clearly in Fig. 1, where the results from the single-site and the cluster mean-field analysis are compared. Furthermore, the cluster method hints at ordering with a non-trivial spatial pattern, a possibility which is not detected within the single-site mean-field Ansatz.
The results that we highlighted here are amenable to an experimental verification. As discussed in Ref. [19], the model considered in this paper can be implemented using trapped ions. Moreover, by changing external controls it is possible to explore the phase diagram, thus allowing to check the results of the present work. Besides the examined system, we think that cluster approaches may be powerful in the general context of driven-dissipative systems, ranging from Rydberg atoms in optical lattices to cavity or opto-mechanical arrays. Our findings point out the importance of the interplay between short-range fluctuations and dissipation in the physics emerging in such devices.
All the present analysis has been performed by considering a static mean field. It would be of great interest to extend these calculations so as to include also self-energy corrections as in the dynamical mean field, already extended to non-equilibrium for the single-site case [83].
Finally we believe that a very interesting development, left for the future, is the determination of the non-Landau critical exponents. When successful, this will be an important step to establish the power of cluster techniques also in many-body open systems. On this perspective, the combination of our approach with the corner space renormalization method developed in Ref. [55] looks promising and some encouraging results have been already obtained [84].