Cluster Gutzwiller Monte Carlo approach for a critical dissipative spin model

TL;DR

This paper investigates the phase diagram of a dissipative XYZ spin model using cluster Gutzwiller Monte Carlo methods, revealing quantum and classical correlations' effects on phase transitions and reciprocity breaking.

Contribution

It introduces a cluster Gutzwiller Monte Carlo approach to analyze short-range quantum correlations in dissipative spin systems, highlighting their impact on phase behavior.

Findings

Emergence of ferromagnetic and paramagnetic phases

Quantum correlations significantly alter the phase diagram

Reciprocity is broken in the susceptibility tensor

Abstract

We study the influence of short-range quantum correlations and classical spatial correlations on the phase diagram of the dissipative XYZ model by using a Gutzwiller Monte carlo method and a cluster Gutzwiller ansatz for the wave function. Considering lattices of finite size we show the emergence of a ferromagnetic phase, two paramagnetic phases and the possible existence of a phase transition which is entirely quantum in nature. The inclusion of short-range quantum correlations has a drastic effect on the phase diagram but our results show the inclusion of long-range quantum correlations or the use of more sophisticated methods are needed to quantitatively match the exact results. A study of the susceptibility tensor shows that reciprocity is broken, a feature not observed in closed quantum systems.

Equations

$$\hat{H}=\sum _{⟨i,j⟩}\left(J_x{\hat{\sigma }}_i^{(x)}{\hat{\sigma }}_j^{(x)}+J_y{\hat{\sigma }}_i^{(y)}{\hat{\sigma }}_j^{(y)}+J_z{\hat{\sigma }}_i^{(z)}{\hat{\sigma }}_j^{(z)}\right),$$

$${\partial }_t\hat{\rho }=-i\left[\hat{H},\hat{\rho }\right]+\frac{\gamma }{2}\sum _j\left(2{\hat{\sigma }}_j^{(-)}\hat{\rho }{\hat{\sigma }}_j^{(+)}-\left\{{\hat{\sigma }}_j^{(+)}{\hat{\sigma }}_j^{(-)},\hat{\rho }\right\}\right).$$

$$\psi (t)=\frac{\exp \left(-iHt\right)\tilde{\psi }}{||\exp \left(-iHt\right)\tilde{\psi }||},$$

$$H=\hat{H}-i\frac{\gamma }{2}\sum _i{\hat{\sigma }}_i^{(+)}{\hat{\sigma }}_i^{(-)},$$

$$\psi \to \frac{{\hat{\sigma }}_i^{(-)}\psi }{||{\hat{\sigma }}_i^{(-)}\psi ||},$$

$${\Psi }_{GW}\left(\left\{\mathcal{C}\right\}\right)=\prod _{\mathcal{C}}{\psi }_{\mathcal{C}},$$

$$S_{SS}^{xx}(\text{k})=\frac{1}{N(N-1)}\sum _{j\ne l}{e}^{i\text{k}.(\text{j}-\text{l})}⟨{\hat{\sigma }}_j^{(x)}{\hat{\sigma }}_l^{(x)}⟩.$$

$${\hat{H}}_B=h\sum _i\left(\cos (\theta ){\hat{\sigma }}_i^{(x)}+\sin (\theta ){\hat{\sigma }}_i^{(y)}\right).$$

