Exotic synchronization in continuous time crystals outside the symmetric subspace

Parvinder Solanki [email protected] Department of Physics, University of Basel, Klingelbergstrasse 82, CH-4056 Basel, Switzerland Midhun Krishna Department of Physics, Indian Institute of Technology-Bombay, Powai, Mumbai 400076, India Michal Hajdušek [email protected] Keio University Shonan Fujisawa Campus, 5322 Endo, Fujisawa, Kanagawa 252-0882, Japan Keio University Quantum Computing Center, 3-14-1 Hiyoshi, Kohoku, Yokohama, Kanagawa 223-8522, Japan Christoph Bruder Department of Physics, University of Basel, Klingelbergstrasse 82, CH-4056 Basel, Switzerland Sai Vinjanampathy [email protected] Department of Physics, Indian Institute of Technology-Bombay, Powai, Mumbai 400076, India Centre of Excellence in Quantum Information, Computation, Science and Technology, Indian Institute of Technology Bombay, Powai, Mumbai 400076, India Centre for Quantum Technologies, National University of Singapore, 3 Science Drive 2, Singapore 117543

Abstract

Exploring continuous time crystals (CTCs) within the symmetric subspace of spin systems has been a subject of intensive research in recent times. Thus far, the stability of the time-crystal phase outside the symmetric subspace in such spin systems has gone largely unexplored. Here, we investigate the effect of including the asymmetric subspaces on the dynamics of CTCs in a driven dissipative spin model. This results in multistability, and the dynamics becomes dependent on the initial state. Remarkably, this multistability leads to exotic synchronization regimes such as chimera states and cluster synchronization in an ensemble of coupled identical CTCs. Interestingly, it leads to other nonlinear phenomena such as oscillation death and signature of chaos.

Introduction.— Crystals have traditionally been characterized by their spatial periodicity, wherein atoms are arranged in a repeating pattern. Time crystals showcase a distinct form of periodicity, namely in the temporal dimension Wilczek (2012); Watanabe and Oshikawa (2015); Zaletel et al. (2023); Sacha and Zakrzewski (2017). Initially, time crystals breaking discrete time-translation symmetry were conceptualized and termed discrete time crystals Sacha (2015); Else et al. (2016); Khemani et al. (2016); Russomanno et al. (2017); Surace et al. (2019); Heugel et al. (2019); Khasseh et al. (2019); Sakurai et al. (2021); Sacha and Zakrzewski (2017); Zaletel et al. (2023). A more recent development in the field has been the emergence of continuous time crystals (CTCs) Iemini et al. (2018); Tucker et al. (2018); Buča et al. (2019); Booker et al. (2020); Buča and Jaksch (2019); Zhu et al. (2019); Lledó et al. (2019); Seibold et al. (2020); Prazeres et al. (2021); Piccitto et al. (2021); Carollo and Lesanovsky (2022); Hajdušek et al. (2022); Krishna et al. (2023a); Hurtado-Gutiérrez et al. (2020); Cosme et al. (2023); Buonaiuto et al. (2021); Cabot et al. (2023a), which break continuous time-translation symmetry. CTCs were initially discussed in the context of a driven dissipative Dicke model Iemini et al. (2018). Their behavior was also studied in other spin systems like collective $d$ -level systems Prazeres et al. (2021), $p$ - $q$ interacting spin models Piccitto et al. (2021), spin star models Krishna et al. (2023a), and spin systems with infinite-range interactions Cosme et al. (2023). In addition to the time-crystal behavior, the thermodynamics Carollo et al. (2024) of driven Dicke models and the occurrence of entanglement phase transitions Passarelli et al. (2022); Cabot et al. (2023b) have also been investigated recently. The application of time-crystal phases in sensing Cabot et al. (2024); Montenegro et al. (2023) and quantum engines Carollo et al. (2020) has been proposed. Furthermore, experimental platforms such as Bose-Einstein condensate systems have been used to realize CTCs Keßler et al. (2021); Kongkhambut et al. (2022).

All of these examples are variants of the Dicke model: the system is a collective spin that is built up of $N$ spins 1/2. These systems are characterized by permutational symmetry which is a strong symmetry of the underlying open quantum system Buča and Prosen (2012); Albert and Jiang (2014a). Previous investigations of both superradiance Dicke (1954); Andreev et al. (1980); Wang and Hioe (1973); Scully and Svidzinsky (2009) and time-crystal behavior have focused on the totally symmetric Dicke subspace, which corresponds to the maximum-excitation subspace. The assumption of permutational symmetry is not always straightforward to enforce under realistic conditions Ferioli et al. (2023). This leads to the important open question of the role played by the lower-excitation subspaces involving asymmetric Dicke states in the dynamics of time crystals. Recent studies have shown that asymmetric subspaces are crucial for understanding subradiant phases Bhatti et al. (2018); Bojer and von Zanthier (2022); Gegg et al. (2018); Wiegner et al. (2011) and the thermodynamics of permutation-invariant quantum many-body systems Latune et al. (2019, 2020); Yadin et al. (2023).

In this Letter, we investigate the existence of time crystals outside the symmetric Dicke subspace in a spin-only description of the driven dissipative Dicke model Iemini et al. (2018); Hajdušek et al. (2022). We show that incorporating the asymmetric subspaces enhances the complexity and richness of CTC dynamics and yields an extended phase diagram. Using mean-field analysis, we find that CTCs exhibit multistability in the asymmetric subspaces, which was not observed for CTCs restricted to the symmetric subspace. The existence of multistability is confirmed by analyzing the eigenspectra of the Liouville superoperator that describes the Markovian dynamics of open quantum systems Buča and Prosen (2012); Albert and Jiang (2014b); Krishna et al. (2023b). Depending on the system parameters, an initial state outside the symmetric subspace can always be found which breaks time-translational symmetry. Interestingly, this initial-state dependence leads to exotic synchronization phenomena in a network of coupled identical CTCs like chimera states that are characterized by the coexistence of synchronized and unsynchronized regions Kuramoto and Battogtokh (2002); Abrams and Strogatz (2004); Yeldesbay et al. (2014); Chandrasekar et al. (2015); Ujjwal et al. (2017a). Furthermore, such a network of CTCs also exhibits cluster synchronization, where identical oscillators arrange themselves into two or more differently synchronized clusters Haugland (2021). Finally, we show that this highly nonlinear coupled system also features chaotic mean-field dynamics and the phenomenon of oscillation death.

