Dynamical properties of dissipative XYZ Heisenberg lattices

We study dynamical properties of dissipative XYZ Heisenberg lattices where anisotropic spin-spin coupling competes with local incoherent spin flip processes. In particular, we explore a region of the parameter space where dissipative magnetic phase transitions for the steady state have been recently predicted by mean-field theories and exact numerical methods. We investigate the asymptotic decay rate towards the steady state both in 1D (up to the thermodynamical limit) and in finite-size 2D lattices, showing that critical dynamics does not occur in 1D, but it can emerge in 2D. We also analyze the behavior of individual homodyne quantum trajectories, which well reveal the nature of the transition.


I. INTRODUCTION
Quantum many-body physics with light has proved to be an extremely rich and interesting field of study, as it combines the complexity of condensed matter with the intrinsically out-ofequilibrium behavior of optical systems [1][2][3][4]. Collective phenomena among photons, such as Bose-Einstein condensation [5][6][7][8] or superfluidity [9][10][11][12][13], have been observed in planar semiconductor microcavities in the strong light-matter coupling regime. In these systems, the optical confinement and the nonlinearity of the media give rise to a weak photon-photon interaction, which allows the many-photon system to behave as a quantum fluid.
The appearance of strongly correlated states of light is even more evident in regimes where the interaction among photons becomes large. When the nonlinearity of the optical cavity is much larger than its dissipation rate, the presence of a single photon inside the cavity is able to effectively block the entrance of a second one. This effect, known as photon-blockade [14,15], has been observed experimentally at first with optical photons using a single atom in a cavity [16] and is particularly strong in circuit quantum electrodynamics systems in the microwave domain [17]. Non-trivial phases can also arise when several cavities are coupled together and form a lattice of resonators [18]. For instance, correlations can lead to a transition from a photonic Mott insulator to a superfluid [19][20][21][22][23], similar to that observed with ultracold atoms confined in optical lattices [24,25]. Interestingly, a system of coupled resonators in the photon-blockade regime arranged according a lattice geometry can be mapped into an effective spin model [21,26,27]. This class of systems can be realized nowadays using different experimental platforms, such as superconducting quantum simulators [28] or Rydberg atoms [29][30][31][32].
Among the collective phenomena appearing in coupled photonic lattices, dissipative phase transitions are nowadays deserving more and more attention. Dissipative processes are usually at odds with the unitary Hamiltonian evolution of the quantum system and the competition between the incoherent and the coherent dynamics can give rise to criticality for the steady state in the thermodynamic limit [33]. Dissipative phase transitions have been discussed theoretically for single cavity photonic systems [34][35][36], as well as for lattices of cavities with mean field methods [37][38][39] or full-size lattice simulations [40,41]. An experimental observation of these critical phenomena seems feasible with state-of-the-art techniques, and some remarkable results have already been obtained [42][43][44].
In this context, the dissipative XYZ Heisenberg model [45] has attracted a considerable attention. It describes a lattice of spins interacting via an anisotropic Heisenberg Hamiltonian coupled to an environment which forces spins to align along the z-axis. The single-site Gutzwiller meanfield theory predicts a rich phase diagram for the magnetic properties of the steady state of this model [45]. More refined calculations [46][47][48][49][50], based on numerical methods including many-body correlations, have confirmed the emergence of a critical behavior in two-dimensional lattices, while the phase transition disappear when the spins are arranged according to a one-dimensional geometry. All these works, however, focussed on the calculation of steady-state properties and a full description of the dynamics of the system is still lacking.
In this work, we explore the dynamical properties of the dissipative XYZ model in the region where a second-order phase transition from a paramagnetic to a ferromagnetic steady state has been predicted. For finite-size 1D arrays and 2D lattices, we have performed an exact integration of the master equation using the whole Hilbert space via the Wave Function Monte Carlo method [51]. Moreover, for 1D arrays of infinite length we have applied the infinite Matrix Product Operator (iMPO) technique [52,53].
This article is organized as follows. In Sec. II we discuss the theoretical framework and describe the methods used for the calculations. In Sec. III we show the main results of the work. In Sec IV we draw our conclusions and present some perspectives.
whereσ α i (α = x, y, z) represent the Pauli matrices acting on the i-th site. The sum runs over the nearest neighbour sites i, j . The dissipative part describes incoherent spin-flip processes which tend to align a single spin towards the negative direction of the z-axis with a rate γ. The density matrixρ(t) dynamics is obtained from the Lindblad master equation whereσ ± j = (σ x j ± iσ y j )/2 are the spin raising and lowering operators acting on the j-th spin and L is the Liouvillian superoperator. The latter is non-Hermitian and has a spectrum of complex eigenvalues, defined by the equation L[ρ r ] = λ rρr . The dissipative XYZ model evolves towards a steady stateρ ss , which depends on the parameters in (2) and corresponds to the zero eigenvalue of L (∂ tρss = L[ρ ss ] = 0). All the other eigenvalues λ r are such that their real part is negative and describe the relaxation dynamics ofρ(t) towards the steady state. Since a dissipative phase transition is expected to be characterized by a critical slowing down in the dynamics of the system, a particular relevance has to be given to the so-called Liouvillian gap λ = min r |Re(λ r )|, which is also called asymptotic decay rate [33]. The emergence of a critical behavior is associated to a closing of the Liouvillian gap in the thermodynamic limit [33,41,54].
The Lindblad master equation (Eq. 2) is invariant under a π-rotation of all the spins around the z-axis (σ x i → −σ x i ,σ y i → −σ y i ∀i). In the thermodynamic limit, the Z 2 symmetry associated to this transformation may spontaneously break, resulting in the appearance of several magnetic phases for the steady state of the model. In this work, we will focus on a particular regime where previous calculations have predicted a transition from a paramagnetic phase with no magnetization in the xy plane ( σ x = Tr(ρ ssσ x j ) = 0 , σ y = Tr(ρ ssσ y j ) = 0) to a ferromagnetic phase with finite magnetization in the xy plane ( σ x = 0 , σ y = 0) [45][46][47][48][49][50].
From a computational point of view, the numerical solution of the master equation (2) is a formidable task when considering extended lattices. The corner-space renormalization method [55], which has shown the criticality of several steady-state observables in 2D lattices [47], does not give access to the dynamic properties of the system. For small systems with a number N < 10 spins, the problem can be solved via a standard Runge-Kutta integration of Eq. (2). For 10 ≤ N ≤ 16, instead, we have solved the master equation stochastically via the Wave Function Monte Carlo method [51]. This method describes the time evolution of the open quantum system in terms of a set of N T pure states |Ψ k (t) (usually called quantum trajectories), obtained independently according to a stochastic evolution protocol [56][57][58][59]. The density matrix is retrieved by averaging over the N T sampled trajectories, according to the formulaρ The computational advantage of this method is clear, as it allows to study the evolution of the open system dealing with pure states (which are vectors of size 2 N ), instead of the density matrix (which has size 2 N × 2 N ).
It is important to notice that quantum trajectories are useful not only to reduce the complexity of the integration of the Lindblad master equation (2), but their analysis of can also provide insightful results about the nature of the dissipative phase transition [40,41]. To this aim, we have investigated the stochastic evolution of individual quantum trajectories for the dissipative XYZ model, obtained according to the homodyne protocol described by the following equation: and dW j are stochastic Wiener increments with zero expectation value, variance equal to √ dt and uncorrelated among the different spins (the detailed derivation can be found, e.g., in [58]). Contrarily to the master equation (2), the stochastic equation in (3) does not conserve the Z 2 symmetry of the Liouvillian superoperator, due to the presence of the terms s j (t). Therefore, by studying the time evolution of the magnetic order parameter over an individual quantum trajectory, it is possible to reveal the emergence of different magnetic phases, when we change the parameters of the system. Nevertheless, the symmetry of the Liouvillian is restored when we consider the density matrix, obtained by averaging over many trajectories.
Alternative approaches for the simulation of 1D arrays are based on tensor networks techniques [60] making use of the Matrix Product Operator (MPO) ansatz for the density matrix [61,62] (see for example Refs. [38,46,[63][64][65][66]). The MPO ansatz for the many-body mixed state can be controlled by changing a single parameter, i.e. the bond-link dimension χ: the more χ increases, the more non-local quantum correlations can be encoded. The dynamics of the open system is obtained via a time-evolving block decimation scheme [67,68]. In the case of translational invariant systems, the MPO ansatz and the time evolution procedure can be further simplified leading to the infinite MPO (iMPO) representation [52,53], which allows to directly access the thermodynamic limit of an infinite number of sites. Very recently, this technique has been extended to the case of 2D lattices [48] although with a very reduced bond dimension.

III. RESULTS AND DISCUSSION
We start our discussion on the dynamics of the dissipative XYZ model by studying the time evolution of the average lattice magnetization M x (t) = i Tr [ρ(t)σ x i ] /N , N being the number of spins in the lattice. In Fig. 1, we plot M x (t) for a fixed choice of the parameters of the Hamiltonian (1) in vicinity of the critical point, for spin systems of different size, both with 1D ( Fig. 1-(a)) and 2D geometry ( Fig. 1-(b)). In all these calculations, the master equation has been solved assuming an initial configuration where all the spins point along the positive direction of the x-axis (therefore M x (t = 0) = 1) and imposing periodic boundary conditions to the finite-size lattice.
For t 5γ, all the curves M x (t) decay exponentially towards the steady-state expectation value M x ss = 0 (notice that we have M x ss = 0 for all the values of the parameters since we do not break explicitly the Z 2 symmetry of the Liouvillian superoperator in our simulations). The presence of an asymptotic exponential behavior for M x (t) indicates that, at large times, the dynamics of the system can be described uniquely in terms of the eigenstate associated to the Liouvillian gap. The density matrix can be approximated asρ(t) =ρ ss + Aρ 1 e −λt , where A is a real number depending on the choice of the initial configuration. From our results, we notice also that the dynamics gets slower when increasing the size of the system, both in 1D arrays and in 2D lattices (respectively Fig.  1(a) and Fig. 1-(b)). In 1D arrays the decay rate saturates when the size of the system increases. For an array with 16 sites the decay curve is nearly indistinguishable from what obtained for an array of infinite length (obtained via the iMPO technique). Instead, in 2D lattices no saturation of the decay rate is observed.
By fitting the curves for M x (t) at large t with a simple exponential, we can extract the value of the Liouvillian gap λ. The results for λ obtained with this procedure have been successfully benchmarked against those calculated with an exact diagonalization of the Liouvillian superoperator in small systems (4 × 1 array and 2 × 2 lattice). In Fig. 1-(c,d) we plot λ as a function of the normalized coupling parameter J y /γ (the other coupling parameters J x /γ and J z /γ are kept fixed). Both in the 1D and in the 2D case, all the curves λ(J y ) present a minimum close to the critical value of J c , indicating a slowing down in the dynamics of the system. Nevertheless, we clearly notice that this slowing down is not critical in 1D systems. Indeed, the results for λ(J y ) in the largest 1D systems (with N ≥ 12) overlap and are in good agreement with the prediction for the infinite array obtained with iMPO [69], showing a finite value of the Liouvillian gap. Instead, in 2D systems, the minimum of λ(J y ) becomes smaller and smaller when the size of the lattice increases. This behavior is consistent with a closure of the Liouvillian gap in the thermodynamic limit.
In order to better characterize the behavior of the 2D system across the critical point, we study the average magnetization of the lattice M x Ψ (t) = Ψ(t)| i σ x i |Ψ(t) along a single trajectory |Ψ(t) . To this extent, we have computed |Ψ(t) following the homodyne protocol in Eq.  In the three panels of Fig. 2, we show the results for M x Ψ (t) in a 3 × 3 lattice for J y = 0.95γ, J y = 1.25γ and J y = 1.8γ. When the steady state presents a paramagnetic phase (J y = 0.95γ, Fig. 2-(a)), the curve for M x Ψ (t) presents only small fluctuations around the zero value for the magnetization. The behavior of the quantum trajectory is strikingly different in the ferromagnetic phase (J y = 1.25γ, Fig. 2-(b)). In this case, we can clearly distinguish intervals of time where the curve for M x Ψ (t) fluctuates around a positive value of the magnetization and others where it fluctuates around the opposite value. The duration of these time intervals is of the order ∆t ∼ λ −1 . Finally, for large values of the coupling parameter J y (J y = 1.8, γ, Fig. 2-(c)), M x Ψ (t) presents yet another different behavior. It is reminiscent of what observed in the paramagnetic phase (see Fig.  2-(a)), since it fluctuates around the zero value of the magnetization, but the amplitude of the fluctuations is much larger than in the regime J x J y . This peculiar behavior can be ascribed to the strongly mixed and correlated nature of the steady state of the model in this regime.
To better understand the nature of those three regimes, we studied the probability distribution of M x Ψ (t) over many trajectories, which we will call p(M x ), defined as follows. We consider a time t s where the density matrix of the system has reached the steady state, and statistically collect all the values of M x Ψ (t) for t > t s over many trajectories. The results for p(M x ) are presented in the top panel of Fig. 3, as a function of the coupling J y . We notice that for small J y the distribution is monomodal around zero. As J y increases, one reaches a point J c 1.05γ where p(M x ) starts to present two distinct peaks, which are symmetric around the value M x = 0. If we continue to increase J y , the two peaks broaden and they move apart, until they reach their maximum distance for J y 1.2γ. Above this value of J y , the peaks continues to broaden and they start to approach one to the other, until they merge again into a single peak for J y 1.6γ. The broadening, the separation and the merging of the peaks in the probability distribution is even more evident in the panels in Fig. 3-(a,f), where we plot the curves for p(M x ) for some values of the coupling parameter J y .
In order to perform a more quantitative analysis of the distribution p(M x ), we compute the bimodality coefficient b [70] that for an even distribution reads: b is an indicator of the bimodal character of the distribution, which in the present study is related to the ferromagnetic nature of the steady state. Indeed, when p(M x ) presents two narrow peaks, then the quantity b approaches its maximum value b max = 1. Instead, unimodal distributions are characterized by smaller values of b (for instance, a Gaussian distribution centered at M x = 0 would have b = 1/3). In Fig. 4, we plot the value of b as a function of J y , for different sizes of the 2D lattice. The emergence of the phase transition at J y /γ 1.05 is signaled by a steep increasing of the ratio b, which is almost independent of the lattice size. Furthermore, the decreasing of b for J y /γ > 1.2 indicates the disappearance of the ferromagnetic order for large anisotropies. In this case, however, the drop of b is not particularly sharp and tends to become smoother and smoother as the size of the lattice increases.
The study of the behavior of b(J y ) is interesting to address the open question about the nature of the steady state of the dissipative XYZ model for large anisotropies. Several works in literature [46,48,49] have predicted a ferromagnetic to paramagnetic phase transition for J y /γ > 1.5. However, the critical value of J y for this second transition depends strongly on the method used and on the size of the cluster considered in the calculation [46,48,49]. Moreover, the behavior of the magnetic susceptibility and of the von-Neumann entropy as a function of J y do not present any feature signaling the emergence of a critical point for J y > 1.2γ [47]. Our results in Fig. 4, showing a smooth decreasing of b at large J y , together with the absence of a slowing down J y > 1.2γ (see Fig. 1-(d)), suggest that the disappearance of the ferromagnetic order for large anisotropies might be due to a crossover and not to another second-order phase transition.

IV. CONCLUSIONS
In this paper we investigated numerically the dynamics of a dissipative spin-1 2 lattice interacting through an XYZ-Heisenberg Hamiltonian. This model is particularly relevant in the context of strongly correlated open quantum systems since it is known to support a second-order dissipative phase transition in two dimensions, associated with the breaking of the Z 2 symmetry.
By performing stochastic quantum trajectories simulations on finite-size systems, we determined the Liouvillian gap from the asymptotic decay rate of the dynamics towards the steady state. When the system is driven across the critical point, we found that the relaxation exhibits a slowing down. For 1D systems, the Liouvillian gap remains finite as the length of the chain is increased up to the thermodynamical limit, thus indicating the absence of a phase transition. Instead, results for 2D lattices do not show a saturation of the Liouvillian gap, which is consistent with the emergence of critical slowing down. By analyzing individual stochastic homodyne trajectories in 2D lattices, we characterized the emergence and disappearance of two metastable states with opposite magnetization. Our predictions might be tested in quantum simulators based on superconducting quantum circuits or Rydberg atoms. As a perspective, the effects of disorder on the dynamics of these systems is a very interesting aspect that needs to be investigated in the future.