Quantum Synchronisation Enabled by Dynamical Symmetries and Dissipation

In nature, instances of synchronisation abound across a diverse range of environments. In the quantum regime, however, synchronisation is typically observed by identifying an appropriate parameter regime in a specific system. In this work we show that this need not be the case, identifying conditions which, when satisfied, guarantee that the individual constituents of a generic open quantum system will undergo completely synchronous limit cycles which are, to first order, robust to symmetry-breaking perturbations. We then describe how these conditions can be satisfied by the interplay between several elements: interactions, local dephasing and the presence of a strong dynamical symmetry - an operator which guarantees long-time non-stationary dynamics. These elements cause the formation of entanglement and off-diagonal long-range order which drive the synchronised response of the system. To illustrate these ideas we present two central examples: a chain of quadratically dephased spin-1s and the many-body charge-dephased Hubbard model. In both cases perfect phase-locking occurs throughout the system, regardless of the specific microscopic parameters or initial states. Furthermore, when these systems are perturbed, their non-linear responses elicit long-lived signatures of both phase and frequency-locking.

In nature, instances of synchronisation abound across a diverse range of environments. In the quantum regime, however, synchronisation is typically observed by identifying an appropriate parameter regime in a specific system. In this work we show that this need not be the case, identifying conditions which, when satisfied, guarantee that the individual constituents of a generic open quantum system will undergo completely synchronous limit cycles which are, to first order, robust to symmetry-breaking perturbations. We then describe how these conditions can be satisfied by the interplay between several elements: interactions, local dephasing and the presence of a strong dynamical symmetry -an operator which guarantees long-time non-stationary dynamics. These elements cause the formation of entanglement and off-diagonal long-range order which drive the synchronised response of the system. To illustrate these ideas we present two central examples: a chain of quadratically dephased spin-1s and the many-body charge-dephased Hubbard model. In both cases perfect phase-locking occurs throughout the system, regardless of the specific microscopic parameters or initial states. Furthermore, when these systems are perturbed, their non-linear responses elicit long-lived signatures of both phase and frequency-locking.

I. INTRODUCTION
Synchronization is a fascinating and multi-disciplinary topic in modern science, focussed on understanding how a collection of individual bodies adjust their natural rhythms and phases through their interactions with each other and the environment [1][2][3][4][5]. In a striking display of cooperative behaviour, this adjustment can lead to a variety of phenomena such as the 'winking' of fireflies, the behavioural synchrony of groups of strangers or the coupling of a pair of pendulums through a mutual support [6][7][8].
The study of synchronisation in quantum systems has attracted significant attention [9][10][11][12][13][14][15]. In this regime, synchronisation takes on a fairly broad definition due to the variety of cooperative, entangled behaviour that can occur [16]. The formation of a Bose-Einstein condensate (BEC), for example, could be considered perfect synchronisation [17] due to the collective condensation of the atoms in the bosonic gas. In closer analogy to classical systems, models of self-sustained quantum oscillators, such as quantum Van der Pol oscillators [13,18,19] or pairs of micromasers [20], have been shown to lock phases and reach coupled limit cycles. Quantum effects play a decisive role in either enhancing [21,22] or hindering [23] this synchronicity. Under the mean-field approximation, these results have been extended to larger systems of oscillators where the underlying mechanism for synchronisation is a reduction in the uncertainty in the phase distribution at the expense of the certainty in the number distribution [18].
Recently, there has been a focus on observing synchronisation in the limit cycles of quantum systems which have no classical analogue [9][10][11]. The qutrit has been proposed as a logical candidate for this and recent work has demonstrated that it can be entrained to an external signal [9,11] or phase-locked and entangled with a second spin [10]. In these single or two qutrit systems, synchronisation emerges due to careful control over the Hamiltonian and dissipation parameters and is witnessed through both the phase space portrait and entanglement profile of the spins.
One of the most remarkable features of synchronisation in the classical regime, however, is that it occurs in such a diverse range of systems -with completely different sizes, structures and microscopic parameters [6][7][8]24]. This diversity, in turn, leads to a rich variety of observable, complex behaviour. Hence, instead of identifying specific quantum systems and regions of parameter space where a synchronised response can be observed, we consider it pertinent to take a different route and determine, more generally, conditions which will ensure synchronisation in a quantum system.
In this work we adopt this approach, identifying these conditions and uncovering a novel mechanism which guarantees synchronisation in a generic open quantum system, independent of its microscopic details. We show how these conditions can be satisfied via the interplay between several elements: interactions, local dephasing and the existence of a strong dynamical symmetry (an operator which guarantees non-stationary dynamics in the long-time limit of the system [25]). The coaction of these elements underpins the formation of a structure to the long-time density matrix which ensures limit cycles describing entangled, cooperative behaviour. These cycles capture the essence of quantum synchronisation, describing oscillations where the constituents of the system are locked to a common phase and frequency whilst also featuring the off-diagonal long-range order present in states such as BECs and superconductors [26,27]. Furthermore, we prove that this mechanism for synchronisation is, to first order, completely robust to the presence of symmetry-breaking perturbations.
We then present several physical examples, which have no classical analogue, where this phenomenon arises -a chain of interacting spin 1s and the many-body charge-dephased Hubbard model. These systems exhibit perfect distance-invariant phase synchronisation for a wide range of parameters and initial states. Moreover, in these examples, we are able to identify analytical expressions for the long-time density matrix -which is typically an unfeasible task in strongly-correlated many-body systems. Finally, we peturb these systems away from the dynamical symmetry regime where the non-linear response facilitates the observation of strong, exceptionally long-lived signatures of both phase and frequency locking.
Quantum synchronisation is sometimes viewed in terms of a locking in phase space of self-sustained oscillators, measured through the Husimi-Q or Wigner phase space distributions [9,10,19]. In the work and examples in this manuscript we have, instead, opted to focus on the explicit limit cycles of the bodies in the system and plot them alongside the various quantum synchronisation measures we use (such as the entanglement or off-diagonal coherences). This is in order to reflect our intuition of quantum synchronisation as an intrinsically rhythmic process underpinned by quantum properties not available in classical systems. Our definition exposes the presence of quantum properties in the synchronised states which are less easy to see in, for example, the Husimi-Q distribution.