Refer to caption — Figure 1: (a) Phase diagram of a continuous time crystal (CTC) characterized using fixed-point analysis of mean-field equations with $m=\langle\hat{S}\rangle/S$ . (b) Schematic representation of block diagonal structure of Liouville superoperator $\mathcal{L}=\mathcal{L}_{S}+\mathcal{L}_{A}$ where $\mathcal{L}_{S(A)}$ corresponds to symmetric (asymmetric) Dicke subspace. Here, dotted lines separate the blocks with $m>\Omega/\kappa$ from those with $m<\Omega/\kappa$ . Only the sub-blocks having rescaled magnetization value $m<\Omega/\kappa$ will have purely imaginary eigenvalues in the thermodynamic limit $N=2S\rightarrow\infty$ .

We begin with an investigation of the mean-field dynamics of a collective spin using fixed-point analysis, which describes the steady-state dynamics of the system. Then, we will use Liouville superoperator theory to validate the mean-field results. Finally, we discuss the existence of a chimera state, cluster synchronization, oscillation death and chaos in coupled CTCs before summarizing our results.

Mean-field analysis and phase diagram.— The master equation governing the dynamics of the CTC Iemini et al. (2018); Hajdušek et al. (2022) is given by

\dot{\rho}=\mathcal{L}_{\Omega,\kappa}(\rho)=-i[\Omega\hat{S}_{x},\rho]+\frac{\kappa}{S}(\hat{S}_{-}\rho\hat{S}_{+}-\frac{1}{2}\{\hat{S}_{+}\hat{S}_{-},\rho\}),

(1)

where $\mathcal{L}_{\Omega,\kappa}$ is the Liouville superoperator, $\hat{S}_{\alpha}=\sum_{i}\sigma^{i}_{\alpha}/2$ are collective spin operators where $\alpha\in\{x,y,z\}$ , $S=N/2$ is the maximum spin value for $N$ spins $1/2$ and $\hat{S}_{\pm}=\hat{S}_{x}\pm i\hat{S}_{y}$ . The dynamics of the system in the thermodynamic limit ( $N\rightarrow\infty$ ) can be well understood using mean-field analysis. The evolution of mean-field expectation values of the rescaled spin operators, $m_{\alpha}=\langle\hat{S}_{\alpha}\rangle/S$ is governed by

$\displaystyle\dot{m}_{x}$	$\displaystyle=\kappa m_{x}m_{z},$
$\displaystyle\dot{m}_{y}$	$\displaystyle=-\Omega m_{z}+\kappa m_{y}m_{z},$	(2)
$\displaystyle\dot{m}_{z}$	$\displaystyle=\Omega m_{y}-\kappa[m_{x}^{2}+m_{y}^{2}].$

The total spin of the system is conserved since $[\hat{S}^{2},S_{x,y,z,\pm}]=0$ and hence is a strong symmetry of the system Buča and Prosen (2012); Albert and Jiang (2014b); Manzano and Hurtado (2018). At the mean-field level, this implies that $m^{2}_{x}+m^{2}_{y}+m^{2}_{z}=m^{2}$ is a conserved quantity since correlations are assumed to vanish Mukherjee et al. (2024).

Using the fixed-point analysis, we can characterize the phase transition of the system as described by Fig. 1(a) (see sup for details). Initial states on the left of the $m=\Omega/\kappa$ (blue) line ( $m>\Omega/\kappa$ ) settle down to a saddle fixed point. In the parameter region $m<\Omega/\kappa$ , the systems exhibit stable oscillations around the center fixed point leading to time-translational symmetry breaking. The frequency of oscillations is determined by the imaginary eigenvalues of the Jacobian matrix $\pm i\kappa\sqrt{(\Omega/\kappa)^{2}-m^{2}}$ . The symmetric Dicke subspace corresponds to $m=1$ , which was shown to lead to the breaking of time-translational symmetry at $\Omega/\kappa=1$ in earlier studies Iemini et al. (2018); Hajdušek et al. (2022). The system exhibits multistability for $\Omega/\kappa<1$ . Depending upon the initial state, the system can exhibit two different phases, namely, a melted and a time-crystal phase.

Role of asymmetric subspaces.— We now go beyond mean-field theory and study the full master equation Eq. (1). The Hilbert space of the ensemble consisting of $N$ spins can be described as a direct sum of irreducible representations of $SU(2)$ using the Clebsch-Gordan decomposition Sakurai and Commins (1995). Since $\hat{S}^{2}$ is a strong symmetry of the dynamics, the Liouville superoperator block-diagonalizes Buča and Prosen (2012); Thingna and Manzano (2021) in respective spin- $J$ eigenbases such that $\mathcal{L}=\oplus\mathcal{L}_{J}$ , as shown in Fig. 1(b). Therefore, the dynamics is constrained within each spin- $J$ subspace (invariant subspaces of dynamics). The maximum angular momentum value $J=S=N/2$ corresponds to the totally symmetric Dicke subspace; lower $J$ values correspond to asymmetric subspaces. The spin- $J$ irreducible representation has dimension $2J+1$ and multiplicity $n_{J}=N!(2J+1)/(N/2+J+1)!(N/2-J)!$ Georgi (2000); Ramadevi and Dubey (2019); Yadin et al. (2023). The total Hilbert space is spanned by the generalized Dicke basis, $\{|J,m_{J},\mu_{J}\rangle\}$ . Here, $J\in\{S,S-1\ldots J_{0}\}$ , $m_{J}\in\{J,J-1\ldots-(J-1),-J\}$ and $1\leq\mu_{J}\leq n_{J}$ , with $J_{0}=0\text{ or }1$ depending on whether $N$ is even or odd. We can associate each of the spin- $J$ subspaces with the mean-field total angular momentum value $m=J/S$ . The states with $m<1$ belong to the asymmetric angular momentum sectors.

The time-crystal phase is characterized by the spectral gap closing of the Liouvillian eigenspectra in the thermodynamic limit, with the leading eigenvalue being purely imaginary. For the symmetric subspace, the spectral gap of $\mathcal{L}_{S}$ closes only for $\Omega/\kappa>1$ . Including the asymmetric subspaces $\mathcal{L}_{A}=\oplus_{J<S}\mathcal{L}_{J}$ introduces new leading-order eigenvalues, as shown in Fig. 2(a). Therefore, the spectral gap closes even for $\Omega/\kappa<1$ , as depicted in Figs. 2(b)–(c). The first dominant eigenvalue in the asymmetric subspace scales differently for even and odd $N$ , but converges to the same value for large values of $N$ , as shown in Fig. 2(b). The real part of the eigenvalue $\lambda_{1}$ for both even and odd values of $N$ scales as $1/N$ with the same slope as shown in Fig. 2(c) and vanishes in the thermodynamic limit. This means that the spectral gap closes even for $\Omega/\kappa<1$ for the asymmetric subspaces. Therefore, time crystal behavior can be observed outside the permutational symmetric subspace even for $\Omega/\kappa<1$ if the initial state belongs to an asymmetric subspace having $m<\Omega/\kappa$ .