$$\{\begin{array}{rl}& -\frac{\gamma }{2}{M}_x+2d(J_y-J_z)M_y{M}_z+2h_y{M}_z=0,\\ & -\frac{\gamma }{2}{M}_y+2d(J_z-J_x)M_x{M}_z-2h_x{M}_z=0,\\ & -\gamma \left(M_z+1\right)+2d(J_x-J_y)M_x{M}_y\\ & +2(h_x{M}_y-h_y{M}_x)=0,\end{array}$$

$$M_x=\frac{1}{\gamma }\frac{4h_y{M}_z}{1-\frac{16d^2}{{\gamma }^2}\left(J_y-J_z\right)\left(J_z-J_x\right)M_z^2},$$

$$M_y=\frac{4d}{\gamma }(J_z-J_x)M_x{M}_z,$$

$$M_z=-\frac{1}{8}\frac{{\gamma }^2}{h_y^2}.$$

$$\{\begin{array}{rl}& M_x=-\frac{h_y}{\frac{1}{4}+2\left(h_x^2+h_y^2\right)},\\ & M_y=\frac{h_x}{\frac{1}{4}+2\left(h_x^2+h_y^2\right)},\\ & M_z=-\frac{1}{1+8\left(h_x^2+h_y^2\right)}.\end{array}$$

$$\vec{M}(h,\theta )=\left(\begin{array}{cc}{\chi }_{xx}& {\chi }_{xy}\\ {\chi }_{yx}& {\chi }_{yy}\end{array}\right).\left(\begin{array}{c}h\cos (\theta )\\ h\sin (\theta )\end{array}\right).$$

$${\chi }_{\alpha \beta }={\frac{\partial M_{\alpha }}{\partial h_{\beta }}|}_{h=0}.$$

$${\chi }_{av}=\frac{1}{2\pi }{\int }_0^{2\pi }d\theta {\frac{\partial |\vec{M}(h,\theta )|}{\partial h}|}_{h=0}.$$

Full text

I Introduction

In recent years the field of dissipative phase transitions has been under intense study due to technological advances which have given access to a wide range of systems where it is possible to study the quantum dissipative behaviour of open quantum systems. The experimental platforms include trapped ions, cold atoms, semiconductor microcavities, cavity and circuit QED Bloch ; Carusotto ; Hartmann ; Noh . Driven-dissipative many-body systems find the origin of their interesting dynamics in the interplay of coherent driving and dissipation. Typically, in the long time limit, the system is driven into a steady state that is determined by the competition between both Hamiltonian and dissipative dynamics, resulting in a wide range of possible phases. Dissipative phase transitions have been observed in a wide range of experiments Rodriguez ; Fitzpatrick ; Baumann1 ; Baumann2 ; Brennecke ; Fink1 ; Fink2 . Consequently, the interest in the study of dissipative phase transitions has spiked Casteels0 ; Hwang ; Tomadin1 ; Leboite ; Diehl1 ; Tomadin2 ; Biella ; Rota ; Overbeck2 ; Lee .

Analytical solutions for open quantum systems are scarce and if one wants to calculate the dynamics of the system, one has to rely on numerical tools. As usual, a numerically exact description is unfeasible due to the exponential scaling of the Hilbert space dimension with system size and approximations are needed. Different approaches have been applied to a range of open quantum systems. Among them are the mean-field approximations Nissen ; Jin1 ; Tomadin1 ; Owen ; Leboite ; Lee ; Jin ; Diehl1 ; Tomadin2 ; Biella , the approaches based on matrix product operators and matrix product states Schroder ; Manzoni ; Orus ; Verstraete1 ; Schollwoeck ; Mascarenhas1 ; Cui ; Verstraete2 ; Zwolak - which have proven to be of great success in one dimension -, the corner space renormalization method Finazzi ; Rota , variational methods for the master equation Weimer ; Overbeck1 ; McCutcheon ; Pollock ; Suri ; Overbeck2 and variational methods at the level of the wave function Casteels1 ; Casteels2 ; Pichler ; Verstraelen ; Mascarenhas2 ; Nagy .

In this work we will study the dissipative XYZ Heisenberg model. Previous works have investigated the steady-state phase diagram by using a mean-field approach Lee , a cluster mean-field approach (CMF) Jin and the corner space renormalization method Rota . The CMF has also been used together with a self-consistent Mori projector approach to investigate the (non)existence of limit cycles Owen . A tensor network algorithm Kshetrimayum and a driven-dissipative quantum Monte Carlo Nagy have also been applied to this model. In Ref. Casteels2 this model and its phase transitions are studied using quantum trajectories and a single site Gutzwiller ansatz by using the Gutzwiller Monte Carlo approach (GMC). In this work we will extend this approach by applying a cluster Gutzwiller ansatz for the wave function. We will refer to this method as the cluster Gutzwiller Monte Carlo approach (CGMC). By including quantum correlations we show the possible existence of a dissipative phase transition entirely quantum in nature, and not predicted by the classical mean-field method, through studying the steady-state spin structure factor and the susceptibility.

We will introduce the model we will be working with in section II. The method of quantum trajectories, the used unravelling and the cluster Gutzwiller ansatz are explained in section III. In section IV we discuss the steady-state spin structure factor. The influence of an applied magnetic field will be studied in section V. In section VI the susceptibility and the possible existence of a new phase transition is looked into. Finally, conclusions are formulated in section VII.

II The model system

We study the anisotropic XYZ Heisenberg Hamiltonian (with $\hbar=1$) on a spin lattice

[TABLE]

where $J_{x}$, $J_{y}$ and $J_{z}$ are the coupling strengths in the $x$, $y$ -and $z$-direction, $\hat{\sigma}_{i}^{\alpha}$ the Pauli matrices ($\alpha=x,y,z$) and the sum goes over the nearest neighbours in the lattice. Periodic boundary conditions are applied at the edge of the lattice. A proposal for an experimental setup based on the optical pumping of two-level atoms has been made in Ref. Lee resulting in the above Hamiltonian with effective spins. This Hamiltonian governs the unitary part of the time evolution of the system. The total time evolution of the system is governed by a Lindblad Equation with dissipation along the z-axis

[TABLE]

With $\gamma$ the decay rate of the spins and $\hat{\sigma}_{i}^{(+)}$ ($\hat{\sigma}_{i}^{(-)}$) the raising (lowering) operators along the z-axis. This driven-dissipative Heisenberg model has recently been subject of several studies attempting to describe the dynamics. In this work we will compare our results with the CMF Jin . We will consider the parameters $J_{x}=0.9\gamma$ and $J_{z}=\gamma$, unless stated differently, and vary $J_{y}$. The mean-field approach Lee predicts a transition from a paramagnetic phase to a ferromagnetic phase for this parameter set at $J_{y}\approx 1.04\gamma$. By including quantum correlations in the CMF Jin or classical spatial correlations in the GMC Casteels2 the existence of another transition from the ferromagnetic phase to the paramagnetic phase is observed.

III Quantum trajectories and the wave function ansatz

The theory of quantum stochastic processes was first introduced by Davies Davies and has been further developed Barchielli ; Dalibard ; Dum ; Molmer ; Carmichael1 ; Carmichael2 into the quantum trajectory formalism, also known as the Monte Carlo wave function method. It offers an alternative to calculating the dynamics of the system with a Lindblad equation for the density matrix of the open system. Instead of solving for the density matrix, the dynamics of the system are calculated on the level of the wave function. Through a stochastic process, many realizations of the wave function, so called quantum trajectories, are simulated and by averaging over these trajectories one recovers the dynamics of the open quantum system. This method is often referred to as the unravelling of the master equation. The stochastic process finds its origin in the continuous measurement of the environment. This continuous measurement results in random changes of the wave function of the system. The unravelling of the master equation is not unique, as there are several possible measurements of the environment. In this work we will unravel the master equation through the process of photon counting Breuer . Usually the environment is under constant observation for emitted photons in this measurement scheme. As we are working with a spin system we will not be counting the excitations of the environment caused by emitted photons, we do however monitor the environment for excitations caused by spin-flips in the system. In between the detection of these excitations the wave function will evolve according to

[TABLE]

with $\tilde{\psi}$ an initial (normalized) wave function. Note that the time evolution generated by

[TABLE]

does not preserve the norm. When an excitation is detected a quantum jump is made in the evolution of the wave function

[TABLE]

after which the wave function continues evolving according to (3).

Numerically solving the master equation for the solution of the density matrix is computationally very demanding. Exact solutions are unfeasible already for small systems. One of the reasons for this inconvenient characteristic of open many-body systems is the exponentially large Hilbert space. This is were one makes some efficiency gain by using quantum trajectories: in the master equation approach one has to work with the quadratically large Hilbert space $\mathcal{H}^{2}$ of the density matrix, whereas in the trajectory approach one works with the Hilbert space $\mathcal{H}$ of the wave function. This advantage comes at the cost of needing to average over multiple realizations to obtain the dynamics of the open system. The number of realizations needed, however, is usually much smaller than the dimension of the Hilbert space.

The quantum trajectory method however does not solve the exponential scaling of computational resources with system size. In order to reduce the dimensionality of the problem, an ansatz for the wave function will be considered. Previous work has investigated the Gutzwiller ansatz Gutzwiller applied to the dissipative XYZ Heisenberg model Casteels2 . We extend this approach by including short-range quantum correlations through the use of the cluster Gutzwiller ansatz. This ansatz is realized by considering a sublattice of clusters, as shown in figure 1, where each cluster contains a set of (nearest neighbour) lattice points

[TABLE]

where the product runs over the different clusters $\mathcal{C}$ with cluster wave function $\psi_{\mathcal{C}}$. This approach allows for the inclusion of both classical and quantum correlations inside the clusters and only classical correlations between the clusters. To keep the dimension of the Hilbert space limited only small clusters will be studied and the importance of short-range quantum correlations will be determined.

IV Steady-state spin structure factor

In order to investigate the dissipative phase transition between a paramagnetic and ferromagnetic state, we will consider the steady-state spin structure factor $S_{SS}^{xx}(\textbf{k}=0)$, where:

[TABLE]

We use this correlation function rather than the spontaneous magnetization itself, because in a finite system, the $\mathbb{Z}_{2}$-symmetry does not spontaneously break. Alternatively, a (small) magnetic field could be applied to break the symmetry, as will be discussed in section V.

A non-zero value of the steady-state spin structure factor indicates a ferromagnetic phase. A zero value has a wider range of possibilities such as a paramagnetic phase, an anti-ferromagnetic phase and spin density waves. To distinguish between these phases different values of the wave vector have to be studied.

We simulate the dynamics of a trajectory over a minimum total time of $10.000/\gamma$ and obtain the steady-state solutions by time averaging over this single trajectory. The results of our numerical simulations are shown in Fig. 2, where the spin structure factor was obtained for a 4x4 lattice with various cluster sizes. It is clear that increased incorporation of quantum correlations present for larger cluster sizes significantly affects the spin correlations.

Both the $1\times 2$ and $2\times 2$ clusters show the existence of the ferromagnetic phase and show qualitatively the same behaviour as predicted by the CMF Jin1 and the GMC Casteels2 . The clusters however, are able to find a non-zero value for $S_{ss}^{xx}(0)$ for values of $J_{y}<0.9\gamma$. This behaviour is not captured by the single-site Gutzwiller ansatz Casteels2 or the mean-field Lee and thus originates from quantum correlations. In the GMC Casteels2 the ferromagnetic region becomes smaller with growing system size and the transition to the paramagnetic phase becomes sharper. The inset of Fig. 2 shows that by including clusters of size $1\times 2$ and $2\times 2$ this behaviour is captured already for smaller lattice sizes. Increasing the cluster size makes the sharpening steeper and occur for smaller values of $J_{y}$. This shows the importance of quantum correlations in the simulation of an open quantum system. This sharpening is what is to be expected when the size of the lattice grows (i.e. when the thermodynamic limit is approached).

By comparing the results for different cluster sizes with the exact solution of this lattice we see that for increasing cluster sizes the exact solution is approached more closely, but differences persist. It has to be noted that short-range quantum correlations are not enough to accurately approximate the exact solution for the $4\times 4$ system. It remains to be seen if this stays true for larger lattices and if longer-range quantum correlations have to be taken into account as well.

As mentioned earlier, for values of the parameter $J_{y}$ smaller than $0.9\gamma$ we find an unexpected buildup of spin-spin correlations. Where the mean-field theory predicts an all-zero steady-state spin structure factor we find a non-zero value by including clusters. This non-zero value is most pronounced in the exact solution. This behaviour is completely neglected by the classical mean-field theory and thus entirely driven by quantum fluctuations. The question remains whether a phase transition is present or not. It is however clear that short-range quantum correlations do not capture the exact behaviour for small lattices and long-range correlations have to be included for a more accurate description.

A comparison with the CMF used in Ref. Jin can show us the importance of spatial correlations between the clusters, as they are not captured by the CMF. Fig. 2 shows the steady-state spin structure factor for several cluster sizes on a $4\times 4$ lattice. We notice two distinct areas, again for $J_{y}<0.9\gamma$ and $J_{y}>0.9\gamma$ (we will resp. call them the left hand and right hand side). On the left hand side both the CMF and the CGMC match closely, giving a strong indication that only quantum correlations are important in this regime. On the right hand side of the figure however, the CMF and CGMC match only qualitatively, in this regime both quantum and classical spatial correlations contribute. This confirms the difference between the buildup of order on the left and right hand sides respectively.

In 1D we can take clusters with larger linear size. In the top panel of Fig. 3 we look at a $1\times 12$ lattice for which we can go up to $1\times 6$ clusters. We find that the behaviour of the steady-state spin structure factor is qualitatively captured by the cluster approach both in the left and right region for $1\times 6$ clusters. It has to be noted that the system has no phase transition in the thermodynamic limit in 1D, which has been shown by using a matrix product operator ansatz for the density matrix Jin . The behavior of the spin structure factor for finite size systems does however gives insight in the importance of the longer-range quantum correlations to describe the exact dynamics of the open quantum system. As one can see in Fig. 3, decreasing the size of the clusters results in a steady-state spin structure factor that differs completely from the exact value, even negative values are found for values of $J_{y}>\gamma$ where they should be positive. For values of $J_{y}<\gamma$ the influence of increasing the cluster size can be clearly observed. By including longer-range quantum correlations the exact behaviour is matched more closely. For $1\times 2$ clusters we find the same linear behaviour for the steady-state spin structure factor as found in the 2D case. By increasing the cluster size we can see a clear convergence to the same behaviour as the exact solution. Short-range quantum correlations are as such not sufficient for the description of the system and longer-range quantum correlations play an important role.

To further confirm the existence of the phase transition driven by quantum correlations we show the steady-state spin structure factor for the exact solution of a $2\times 2$, $3\times 3$ and $4\times 4$ lattice on the top panel of figure 4. As can be seen from this figure, and as is expected for finite size systems, the region where the phase transition occurs is smoothed out and one could suspect that there is only a continuous change of the order parameter rather than an actual phase transition. In the parameter region $J_{y}/\gamma\in]0.9,1[$ however the spin structure factor does decrease when the lattice size is increased ( Note that the steady-state spin structure factor is always zero for the values $J_{y}=J_{x}$ and $J_{y}=J_{z}$ due to the unitary dynamics conserving respectively the magnetization in the z-direction and x-direction Jin ). In order to check for the convergence as a function of increasing system size, we show in the bottom panel of Fig. 4 the behavior as a function of the system size, together with a fit to the power law dependence $S_{ss}^{xx}(0,L)=aL^{b}$, where $L$ is the number of points along one dimension of the $L\times L$ lattice. In all cases, we find a negative exponent $b$, which is compatible with a vanishing of the spin structure factor in the thermodynamic limit. For $J_{y}=0.95\gamma$ (in the middle of the interval), we find $b\approx-2.1$, close to the value $b=-2$, that is expected for a two-dimensional system with a finite correlation length.

V An applied magnetic field

The $\mathbb{Z}_{2}$-symmetry can be explicitly broken by applying a small magnetic field. In this section we will study the behaviour of the magnetization of the system as a function of the applied field in the $x$-direction and $y$-direction. An applied field $\vec{h}=h_{x}\vec{e}_{x}+h_{y}\vec{e}_{y}$ translates in adding a term $\hat{H}_{B}$ to the Hamiltonian $\hat{H}$ from (1)

[TABLE]

In Fig. 5 the magnetization in the $y$-direction is shown for a $4\times 4$ lattice with $1\times 2$ and $2\times 2$ clusters as a function of $h_{y}$ $(h_{x}=0)$ in the CMF. From the theory of closed quantum systems one would expect the magnetization to saturate at $\pm 1$. This however is not the case as can be seen in the figure, both the $x$ and $y$ magnetization go to zero when the field is increased. To obtain a closer understanding of this behaviour we note that it is also present in the mean-field approximation, for which we can find analytic expressions. These expressions will enable us to shed light on this behaviour.

The system of mean-field equations in the steady-state can be written as

[TABLE]

with $d$ the dimensionality, $\gamma$ the dissipation rate and $M_{x}$, $M_{y}$ and $M_{z}$ resp. the magnetization in the $x$, $y$ and $z$ direction. We look at the case where $h_{y}\neq 0$ and $h_{x}=0$ (the reverse situation is analogue). With these parameters one can rewrite the system of equations as an expression for $M_{x}$ and $M_{y}$ in terms of $M_{z}$

[TABLE]

[TABLE]

and an equation for $M_{z}$ of which the solution remains to be found by substituting the above solutions for $M_{x}$ and $M_{y}$ into the last equation in (9). This equation has no analytic solution and has to be solved numerically. We can however use the knowledge that $M_{z}\rightarrow 0$ as $h_{y}\rightarrow\infty$. A more close study shows that for a growing field $h_{y}$, $h_{y}M_{z}\rightarrow 0$. With these conditions we can approximate the third equation of (9) up till order $M_{z}$. We then find for large $h_{y}$

[TABLE]

This relation explains why $M_{x}$ and $M_{y}$ go to zero for big applied fields, rather than $\pm 1$. This is shown in Fig. 5 as the magenta line. Unlike in thermal equilibrium, the magnetization goes to zero for large magnetic field. The reason is that in the limit $|h_{y}|\rightarrow\infty$, the Zeeman term dominates the Hamiltonian, so that the eigenstates are simply the eigenstates of $\sigma_{y}$. The dissipation being in the orthogonal direction does not drive the system to the ground state, but rather destroys the coherence between the eigenstates. The resulting steady state is then the unit matrix, from which the zero magnetization follows.

For the special case of $J_{x}=J_{y}=J_{z}$ ($\gamma=1$) the system of mean-field equations is analytically solvable without any approximations:

[TABLE]

This is in agreement with the above result and is also true for small $h_{x}$ and $h_{y}$. We can conclude that a large magnetic field will cause the system to have no magnetization at all.

Note the difference in sign between $M_{x}$ and $M_{y}$ in this special case. This sign difference implicates that we cannot interchange $x$ and $y$ without introducing a sign change. At first sight this might appear confusing because of the identical parameters $J_{x}=J_{y}=J_{z}$. Only the $z$-direction is fixed by the dissipation and so there appears to be no clear reason for a distinction between $x$ and $y$. A closer look shows that one cannot interchange $x$ and $y$ because this changes the handedness of our coordinate system. This is reflected in the commutation relations of the Pauli matrices requiring that $\left[\hat{\sigma}^{(a)},\hat{\sigma}^{(b)}\right]=2i\epsilon_{abc}\hat{\sigma}^{(c)}$. To interchange $x\rightarrow y$ one could do $\hat{\sigma}^{(x)}\rightarrow\hat{\sigma}^{(x)}$ and $\hat{\sigma}^{(y)}\rightarrow-\hat{\sigma}^{(y)}$. This would however result in $\hat{\sigma}^{(z)}\rightarrow-\hat{\sigma}^{(z)}$ in order to satisfy the Pauli commutation relations. So $x$ and $y$ cannot be interchanged without changing the sign of $z$. This also results in the same coordinate system and thus no interchange was made in the end. As such there is no symmetry to transform $x$ into $y$ explaining why a sign difference can be present.

VI Susceptibility

A divergence in the susceptibility indicates the existence of a phase transition. To calculate the susceptibility we use the scheme presented in Rota and apply a field in the $xy$-plane which corresponds to adding a term (8) to the system Hamiltonian. The resulting magnetization is given by

[TABLE]

The susceptibility tensor can be extracted through

[TABLE]

With $M_{\alpha}$ the magnetization in the $x$ or $y$ direction. To obtain a scalar value the average angular susceptibility can be calculated as follows

[TABLE]

For a more elaborate discussion we refer to Rota .

When we calculate the susceptibility tensor for different cluster and lattice sizes we find that $\chi_{xy}\neq\chi_{yx}$. This is strikingly different from the case of closed systems, where the susceptibility is found from the free energy $F$, $\chi_{xy}=\frac{\partial^{2}F}{\partial h_{x}\partial h_{y}}=\frac{\partial^{2}F}{\partial h_{y}\partial h_{x}}=\chi_{yx}$.

We find this result even in the mean-field approximation. If we take the parameter values $J_{x}=J_{y}=J_{z}$ we find that $\chi_{xy}=-\chi_{yx}$ from (13). Numerical results show that for general coupling parameters the magnitudes differ and in general we find $|\chi_{xy}|\neq|\chi_{yx}|$ and reciprocity is broken.

For the 2D lattices we have two regions of interest. The earlier mentioned right hand region, studied in Casteels2 ; Lee ; Jin , and secondly the left hand region, discussed in section IV. In Fig. 6 the exact solution of the susceptibility for a $2\times 2$ and $3\times 3$ lattice is shown. Larger lattices are computationally not feasible and a more sophisticated method would have to be used, such as the Corner Space method Rota . The right hand region of the susceptibility has been studied in Rota . The presence of a peak in the susceptibility could indicate a phase transition. Note the ‘shoulder’ which is present on the left side. This could indicate the presence of a second peak, partially masked by the higher peak on the right hand side. These two peaks move away from each other when the lattice size is increased. The peak on the left however is not sharp and it remains to be seen whether it diverges for larger lattice sizes and whether a true phase transition is present. It should be noted that the positions of the shoulder do coincide with a sharp decrease of the steady-state spin structure factor. These findings show the possible presence of a phase transition which is entirely quantum in nature, which is completely missed when the quantum correlations are neglected.

VII Conclusions

We studied the dissipative XYZ Heisenberg model with the cluster Gutzwiller Monte Carlo. This method allows for the inclusion of short-range quantum correlations as well as classical spatial correlations. Calculation of the steady-state spin structure factor shows the appearance of a ferromagnetic phase and two paramagnetic phases also found in Ref. Jin1 and Ref. Casteels2 . We show the possible existence of another phase transition which is entirely quantum in nature and which is completely missed when quantum correlations are neglected. A calculation of the susceptibility tensor shows how reciprocity is broken, a feature not observed in closed quantum systems. We show that the inclusion of short-range quantum correlations causes the sharpening of the phase transition to occur for smaller lattice sizes. Even though short-range quantum correlations have a big impact on the phase diagram of the system it has to be noted that they only qualitatively match the exact results and long-range quantum correlations play an important role in the dynamics of the system.

Acknowledgements.

Discussions with C. Ciuti, F. Storme, W. Verstraelen and M. Van Regemortel are greatfully acknowledged. This work is supported by UAntwerpen/DOCPRO/34878. Part of the computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and the Flemish Government department EWI.