A. Strong Dynamical Symmetries
Firstly, we introduce the concept of a strong dynamical symmetry by providing a brief summary of the work in [25], we restrict ourselves to Markovian dynamics for simplicity. Consider the time evolution of the density matrix of an open quantum system via the Lindblad equation (here, and in the remainder of this work, we set = 1) where H is the Hamiltonian of the system and {L j } are a set of 'jump' operators which model the interaction between the system and the environment with associated coupling strengths {γ j }. The jump operators are used to form the dissipator D[ρ] which competes with the coherent evolution due to the Hamiltonian H. We denote the Liouvillian superoperator with L and the steady state(s) as ρ ss , which satisfy Lρ ss = 0. If we can identify an operator A which satisfies then we say that the system posseses a 'strong dynamical symmetry'. The relation [H, A] = ωA describes the presence of a dynamical symmetry operator. We then refer to this as a strong dynamical symmetry operator because it commutes with all the jump operators and their conjugates [28]. It is then straightforward to prove from these relations that there exists a series of eigenmodes of L of the form where the corresponding imaginary eigenvalues indicate the presence of non-stationary dynamics in the long-time limit of the system. The operator A acts as a raising/lowering operator, generating a ladder of equidistant mixed states within the kernel of the Liouvillian. These results extend beyond that of a decoherence free subspace [29,30] as the imaginary modes are, in general, mixed and cannot be written as a convex superposition of pure states |φ φ| which are immune to the dissipation L j |φ = 0 ∀j.

B. Quantum Synchronisation via Dynamical Symmetries, Interactions and Dephasing
We now show how a generic open quantum system can provide a natural environment for observing quantum synchronisation. Consider an open quantum system where the Hilbert space is constructed from a series of N identical, local spaces or 'bodies' H = ⊗ j H j .
We now imagine the system has a strong dynamical symmetry operator satisfying Eq. (2) and assume the imaginary modes from Eq. (3) form a complete basis for the long-time density matrix of the system. Conseqently, we can write this state as lim t→∞ ρ(t) = ρ ∞ (t) = n,m; n≥m C n,m e iω(m−n)t ρ n,m + h.c. , (4) where the C n,m are a set of real coefficients associated with the overlap between the initial state and the ρ n,m . Now, consider the expectation value of some M -point observable X = j∈B X j , where B = {a, b, c, ...} is a set of M local spaces containing no duplicates and X j is some hermitian local operator acting on site j. It follows from Eq. (4) that lim t→∞ X (t) = n,m; n>m D n,m cos (ω(m − n)t) + const., (5) with D n,m = 2Tr(Xρ n,m )C n,m . Provided that D n,m = 0 for at least one n, m where |n − m| = 0 then Eq. (5) describes coherent, non-decaying limit cycles in the associated observable. The equidistance of the imaginary eigenspectrum is crucial and ensures the frequencies involved are commensurate and do not destructively interfere with each other. These limit cycles, along with the well-defined, coherent phase-evolution described in Eq. (4) are some of the hallmark features of temporal synchronisation. Importantly, to have full synchronisation we need each of the bodies to undergo the same coherent phase evolution. Whilst the existence of a strong dynamical symmetry ensures non-stationarity, it does not mean that the bodies will lock together in phase space and undergo identical limit cycles. The fundamental requirements for this to happen are that the steady state and strong dynamical symmetry operator(s) are translationally invariant and, as assumed earlier, form a complete basis for the long-time density matrix ρ ∞ (t) of the system. If these requirements are met then, through Eq. (3), the ρ n,m and thus ρ ∞ (t) will inherit the translational symmetry of these operators and the limit cycles described in Eq. (5) will be independent of the specific bodies in the set B = {a, b, c, ...} (only the cardinality of the set matters). Consequently, the system will be perfectly synchronised as all the bodies in the system will be locked to the same frequency and phase -independent of the specific value of any microscopic parameters.
Recent work has shown that the interplay between interactions and local, homogeneous dephasing in an open quantum system can 'wash' out any geometry associated with the system: creating steady states with off-diagonal long-range order and ensuring they, along with any strong dynamical symmetry operators, are completely translationally symmetric (see [25,31]). These states and operators form a complete basis for the longtime density matrix of the system. Hence, following the discussion in the previous paragraph, we identify interactions between the bodies in our system as well as local, homogeneous dephasing [? ] as elements which, when combined with the existence of a strong dynamical symmetry, can ensure the system reaches a completely quantum synchronised state: i.e. with locked limit cycles underpinned by intrinsically quantum properties such as entanglement and off-diagonal long-range order. We will illuminate these ideas with a pair of examples in Sec. 3 and explicitly show these fully synchronised cycles alongside their intrinsically quantum behaviour.
We anticipate that for synchronisation to occur in our framework via homogeneous, local dephasing the local Hilbert space dimension should satisfy Dim(H j ) > 2. This is because there must be local coherences available in the long-time limit where a valid phase relationship can be established and the system can undergo long-time oscillations. Any non-trivial local dephasing in an array of 2-level systems will destroy the available coherences and prevent the qubits from undergoing a valid limitcycle, a pre-requisite for synchronisation. This argument does not apply to arrays of 2-level systems under more general non-local dissipation.
In this work we consider synchronisation under the Markov approximation and so are limited to weakly interacting systems. However, in our examples synchronisation occurs for any finite interaction strength -its amplitude only sets the timescale on which a synchronised state is reached. Hence, even if the interaction strength is small, the system will eventually reach a synchronised state (in experimental setups care should also be taken that the timescale on which synchronisation occurs is shorter than the coherence time of the system). Moreover, we emphasize that highly controllable quantum systems such as lattices of ultracold atoms immersed in a Bose-Einstein Condensate [25,31], can be engineered to accurately implement dynamics described by local master equations.