Chimera state and cluster synchronization.— We now discuss the effect of including asymmetric subspaces on the dynamics of $n$ coherently coupled CTCs described by the following master equation,

\dot{\rho}=\sum_{\alpha}^{n}\mathcal{L}_{\Omega_{\alpha},\kappa_{\alpha}}(\rho)-\frac{i}{n}[\sum_{\beta\neq\alpha}^{n}\Gamma_{\alpha,\beta}(\hat{S}^{\alpha}_{+}\hat{S}^{\beta}_{-}+\hat{S}^{\alpha}_{-}\hat{S}^{\beta}_{+}),\rho].

(3)

Here, $\mathcal{L}_{\Omega_{\alpha},\kappa_{\alpha}}(\rho)$ describes the dynamics of uncoupled CTCs, and $\Gamma_{\alpha,\beta}$ is the coupling strength between two CTCs $\alpha$ and $\beta$ . Such an interaction resulted in the observation of the seeding effect and the synchronization of detuned CTCs considering only the symmetric subspace Hajdušek et al. (2022). We now also take into account asymmetric subspaces and assume that all the CTCs are identical, $\mathcal{L}_{\Omega_{\alpha},\kappa_{\alpha}}=\mathcal{L}_{\Omega,\kappa},~{}\forall\alpha$ . Even if they are identical, different choices of initial states can lead to surprising disparate effects in the dynamics. We investigate the mean-field dynamics of $n=20$ coupled CTCs in Fig. 3, where all CTCs are initialized in asymmetric subspaces belonging to different values of $m$ . The corresponding frequency and the dynamics of $m_{z}^{\alpha}$ of uncoupled CTCs depending on the value $m$ are shown in Fig. 3(a) and Fig. 3(d), respectively. For simplicity, we consider two subgroups based on the separation between the frequency distribution of uncoupled CTCs, though a general distribution could result in more than two subgroups. In Fig. 3(a), the system is split into two ensembles, ${\varepsilon_{1}}$ (black) and ${\varepsilon_{2}}$ (red), with uniform frequency distributions: ${\varepsilon_{1}}$ ranges from 0.2 to 0.25, and ${\varepsilon_{2}}$ from 0.75 to 0.85.

We use Pearson’s correlation coefficient as an indicator of synchronization. For any two coupled CTCs $\alpha$ and $\beta$ , it can be defined as follows

\mathcal{C}_{\alpha\beta}=\frac{\langle{m^{z}_{\alpha}m^{z}_{\beta}}\rangle-\langle{m^{z}_{\alpha}}\rangle\langle{m^{z}_{\beta}}\rangle}{\sqrt{(\langle{m^{z}_{\alpha})^{2}}\rangle-\langle{m^{z}_{\alpha}}\rangle^{2})(\langle{(m^{z}_{\beta})^{2}}\rangle-\langle{m^{z}_{\beta}}\rangle^{2})}},

(4)

where $\langle f\rangle=\frac{1}{T}\int_{0}^{T}f(t)dt$ represents the time average. The off-diagonal elements $\mathcal{C}_{\alpha\neq\beta}$ of the correlation matrix $\mathcal{C}_{n\times n}$ indicates the amount of correlation between the CTCs $\alpha$ and $\beta$ . Therefore, $\mathcal{C}_{\alpha,\beta}\rightarrow 1$ corresponds to complete synchronization and $\mathcal{C}_{\alpha,\beta}\rightarrow 0$ indicates no synchronization. E.g., Fig. 3(g) shows that this matrix is diagonal for uncoupled CTCs.

First, we consider the case where CTCs within the same ensemble are coupled with a finite interaction strength but there is no coupling between CTCs of different ensembles such that $\Gamma_{\alpha,\beta}=\Gamma\delta_{\varepsilon_{i}\varepsilon_{j}}~{}\forall~{}(\alpha\in\varepsilon_{i},\beta\in\varepsilon_{j})$ . This results in the synchronization of CTCs within the same ensemble as shown by Fig. 3(b,e,h). Surprisingly, considering all-to-all coupled CTCs ( $\Gamma_{\alpha,\beta}/\kappa=\Gamma/\kappa=0.35~{}\forall~{}(\alpha,\beta)$ ) leads to the synchronization of the CTCs only in one ensemble (here $\varepsilon_{1}$ ), while the CTCs in the other ensemble (here $\varepsilon_{2}$ ) exhibit unsynchronized dynamics as shown in Fig. 3(c,f,i). Such surprising patterns of partial synchronization in systems of coupled identical nonlinear oscillators were first observed in Kuramoto and Battogtokh (2002) and were called chimera states Abrams and Strogatz (2004); Zakharova (2020). First strides towards chimera states in quantum systems involved exploring one-dimensional arrays of quantum oscillators Bastidas et al. (2015) and periodically driven networks of spins Sakurai et al. (2021). Our work demonstrates a chimera state in a system of coupled identical CTCs. Here, the chimera state results from including the asymmetric subspaces into the dynamics of identical CTCs, allowing us to choose an initial state with different frequencies in different subspaces Ujjwal et al. (2017b). Increasing the coupling strength results in the synchronization of all CTCs, see sup .

Another example of an exotic synchronization regime that can be realized in coupled identical CTCs is cluster synchronization, where the system separates into two differently synchronized clusters. Such a state is shown in Fig. 4(b), where the group $\varepsilon_{1}$ (black) and the group $\varepsilon_{2}$ (red) of CTCs get synchronized separately at different frequencies. This results from the choice of the initial state shown in Fig. 4(a), where the frequency distribution of the ensemble $\varepsilon_{1}$ is the same as in Fig. 3(a), while it ranges from 0.75 to 0.8 for $\varepsilon_{2}$ .

