Quantum synchronisation enabled by dynamical symmetries and dissipation

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 ${\hslash }=1$)

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}={ \mathcal L }\rho =-{\rm{i}}[H,\rho ]+\displaystyle \sum _{j}{\gamma }_{j}\left({L}_{j}\rho {L}_{j}^{\dagger }-\displaystyle \frac{1}{2}\{{L}_{j}^{\dagger }{L}_{j},\rho \}\right)=-{\rm{i}}[H,\rho ]+D[\rho ],\end{eqnarray} \tag{ 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 ${ \mathcal L }$ and the steady state(s) as ρ_ss, which satisfy ${ \mathcal L }{\rho }_{{ss}}=0$.

If we can identify an operator A which satisfies

$$\begin{eqnarray}&&[H,A]=\omega A,\quad [{L}_{j},A]=[{L}_{j}^{\dagger },A]=0\quad \forall j,\omega \in {\mathbb{R}}\end{eqnarray} \tag{ 2 }$$

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 ${ \mathcal L }$ of the form

$$\begin{eqnarray}&&{\rho }_{n,m}\propto {\left(A\right)}^{n}{\rho }_{{ss}}{\left({A}^{\dagger }\right)}^{m},\quad { \mathcal L }{\rho }_{n,m}={\rm{i}}\omega (m-n){\rho }_{n,m},\end{eqnarray} \tag{ 3 }$$

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 $| \phi \rangle \langle \phi |$ which are immune to the dissipation ${L}_{j}| \phi \rangle =0\ \forall j$.

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' ${ \mathcal H }={\otimes }_{j}{{ \mathcal H }}_{j}$.

We now imagine the system has a strong dynamical symmetry operator satisfying equation (2) and assume the imaginary modes from equation (3) form a complete basis for the long-time density matrix of the system. Conseqently, we can write this state as

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{t\to \infty }\rho (t)={\rho }_{\infty }(t)=\displaystyle \sum _{n,m;n\geqslant m}\left({C}_{n,m}{{\rm{e}}}^{{\rm{i}}\omega (m-n)t}{\rho }_{n,m}+{\rm{h}}.{\rm{c}}.\right),\end{eqnarray} \tag{ 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={\prod }_{j\in 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 equation (4) that

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{t\to \infty }\langle X\rangle (t)=\displaystyle \sum _{n,m;n\gt m}{D}_{n,m}\cos \left(\omega (m-n)t\right)+\mathrm{const}.,\end{eqnarray} \tag{ 5 }$$

with D_n,m = 2Tr(Xρ_n,m)C_n,m. Provided that ${D}_{n,m}\ne 0$ for at least one n, m where $| n-m| =0$ then equation (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 equation (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 ${\rho }_{\infty }(t)$ of the system. If these requirements are met then, through equation (3), the ρ_n,m and thus ${\rho }_{\infty }(t)$ will inherit the translational symmetry of these operators and the limit cycles described in equation (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 long-time 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 section 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 limit-cycle, 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 emphasise that highly controllable quantum systems such as lattices of ultracold atoms immersed in a BEC [25, 31] can be engineered to accurately implement dynamics described by local master equations.

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 equation (2) arises due to some homogeneous field $F=\omega {\sum }_{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={\sum }_{j}{\omega }_{j}{f}_{j}$, then the Hamiltonian H can be split into two terms. The first contains a homogeneous term ${\sum }_{j}{\bar{\omega }}_{j}{f}_{j}$ and all other non-field terms, the second contains only the inhomogeneous part ${\sum }_{j}{\delta }_{j}{f}_{j}$. We have parametrised ω_j as ${\omega }_{j}={\bar{\omega }}_{j}+{\delta }_{j}$ where ${\bar{\omega }}_{j}$ is the average of the set {ω_j}. We then, correspondingly, split the Liouvillian into two parts and scale by $1/{\bar{\omega }}_{j}$

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal L } & = & {{ \mathcal L }}^{(0)}+\epsilon {{ \mathcal L }}^{(1)},\\ {{ \mathcal L }}^{(0)} & = & -\displaystyle \frac{{\rm{i}}}{{\bar{\omega }}_{j}}[H-\displaystyle \sum _{j}{\delta }_{j}{f}_{j},\cdot ]+\displaystyle \frac{1}{{\bar{\omega }}_{j}}D[\cdot ],\\ {{ \mathcal L }}^{(1)} & = & -{\rm{i}}\left[\displaystyle \sum _{j}\displaystyle \frac{{\delta }_{j}}{{\bar{\delta }}_{j}}{f}_{j},\cdot \right],\end{array}\end{eqnarray} \tag{ 6 }$$

where $\epsilon ={\bar{\delta }}_{j}/{\bar{\omega }}_{j}$ is the, small, perturbation parameter and ${\bar{\delta }}_{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 ${{ \mathcal L }}^{(0)}$

$$\begin{eqnarray}\begin{array}{rcl}\rho & = & {\rho }^{(0)}+\epsilon {\rho }^{(1)}+{\epsilon }^{2}{\rho }^{(2)}+\ldots ,\\ \lambda & = & {\lambda }^{(0)}+\epsilon {\lambda }^{(1)}+{\epsilon }^{2}{\lambda }^{(2)}+\ldots .\end{array}\end{eqnarray} \tag{ 7 }$$

It then follows (see supplemental material, SM is available online at stacks.iop.org/NJP/22/013026/mmedia) that the first order eigenvalue shift is ${\lambda }^{(1)}=\mathrm{Tr}\left[{\left({\rho }^{(0)}\right)}^{\dagger }{{ \mathcal L }}^{(1)}{\rho }^{(0)}\right]$, 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 equation (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 λ⁽¹⁾ = 0 as $\mathrm{Tr}\left[{\left({\rho }^{(0)}\right)}^{\dagger }{f}_{j}{\rho }^{(0)}\right]$ is independent of j. As discussed in section 2.2 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 nonlinear response to peturbations away from the dynamical symmetry regime. To first order, the solutions in equation (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].

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 $\{| \downarrow \rangle ,| 0\rangle ,| \uparrow \rangle \}$. 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: ${S}_{j}^{x}=(1/2)({S}_{j}^{+}+{S}_{j}^{-})$, ${S}_{j}^{y}=(i/2)({S}_{j}^{-}-{S}_{j}^{+})$. By dropping the local subscript we denote the total of an operator, e.g. ${S}^{z}={\sum }_{j}{S}_{j}^{z}$.

We take a spin-1 anisotropic Heisenberg model in a chain geometry [33] (see figure 1)

$$\begin{eqnarray}&&H=\displaystyle \sum _{j=1}^{N}{\omega }_{j}{S}_{j}^{z}+\displaystyle \sum _{j=1}^{N-1}J\left({S}_{j}^{+}{S}_{j+1}^{-}+{S}_{j}^{-}{S}_{j+1}^{+}\right)+{\rm{\Delta }}{S}_{j}^{z}{S}_{j+1}^{z},\end{eqnarray} \tag{ 8 }$$

where the spins each have natural frequency ω_j and nearest-neighbour coupling strengths J and ${\rm{\Delta }}\ne 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

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}={ \mathcal L }\rho =-{\rm{i}}[H,\rho ]+\gamma \displaystyle \sum _{j=1}^{N}{\left({S}_{j}^{z}\right)}^{2}\rho {\left({S}_{j}^{z}\right)}^{2}-\displaystyle \frac{1}{2}\{{\left({S}_{j}^{z}\right)}^{4},\rho \},\end{eqnarray} \tag{ 9 }$$

which is the master equation in equation (1) with jump operators ${L}_{j}={\left({S}_{j}^{z}\right)}^{2}$ applied to each spin at a rate $\gamma \ \forall j$.

Zoom In Zoom Out Reset image size

Initially, in order to derive an analytical solution to the long-time dynamics, we focus on the 'frequency-matched' case ${\omega }_{j}=\omega ,\forall j$. In the SM we prove that the 2N + 2 steady states alway take the form

$$\begin{eqnarray}&&{\rho }_{{ss}}=\displaystyle \sum _{m=-N}^{N}{\lambda }_{m}\left(\displaystyle \sum _{i}| {m}_{i}\rangle \langle {m}_{i}| \right)+{\lambda }_{0}^{\prime} \displaystyle \sum _{i}| {0}_{i}\rangle \langle 0{{\prime} }_{i}| ,\end{eqnarray} \tag{ 10 }$$

where $| {m}_{i}\rangle$ is an eigenvector of S^z with eigenvalue m: ${S}^{z}| {m}_{i}\rangle =m| {m}_{i}\rangle$, and i indexes the possible eigenvectors for each m. We have also defined $| -m{{\prime} }_{i}\rangle =\mathrm{SF}| {m}_{i}\rangle$, where $\mathrm{SF}={\otimes }_{j=1}^{N}\left(| \uparrow \rangle \langle \downarrow | +| \downarrow \rangle \langle \uparrow | +| 0\rangle \langle 0| \right)$ is the spin-flip operator. For example if $| {2}_{1}\rangle =| 0\uparrow \uparrow \rangle$ then $| -{2}_{1}^{\prime} \rangle =| 0\downarrow \downarrow \rangle$, or if $| {0}_{2}\rangle =| 0\uparrow \downarrow \rangle$ then $| {0}_{2}^{\prime} \rangle =| 0\downarrow \uparrow \rangle$. In order for Tr(ρ_ss) = 1 the elements {λ_m} and ${\lambda }_{0}^{\prime}$ must satisfy the equation

$$\begin{eqnarray}&&{\lambda }_{0}^{\prime} +\displaystyle \sum _{m=-N}^{N}{\lambda }_{m}\displaystyle \sum _{s=0}^{N}\displaystyle \left(\genfrac{}{}{0em}{}{N}{s}\right)\displaystyle \left(\genfrac{}{}{0em}{}{N-s}{(N-s+m)/2}\right)=1,\end{eqnarray} \tag{ 11 }$$

where the terms in the second summation are skipped if (N − s + m)/2 is not an integer.

The long-time dynamics of ${ \mathcal L }$ is not, however, solely governed by this steady state. We identify (see SM) multiple strong dynamical symmetry operators of the form

$$\begin{eqnarray}&&{A}_{m}=\displaystyle \sum _{i}| {m}_{i}\rangle \langle -{m}_{i}^{\prime} | ,\quad m\ne 0,m=-N,...,N,\end{eqnarray} \tag{ 12 }$$

which each satisfy equation (2)

$$\begin{eqnarray}&&[H,{A}_{m}]=2m\omega {A}_{m},\quad [{L}_{j},{A}_{m}]=[{L}_{j}^{\dagger },{A}_{m}]=0,\quad \forall j.\end{eqnarray} \tag{ 13 }$$

Following this we can determine, see equation (3), the imaginary eigenmodes of the Liouvillian through the action of these operators on the steady state. Explicitly, we have

$$\begin{eqnarray}&&{\rho }_{1,0}^{m}\propto {A}_{m}{\rho }_{{ss}}\propto {A}_{m},\qquad { \mathcal L }{\rho }_{1,0}^{m}=-2{\rm{i}}{m}\omega ,\end{eqnarray} \tag{ 14 }$$

which is a non-trivial result as the steady-state ρ_ss is inherently singular. Further application of A_m is redundant as A_mA_m = 0 and thus, each A_m generates a unique imaginary eigenmode ${\rho }_{1,0}^{m}$ 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 equations (10) and (12) form a complete basis for the long-time dynamics of equation (9) and so, similarly to equation (4), the density matrix can be expressed as a superposition of these modes in the limit $t\to \infty$. 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 section 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 equation (12), which occur at all length-scales of the chain and are completely uniform with respect to distance.

In figure 2 we visualise the analytical results in equations (10), (12) and (14). We present a plot of the eigenspectrum of ${ \mathcal L }$ (figure 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 equation (9) (figure 2(b)). The system is in a superposition of the steady state in equation (10) and the imaginary modes in equation (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 $| {m}_{i}\rangle \langle -m{{\prime} }_{i}|$.

Zoom In Zoom Out Reset image size

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 $X={\prod }_{j\in B}{({S}_{j}^{x})}^{2}$, where B = {a, b, c, ...} 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 $| \downarrow \rangle$ and $| \uparrow \rangle$—any operator formed solely from linear local operators will relax to stationarity. We prove (see SM) that in the limit $t\to \infty$

$$\begin{eqnarray}&&\langle X\rangle ={D}_{0}\mathrm{Tr}({\rho }_{{ss}}X)+\displaystyle \sum _{m=1}^{M}{D}_{m}\cos (2m\omega t)\mathrm{Tr}({\rho }_{1,0}^{m}X),\end{eqnarray} \tag{ 15 }$$

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 $\mathrm{Tr}({\rho }_{1,0}^{m}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 $\langle X\rangle$ 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-cycles - regardless of the specific values of the initial state or the Liouvillian parameters.

Moreover, the modes in equation (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 $\langle X\rangle -{\prod }_{j\in B}\langle {X}_{j}\rangle$ is distance-invariant and non-zero. Hence, we consider the synchronisation observed in equation (15) to be inherently quantum—underpinned by long-range correlations, which are a result of entanglement between the bodies in the system.

In equation (8), the ${S}_{i}^{z}{S}_{i+1}^{z}$ (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 synchronisation (see section 2.2). Provided ${\rm{\Delta }}\ne 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 equation (9). 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}| \gg | {D}_{m\ne 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 equation (9) through numerical exponentiation of the Liouvillian superoperator ${ \mathcal 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 equation (9) 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 Pearson-correlation 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]

$$\begin{eqnarray}&&{C}_{f,g}^{{\rm{\Delta }}t}(t)=\displaystyle \frac{{\int }_{t}^{t+{\rm{\Delta }}t}(f-\bar{f})(g-\bar{g}){\rm{d}}{t}}{\sqrt{{\int }_{t}^{t+{\rm{\Delta }}t}{\left(f-\bar{f}\right)}^{2}{\rm{d}}{t}{\int }_{t}^{t+{\rm{\Delta }}t}{\left(g-\bar{g}\right)}^{2}{\rm{d}}{t}}},\end{eqnarray} \tag{ 16 }$$

with the function average $\bar{f}=\tfrac{1}{{\rm{\Delta }}t}{\int }_{t}^{t+{\rm{\Delta }}t}{f}{\rm{d}}{t}$. This correlation factor is maximal (minimal), 1 (−1) when the two signals f and g are perfectly synchronised (anti-synchronised), and 0 when they display no correlations.

In figure 3 we set $f=\langle {({S}_{j}^{x})}^{2}\rangle$ and $g=\langle {({S}_{l}^{x})}^{2}\rangle$ 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 $\langle {({S}_{j}^{x})}^{2}\rangle$ over time for each spin. In agreement with equation (15), when both the environment and interactions are present ($\gamma ,{\rm{\Delta }}\ne 0$) the dynamics causes the spins to synchronise perfectly (figures 3(a) and (d)) to the same frequency and phase, despite being initialised with random phases. The frequency of the oscillations is directly determined by the equidistant spacing of the imaginary eigenvalues. For comparison (figure 3(b)) we show the case when the environment is present but there are no interactions ($\gamma \ne 0,{\rm{\Delta }}=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 (figures 3(c) and (f)), the dynamics are completely chaotic and unsynchronised due to the multitude of incommensurate frequencies in the eigenvalues of the Hamiltonian.

Zoom In Zoom Out Reset image size

Whilst the plots in figure 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 figure 4, we plot the reduced correlator ${ \mathcal A }\,=\langle {({S}_{j}^{x})}^{2}{({S}_{l}^{x})}^{2}\rangle -\langle {({S}_{j}^{x})}^{2}\rangle \langle {({S}_{l}^{x})}^{2}\rangle$, which is non-zero and identical for any choice of spins j and l. The observed oscillations contain several frequencies (figure 4(b)) due to the excitement of multiple imaginary modes in equation (12). These imaginary modes all contain coherences between different spins, the synchronisation observed in figure 3 is dependent on the existence of these inter-spin coherences and figure 4 shows that they give rise to perfect, distance-invariant correlations throughout the system.

Zoom In Zoom Out Reset image size

So far we have considered the 'homogeneous' case ${\omega }_{j}=\omega \ \forall 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 equation (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 equation (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 figures 3 and 4, we introduce two common witnesses for quantum synchronisation: the negativity ${ \mathcal N }$ [35] and off-diagonal coherences ${ \mathcal C }$ [11]

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal N }}_{j}(\rho ) & = & \displaystyle \frac{| | {\rho }^{{T}_{j}}| | -1}{2},\\ { \mathcal C } & = & \displaystyle \sum _{i\ne j}| {\rho }_{{ij}}| ,\end{array}\end{eqnarray} \tag{ 17 }$$

with T_j indicating the partial transpose with respect to site j and $| | X| | =\mathrm{Tr}\sqrt{{X}^{\dagger }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 ${ \mathcal 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 equation (3) quantities such as these will remain finite indefinitely.

In figure 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 ${ \mathcal 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 (figures 5(c), (d)): 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 ${\bar{\omega }}_{j}$ (due to the factor of 2 in equation (13)). The life-time of this cycle is large and as $\delta \to 0$ diverges to infinity, evidenced by the tongue-like behaviour seen in figures 5(a), (b). These figures show how the corresponding measures act as strong witnessess to the synchronisation in the system—emphasising its quantum nature.

Zoom In Zoom Out Reset image size

The imaginary eigenmodes in equation (3) are translationally-invariant and hence the robust, synchronised behaviour observed (figure 5) is the result of a second-order response to the detuning. The cross-sections included in figure 5(b) 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 nonlinear scaling in figure 6. We calculate the shift in the imaginary eigenvalues, $| \lambda -{\lambda }^{(0)}|$ (see section 2.3) from their original value λ⁽⁰⁾ at δ = 0 as a function of the detuning δ. There is no noticeable shift to first order in δ—the fitted curve is proportional to δ². Notably, the highest imaginary eigenmode ${\rho }_{1,0}^{N}=| \uparrow \uparrow ...\rangle \langle \downarrow \downarrow ...|$ is always unshifted, and remains imaginary regardless of the distribution of natural frequencies.

Zoom In Zoom Out Reset image size

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 emphasised that these results are solely based on symmetry and thus can be observed in any bi-partite d-dimensional realisation of the Hubbard model. The Hamiltonian reads

$$\begin{eqnarray}&&H=-\tau \displaystyle \sum _{\langle {jl}\rangle ,\sigma }({c}_{\sigma ,j}^{\dagger }{c}_{\sigma ,l}+{\rm{h}}.{\rm{c}})+U\displaystyle \sum _{j}{n}_{\uparrow ,j}{n}_{\downarrow ,j}+\displaystyle \frac{1}{2}\displaystyle \sum _{j}{\omega }_{j}({n}_{\uparrow ,j}-{n}_{\downarrow ,j})+\displaystyle \sum _{j}{\mu }_{j}({n}_{\uparrow ,j}+{n}_{\downarrow ,j}),\end{eqnarray} \tag{ 18 }$$

where ${c}_{\sigma ,j}^{\dagger }$ and its adjoint are the usual creation and annihilation operators for a fermion of spin $\sigma \in \{\uparrow ,\downarrow \}$ on site j. Additionally, ${n}_{\sigma ,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:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}={ \mathcal L }\rho =-{\rm{i}}[H,\rho ]+\gamma \displaystyle \sum _{j=1}^{N}{n}_{j}\rho {n}_{j}-\displaystyle \frac{1}{2}\{{\left({n}_{j}\right)}^{2},\rho \},\qquad {n}_{j}={n}_{j,\uparrow }+{n}_{j,\downarrow }.\end{eqnarray} \tag{ 19 }$$

The case ${\omega }_{j}=\omega \ \forall j$ of this model was originally studied in [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 long-time dynamics of this model: alongside the strong dynamical symmetry there is an inter-site coupling and homogeneous local dephasing which ensure the long-time dynamics is completely cooperative and translationally invariant. Moreover, later in the section we break the homogeneity of the magnetic field and demonstrate the robustness of this synchronisation to perturbations away from the dynamical symmetry regime.

Firstly, when ${\omega }_{j}=\omega \ \forall j$, the magnetic field in equation (18) breaks the spin SU(2) symmetry of the model: [H, S^± ] = ωS^±, $[{L}_{j},{S}^{\pm }]=0\ \forall j$, where ${S}^{+}={\sum }_{j}{c}_{j,\uparrow }^{\dagger }{c}_{j,\downarrow }$ is the global magnetic raising operator. Hence, there is a single strong dynamical symmetry operator which can be used to form the set of equidistant imaginary eigenmodes

$$\begin{eqnarray}&&{\rho }_{{nm}}\propto {\left({S}^{+}\right)}^{n}{\rho }_{{ss}}{\left({S}^{-}\right)}^{m},{ \mathcal L }{\rho }_{{nm}}={\rm{i}}\omega (m-n){\rho }_{{nm}}.\end{eqnarray} \tag{ 20 }$$

These modes cause the existence of a persistent limit cycles in the magnetisation (in the x and y directions) of the system. For example, by defining the operator $X={\prod }_{j\in B}{S}_{j}^{x}$ (where B = {a, b, c ..} is a set of M sites containing no duplicates) we can use equation (5) to prove (see SM):

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{t\to \infty }\langle X\rangle (t)=\displaystyle \sum _{i=0}^{\lfloor M/2\rfloor }{D}_{i}\cos \left(\left(2i+d\right)\omega t\right),d=M\ \mathrm{mod}\ 2,\end{eqnarray} \tag{ 21 }$$

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 equation (20) are completely translationally symmetric (see [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 ${\mathrm{lim}}_{N\to \infty }\langle {S}^{x}\rangle /N\ne 0$ will ensure the long-time synchronised oscillations in $\langle {S}_{j}^{x}\rangle$ or $\langle {S}_{{jS}}^{x}{S}_{j+1}^{x}\rangle$ 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 equation (19) and simulate the dynamics at the level of an ensemble of pure wavefunctions. Furthermore in figure 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].

Zoom In Zoom Out Reset image size

In figure 7 we demonstrate the synchronicity which results from the eigenmodes in equation (20). We initialise the system in a product state and, after quenching under the master equation in equation (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 (figure 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 − ω₁ = 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 equation (19). In figure 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 (figures 8(a) 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 figure 8(d)). Remarkably, this harmonised response is occuring even in the presence of both magnetic and chemical disorder—emphasising the robustness of a symmetry-based approach to observing quantum synchronisation. The imaginary eigenmodes in equation (20) are translationally-invariant and hence λ⁽¹⁾ evaluates to 0 (see section 2.3). As with the previous example, this robust, synchronised behaviour is a result of a second-order response to the detuning.

Zoom In Zoom Out Reset image size