C. Perturbations away from the Dynamical Symmetry Regime
We now show how the synchronisation discussed in the previous section is robust to perturbations away from the dynamical symmetry regime. Typically, we can imagine that the relation [H, A] = ωA from Eq. (2) arises due to some homogeneous field F = ω j f j in the Hamiltonian for which A is a raising/ lowering operator. The synchronisation described in the previous section is then a consequence of perfect phase locking of the individual constituents of the systems, frequency locking will occur because the individual bodies j share the same natural frequency ω.
If the field is not homogeneous, i.e. F = j ω j f j , then the Hamiltonian H can be split into two terms. The first contains a homogeneous term jω j f j and all other nonfield terms, the second contains only the inhomogeneous part j δ j f j . We have parametrised ω j as ω j =ω j + δ j whereω j is the average of the set {ω j }. We then, correspondingly, split the Liouvillian into two parts and scale by 1/ω j : where =δ j /ω j is the, small, perturbation parameter andδ j is the average over the set of detunings {δ j }. We then assume that the eigenvectors and eigenvalues of this new Liouvillian are perturbations on those for L (0) It then follows (see Supplemental Material, SM) that the first order eigenvalue shift is λ (1) = Tr ρ (0) † L (1) ρ (0) , which is purely imaginary as it can be rearranged to be the trace of a skew-hermitian matrix. As a result, to first order, the eigenvalues of the eigenmodes in Eq. (3) remain imaginary and thus there is no decay in the system's long-time dynamics when perturbed away from the dynamical symmetry regime. Moreover, if the imaginary eigenmodes are unchanged under a swap between two bodies j and l then we have that λ (1) = 0 as Tr ρ (0) † f j ρ (0) is independent of j.
As discussed in Sec. II B this symmetry in the imaginary modes is seen for local, translationally-invariant dephasing in an interacting system [25,31]. Hence, the system will undergo completely non-linear response to peturbations away from the dynamical symmetry regime. To first order, the solutions in Eq. (3) are still eigenmodes and so will decay with a rate that scales at least quadratically with the perturbation parameter . Hence they correspond to 'slow' modes which will, in general, decay much slower than rest of the eigenmodes of the Liouvillian. We can therefore expect to transiently observe the corresponding synchronised features these eigenmodes possess, with a lifetime that scales at least quadratically with the perturbation parameter . The individual constituents will be locked in both phase and frequency, despite having different natural frequencies.
This appearance of synchronisation due to the formation of 'slow' decay modes in the Liouvillian which lock the system to specific frequencies is consistent with the mechanism for transient synchronisation discussed in [16] and observed in [32].

III. EXAMPLES
We have shown how, in a generic interacting open quantum system, the combination of interactions, local dephasing and a strong dynamical symmetry can underpin a coherent, distance-invariant, synchronised structure to the long-time density matrix. Furthermore, the system is robust to perturbations away from the dynamical symmetry regime. In order to elucidate these results we present a pair of paradigmatic examples where they can be observed.