Characterization of different phases.— We now extend our analysis to investigate the effect of choosing different initial states and coupling strengths $\Gamma$ on the synchronization properties of $n$ coupled CTCs. These CTCs are initialized in different $m$ subspaces such that the frequency of CTCs in each subgroup $\varepsilon_{1,2}$ are sampled from Gaussian distributions with distinct mean frequencies $\bar{\omega}_{1,2}$ and standard deviations $\sigma_{1,2}$ . We investigate the synchronization and dynamical properties as a function of coupling strength $\Gamma$ and detuning $\Delta=\bar{\omega}_{2}-\bar{\omega}_{1}$ for fixed values of $\sigma_{1,2}$ , where $\Delta$ and $\sigma_{1,2}$ depend on the choice of initial states. The mean value of Pearson’s correlation $\bar{\mathcal{C}}=\sum_{\alpha<\beta}\mathcal{C}_{\alpha\beta}/\mathcal{N}$ will be used to quantify the amount of synchronization, where $\mathcal{N}=n(n-1)/2$ . This measure takes the value $\bar{\mathcal{C}}=0$ for no synchronization, $\bar{\mathcal{C}}=1$ for complete synchronization, and $0<\bar{\mathcal{C}}<1$ in the regime of partial synchronization, such as the chimera and cluster-synchronized states.

The synchronization measure alone is not sufficient to understand the complex Arnold-tongue structure observed in Fig. 4(e). Therefore, we also investigate the largest Lyapunov exponent $\lambda_{L}$ (Fig. 4(f)), which measures the sensitivity of the dynamics of a nonlinear system to its initial conditions Datseris and Parlitz (2022); Skokos et al. (2016). A positive value of $\lambda_{L}$ corresponds to chaotic dynamics where two nearby trajectories diverge exponentially over time, whereas the system settles down to a stable fixed point for $\lambda_{L}<0$ . Finally, $\lambda_{L}\approx 0$ represents limit cycles and quasi-periodic dynamics.

In Figs. 4(e,f), we characterize different regimes making use of both the synchronization measure $\bar{\mathcal{C}}$ and the maximum Lyapunov exponent $\lambda_{L}$ . For a finite detuning $\Delta$ and small coupling strength $\Gamma$ , coupled CTCs exhibit unsynchronized dynamics indicated by $\bar{\mathcal{C}}\approx\lambda_{L}\approx 0$ . Increasing $\Gamma$ leads to a regime of partial synchronization with $\bar{\mathcal{C}}>0$ and $\lambda_{L}\approx 0$ , which includes chimera and cluster synchronization states. These two states can be distinguished by the correlation matrix $\mathcal{C}$ (see sup for more details). Further increasing $\Gamma$ leads to $\lambda_{L}<0$ and $\bar{\mathcal{C}}\approx 1$ , causing the system to settle down into a highly-correlated fixed point. This sudden termination of oscillations is termed oscillation death. Higher values of $\Gamma$ lead to a less correlated chaotic mean-field dynamics, which is characterized by $\lambda_{L}>0$ and $\bar{\mathcal{C}}<1$ . The reduced synchronization measure results from the chaotic nature of the dynamics which makes the system weakly correlated, see Fig. 4(e). Finally, even larger values of the coupling strength lead to complete synchronization such that $\bar{\mathcal{C}}=1$ . Thus, depending on the parameters and initial state configuration, the system features various exotic synchronization regimes as shown in Fig. 4(e,f).

Conclusions.— In this work, we focused on the fundamental question of whether continuous time crystals in spin models only exist if initialized in the symmetric subspace. Our first main result is that the inclusion of asymmetric subspaces leads to time-crystal behavior even if the ratio of drive to dissipation strength $\Omega/\kappa$ is less than one. A second result is that taking into account asymmetric subspaces leads to multistability: depending on the parameter values, the system can exhibit a time crystal or melted phase for different initial states. For an ensemble of coupled CTCs, this leads to our third result, the existence of chimera states and cluster synchronization, which is not possible for CTCs realized in the symmetric subspace. Finally, our system can also exhibit chaotic mean-field dynamics and oscillations death. Our results provide new insights into the dynamics of time-crystal phases and point out strategies in the design of quantum networks that may lead to the observation of exotic synchronization phenomena.

An interesting future direction is to study the effect of weak permutational symmetry-breaking interactions, which allow the system to move between various spin sectors. Such interactions are inevitable in various experimental settings and will offer further valuable insights into the fundamental characteristics of time crystals.

Acknowledgements.

Acknowledgments.— P.S. and S.V. thank R. Fazio for discussions. Numerical calculations were performed using QuTiP Johansson et al. (2012, 2013) and PIQS Shammah et al. (2018). M.H. was supported by [MEXT Quantum Leap Flagship Program] Grants No. JPMXS0118067285 and No. JPMXS0120319794. P.S. and C.B. acknowledge financial support from the Swiss National Science Foundation individual grant (grant no. 200020 200481). S.V. acknowledges support from DST-QUEST grant number DST/ICPS/QuST/Theme-4/2019. While preparing the manuscript, we became aware of related work Iemini et al. (2024).

References

