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

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–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 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 self-consistent 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–24] demonstrating that higher-order 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.

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 ω₀ 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 = ω₀, make the rotating-wave approximation and move to a frame that rotates at frequency ω₀, such that the uncoupled Hamiltonian is given by

$$\begin{eqnarray}&&{H}_{0}=\displaystyle \frac{{\rm{\Omega }}}{2}\displaystyle \sum _{i}{\sigma }_{i}^{x},\end{eqnarray} \tag{ 1 }$$

where ${\sigma }_{i}^{\alpha }=\{{\sigma }_{i}^{x},{\sigma }_{i}^{y},{\sigma }_{i}^{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

$$\begin{eqnarray}&&H={H}_{0}+\displaystyle \sum _{\langle i,j\rangle }\displaystyle \sum _{\alpha }\displaystyle \frac{{J}_{\alpha }}{z}{\sigma }_{i}^{\alpha }{\sigma }_{j}^{\alpha },\end{eqnarray} \tag{ 2 }$$

where z is the coordination of the lattice and $\langle i,j\rangle$ denotes nearest-neighbour interactions. The factor of z⁻¹ 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

$$\begin{eqnarray}&&\dot{R}(t)=-{\rm{i}}[H,R(t)]+\displaystyle \frac{\gamma }{2}\displaystyle \sum _{i}(2{\sigma }_{i}^{-}R(t){\sigma }_{i}^{+}-\{{\sigma }_{i}^{+}{\sigma }_{i}^{-},R(t)\}),\end{eqnarray} \tag{ 3 }$$

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].

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 mean-field 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 self-consistent classical field generated by the site's z nearest neighbours

$$\begin{eqnarray}&&\displaystyle \sum _{\langle i,j\rangle }\displaystyle \sum _{\alpha }\displaystyle \frac{{J}_{\alpha }}{z}{\sigma }_{i}^{\alpha }{\sigma }_{j}^{\alpha }\to \displaystyle \sum _{\alpha }{J}_{\alpha }{\sigma }_{A}^{\alpha }\langle {\sigma }_{B}^{\alpha }\rangle \end{eqnarray} \tag{ 4 }$$

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.

Zoom In Zoom Out Reset image size

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.

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

4.1. 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

$$\begin{eqnarray}&&\dot{R}(t)={{ \mathcal L }}_{0}R(t)+{{ \mathcal L }}_{I}R(t),\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}&&{{ \mathcal L }}_{0}R(t)=\displaystyle \sum _{j}-\displaystyle \frac{{\rm{i}}{\rm{\Omega }}}{2}[{\sigma }_{j}^{x},R(t)]+\displaystyle \frac{\gamma }{2}\displaystyle \sum _{j}(2{\sigma }_{j}^{-}R(t){\sigma }_{j}^{+}-\{{\sigma }_{j}^{+}{\sigma }_{j}^{-},R(t)\})\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\equiv \displaystyle \sum _{j}{{ \mathcal L }}_{j}R(t)\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{{ \mathcal L }}_{I}R(t)=\displaystyle \sum _{\langle i,j\rangle }\displaystyle \sum _{\alpha }\displaystyle \frac{{J}_{\alpha }}{z}[{\sigma }_{i}^{\alpha }{\sigma }_{j}^{\alpha },R(t)].\end{eqnarray} \tag{ 8 }$$

The equation of motion for the reduced density matrix ρ_n of site n, derived in [20], can be expressed as a Dyson series in ${{ \mathcal L }}_{I}$ and is given by

$$\begin{eqnarray}{\dot{\rho }}_{n}(t) & = & {{ \mathcal L }}_{n}{\rho }_{n}(t)+\displaystyle \sum _{m=1}^{z}\displaystyle \sum _{\alpha }\displaystyle \frac{{J}_{\alpha }}{z}{\mathrm{Tr}}_{\not\hspace{-2pt}{n}}\left\{\left[{\sigma }_{n}^{\alpha }{\sigma }_{m}^{\alpha },\mathop{\displaystyle \bigotimes }\limits_{k}{\rho }_{k}(t)\right]\right\}\\ & & +{\displaystyle \int }_{{t}_{0}}^{t}{\rm{d}}{t}^{\prime} {\mathrm{Tr}}_{\not\hspace{-2pt}{n}}\left\{{{ \mathcal L }}_{I}{{\rm{e}}}^{{{ \mathcal L }}_{n}(t-t^{\prime} )}{{ \mathcal C }}_{t^{\prime} }{{ \mathcal L }}_{I}\mathop{\displaystyle \bigotimes }\limits_{k}{\rho }_{k}(t^{\prime} )\right\}\\ & & +\displaystyle \sum _{m=3}^{\infty }{{ \mathcal Y }}_{m}^{n}(t),\end{eqnarray} \tag{ 9 }$$