A. Synchronisation in a chain of Spin-1s
For our first example, we take a system formed from a series of spin-1s or qutrits. The local basis for each spin-1 is spanned by the three states {|↓ , |0 , |↑ }. The key operators are S + j , S − j and S z j which are, respectively, the spin-1 raising, lowering and magnetisation operators for spin j. The x and y components of the spin-1 operator can be formed from the raising and lowering operators: . By dropping the local subscript we denote the total of an operator, e.g. S z = j S z j . We take a spin-1 anisotropic Heisenberg model in a chain geometry [33] (see Fig. 1) where the spins each have natural frequency ω j and nearest-neighbour coupling strengths J and ∆ = 0. The system is then immersed in a bath which induces local, quadratic dephasing in a spin-agnostic manner. The ensuing dynamics is modelled via the equation which is the master equation in Eq. (1) with jump operators L j = (S z j ) 2 applied to each spin at a rate γ ∀j. Initially, in order to derive an analytical solution to the long-time dynamics, we focus on the 'frequency-matched' case ω j = ω, ∀j. In the SM we prove that the 2N + 2 steady states alway take the form where |m i is an eigenvector of S z with eigenvalue m: S z |m i = m |m i , and i indexes the possible eigenvectors for each m. We have also defined |−m i = SF |m i , where SF = ⊗ N j=1 |↑ ↓|+|↓ ↑|+|0 0| is the spin-flip operator. For example if |2 1 = |0 ↑↑ then |−2 1 = |0 ↓↓ , or if |0 2 = |0 ↑↓ then |0 2 = |0 ↓↑ . In order for Tr(ρ ss ) = 1 the elements {λ m } and λ 0 must satisfy the equation where the terms in the second summation are skipped if (N − s + m)/2 is not an integer. The long-time dynamics of L is not, however, solely governed by this steady state. We identify (see SM) multiple strong dynamical symmetry operators of the form Following this we can determine, see Eq. (3), the imaginary eigenmodes of the Liouvillian through the action of these operators on the steady state. Explicitly, we have, which is a non-trivial result as the steady-state ρ ss is inherently singular. Further application of A m is redundant as A m A m = 0 and thus, each A m generates a unique imaginary eigenmode ρ m 1,0 via left-multiplication of the steady-state. Crucially, however, the eigenspectrum is still equidistant as the eigenvalues of the different modes form a ladder with a spacing of 2ω. Hence, the structure of the long-time eigenspace is analogous to a system with a single strong dynamical symmetry operator.
The steady state and imaginary eigenmodes in Eqs. (30) and (12) form a complete basis for the long-time dynamics of Eq. (27) and so, similarly to Eq. (4), the density matrix can be expressed as a superposition of these modes in the limit t → ∞. The imaginary modes describe coherences between sectors of opposite magnetisation, their excitement will ensure the system reaches a limit cycle in the long-time limit. Moreover, the density matrix is completely translationally invariant; as described in Sec. 2B the dephasing and interactions have washed out any geometry in the system, which now has no characteristic length-scale. This invariance can be seen in the off-diagonal coherences described in Eq. (3), which occur at all length-scales of the chain and are completely uniform with respect to distance.
In Fig. 2 we visualise the analytical results in Eqs. (30), (12) and (14). We present a plot of the eigenspec-trum of L [ Fig. 2(a)], the formation of these imaginary eigenmodes is clear and their spacing is set by the value of ω. We also show the structure of the density matrix in the long-time limit of Eq. (27) [Fig. 2(b)]. The system is in a superposition of the steady state in Eq. (30) and the imaginary modes in Eq. (12), hence it only has elements along the diagonal and anti-diagonal in the configuration basis. The magnitude of each of these matrix elements is constant in time. The phase of the elements along the anti-diagonal is well-defined and evolves in time at a frequency f = 2mω for the corresponding matrix element We can explicitly prove that this coherent density matrix structure leads to observable synchronisation in the long-time limit of the system. Specifically, consider the operator ..} is a set of M sites containing no duplicates. The operator is formed from quadratic, local operators which measure fluctuations in the magnetisation. The quadratic form is necessary in order to be able to measure the coherences between the basis states |↓ and |↑ -any operator formed solely from linear local operators will relax to stationarity. We prove (see SM) that in the limit t → ∞ where the real coefficients D m are those associated with the overlap between the initial state ρ(0) and either the imaginary eigenmodes or the steady state. The trace overlap Tr(ρ m 1,0 X) only depends on the cardinality of B, not the specific sites within the set; a direct consequence of the complete translational invariance of the modes spanning the kernel. Consequently X is also independent of the specific choice of sites over which we measure the correlator X, only the number of sites matters. Thus we see that the long-time dynamics will perfectly synchronise the spin-1s to clean, coherent limit-cyclesregardless of the specific values of the initial state or the Liouvillian parameters.
Moreover, the modes in Eq. (12) are entangled and cannot be written as a superposition of separable states -which we explicitly show in the following numerics. We also demonstrate that the reduced correlator X − j∈B X j is distance-invariant and non-zero. Hence, we consider the synchronisation observed in Eq. (15) to be inherently quantum -underpinned by long-range correlations, which are a result of entanglement between the bodies in the system.
In Eq. (28), the S z i S z i+1 (ZZ) term is an interaction term and hence the parameter ∆ sets the strength of the interactions[? ] and plays a critical role in the formation of synchronisation. Specifically, the ZZ interaction ensures that only translationally invariant strong dynamical symmetries and steady states are present. Therefore, the asymptotic time-dependent density matrix also possesses this symmetry which, in turn, implies perfect synchronization (see Section 2B). Provided ∆ = 0 this will always be case, with the explicit value of ∆ only effecting the time-scale on which synchronisation occurs. When ∆ = 0 there are additional steady states and strong dynamical symmetries of the Liouvillian (see SM for an example) which are not translationally invariant and interfere with the symmetry of the known solutions described earlier, disrupting the synchronicity of the system.
All of these results are valid for any arbitrary length chain of N spin-1s under the Liouvillian in Eq. (27). In the subsequent numerics we focus on a small series of spin-1s, i.e. N = 3 or N = 4. This is because as the system sizes increases, for generic initial states (product states for example), the diagonal correlations become increasingly dominant in the long-time limit compared to the off-diagonal coherences (|D 0 | |D m =0 |), reducing the amplitude of the synchronisation measures in the system (this amplitude will, however, remain finite for any finite-size system). Consequently, by focussing on a small chain, we can readily resolve the features of quantum synchronisation and directly witness our analytical calculations by solving the master equation in Eq. (27) through numerical exponentiation of the Liouvillian superoperator L. Later in the text we will present our second example where synchronisation is induced via our mechanism and is observable even in the thermodynamic limit. The combination of these two examples emphasizes how this symmetry-induced synchronisation occurs in systems of varying size and structure.
For the following results, we start in a specified initial state and then time-evolve it under the Liouvillian in Eq. (27) measuring various time-dependent quantities in order to observe the formation of synchronisation. As a first synchronisation measure for the local observables in our model we consider the time-dependent Pearsoncorrelation factor [14,34]. It can be used to measure the correlation over time for two functions f, g defined on a domain [t, t + ∆t] with the function averagef = 1 ∆t t+∆t t f dt. This correlation factor is maximal (minimal), 1 (−1) when the two signals f and g are perfectly synchronised (antisynchronised), and 0 when they display no correlations.
In Fig. 3 we set f = (S x j ) 2 and g = (S x l ) 2 in order to measure the synchronisation over time between two of the spin-1s j and l -we start from a completely random product state. We also include the individual functions (S x j ) 2 over time for each spin. In agreement with Eq.  tions (γ = 0, ∆ = 0). The presence of dissipation causes the system to converge to the expected clean coherent limit cycles [25] but the limit cycles for each spin are out of phase. Despite the fact the system still has a strong dynamical symmetry, the absence of interactions prevents synchronisation as the kernel of the Liouvillian still has a memory of the initial geometry of the system. We also show the closed case when γ = 0 [Figs. 3(c) and (f)], the dynamics are completely chaotic and unsynchronised due to the multitude of incommensurate frequencies in the eigenvalues of the Hamiltonian.