Wilczek (2012) F. Wilczek, Phys. Rev. Lett. 109, 160401 (2012).
Watanabe and Oshikawa (2015) H. Watanabe and M. Oshikawa, Phys. Rev. Lett. 114, 251603 (2015).
Zaletel et al. (2023) M. P. Zaletel, M. Lukin, C. Monroe, C. Nayak, F. Wilczek, and N. Y. Yao, Rev. Mod. Phys. 95, 031001 (2023).
Sacha and Zakrzewski (2017) K. Sacha and J. Zakrzewski, Reports on Progress in Physics 81, 016401 (2017).
Sacha (2015) K. Sacha, Phys. Rev. A 91, 033617 (2015).
Else et al. (2016) D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. Lett. 117, 090402 (2016).
Khemani et al. (2016) V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 116, 250401 (2016).
Russomanno et al. (2017) A. Russomanno, F. Iemini, M. Dalmonte, and R. Fazio, Phys. Rev. B 95, 214307 (2017).
Surace et al. (2019) F. M. Surace, A. Russomanno, M. Dalmonte, A. Silva, R. Fazio, and F. Iemini, Phys. Rev. B 99, 104303 (2019).
Heugel et al. (2019) T. L. Heugel, M. Oscity, A. Eichler, O. Zilberberg, and R. Chitra, Phys. Rev. Lett. 123, 124301 (2019).
Khasseh et al. (2019) R. Khasseh, R. Fazio, S. Ruffo, and A. Russomanno, Phys. Rev. Lett. 123, 184301 (2019).
Sakurai et al. (2021) A. Sakurai, V. M. Bastidas, W. J. Munro, and K. Nemoto, Phys. Rev. Lett. 126, 120606 (2021).
Iemini et al. (2018) F. Iemini, A. Russomanno, J. Keeling, M. Schirò, M. Dalmonte, and R. Fazio, Phys. Rev. Lett. 121, 035301 (2018).
Tucker et al. (2018) K. Tucker, B. Zhu, R. J. Lewis-Swan, J. Marino, F. Jimenez, J. G. Restrepo, and A. M. Rey, New Journal of Physics 20, 123003 (2018).
Buča et al. (2019) B. Buča, J. Tindall, and D. Jaksch, Nature Communications 10, 1 (2019).
Booker et al. (2020) C. Booker, B. Buča, and D. Jaksch, New Journal of Physics 22, 085007 (2020).
Buča and Jaksch (2019) B. Buča and D. Jaksch, Phys. Rev. Lett. 123, 260401 (2019).
Zhu et al. (2019) B. Zhu, J. Marino, N. Y. Yao, M. D. Lukin, and E. A. Demler, New Journal of Physics 21, 073028 (2019).
Lledó et al. (2019) C. Lledó, T. K. Mavrogordatos, and M. H. Szymańska, Phys. Rev. B 100, 054303 (2019).
Seibold et al. (2020) K. Seibold, R. Rota, and V. Savona, Phys. Rev. A 101, 033839 (2020).
Prazeres et al. (2021) L. F. d. Prazeres, L. d. S. Souza, and F. Iemini, Phys. Rev. B 103, 184308 (2021).
Piccitto et al. (2021) G. Piccitto, M. Wauters, F. Nori, and N. Shammah, Phys. Rev. B 104, 014307 (2021).
Carollo and Lesanovsky (2022) F. Carollo and I. Lesanovsky, Phys. Rev. A 105, L040202 (2022).
Hajdušek et al. (2022) M. Hajdušek, P. Solanki, R. Fazio, and S. Vinjanampathy, Phys. Rev. Lett. 128, 080603 (2022).
Krishna et al. (2023a) M. Krishna, P. Solanki, M. Hajdušek, and S. Vinjanampathy, Phys. Rev. Lett. 130, 150401 (2023a).
Hurtado-Gutiérrez et al. (2020) R. Hurtado-Gutiérrez, F. Carollo, C. Pérez-Espigares, and P. I. Hurtado, Phys. Rev. Lett. 125, 160601 (2020).
Cosme et al. (2023) J. G. Cosme, J. Skulte, and L. Mathey, Phys. Rev. B 108, 024302 (2023).
Buonaiuto et al. (2021) G. Buonaiuto, F. Carollo, B. Olmos, and I. Lesanovsky, Phys. Rev. Lett. 127, 133601 (2021).
Cabot et al. (2023a) A. Cabot, G. Giorgi, and R. Zambrini, (2023a), arXiv:2308.12080 [quant-ph] .
Carollo et al. (2024) F. Carollo, I. Lesanovsky, M. Antezza, and G. D. Chiara, Quantum Science and Technology 9, 035024 (2024).
Passarelli et al. (2022) G. Passarelli, P. Lucignano, R. Fazio, and A. Russomanno, Phys. Rev. B 106, 224308 (2022).
Cabot et al. (2023b) A. Cabot, L. S. Muhle, F. Carollo, and I. Lesanovsky, Phys. Rev. A 108, L041303 (2023b).
Cabot et al. (2024) A. Cabot, F. Carollo, and I. Lesanovsky, Phys. Rev. Lett. 132, 050801 (2024).
Montenegro et al. (2023) V. Montenegro, M. G. Genoni, A. Bayat, and M. G. Paris, Communications Physics 6, 304 (2023).
Carollo et al. (2020) F. Carollo, K. Brandner, and I. Lesanovsky, Phys. Rev. Lett. 125, 240602 (2020).
Keßler et al. (2021) H. Keßler, P. Kongkhambut, C. Georges, L. Mathey, J. G. Cosme, and A. Hemmerich, Phys. Rev. Lett. 127, 043602 (2021).
Kongkhambut et al. (2022) P. Kongkhambut, J. Skulte, L. Mathey, J. G. Cosme, A. Hemmerich, and H. Keßler, Science 377, 670 (2022).
Buča and Prosen (2012) B. Buča and T. Prosen, New Journal of Physics 14, 073007 (2012).
Albert and Jiang (2014a) V. V. Albert and L. Jiang, Phys. Rev. A 89, 022118 (2014a).
Dicke (1954) R. H. Dicke, Phys. Rev. 93, 99 (1954).
Andreev et al. (1980) A. V. Andreev, V. I. Emel’yanov, and Y. A. Il’inskiĭ, Soviet Physics Uspekhi 23, 493 (1980).
Wang and Hioe (1973) Y. K. Wang and F. Hioe, Phys. Rev. A 7, 831 (1973).
Scully and Svidzinsky (2009) M. O. Scully and A. A. Svidzinsky, Science 325, 1510 (2009).
Ferioli et al. (2023) G. Ferioli, A. Glicenstein, I. Ferrier-Barbut, and A. Browaeys, Nature Physics 19, 1345 (2023).
Bhatti et al. (2018) D. Bhatti, R. Schneider, S. Oppel, and J. von Zanthier, Phys. Rev. Lett. 120, 113603 (2018).
Bojer and von Zanthier (2022) M. Bojer and J. von Zanthier, Phys. Rev. A 106, 053712 (2022).
Gegg et al. (2018) M. Gegg, A. Carmele, A. Knorr, and M. Richter, New Journal of Physics 20, 013006 (2018).
Wiegner et al. (2011) R. Wiegner, J. von Zanthier, and G. S. Agarwal, Phys. Rev. A 84, 023805 (2011).
Latune et al. (2019) C. L. Latune, I. Sinayskiy, and F. Petruccione, Phys. Rev. Res. 1, 033192 (2019).
Latune et al. (2020) C. L. Latune, I. Sinayskiy, and F. Petruccione, New Journal of Physics 22, 083049 (2020).
Yadin et al. (2023) B. Yadin, B. Morris, and K. Brandner, Phys. Rev. Res. 5, 033018 (2023).
Albert and Jiang (2014b) V. V. Albert and L. Jiang, Phys. Rev. A 89, 022118 (2014b).
Krishna et al. (2023b) M. Krishna, P. Solanki, and S. Vinjanampathy, Journal of the Indian Institute of Science 103, 513 (2023b).
Kuramoto and Battogtokh (2002) Y. Kuramoto and D. Battogtokh, Nonlinear Phenom. Complex Syst. 5, 380 (2002).
Abrams and Strogatz (2004) D. M. Abrams and S. H. Strogatz, Phys. Rev. Lett. 93, 174102 (2004).
Yeldesbay et al. (2014) A. Yeldesbay, A. Pikovsky, and M. Rosenblum, Phys. Rev. Lett. 112, 144103 (2014).
Chandrasekar et al. (2015) V. Chandrasekar, R. Suresh, D. Senthilkumar, and M. Lakshmanan, Europhysics Letters 111, 60008 (2015).
Ujjwal et al. (2017a) S. R. Ujjwal, N. Punetha, A. Prasad, and R. Ramaswamy, Phys. Rev. E 95, 032203 (2017a).
Haugland (2021) S. W. Haugland, Journal of Physics: Complexity 2, 032001 (2021).
Manzano and Hurtado (2018) D. Manzano and P. Hurtado, Advances in Physics 67, 1 (2018).
Mukherjee et al. (2024) A. Mukherjee, Y. Ibrahim, M. Hajdušek, and S. Vinjanampathy, Phys. Rev. A 110, 012220 (2024).
(62) See Supplemental Material, which includes the mean-field equations for coupled continuous time crystals as well as additional details on the formation of chimera states and cluster synchronization. It includes Refs. [13,24].
Sakurai and Commins (1995) J. J. Sakurai and E. D. Commins, Modern quantum mechanics, revised edition (American Association of Physics Teachers, 1995).
Thingna and Manzano (2021) J. Thingna and D. Manzano, Chaos: An Interdisciplinary Journal of Nonlinear Science 31, 073114 (2021).
Georgi (2000) H. Georgi, Lie algebras in particle physics: from isospin to unified theories (Taylor & Francis, 2000).
Ramadevi and Dubey (2019) P. Ramadevi and V. Dubey, Group Theory for Physicists: With Applications (Cambridge University Press, 2019).
Zakharova (2020) A. Zakharova, Chimera patterns in networks (Springer, 2020).
Bastidas et al. (2015) V. M. Bastidas, I. Omelchenko, A. Zakharova, E. Schöll, and T. Brandes, Phys. Rev. E 92, 062924 (2015).
Ujjwal et al. (2017b) S. R. Ujjwal, N. Punetha, A. Prasad, and R. Ramaswamy, Phys. Rev. E 95, 032203 (2017b).
Datseris and Parlitz (2022) G. Datseris and U. Parlitz, Nonlinear dynamics: a concise introduction interlaced with code (Springer Nature, 2022).
Skokos et al. (2016) C. H. Skokos, G. A. Gottwald, and J. Laskar, Chaos detection and predictability, Vol. 1 (Springer, 2016).
Johansson et al. (2012) J. Johansson, P. Nation, and F. Nori, Computer Physics Communications 183, 1760 (2012).
Johansson et al. (2013) J. Johansson, P. Nation, and F. Nori, Computer Physics Communications 184, 1234 (2013).
Shammah et al. (2018) N. Shammah, S. Ahmed, N. Lambert, S. De Liberato, and F. Nori, Phys. Rev. A 98, 063815 (2018).
Iemini et al. (2024) F. Iemini, D. Chang, and J. Marino, Phys. Rev. A 109, 032204 (2024).