where ${\mathrm{Tr}}_{\not\hspace{-2pt}{n}}\{\cdot \}$ denotes the trace over all degrees of freedom except for those of site n, ${P}_{t}^{n}(\cdot )={\mathrm{Tr}}_{\not\hspace{-2pt}{n}}\{\cdot \}$ $\otimes \,{\bigotimes }_{m\ne n}{\rho }_{m}(t)$ are the time-dependent Mori projectors and ${{ \mathcal C }}_{t}={\mathbb{I}}-{\sum }_{n}{P}_{t}^{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 ${{ \mathcal Y }}_{m}^{n}(t)=0$ 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 ${{ \mathcal L }}_{I}$, we partition the sites onto A and B sublattices such that

$$\begin{eqnarray}{\dot{\rho }}_{A}(t) & = & {{ \mathcal L }}_{A}{\rho }_{A}(t)-{\rm{i}}\displaystyle \sum _{\alpha }{J}_{\alpha }\mathrm{Tr}\{{\sigma }^{\alpha }{\rho }_{B}(t)\}[{\sigma }^{\alpha },{\rho }_{A}(t)]\\ & & -\displaystyle \sum _{\alpha \beta }\displaystyle \frac{{J}_{\alpha }{J}_{\beta }}{z}\left[{\sigma }^{\alpha },{\displaystyle \int }_{0}^{t}{\rm{d}}{t}^{\prime} ({d}_{\alpha \beta }(t,t^{\prime} ){{\rm{e}}}^{{{ \mathcal L }}_{A}(t-t^{\prime} )}\delta {\sigma }_{A}^{\beta }(t^{\prime} ){\rho }_{A}(t^{\prime} )\right.\\ & & -{s}_{\alpha \beta }(t,t^{\prime} ){{\rm{e}}}^{{{ \mathcal L }}_{A}(t-t^{\prime} )}{\rho }_{A}(t^{\prime} )\delta {\sigma }_{A}^{\beta }(t^{\prime} ))],\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&{d}_{\alpha \beta }(t,t^{\prime} )=\mathrm{Tr}\{{\sigma }^{\alpha }{{\rm{e}}}^{{{ \mathcal L }}_{B}(t-t^{\prime} )}\delta {\sigma }_{B}^{\beta }(t^{\prime} ){\rho }_{B}(t^{\prime} )\},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{s}_{\alpha \beta }(t,t^{\prime} )=\mathrm{Tr}\{{\sigma }^{\alpha }{{\rm{e}}}^{{{ \mathcal L }}_{B}(t-t^{\prime} )}{\rho }_{B}(t^{\prime} )\delta {\sigma }_{B}^{\beta }(t^{\prime} )\},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\delta {\sigma }_{A}^{\beta }(t)={\sigma }^{\beta }-\mathrm{Tr}\{{\sigma }^{\beta }{\rho }_{A}(t)\}{{\mathbb{I}}}_{2},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\delta {\sigma }_{B}^{\beta }(t)={\sigma }^{\beta }-\mathrm{Tr}\{{\sigma }^{\beta }{\rho }_{B}(t)\}{{\mathbb{I}}}_{2}\end{eqnarray} \tag{ 14 }$$

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.

4.2. 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 ${ \mathcal C }$ which are identical due to the translational symmetry of the lattice

$$\begin{eqnarray}&&R(t)\approx \mathop{\displaystyle \bigotimes }\limits_{i}{\rho }_{{ \mathcal C }}.\end{eqnarray} \tag{ 15 }$$

This density matrix then evolves according to the decoupled cluster mean-field Liouvillian which can be written as

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{CMF}}={{ \mathcal L }}_{{ \mathcal C }}+{{ \mathcal L }}_{{ \mathcal B }({ \mathcal C })},\end{eqnarray} \tag{ 16 }$$