FIG. 4. Dynamics of
Whilst the plots in Fig. 3 demonstrate that the spins are able to perfecly lock phases they do not capture the collective origin of this synchronisation. In this vein, in Fig. 4, we plot the reduced correlator which is non-zero and identical for any choice of spins j and l. The observed oscillations contain several frequencies [ Fig. 4(b)] due to the excitement of multiple imaginary modes in Eq. (12). These imaginary modes all contain coherences between different spins, the synchronisation observed in Fig. 3 is dependent on the existence of these inter-spin coherences and Fig. 4 shows that they give rise to perfect, distance-invariant correlations throughout the system. So far we have considered the 'homogeneous' case ω j = ω ∀j. The spins share the same natural frequency and we have shown how, under dephasing and interactions, their phases will align perfectly. In order to discuss synchronisation in full we now set the frequencies of the spins to be mismatched. In this case, the imaginary modes in Eq. (3) are no longer exact eigenvectors of the Liouvillian and, in the long-time limit, the system will decay to an ensemble which is diagonal in the configuration basis -where no synchronisation can occur. However, as was shown in Eq. (6), for sufficiently small values of , the system is only slightly perturbed from a 'dynamical symmetry' regime defined by the spins having a common frequency which is the average of their natural frequencies. Furthermore due to the translational symmetry of the imaginary eigenmodes, we know the system is, to first order, completely robust to this kind of perturbation. Hence, it is interesting to observe whether the spins are able to synchronise to the dynamical symmetry regime on an intermediate time-scale and, if so, how long the system takes to desynchronise and reach a diagonal ensemble.
As measures to track this, and to highlight the quantum nature of the synchronisation in Figs. 3 and 4, we introduce two common witnesses for quantum synchronisation: the negativity N [35] and off-diagonal coherences with T j indicating the partial transpose with respect to site j and ||X|| = Tr √ X † X denoting the trace norm of an operator. The negativity can be seen as a measure of the degree to which spin j is entangled with the rest of the system whilst the coherence quantifier C describes the total magnitude of the off-diagonal elements in the density matrix. When the frequencies are matched the system is synchronised, and due to the off-diagonal, entangled nature of the modes in Eq. (3) quantities such as these will remain finite indefinitely.
In Fig. 5 we show how these synchronisation witnesses evolve in time when the system is perturbed from the dynamical symetry regime. We use the detuning strength δ to characterise the range of the natural frequencies. The explicit distribution of natural frequencies is not important, the key parameter is its width and in the SM we obtain similar results when the natural frequencies are drawn from a uniform random distribution. Initially, the system is in a product state where N = 0, the transient dynamics then causes the formation of entanglement and anti-diagonal coherences which decay away at a rate set by δ. We show how this entanglement forms [Figs. 5(cd)]: despite having mismatched frequencies and phases the spins lock to an intermediate limit cycle with identical phase and frequency -which is twice the average of the natural frequenciesω j (due to the factor of 2 in Eq. (13)). The life-time of this cycle is large and as δ → 0 diverges to infinity, evidenced by the tongue-like behaviour seen in Figs. 5(a-b). These figures show how the corresponding measures act as strong witnessess to the synchronisation in the system -emphasizing its quantum nature.
The imaginary eigenmodes in Eq.
(3) are translationally-invariant and hence the robust, synchronised behaviour observed [ Fig. 5] is the result of a secondorder response to the detuning. The cross-sections included in Fig. 5b are evidence of this. At a given time the coherences are well-approximated by a gaussian profile (see SM) as a function of the detuning. Meanwhile at a given detuning the coherences decay away exponentially as a function of time, the decay rate d is proportional to the square of the detuning (see SM for numerical evidence of this). Furthermore, we explicitly show this non-linear scaling in Fig. 6. We calculate the shift in the imaginary eigenvalues, |λ − λ (0) | (see Sec. (II C)) from their original value λ (0) at δ = 0 as a function of the detuning δ. There is no noticeable shift to first order in δ -the fitted curve is proportional to δ 2 . Notably, the highest imaginary eigenmode ρ N 1,0 = |↑↑ ... ↓↓ ...| is always unshifted, and remains imaginary regardless of the distribution of natural frequencies.