Appendix A Mean-field and fixed-point analysis of a single continuous time crystal

In this section, we extend the mean-field analysis of continuous time crystals (CTCs) outside the symmetric subspace discussed in the main text. The mean-field dynamics of the system described by Eq. (1) of the main text is governed by the following set of coupled differential equations

$\displaystyle\dot{m}_{x}$	$\displaystyle=\kappa m_{x}m_{z},$
$\displaystyle\dot{m}_{y}$	$\displaystyle=-\Omega m_{z}+\kappa m_{y}m_{z},$	(5)
$\displaystyle\dot{m}_{z}$	$\displaystyle=\Omega m_{y}-\kappa[m_{x}^{2}+m_{y}^{2}],$

where $m_{\alpha}=\langle S_{\alpha}\rangle/S$ are rescaled spin operators with $\alpha\in\{x,y,z\}$ . For the symmetric Dicke subspace, $m$ is constrained to take the value $1$ , which led to the breaking of time-translational symmetry at $\Omega/\kappa=1$ in earlier studies Iemini et al. (2018); Hajdušek et al. (2022). We will now consider the more general case $0\leq m<1$ . The fixed points, $(m_{x}^{*},m_{y}^{*},m_{z}^{*})$ , of the coupled differential Eqs. (B) for this case are given as $\mathcal{M}_{1}=(\pm m\sqrt{1-(m\kappa/\Omega)^{2}},m^{2}\kappa/\Omega,0)$ and $\mathcal{M}_{2}=(0,\Omega/\kappa,\pm\sqrt{m^{2}-(\Omega/\kappa)^{2}})$ . The stability of these fixed points is given by the eigenvalues of the Jacobian $\mathcal{J}$ , whose matrix elements are defined as $\mathcal{J}_{\alpha\beta}=d\dot{m}_{\alpha}/dm_{\beta}$ . The eigenvalues of $\mathcal{J}$ corresponding to fixed points $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ are $\{\Lambda_{1}\}=\{0,\pm\kappa\sqrt{m^{2}-(\Omega/\kappa)^{2}}\}$ and $\{\Lambda_{2}\}=\{0,\pm\kappa\sqrt{m^{2}-(\Omega/\kappa)^{2}},\pm\kappa\sqrt{m^{2}-(\Omega/\kappa)^{2}}\}$ respectively, where both the upper or both the lower signs have to be chosen together. Note that the dynamics of the system are now dependent on both the system parameters $\Omega/\kappa$ and rescaled total magnetization $m$ of the given spin-subspace. For $\Omega/\kappa<m$ , $\mathcal{M}_{2}$ is the only physical state of the system since $\mathcal{M}_{1}$ has unphysical imaginary values for the expectation of spin components. The corresponding eigenvalues $\{\Lambda_{2}\}$ for the fixed point $\mathcal{M}_{2}$ are purely real, and hence it is a saddle fixed point of the system. Therefore, the system relaxes to a time-invariant state given by the fixed point $\mathcal{M}_{2}$ for initial states having $m>\Omega/\kappa$ , and no time-crystal behavior is observed. For states with $m<\Omega/\kappa$ , $\mathcal{M}_{1}$ is the only physical fixed point of the system. The eigenvalues $\{\Lambda_{1}\}$ for the fixed point $\mathcal{M}_{1}$ are purely imaginary, and hence, the corresponding fixed point is a center. Therefore, the system oscillates in time for initial states having $m<\Omega/\kappa$ , and time-translational symmetry is broken even for $\Omega/\kappa<1$ .