where ${{ \mathcal L }}_{{ \mathcal C }}={\sum }_{j\in { \mathcal C }}{{ \mathcal L }}_{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 ${{ \mathcal L }}_{{ \mathcal B }}({ \mathcal C })$ which acts only on sites at the boundary of the cluster. For the driven-dissipative Heisenberg model, the boundary Liouvillian is given by

$$\begin{eqnarray}&&{{ \mathcal L }}_{{ \mathcal B }}({ \mathcal C })=\displaystyle \sum _{j\in { \mathcal B }({ \mathcal C })}-{\rm{i}}[{{\bf{B}}}_{j}^{\mathrm{eff}}(t)\cdot {{\boldsymbol{\sigma }}}_{j},{\rho }_{{ \mathcal C }}],\end{eqnarray} \tag{ 17 }$$

where ${{\bf{B}}}_{j}^{\mathrm{eff}}(t)$ is related to the average polarisation of the sites adjacent to the boundary ${ \mathcal B }({ \mathcal C })$ 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 ${{ \mathcal C }}_{A}$ and ${{ \mathcal C }}_{B}$ where the four sites in the corners and the central site of ${{ \mathcal C }}_{A}$ are assigned to the A sublattice and ${{ \mathcal C }}_{B}$ is the complement of ${{ \mathcal C }}_{A}$.

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

$$\begin{eqnarray}&&R(0)=\mathop{\displaystyle \bigotimes }\limits_{i\in A}\displaystyle \frac{1}{2}{\left({{\mathbb{I}}}_{2}+\displaystyle \sum _{\alpha }{n}_{A}^{\alpha }{\sigma }^{\alpha }\right)}_{i}\ \mathop{\displaystyle \bigotimes }\limits_{j\in B}\displaystyle \frac{1}{2}{\left({{\mathbb{I}}}_{2}+\displaystyle \sum _{\alpha }{n}_{B}^{\alpha }{\sigma }^{\alpha }\right)}_{j},\end{eqnarray} \tag{ 18 }$$

where ${n}_{A}^{\alpha }$ and ${n}_{B}^{\alpha }$ are the components of two vectors within the unit sphere. For this initial state, we time-integrate 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 ≫ γ⁻¹, 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 time-independent steady state. We performed the simulations for two sets of parameters: $\{{J}_{x},{J}_{y},{J}_{z},{\rm{\Omega }}\}$ $=\,\{-7.0,6.0,2.0,1.0\}$ and $\{{J}_{x},{J}_{y},{J}_{z},{\rm{\Omega }}\}=\{-6.4,3.0,6.0,2.25\}$ 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.

5.1. 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 $\mathrm{Tr}\{{\sigma }^{z}{\rho }_{A}\}$ and $\mathrm{Tr}\{{\sigma }^{z}{\rho }_{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.

Zoom In Zoom Out Reset image size

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 ${n}_{A}^{\alpha }$ and ${n}_{B}^{\alpha }$. In the limit $z\to \infty$, 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 ${{ \mathcal Y }}_{m}^{n}(t)$ 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}^{* }\gt 10$ in both cases. In contrast to $\{{J}_{x},{J}_{y},{J}_{z},{\rm{\Omega }}\}=\{-7.0,6.0,2.0,1.0\}$, for $\{{J}_{x},{J}_{y},{J}_{z},{\rm{\Omega }}\}=\{-6.4,3.0,6.0,2.25\}$ all initial states enter limit cycles for $z\gt {z}^{* }$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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\lt {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 parameters considered here, the truncated terms in the Dyson series given in equation (9) scale as ${{ \mathcal Y }}_{m}^{n}(t)\sim {({J}_{\alpha }/\gamma )}^{m}$ and therefore convergence is not guaranteed as $| {J}_{\alpha }/\gamma | \geqslant 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.

5.2. 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 two-dimensional 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) for $\{{J}_{x},{J}_{y},{J}_{z}\}=\{-7.0,6.0,2.0,1.0\}$. 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For $\{{J}_{x},{J}_{y},{J}_{z},{\rm{\Omega }}\}=\{-6.4,3.0,6.0,2.25\}$, 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.

Zoom In Zoom Out Reset image size

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}^{* }\gtrsim 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.

ETO and MJH were supported by the EPSRC via EP/N009428/1. ETO acknowledges funding from the SFC PECRE scheme coordinated by SUPA. JJ acknowledges support from the National Natural Science Foundation of China No. 11605022 and No. 11547119, Natural Science Foundation of Liaoning Province No. 2015020110, and the Xinghai Scholar Cultivation Plan and the Fundamental Research Funds for the Central Universities. RF acknowledges EU-QUIC and MIUR-QUANTRA.