B. Many-body synchronisation in the Hubbard model
As our second example, we take the 1D N -site Hubbard model [36] in, potentially, disordered magnetic and chemical fields. We focus on 1D lattices for numerical tractability, nonetheless it should be emphasized that these results are solely based on symmetry and thus can be observed in any bi-partite d-dimensional realisation of where c † σ,j and its adjoint are the usual creation and annihilation operators for a fermion of spin σ ∈ {↑, ↓} on site j. Additionally, n σ,j is the number operator for a particle of spin σ on site j and τ , U , ω j and µ j play the role of kinetic, interaction, magnetic and chemical energy scales respectively.
We then couple the system to a bath which induces spin-agnostic dephasing on each site. Hence, the system's time evolution can be described by the Lindblad Equation: The case ω j = ω ∀j of this model was originally studied in Ref. [25], where the existence of a strong dynamical symmetry was shown to ensure long-time non-stationary dynamics in this strongly-correlated system. A possible experimental realisation of the system is also described in this reference. In this section we show how, further to this, perfect synchronisation is induced in the longtime dynamics of this model: alongside the strong dynamical symmetry there is an inter-site coupling and homogeneous local dephasing which ensure the long-time where the coefficients D i are set by the initial state of the system. Due to the inter-site coupling and local dephasing the imaginary modes in Eq. (20) are completely translationally symmetric (see Ref. [25] for the explicit form of the steady state) and thus this observable is depent only on the cardinality M of the set B, not the specific sites within the set.
As a result, even in the presence of disorder in the chemical potential, the system displays perfectly synchronised magnetic oscillations in the long-time limit. This will occur for a wide range of specific parameters and initial states of the system, the only requirements are that the appropriate coefficients, D i , are finite and the hopping amplitude τ (which couples the different sites together) is non-zero. Moreover, in the thermodynamic limit, initial states which have lim N →∞ S x /N = 0 will ensure the long-time synchronised oscillations in S x j or S x jS S x j+1 have a finite amplitude [25]. Similarly to the previous example these oscillations are underpinned by long-range correlations in the system which arise due to the entangled nature of the long-time density matrix. We demonstrate these features in the following numerics, showing synchronisation in a fully many-body quantum system.
In order to increase the system size accessible to our numerical calculations, we have used a 'quantum trajectories' approach [37] to perform a stochastic unravelling of Eq. (19) and simulate the dynamics at the level of an ensemble of pure wavefunctions. Furthermore in Fig. 7, as the simulation is only on a short time-scale, we were able to use the time-evolving block decimation [38] algorithm on a Matrix Product State [39] decomposition of the trajectory wavefunctions, further increasing the available system size. These simulations were performed with the aid of the Tensor Network Theory library [40].
In Fig. 7 we demonstrate the synchronicity which results from the eigenmodes in Eq. (20). We initialise the system in a product state and, after quenching under the master Equation in Eq. (19) observe how the x-magnetisation on each site synchronises perfectly, oscillating at the anticipated frequency. The Pearson coefficient for the magnetisation on any two sites saturates to 1 in the long-time limit [ Fig. 7(b)], with the dip at tτ ≈ 2 being a transient effect which occurs at the first turning point in the magnetisation.
We now perturb the system from the dynamical symmetry regime by setting the natural frequencies ω j to be inhomogeneous, here we draw them from an evenly spaced distribution, i.e. ω j+1 −ω j = const and ω N −ω 1 = 2δ. Again, as for the spin-1 case, we choose this distribution for simplicity, our observations are independent of the explicit distribution -the key parameter is its width δ. We initialise the system in a specified state and time-evolve under the Liouvillian in Eq. (19). In Fig. 8 we show how, similarly to the previous spin-1 example, the system is still attracted to the synchronised state, in both phase and frequency, on an intermediate time-scale. There is a significant band of detunings where the system stays in this long-lived synchronisation phase (Figs. 8a and b) and the spin on each site locks to the same phase and frequency (which is set by the average frequency of the individual sites, see Fig. 8d). Remarkably, this harmonized response is occuring even in the presence of both magnetic and chemical disorder -emphasizing the robustness of a symmetry-based approach to observing quantum synchronisation. The imaginary eigenmodes in Eq. (20) are translationally-invariant and hence λ (1) evaluates to 0 (see Sec. II C). As with the previous example, this robust, synchronised behaviour is a result of a second-order response to the detuning.

IV. CONCLUSION
We have provided condtions which, when satisfied, guarantee synchronisation in a generic qopen quatnum system. We have then shown how, using a combination of analytics and numerics, the interplay between interactions, local dephasing and a strong dynamical symmetry can satisfy these conditions and facilitate the combination of entanglement and perfect phase synchronisation between the individual constituents of the system. This is a direct result of the formation, in the long-time limit, of a well-defined phase relationship in the off-diagonal coherences of the density matrix, at all length-scales of the system. Furthermore, when perturbed from the dynamical symmetry regime these systems exhibit a secondorder response which results in both phase and frequency locking throughout the system. These observations orginate at the level of the symmetries of the system. Thus, we believe, this work marks an important step in understanding how fully-quantum synchronisation can originate in a wide range of generic physical systems -as opposed to in a single delicately controlled setup. We anticipate further examples of complex quantum networks where symmetry can guide the individual nodes into an entangled, fully synchronised state.
We highlight the potential role such a harmonised response can play in developing quantum technologies such as atomic clocks and other metrological instrumentswhich rely on quantum-enhanced synchronicity and cooperative behaviour in order to outperform their classical counterparts [41,42].
We also consider it pertinent to explore the role of dynamical symmetries in synchronising closed, stronglycorrelated systems -where dissipation is absent. Recently, it was shown how the presence of a quasi-local dynamical symmetry can prevent stationarity and guarantee oscillatory dynamics in the XXZ model [43].
Finally we note that strong dynamical symmetries are also defined outside of this approximation [25] and so we anticipate our results to be extendible beyond the Markov regime. This will be particularly important for understanding whether our results can be observed for non-local dissipation: strongly interacting subsystems are often not described by local master equations [44].