Using this fixed-point analysis, we can characterize the phase transition of the system as described by Fig. 1(a) of the main text. The system exhibits multistability for $\Omega/\kappa<1$ . The two fixed points $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ have different basins of attraction: depending upon the initial state, the system can exhibit two different phases, namely, a melted and a time-crystal phase. For the Dicke subspace, this bistability does not exist, as depicted in Fig. 1(a) of the main text. The Dicke subspace corresponds to initial states having $m=1$ , which shows that for $\Omega/\kappa<1$ , the system always settles down to a saddle fixed point (dashed red line), and for $\Omega/\kappa>1$ , it shows oscillatory behavior (solid red line). However, for $m<1$ , both time-translation invariant and symmetry-broken phases coexist and are separated by the blue line corresponding to $m=\Omega/\kappa$ . Initial states having $m>\Omega/\kappa$ on the left of $m=\Omega/\kappa$ (blue) line settle down to the saddle fixed point $\mathcal{M}_{2}$ . The parameter region right to the $m=\Omega/\kappa$ (blue) line corresponds to $m<\Omega/\kappa$ and gives stable oscillations around the fixed point $\mathcal{M}_{1}$ leading to time-translational symmetry breaking.

To illustrate the phase-space portrait of the fixed points discussed above, we consider the following coordinate transformation $P=\tan^{-1}\left[m_{y}/m_{x}\right],~{}Q=m_{z}$ Iemini et al. (2018). In terms of the new co-ordinates $P$ and $Q$ , the fixed points $\{P^{*},Q^{*}\}$ are given as $\mathcal{M}^{P,Q}_{1}=(\tan^{-1}[\pm 1/\sqrt{(\Omega/m\kappa)^{2}-1}],0)$ and $\mathcal{M}^{P,Q}_{2}=(\pm\pi/2,\pm\sqrt{m^{2}-(\Omega/\kappa)^{2}})$ . The phase-space portrait of the CTC for different $m$ values is depicted in Fig. S1(a) for $\Omega/\kappa=0.75$ . Initial states having $m<\Omega/\kappa$ keep oscillating in the vicinity of the fixed point $\mathcal{M}^{P,Q}_{1}$ . As soon as we choose an initial state with $m>\Omega/\kappa$ , the corresponding trajectory in Fig. S1(a) settles down to the saddle fixed point $\mathcal{M}^{P,Q}_{2}$ . For $\Omega/\kappa>1$ , there exists no physical state of the system that can have $m>\Omega/\kappa$ since $m\leq 1$ for the rescaled spin operators. The stability of the oscillations around the fixed point $\mathcal{M}_{1}^{P,Q}$ can be verified by observing Fig. S1(a), where each limit cycle trajectory consists of around 4000 cycles of oscillations. These trajectories are obtained by the time evolution of the exact mean-field equations (B) (including non-linear terms) and show no deviation even at such long times.

Figure S1(b) depicts the time evolution of the master equation (Eq. (1) of the main text) with the initial state belonging to the $m=0.1$ subspace, with $\{m_{x},m_{y},m_{z}\}=\{0,0,0.1\}$ for $\Omega/\kappa=0.9$ . As the system size increases, the oscillation in $m_{z}$ becomes persistent and approaches the mean-field limit. Any initial state belonging to subspaces with $m>0.9$ will decay and approach a time-invariant steady state irrespective of system size. Thus, the exact treatment confirms the mean-field analysis for $\Omega/\kappa<1$ .

Appendix B Mean-field dynamics of coupled continuous time crystals

Here we discuss the mean-field dynamics of the coherently coupled CTCs described by Eq. (3) of the main text. These equations are formulated in terms of the rescaled spin operators $m^{\alpha}_{\beta}=\langle\hat{S}^{\alpha}_{\beta}\rangle/S$ with $\alpha:\{x,y,z\}$ and $\beta:\{A,B\}$ , and are defined as follows,

$\displaystyle\frac{dm^{x}_{\alpha}}{dt}$	$\displaystyle=\kappa m^{x}_{\alpha}m^{z}_{\alpha}+\frac{\Gamma}{n}m^{z}_{\alpha}\sum_{\beta\neq\alpha}m^{y}_{\beta},$
$\displaystyle\frac{dm^{y}_{\alpha}}{dt}$	$\displaystyle=-\Omega_{\alpha}m^{z}_{\alpha}+\kappa m^{y}_{\alpha}m^{z}_{\alpha}-\frac{\Gamma}{n}m^{z}_{\alpha}\sum_{\beta\neq\alpha}m^{x}_{\beta},$	(6)
$\displaystyle\frac{dm^{z}_{\alpha}}{dt}$	$\displaystyle=\Omega_{\alpha}m^{y}_{\alpha}-\kappa\left[(m^{x}_{\alpha})^{2}+(m^{y}_{\alpha})^{2}\right]+\frac{\Gamma}{n}\left[m^{y}_{\alpha}\sum_{\beta\neq\alpha}m^{x}_{\beta}-m^{x}_{\alpha}\sum_{\beta\neq\alpha}m^{y}_{\beta}\right].$

We numerically solve these coupled differential equations to study the synchronization of coupled CTCs.

Appendix C Chimera state and synchronization via seeding

The initial state leading to the chimera state shown in Fig. 3 of the main text is chosen such that the absolute values of the corresponding Jacobian eigenvalues $\{|\Lambda_{1}^{\alpha}|;\alpha=1,2,\ldots 20\}$ for the fixed point $\mathcal{M}_{1}$ are uniformly sampled in the window from $0.2$ to $0.25$ for CTCs in ensemble $\varepsilon_{1}$ (represented by black color) and from $0.75$ to $0.85$ for CTCs in ensemble $\varepsilon_{2}$ (represented by red color). For the given initial state, we investigate the effect of increasing the coupling strength between the all-to-all coupled CTCs on the synchronization properties. For coupling strength $\Gamma/\kappa=0.5$ , all the CTCs in ensemble $\varepsilon_{1}$ are synchronized while CTCs in ensemble $\varepsilon_{2}$ only show partial synchronization, as shown in Figs. S2(a,d). Increasing $\Gamma$ we see more complex dynamics (not shown here), but the system continues to exhibit a chimera state. At $\Gamma/\kappa\approx 0.605$ , a transition occurs: the strength of the coupling leads to the melting of the CTCs of ensemble $\varepsilon_{1}$ as depicted in Figs. S2(b,e). Further increasing the coupling strength leads to the synchronization of the CTCs in ensemble $\varepsilon_{2}$ which seeds oscillations to ensemble $\varepsilon_{1}$ . As a result, at coupling strength $\Gamma/\kappa=1.2$ , all the CTCs oscillate with a common frequency as shown in Figs. S2(c,f). Such a synchronization via seeding is consistent with what was observed in Hajdušek et al. (2022).

