Quantum correlations and limit cycles in the driven-dissipative Heisenberg lattice

Driven-dissipative quantum many-body systems have attracted increasing interest in recent years as they lead to novel classes of quantum many-body phenomena. In particular, mean-field calculations predict limit cycle phases, slow oscillations instead of stationary states, in the long-time limit for a number of driven-dissipative quantum many-body systems. Using a cluster mean-field and a self-consistent Mori projector approach, we explore the persistence of such limit cycles as short range quantum correlations are taken into account in a driven-dissipative Heisenberg model.


Introduction
Understanding the phases of quantum many-body systems is one of the central goals of modern physics. Phases of matter emerging from cooperative behaviour in equilibrium systems have proven to be of fundamental and technological importance, with notable examples including superconductors [1] and topological materials [2]. Recent advances have provided the opportunity to extend this field into the exploration of the phase diagrams of non-equilibrium quantum systems where excitations which dissipate from the system are replenished using an external driving field [3][4][5]. Experimental platforms, such as cavity arrays, superconducting circuits and polariton waveguides, have introduced a new class of systems where the interplay between coherent driving and incoherent dissipation has led to the discovery of novel phenomena. Bistability [6,7] and crystallisation [8] in the driven-dissipative nonlinear resonator arrays and synchronised switching in an array of coupled Josephson junctions [9] provide a couple of examples where non-equilibrium phenomena are essential for the understanding of quantum photonic systems.
An intriguing possibility of a non-equilibrium phase that appears in driven-dissipative systems are limit cycles, whereby the system enters a periodic trajectory which breaks the time-translation symmetry of the master equation. Mean-field studies suggest that limit cycles could exist in a range of non-equilibrium systems, including optomechanical arrays [10], anisotropic Heisenberg lattices [11,12], Rydberg lattices [13,14] and Bose-Hubbard lattices with cross-Kerr interactions [15,16]. Experimentally realising limit cycles would not only be the discovery of a new class of phases in driven dissipative quantum many-body systems but could also have important technological applications, for example in synchronising quantum many-body devices [10,17].
Existing predictions of the occurrence of limit cycles are almost exclusively based on Gutzwiller mean-field approaches, which assume a factorised density matrix and ignore quantum fluctuations. It is therefore important to investigate to what extent these limit cycles are affected by quantum fluctuations or correlations [18], which often play a significant role in determining the structure of exotic materials [8,19]. Recently, Chan Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. et al have shown that limit cycles persist in the anisotropic Heisenberg lattice even when Gaussian fluctuations are taken into account [12].
In this work, we present simulations of the driven-dissipative anisotropic Heisenberg lattice using the selfconsistent Mori projector [20] and cluster mean-field methods [21] in order to explore the role of short-range fluctuations beyond the Gaussian approximation. Within the limits of our approximation limit cycles in the Heisenberg model disappear for reduced dimensionality and increasing cluster size respectively. Both methods implicitly include fluctuations beyond the Gaussian approximation [20,[22][23][24] demonstrating that higherorder correlation functions have an important influence on the existence of limit cycles. The distinct approaches of these two numerical methods in simulating local quantum correlations complement each other and both show a disappearance of the limit cycle phase of the driven-dissipative anisotropic Heisenberg model at low coordination numbers.

Model
We investigate the long-time behaviour of the driven-dissipative anisotropic Heisenberg lattice, the phase diagram of which was first studied in [11]. The system consists of a regular d-dimensional lattice of two-level sites with an energy splitting ω 0 which are coherently driven by an external drive field of strength Ω and frequency ω D . We consider such frequency to be resonant with the energy splitting of the two level system, that is ω D =ω 0 , make the rotating-wave approximation and move to a frame that rotates at frequency ω 0 , such that the uncoupled Hamiltonian is given by z are the Pauli matrices acting on site i and we have set ÿ=1. Each site has z=2d nearest neighbours which are coupled by an anisotropic Heisenberg term such that the Hamiltonian of the full system is given by where z is the coordination of the lattice and á ñ i j , denotes nearest-neighbour interactions. The factor of z −1 in the coupling term is required in order to ensure that the energy of the system is extensive. Additionally, individual sites can spontaneously decay from the excited state to the ground state at a rate γ. This gives rise to a Markovian dynamics that is ruled by the following master equation in the Lindblad form å g s s s s = -+ - where R(t) is the density matrix of the whole system and σ ± i is the annihilation/creation operator for site i. Avenues to realising this model in experiments with Rydberg atoms and trapped ions have been discussed in [11].

Mean-field phase diagram
Obtaining a solution for equation (3) is impractical under almost all circumstances and approximations must be made in order to determine the phase diagram of the model. The simplest approximation is to ignore quantum correlations between the individual two-level systems and treat interactions as though each site is coupled to the mean-field generated by its nearest neighbours. Linear stability analysis of the fixed points of the resultant meanfield master equation was performed by Chan et al [12] and, for completeness, we summarise some of their results here.
The mean-field phase diagram indicates that an antiferromagnetic phase can be realised in this system. Therefore, in order to allow the density matrix to break the translational symmetry of the model, we divide the sites into an A and B sublattice where all sites on the A sublattice interact only with the B sublattice and vice versa. In the mean-field approximation, the equations of motion for the reduced density matrices of the A and B sublattice can then be simplified by replacing the anisotropic Heisenberg coupling with an effective selfconsistent classical field generated by the site's z nearest neighbours thus excluding quantum correlations. In figure 1, we reproduce a part of the mean-field phase diagram. Within this approximation, the model supports uniform and antiferromagnetic phases which would also be expected in an equilibrium system. However, the inclusion of driving and dissipation allows for regions of the phase diagram where the system enters a limit cycle [11,12], see region LC in figure 1.
In the limit cycle phase, the sublattice symmetry of the system is broken and the local system observables of the two sublattices oscillate periodically with a relative phase π. The limit cycles exist even though the mean-field Liouvillian is time-independent. Note that the breaking of the sublattice symmetry is due to the instability of the mean-field solution and is not an artefact of assuming a bipartite lattice [11,21]. In this paper, we explore how the limit cycles predicted by these mean-field calculations are affected by quantum fluctuations and correlations.

Methods
We explore the existence of limit cycle phases for the model specified in equation (3) via two methods, selfconsistent Mori projectors and cluster mean field.

Self-consistent Mori projector method
The self-consistent Mori projector approach solves for the reduced density matrices of individual lattice sites by integrating out the correlations to give non-Markovian equations of motion for the reduced density matrices of the system [20]. For the system described by equation (3), we start by partitioning the Liouvillian into local and interaction parts The equation of motion for the reduced density matrix ρ n of site n, derived in [20], can be expressed as a Dyson series in  I and is given by . The long-time behaviour of the system is separable into antiferromagnetic (AFM), various uniform (U 1 and U 2 ) and limit cycle (LC) phases. Bistable regimes where both phases can be reached depending on the system's initial state are also present. n . The first term in equation (9) describes the free evolution of the nth site whilst the second term is the interaction of the site with the mean field. The third term is referred to as the Born term and is the first-order quantum correction to the mean-field prediction for the dynamics of ρ n (t). The final term is the sum over the remaining terms of the self-consistent Mori projector Dyson series. Its explicit form can be found in appendix C of [20]. In order to make simulation of equation (9) tractable, we make a truncation of the Dyson series expansion by setting  = ( ) t 0 m n for m3, corresponding to a Born approximation. Note that even with this truncation, equation (9) was found to give more accurate results than standard perturbation theory to second order, see figure 3(d) in [20] where a second order expansion in correlators [22,25] is compared to self-consistent Mori projector results in the Born approximation.
With the truncation of equation (9) at second order in  I , we partition the sites onto A and B sublattices such that Tr e , 12 and similarly for ρ B (t). For simplicity, we have omitted the site index for the Pauli operators, which operate on the relevant sublattice. The long-time behaviour of the system can then be calculated by time-integrating the equations of motion over a sufficiently long time such that the transient behaviour has disappeared.

Cluster mean-field theory
Approximate solutions to the master equation become more accurate when reduced density matrices for clusters consisting of multiple lattice sites are considered [20,21]. Quantum correlations within these clusters are then calculated exactly and inaccuracies are limited to interactions between clusters. The exponential scaling of the dimension of the Hilbert space with increasing coordination number unfortunately makes it impossible to consider clusters for lattices with large z but, for low coordination number, simulating the reduced density matrix of a cluster becomes more manageable as the size of the cluster's Hilbert space is reduced. However, when using the self-consistent Mori projector method for this model, evaluating the Born term in equation (9) for a cluster is computationally difficult due to the large number of terms in the interaction Liouvillian. Nevertheless, even in a mean-field calculation, all correlations are taken into account for the internal quantum dynamics of the cluster which interacts with the mean-field exerted by its neighbouring clusters.
In the cluster mean-field approximation, the density matrix of the lattice is divided into a product state of contiguous clusters of sites  which are identical due to the translational symmetry of the lattice This density matrix then evolves according to the decoupled cluster mean-field Liouvillian which can be written as where     = å Î j j describes the evolution of the isolated cluster and the interaction with the mean-field of the neighbouring clusters is described by the nonlinear Liouvillian    ( ) which acts only on sites at the boundary of the cluster. For the driven-dissipative Heisenberg model, the boundary Liouvillian is given by is related to the average polarisation of the sites adjacent to the boundary   ( ) at time t. Clusters containing an odd number of sites break the bipartite symmetry of the lattice and a pair of complementary coupled clusters must be used. For example, in order to perform calculations using a cluster size of 3 × 3, the lattice can be divided into two subclusters  A and  B where the four sites in the corners and the central site of  A are assigned to the A sublattice and  B is the complement of  A .

Results
We apply the methods described in the previous section to find the long-time behaviour of the driven-dissipative Heisenberg lattice. At t = 0, the system is prepared in a product state, that we parametrise as  where a n A and a n B are the components of two vectors within the unit sphere. For this initial state, we timeintegrate the equations of motion of the respective reduced density matrices using both the self-consistent Mori projector method and cluster mean-field theory until t?γ −1 , where the transient behaviour due to the product state initialisation has decayed and the system has either entered into a limit cycle or has reached a timeindependent steady state. We performed the simulations for two sets of parameters: x y z with γ=1 which both correspond to points in the parameter space where mean-field calculations predict the steady state to be composed of limit cycles.

Self-consistent Mori projector method
We start by presenting the results calculated using the self-consistent Mori projector method. The system exhibits two distinct long-time behaviours which we show in figure 2, where we characterise the limit cycle using the average polarisation of the sublattices s r { } Tr z A and s r { } Tr z B . For some initial states, after a transient behaviour, the sublattice symmetry is restored and the system relaxes to a ferromagnetic stationary state. However, for the other initial states, the system enters a limit cycle where the sublattice symmetry is broken and the system oscillates anharmonically.
Whether a random initial state will enter the limit cycle phase depends strongly both on the system parameters and the coordination number of the lattice. Figures 3 and 4 show the proportion of states which do not enter the limit cycle phase as a function of the coordination number z from a sample of 50 initial states with randomised a n A and a n B . In the limit  ¥ z , the equations of motion become equivalent to the mean-field approximation as the prefactor of the Born term in equation (10) vanishes so all initial states enter into limit cycles. Whilst this limit is unlikely to be experimentally practical, it does allow us to connect the mean-field result to more accurate investigations as the coordination number of the system controls the magnitude of the Born term and all higher order terms  ( ) t m n with m3 in equation (9). As the coordination number of the lattice decreases, quantum correlations included within the Mori projector expansion become relevant and certain initial states will evolve towards a stationary ferromagnetic state. For the above parameters, we find that below a critical coordination number z * , the limit cycles are absent and there is no evidence of long-time oscillatory behaviour for all initial states. The critical coordination number differs depending on the system parameters but * > z 10 in both cases. In contrast to x y z all initial states enter limit cycles for * > z z . In the mean-field limit, we see that for these parameter sets, every initial state is attracted to a limit cycle. As the coordination number of the system decreases, a proportion of these initial states instead are attracted to a time-independent steady state. For * < z z , our numerical results indicate that the size of the attractor basin for the limit cycle phase disappears and all initial states converge on a stationary solution. This result indicates that it is not possible to enter into a limit cycle phase for experimentally realistic coordination numbers.
Whilst this transition is accompanied by a shift in the frequency of the limit cycle, it is not possible to perform a quantitative analysis of this shift as higher-order terms in equation (9) may become relevant. For the . Below a critical dimension * » z 100, the limit cycle phase disappears and all initial states relax to a paramagnetic stationary state. This transition is accompanied by a shift in the frequency of the limit cycle. The error in the frequency is due to computational limitations which restricted the period of time over which the density matrix could be evolved whilst the uncertainty in the proportion of initial states not entering a limit cycle is taken from the standard error of a Bernoulli process. . The critical dimension is lower than for the parameters used in figure 3 with * » z 50, but below this value, all initial states still relax to a paramagnetic stationary state. For * > z z , all initial states enter into a limit cycle. parameters considered here, the truncated terms in the Dyson series given in equation (9)  1 for α=x, y and z. Hence, the difference in behaviour of the limit cycle frequency between the two parameter sets presented here for coordination numbers close to the critical value z * cannot necessarily be expected to be a quantitatively reliable prediction. While we cannot exclude the possibility that higher order terms would restore the limit cycle behaviour, we do not expect that expanding equation (9) to first order in the couplings-which is the mean-field approximation-will be more accurate than the expansion to second order that we consider here.

Cluster mean-field theory
Next we investigate the existence of limit cycle phases in the driven-dissipative anisotropic Heisenberg model via cluster mean-field theory. Figure 5 shows the average polarisation of the A and B sublattices for a twodimensional lattice calculated using cluster mean-field simulations of 2×2 and 3×3 clusters. The wave function was initiated in the product state given by equation (18) Figure 5 shows that the limit cycles are not observed in the cluster mean-field simulations. Such conclusions extend to cluster mean-field simulations for a three-dimensional lattice. In figure 6, we show the average polarisation for a 2×2×2 cluster. Once again, the system relaxes into a stationary steady state after an initial transient and limit cycles are not observed. For , , , 6.4, 3.0, 6.0, 2.25 x y z , cluster mean-field theory shows that the limit cycle phase is more robust as the system still exhibits periodic behaviour for a 2×2 cluster (see figure 7). However, for a 3×3 cluster, once again the system relaxes to a stationary steady state.

Conclusions
We have used the self-consistent Mori projector and cluster mean-field methods to simulate the evolution of a driven-dissipative anisotropic Heisenberg model within a regime where mean-field results predict that this system should exhibit limit cycles. The approximations required for these two methods are limited to calculating local properties of the system but have the advantage of including higher-order quantum correlation effects. For a high coordination number, the self-consistent Mori projector method agrees with the mean-field prediction and limit cycles are observed for any initial state. However, as the coordination number of the lattice is reduced, a proportion of the initial states no longer enter the limit cycle phase and relax towards a ferromagnetic steady state. Below a critical coordination number, for which *  z 10 for the parameter sets which we have studied here, the limit cycle phase disappears and all initial states relax towards the same steady state. Recently limit cycles in many-body systems were associated with the appearance of time-crystalline behaviour [26]. Our work seems to indicate that these effects may appear only in long-range systems or for high dimensions.   x y z for which the system was initiated in the product state R II (0) as defined in figure 2. The limit cycles are more robust than in figure 5, as they persist for 2×2 clusters, but the system relaxes to a stationary steady state for 3×3 clusters.