V. ACKNOWLEGDMENTS
We would like to thank J. Mur-Petit and J. Coulthard for useful discussions.
This work has been supported by EPSRC grants No. . The Hamiltonian contains, amongst other terms, an inhomogeneous field j ω j f j where f j is some local field operator. We can split the Hamiltonian into two terms, a homogeneous part: jω j f j and an inhomogeneous part: j δ j f j where ω j =ω j + δ j andω j is the average natural frequency. We then, correspondingly, split the Liouvillian into a perturbed and an unperturbed part, scaling by 1/ω j : withδ j is the average over the set of detunings {δ j }. We work in superket and superbra form and assume that the unperturbed Liouvillian L (0) contains a series of imaginary eigenvectors and values L (0) |ρ For small << 1 we expand the eigenvectors and values of the new Liouvillian L as a peturbative power series on the previous, i.e.
We also know that the orthonormality condition σ i |ρ j = Tr(σ † i ρ j ) = δ i,j must hold ∀ -where we have defined σ i and ρ j as the matrix forms of the corresponding superket and superbras. Using this condition, to 0th and 1st order, we have We now simplify the known expression σ i |L|ρ i = λ i σ i |ρ i by subsitituting the expansions in Eq. (24) and converting to matrix form. As a result we find the first order correction to the imaginary eigenvalue For imaginary eigenmodes formed from a strong dynamical symmetry [25] it can be proved that the left and right eigenmodes are the same, i.e. σ (0) If the mode ρ (0) i is translationally invariant (i.e. it is unchanged under a permutation of any pair of sites) then we notice that the trace in Eq. (26) is independent of j. Using the fact j δ j = 0 we then have λ (1) i = 0. Hence, for translationally-invariant imaginary eigenmodes formed from a strong dynamical symmetry we find that the system exhibits a non-linear response to perturbations in the homogeneity of the natural frequencies. This underpins the strong-synchronised response of the two systems considered in the main text.

VII. IMAGINARY MODES AND STEADY STATES OF A SPIN 1 CHAIN
Here we prove the existence of certain imaginary modes and steady states of a dephased XXZ spin-1 chain of length N . In the main text, we consider the Lindblad equation with the Hamiltonian H We start by proving that any state ρ = Gm i=1 |m i m i |, where |m i is one of the G m eigenvectors satisfying S z |m i = m |m i , is a steady state: Lρ = 0. Firstly we substitute ρ ss into Eq. (27) where it is easy to show that This can be done by considering the two spin-1s on the j and j + 1 positions for a given |m i . Then, we have that S + j S − j+1 |m i is a non-zero vector only if the two spin-1s are in one of the configurations |0 ↑ , |00 , |↓↑ , |↓ 0 . Because for these configurations swapping spins j + 1 and j doesn't change the magnetisation m, we can always find the term |m i m i | in the steady state where |m i is just |m i with spins j and j + 1 swapped. Equation (29) then follows from the fact for every term we can find a corresponding term to cancel it with.
We can also show that L Firstly, it is clear that which is 2N + 2 fold degenerate. The coefficients {λ m } and λ 0 must satisfy in order for Tr(ρ ss ) = 1. The terms in the second summation are skipped if (N − s + m)/2 is not an integer. Furthermore, the imaginary eigenmodes ρ m .., N , which we proved satisfy Lρ m 1,0 = 2inωρ m 1,0 , originate as a series of strong dynamical symmetries [25] of the model because: I.e. for this system the strong dynamical symmetry operators are the imaginary modes because they return themselves upon application to the steady state ρ m 1,0 ρ ss ∝ ρ m 1,0 (the steady state is singular). Further application of the strong dynamical symmetry operators is redundant as ρ m 1,0 ρ m 1,0 = 0. For the case when N > 2 and ∆ = 0 numerical calculations show that these 2N +2 steady states and 2N imaginary modes completely span the kernel of L and thus form a complete description of the system's dynamics in the limit t → ∞.
Noninteracting Solutions -The steady states and strong dynamical symmetries described above provide a complete basis for the long-time density matrix of the system when ∆ = 0. When ∆ = 0 these solutions are still valid (the derivations of the previous section are true ∀∆), however there exist additional solutions which are not translationally invariant and therefore disrupt the synchronicity of the system. This is because the system is no longer interacting due to the absence of the S z i S z i+1 term (in the z-basis the S + i S − i+1 terms only describe the exchange of excitations through the lattice and do not constitute interaction terms).
We now provide an example of one of these translationally invariant solutions. Specifically one can define the operator and easily prove that where H(∆) is the Hamiltonian in Eq. (28) as a function of ∆. The first 2 relations in Eq. (34) are because the interaction term in the Hamiltonian j S z i S z i+1 does not commute with B whilst the hopping term does. Hence B is only a valid dynamical symmetry operator when ∆ = 0. The operator B is clearly not translationally invariant (even and odd sites are distinct) and so this interferes with the perfect translational invariance of the solutions derived earlier. Moreover, there are additional solutions which break the translational invariance further. When these solutions are excited by the initial state then the synchronicity of the system is disrupted (see Fig. 3 of the main text).

VIII. LONG-TIME DYNAMICS OF THE SPIN 1 CHAIN
As the imaginary modes ρ m 1,0 contain coherences between the states |↑ and |↓ then they will only affect the dynamics of quadratic observables such as (S x j ) 2 and (S y j ) 2 . We can always write the long-time density matrix as where C m are a series of real coefficients (to ensure hermicity) associated with the overlap between the initial state ρ(t = 0) and either the steady state ρ ss or the imaginary modes ρ m 1,0 . We also have (ρ m 1,0 ) † = ρ −m 1,0 . We consider the expectation value of the operator (S x j ) 2 = (1/4)(S + j + S − j ) 2 . As the imaginary modes for which |m| ≥ 2 must contain at least two flipped spins between the states |m i and −m i | then we immediately have Tr(ρ m 1,0 (S x j ) 2 ) = 0, |m| ≥ 2 ∀j. Hence, we get: where we have used the fact Tr(ρ −1 (S x j ) 2 ) = Tr(ρ 1 (S x j ) 2 ). Equation (36) proves the formation of clean, single frequency oscillations in the associated observable. Furthermore, the modes ρ ss and ρ m 1,0 are all translationally invariant and so the oscillations are identical for all spins: ensuring perfect phase synchronisation.
In order to observe the excitement of higher order modes we must measure higher order correlators. Specifically consider the operator where the set of M sites {a, b, c, ..., } contains no duplicates. Because we now have Tr(ρ m 1,0 X) = 0, m ≤ M and Tr(ρ m 1,0 X) = Tr(ρ −m 1,0 X) then we find and see the appearance of higher order frequencies due to the excitement of higher order imaginary modes. This explains the Fourier Spectrum observed in Fig. 4b) in the main text.