Appendix D Cluster synchronization

We will now discuss the effect of small changes in the initial state on the synchronization of CTCs. As an example, we consider a change in the initial states for CTCs belonging to ensemble $\varepsilon_{2}$ such that $\{|\Lambda_{1}^{\alpha}|\}$ are uniformly sampled in the interval from $0.75$ to $0.8$ . The choice of initial states for the uncoupled CTCs is shown in Figs. S3(a,d). Now we study the effect of choosing this distribution on the synchronization of CTCs and we only consider all-to-all coupling in this case. This small change results in the emergence of a new synchronization phase, termed cluster synchronization as shown in Figs. S3(b,e). As the name suggests, there exist two clusters of synchronized CTCs with different frequencies. Further increasing the coupling strength to $\Gamma/\kappa=1.25$ results in complete synchronization and these two clusters oscillate at the same frequency as shown in Figs. S3(c,f). These findings demonstrate that different choices of initial states can result in different regimes of synchronization.

Appendix E Characterization of different phases using Pearson coefficient and maximum Lyapunov exponent

Here, we provide a detailed characterization of the various synchronization regimes discussed in the main text. We consider an ensemble of $n=100$ CTCs that are divided into two subgroups $\varepsilon_{1,2}$ containing $50$ CTCs each. The frequencies of the CTCs in each subgroup $\varepsilon_{1,2}$ are sampled from Gaussian distributions with distinct mean frequencies $\bar{\omega}_{1,2}$ and standard deviations $\sigma_{1,2}$ . As in Fig. 4(e,f) of the main text, we use $\sigma_{1}=\sigma_{2}=0.1$ , and vary $\Delta=\bar{\omega}_{1}-\bar{\omega}_{2}$ and $\Gamma$ to explore six regimes of synchronization marked by different symbols in Figs. S4 and S5. We study the average value of the Pearson coefficient $\bar{\mathcal{C}}$ and the maximum Lyapunov exponent $\lambda_{L}$ defined in the main text to characterize these regimes.

Partial Synchronization —

We first consider the three different cases marked by the symbols ‘ $\scriptstyle\bigstar$ ’( $\Gamma=0.8,\Delta=0.7$ ), ‘’( $\Gamma=0.35,\Delta=0.6$ ) and ‘ $\bullet$ ’( $\Gamma=0.8,\Delta=0.6$ ) in Figs. S4 and S5. These three cases exhibit $\bar{\mathcal{C}}<1$ as depicted in Fig. S4(a). However, they possess distinct dynamical properties, as evidenced by the maximum Lyapunov exponent shown in Fig. S5(a). Therefore, we examine the Pearson correlation matrix $\mathcal{C}$ along with the Lyapunov exponent to distinguish between the following different cases of partial synchronization:

•

Coexistence of chimera and cluster synchronization: For $\Gamma=0.8$ and $\Delta=0.7$ (denoted by ‘ $\scriptstyle\bigstar$ ’), the CTCs are in a stable limit-cycle state since $\lambda_{L}=0$ (as shown in Fig. S5(a)). The correlation matrix depicted in Fig. S4(b) reveals two prominent clusters of synchronized CTCs, while the remainder of the matrix is diagonal, indicating unsynchronized CTCs. A similar behavior can be observed in Fig. S5(b), where two synchronized subsets of CTCs have equal frequencies, while the rest of the unsynchronized CTCs oscillate with different frequencies. Thus, cluster synchronization and a chimera state can coexist.
•

Chimera state: For $\Gamma=0.35$ and $\Delta=0.6$ (denoted by ‘’), there is only one prominent cluster of synchronized CTCs while the other CTCs remain unsynchronized. This can be understood from Fig. S4(c), where the CTCs of subgroup $\varepsilon_{1}$ show a high amount of correlation while the remainder of the CTCs are almost uncorrelated. Therefore, the CTCs of subgroup $\varepsilon_{1}$ have a common frequency, while the other CTCs oscillate with different frequencies as shown in Fig. S5(c). Thus, the system exhibits a chimera state.
•

Chimera state with partial oscillation death: We next discuss a case in which the Pearson correlation matrix looks qualitatively similar to that of a chimera state (see Fig. S4(c)). This can be observed for $\Gamma=0.8,\Delta=0.6$ (denoted by ‘ $\bullet$ ’), as shown in Fig. S4(d). Note that $\lambda_{L}<0$ , see Fig. S5(a). Thus, some of the CTCs cease to oscillate leading to partial oscillation death. The same conclusion follows from Fig. S5(d) where all the CTCs of subgroup $\varepsilon_{1}$ have zero frequency and hence melt down. The steady state of the melted CTCs corresponds to highly correlated fixed points, and the correlation matrix exhibits a block diagonal structure similar to what was observed for chimera states.

As we have seen, $\bar{\mathcal{C}}$ alone does not provide sufficient information to distinguish different cases of partial synchronizations. Therefore, we investigated the correlation matrix and the maximum Lyapunov exponent to differentiate between these cases.

Oscillation death —

We now discuss the case of complete oscillation death, where all the CTCs melt down and stop oscillating. This can be observed for $\Gamma=0.8$ and $\Delta=0.3$ (represented by ‘ ’), where $\lambda_{L}\ll 0$ , as shown in Fig. S5(a). Indeed, from Fig. S5(e), we can observe that all the CTCs have zero frequency. Therefore the system settles down to a fixed point corresponding to complete oscillation death. The Pearson coefficient matrix shows a large amount of correlation between all the CTCs.

Mean-field chaos —

By further reducing the detuning to $\Delta=0.2$ while keeping $\Gamma=0.8$ (represented by ‘’), the mean-field dynamics becomes chaotic, indicated by a positive value of $\lambda_{L}$ in Fig. S5(a). This leads to a lower value of $\bar{\mathcal{C}}$ as shown in Fig. S4(a). Although the correlation matrix exhibits a reduced amount of correlations between the CTCs due to the chaotic nature of the mean-field dynamics, the finite value of its off-diagonal elements indicates the presence of chaotic synchronization.

Complete synchronization —

Finally, for small detuning $\Delta=0.1$ and $\Gamma=0.8$ (represented by ‘ $\blacktriangle$ ’), the system becomes completely synchronized such that $\bar{\mathcal{C}}\approx 1$ and $\lambda_{L}\approx 0$ . This regime is indicated by Fig. S4(g), where all the off-diagonal elements of $\mathcal{C}$ are non-zero, representing a large amount of correlation between all CTCs. The complete synchronization can also be observed from Fig. S5(g), where all CTCs oscillate with a common frequency.