IX. PARAMETRISING THE CROSS-SECTIONS OF THE FREQUENCY-DETUNED SPIN-1 CHAIN
In the main text we considered the response of the system when the magnetic field is inhomogeneous, i.e. the system's dynamics is modelled by Eq. (27) with the Hamiltonian now of the form (39) where the ω j are a series of natural frequencies associated with each spin j. We then considered how, for a given range of natural frequencies, synchronisation witnesses such as the negativity N j (ρ) = ||ρ T j ||−1 2 [35] or off-diagonal coherences C = i =j |ρ ij | [11] evolve in time. For the example in the main text we considered N = 3 with the natural frequencies equally spaced {ω 1 , ω 2 , ω 3 } = {1−δ, 1, 1+δ}J, which produced the maps in Fig. 9(a-b), showing the witnesses as a function of time and detuning.
We now parametrise the cross-sections in Fig. 9(b). In Fig. 9(c) we show how, at a given time, and for sufficiently small detunings, the off-diagonal coherences C are well-described by a Gaussian curve as a function of the detuning. Furthermore, in Fig. 9(d) we calculate the decay coefficient (for the exponential decay of the off-diagonal coherences C versus time) versus δ and show how, for small detunings, d ∝ δ 2 . In Figs. 9(c-d) the tails of the distribution aren't captured by this parametrisation due to numerical precision (both synchronisation quantities are very close to 0 for large detunings). This parametrisation also holds for the cross-sections of the average negativityN j in Fig. 9(b). where [0, δ] is a uniform random distribution over the specified interval. The Synchronisation measures are then averaged over 100 instances of disorder associated with the detuning. Insets) top-left, orange, shows the witness versus time at detuning δ = −0.25. Top-right, white, shows the witness versus detuning at time tJ = 55.0. a) Synchronisation is measured by the average negativity,Nj, for each site. b) Synchronisation is measured as the total magnitude of the off-diagonal coherences C. c) Example dynamics of (S x j ) 2 for the same system except with specific natural frequencies {ω1, ω2, ω3} = {1.258, 1.210, 1.160}J. d) Prevalence of angular frequencies, extracted from the distribution of angular frequencies created using the the time-periods between successive turning points for the oscillations in c) but up to tJ = 100.0. The central line is twice the average of the spin's natural frequenciesωj.

X. FURTHER PLOTS OF THE FREQUENCY-DETUNED SPIN-1 CHAIN
In Figure 5 of the main text we showed how, when the natural frequencies of the spins in the chain are inhomogeneous, the system still locks to a long-lived, synchronised cycle with a frequency which is the average of their natural frequencies. This response emerges as a tongue-like profile in the witnesses C andN j as a function of detuning δ and time. In the main text, for simplicity, we considered the case where {ω 1 , ω 2 , ω 3 } = {1.0 − δ, 1.0, 1.0 + δ}, i.e. the natural frequencies form a uniform sequence. Here, in Fig. 10, we show that this distribution is arbitrary, showing how similar tongues and cross-sections emerge when the natural frequencies are drawn from a uniform random distribution of width δ. The spins are able to lock to an intermediate cycle with a frequency which is twice the average of the natural frequenciesω j .
It is known [25] that the imaginary eigenmodes of this Liouvillian are ρ nm = (A) n ρ ss (A † ) m , Lρ nm = i(m − n)ωρ nm , where ρ ss is a grand-canonical-like state containing the strong-symmetries of the system [25,28,31] and ρ † n,m = ρ m,n . Thus in the long-time limit the state of the system can be written as lim t→∞ ρ(t) = n,m n≥m C n,m e i(m−n)ωt ρ n,m + h.c., (43) where the C n,m are a series of real coefficients associated with the overlap between the initial state and the modes ρ n,m . We calculate the expectation value of the operator S x j lim t→∞ S x j (t) = n,m n≥m C n,m e i(m−n)ωt Tr S x j ρ n,m + h.c..
By expressing S x j in terms of raising and lowering operators and using the fact that a) Tr S x j ρ n,m = Tr S x j ρ m,n and b) the trace vanishes unless |m − n| = 1 (as the operator S x j ρ n,m will have no diagonal elements in the eigenbasis of ρ ss ) we get lim t→∞ S x j (t) = 2 cos(ωt) n=1 C n,n−1 Y n,n−1 , with Y n,n−1 = Tr S x j ρ n,n−1 = Tr S x j ρ n−1,n . Hence, we see persistent oscillations in S x j , which are centred around the x-axis. The modes ρ n,m are completely translationally invariant and thus the spins on each site will synchronise to limit cycles perfectly in phase, regardless of the initial state. We can immediately treat higher order modes through the operator where the set of M sites {a, b, c, ..., } contains no duplicates. We can calculate the expectation value of this operator by expanding it in terms of raising and lowering operators and using the fact the trace of each term is only non-vanishing if the difference between the number of raising and lowering operators is equal to |m − n|.
where the D i 's are a series of coefficients based on the initial state and the various traces between the ρ n,m and products of local spin-raising and lowering operators.