Simulating complex quantum networks with time crystals

TL;DR

This paper introduces a graph-theoretic approach to analyze and simulate complex quantum networks using time crystals, revealing emergent phenomena like scale-free networks during phase transitions.

Contribution

It presents a novel method to characterize time crystals through graph analysis, enabling simulation of complex quantum networks with new insights into their melting processes.

Findings

Graph analysis visualizes time-crystalline order.

Melting process shows emergent scale-free network behavior.

Time crystals can simulate complex quantum networks.

Abstract

Crystals arise as the result of the breaking of a spatial translation symmetry. Similarly, translation symmetries can also be broken in time so that discrete time crystals appear. Here, we introduce a method to describe, characterize, and explore the physical phenomena related to this phase of matter using tools from graph theory. The analysis of the graphs allows to visualizing time-crystalline order and to analyze features of the quantum system. For example, we explore in detail the melting process of a minimal model of a period-2 discrete time crystal and describe it in terms of the evolution of the associated graph structure. We show that during the melting process, the network evolution exhibits an emergent preferential attachment mechanism, directly associated with the existence of scale-free networks. Thus, our strategy allows us to propose a previously unexplored far-reaching application of time crystals as a quantum simulator of complex quantum networks.

Equations

$$\hat{H}(t)=\{\begin{array}{ll}{\hat{H}}_1\equiv \hslash g\left(1-ϵ\right){\sum }_l{\sigma }_l^x& 0<t<T_1\\ {\hat{H}}_2\equiv \hslash {\sum }_{lm}{J}_{lm}^z{\sigma }_l^z{\sigma }_m^z+\hslash {\sum }_l{B}_l^z{\sigma }_l^z& T_1<t<T.\end{array}$$

$${\hat{\mathcal{F}}}_{ϵ}={\hat{U}}_{ϵ}(T)=\exp \left(-\frac{\mathrm{i}}{\hslash }{\hat{H}}_2{T}_2\right)\exp \left(-\frac{\mathrm{i}}{\hslash }{\hat{H}}_1{T}_1\right).$$

$${\hat{H}}_{ϵ,T}^{\text{eff}}=\sum _i{\mathcal{E}}_i|i⟩⟨i|+\sum _{i,j}{K}_{ij}|i⟩⟨j|=\sum _s{\lambda }_s|{\Phi }_s⟩⟨{\Phi }_s|,$$

$$|{\mathcal{E}}_j-{\mathcal{E}}_i|<|K_{ij}|$$

$${\hat{\mathcal{F}}}_{ϵ=0}^2=\exp \left(-2i{T}_2\sum _{lm}{J}_{lm}^z{\sigma }_l^z{\sigma }_m^z\right).$$

$${\mathcal{H}}_{ϵ=0,2T}^{\text{eff}}({\theta }_1,{\theta }_2,\dots ,{\theta }_n)=\frac{\hslash }{2}\sum _{lm}{J}_{lm}^z\cos {\theta }_l\cos {\theta }_m,$$

$${\hat{H}}_{ϵ,2T}^{\text{eff}}$$

$$+\sin (B_{l}2T_2){\sigma }_l^y].$$

$$\hat{H}(t)=\{\begin{array}{ll}{\hat{H}}_1\equiv \hslash g\left(1-ϵ\right)\left(|i⟩⟨j|+|j⟩⟨i|\right)& 0<t<T_1\\ {\hat{H}}_2\equiv E_i|i⟩⟨i|+E_j|j⟩⟨j|& T_1<t<T_2.\end{array}$$

$${\hat{H}}_{ϵ,T}^{\text{eff}}={\mathcal{E}}_i|i⟩⟨i|+{\mathcal{E}}_j|j⟩⟨j|+K_{ij}|i⟩⟨j|+K_{ji}|j⟩⟨i|=$$

$${\lambda }_{s_1}|{\Phi }_{s_1}⟩⟨{\Phi }_{s_1}|+{\lambda }_{s_2}|{\Phi }_{s_2}⟩⟨{\Phi }_{s_2}|.$$

$${\hat{G}}_{\text{F}}(z)=\sum _{l=0}^{\infty }{\hat{\mathcal{F}}}^n{e}^{\mathrm{i}nz}=\frac{1}{1-\hat{\mathcal{F}}{e}^{\mathrm{i}z}}.$$

$${\hat{G}}_{\text{F}}(z)$$

$$\approx \frac{\mathrm{i}}{T(z-{\hat{H}}_{ϵ,T}^{\text{eff}})}+B_1+\cdots$$

$$\frac{1}{z-{\hat{H}}_{ϵ,T}^{\text{eff}}}=\frac{1}{z-{\hat{H}}_0}+\left(\frac{1}{z-{\hat{H}}_0}\right)\hat{V}\left(\frac{1}{z-{\hat{H}}_0}\right)+\cdots .$$

$$⟨i|{\hat{G}}_{\text{F}}(z)|j⟩\approx \frac{\mathrm{i}{\delta }_{i,j}}{T(z-{\mathcal{E}}_i)}+\frac{\mathrm{i}{K}_{ij}}{T(z-{\mathcal{E}}_i)(z-{\mathcal{E}}_j)}+B_1{\delta }_{i,j}+\cdots ,$$

$$P_{\text{Poisson}}(r)=\frac{2}{(1+r{)}^2}.$$

$$P_{\text{GOE}}(r)=\frac{27}{4}\frac{r+r^2}{(1+r+r^2{)}^{5/2}}$$

$$P_{\text{COE}}(r)=\frac{2}{3}\left\{\left[\frac{\sin \left(\frac{2\pi r}{r+1}\right)}{2\pi r^2}\right]+\frac{1}{(1+r{)}^2}+\left[\frac{\sin \left(\frac{2\pi }{r+1}\right)}{2\pi }\right]\right\}-\frac{2}{3}\left\{\left[\frac{\cos \left(\frac{2\pi }{r+1}\right)}{2\pi r^2}\right]+\left[\frac{\cos \left(\frac{2\pi r}{r+1}\right)}{r(r+1)}\right]\right\}$$

$$M_{i,Z}(mT)=-\frac{1}{n}\sum _{r=1}^n\sum _{j=1}^{2^n}(-1{)}^{j_r}|A_j(mT){|}^2,$$

$$\{M_{i,Z}(0),M_{i,Z}(T),\dots ,M_{i,Z}(NT)\},$$

$${\mathcal{M}}_{i,Z}(k)=\frac{1}{N}\sum _{m=1}^N{e}^{-\frac{\mathrm{i}2\pi k}{N}m}{M}_{i,Z}(mT)=\frac{1}{N}\sum _{m=1}^N{e}^{-\mathrm{i}{\omega }_{k}mT}{M}_{i,Z}(mT),$$

$${\mathbf{V}}_{i,ϵ}=\{|{\mathcal{M}}_{i,Z}(0){|}^2,|{\mathcal{M}}_{i,Z}(1){|}^2\dots ,|{\mathcal{M}}_{i,Z}(N-1){|}^2\}$$

$$F_i^{\text{ps}}(ϵ)=\sqrt{\frac{{\mathbf{V}}_{i,0}\cdot {\mathbf{V}}_{i,ϵ}}{\Vert {\mathbf{V}}_{i,0}\Vert \Vert {\mathbf{V}}_{i,ϵ}\Vert }}.$$

$${\hat{H}}_{ϵ,2T}^{\text{eff}}$$

$${\mathcal{H}}_{0,2T}({\theta }_1,\dots ,{\theta }_n)=\frac{\hslash T_2}{T}\sum _{l,m}{J}_{l,m}^z\cos {\theta }_l\cos {\theta }_m.$$

$${\mathcal{J}}_{i,i}$$

$${\mathcal{J}}_{i,j}$$

$${\hat{H}}_{ϵ,T}^{\text{eff}}=\sum _i{\mathcal{E}}_i|i⟩⟨i|+\sum _{i,j}{K}_{ij}|i⟩⟨j|=\sum _s{\lambda }_s|{\Phi }_s⟩⟨{\Phi }_s|$$

$$\hat{\rho }(nT)=|\psi (nT)⟩⟨\psi (nT)|=\sum _{i,j}{A}_i(nT)A_j^{\ast }(nT)|i⟩⟨j|,$$

$$C_{i,Z}(nT)$$

$$=\frac{1}{2^n}\text{tr}\left[\hat{\rho }(nT)({\sigma }_1^z-(-1{)}^{i_1}\hat{I})({\sigma }_2^z-(-1{)}^{i_2}\hat{I})\cdots ({\sigma }_n^z-(-1{)}^{i_n}\hat{I})\right],$$

$$C_{2^{n-1},Z}(nT)=\frac{1}{2^n}\text{tr}\left[\hat{\rho }(nT)({\sigma }_1^z+\hat{I})({\sigma }_2^z+\hat{I})\cdots ({\sigma }_n^z+\hat{I})\right],$$

Full text

INTRODUCTION

Symmetries are of utmost importance in condensed matter physics and statistical mechanics due to the strong relation between quantum phases of matter and symmetry breaking Shahar1997 ; Higgs1964 ; BCS1957 ; Luke1998 . Among the zoo of symmetry broken phases, discrete time crystals (DTC) play a fundamental role on their own due to the type of symmetry involved Sacha2018 . The time-crystalline phase, analogous to “space” crystals, arises when time (instead of space) translation symmetry is broken. The existence of quantum time crystals was originally proposed by Wilczek Wilczek2012 , and a discrete version of time crystals can be realized periodically-driven quantum systems Sacha2015 ; Else2016 ; Khemani2016 . It is in such cases that the system dynamics exhibits a subharmonic response with respect to the characteristic period of the drive caused by the synchronization in time of the particles of a many-body system Sacha2018 .

The exploration of time crystals is a very active field of research and several experimental realizations with trapped ions Zhang2017 , dipolar spin impurities in diamond Choi2017 , ordered dipolar many-body systems Barrett2018 , ultracold atoms Smits2018 and nuclear spin-1/2 moments in molecules Pal2018 have been achieved. Yet an intuitive and complete insight of the nature of time-crystals and its characterization, as well as a set of proposed applications, is lacking. We here provide new tools based on graph theory and statistical mechanics to fill this gap. We propose a new strategy for the study and understanding of time symmetry-broken phases and their related phenomena Yao2017 ; Sacha2018-2 . As an example, we characterize the time-crystalline order and its melting.

Our approach allows us to identify an instance of a perturbed time-crystal as a novel physical platform where to simulate complex networks whose evolution is governed by preferential attachment mechanisms Barabasi1999 ; Barabasi2002 . This type of networks, far from being regular or random, contains non-trivial topological structures present in many biological, social and technological systems. Small-world and scale-free networks are two of the most popular examples, the later being commonly characterized through power-law degree distributions that can be explained from the presence of a preferential attachment mechanism. The simulation of such networks has wide applicability, ranging from the study and understanding of behaviors present in communication or internet networks Yuan2012 , the development of new algorithms in deep learning Pal2018a , or the analysis of genetic and neural structures in biological systems Dorogovtsev2003 .

**RESULTS

Floquet theory in a nutshell

**The fundamentals of this work rely on the exploration of the dynamics of driven many-body quantum systems described by time-periodic Hamiltonians $\hat{H}(t+T)=\hat{H}(t)$, with $T=2\pi/\omega$ being the period of that drive. Our approach builds from the calculation of the Floquet operator $\hat{\mathcal{F}}=\hat{U}(T)=\mathcal{T}\exp(-\mathrm{i}\int_{0}^{T}\hat{H}(\tau)d\tau/\hbar)$, which is the evolution operator within one period from time $t=0$ to time $t=T$ Haenggi1998 ; Polkovnikov2015 ; Anisimovas2015 with $\mathcal{T}$ being the time ordering operator. In Floquet theory, one wishes to solve the eigenvalue problem $\hat{\mathcal{F}}|\Phi_{s}\rangle=e^{-\mathrm{i}\lambda_{s}T/\hbar}|\Phi_{s}\rangle$. The eigenvectors $|\Phi_{s}\rangle$ are known as Floquet states and $-\hbar\omega/2\leq\lambda_{s}\leq\hbar\omega/2$ are the quasienergies. At discrete times $t_{n}=nT$, the evolution can be understood in terms of an effective “time-independent” Hamiltonian $\hat{H}^{\text{eff}}$ such that $\hat{\mathcal{F}}=e^{-\mathrm{i}\hat{H}^{\text{eff}}T/\hbar}$. This Hamiltonian is very relevant in the context of Floquet engineering as it contains effective interactions that are absent in equilibrium systems, allowing it to be applied to quantum simulation problems Polkovnikov2015 ; Brandes2016 ; Eckardt2017 .

**Discrete time crystals of period 2T

**We now move to a well-known Period-$2$ Discrete Time Crystal (2T-DTC) model Yao2017 ; Zhang2017 consisting of a one dimensional chain with $n$ spin-$1/2$ particles and governed by a time-dependent Hamiltonian with period $T=T_{1}+T_{2}$,

[TABLE]

Here $\{\sigma_{l}^{x},\sigma_{l}^{y},\sigma_{l}^{z}\}$ are the usual Pauli operators on the $l$-th spin, $J_{lm}^{z}\equiv J_{0}/|l-m|^{\alpha}$ is the long-range interaction between spins $l,m$ that takes the form of an approximate power-law decay with a constant exponent $\alpha$, and $B_{l}^{z}\in\left[0,W\right]$ is a random longitudinal field. From the interplay between the driving, interactions and disorder, an MBL phase arises, being a key ingredient to allow the existence of DTC. Such a phase prevents the system from heating up to infinite temperature due to the effect of the drive and thus destroying the time-crystalline phase (see Sec. II of the Supplementary Material for a more detailed description). The parameter $g$ satisfies the condition $2gT_{1}=\pi$ such that $\hat{U}_{1}=\exp{\left(-\mathrm{i}\hat{H}_{1}T_{1}/\hbar\right)}$ becomes a global $\pi$ pulse around the $x$-axis, with a rotation error $\epsilon$ also introduced. The Floquet operator is then given by:

[TABLE]

Crucially, the Floquet operator depends on the error $\epsilon$, which determines how the time crystal will melt. In our work, we choose $\alpha=1.51$, $J_{0}T_{2}=0.06$ with a disorder strength $WT_{2}=\pi$, which are similar values used in the recent experiment Zhang2017 . The key feature of this unitary evolution is that at each period $T$, as shown in Fig. 1.a, the system evolves between the states $|\uparrow\uparrow\cdots\uparrow\rangle$ and $|\downarrow\downarrow\cdots\downarrow\rangle$, when the state is initialized at one of these two. This periodic evolution is robust against moderate errors $\epsilon$ in $\hat{H_{1}}$.

**Floquet graphs

**Let us now introduce the concept of a Floquet graph as the graphical representation of the effective Hamiltonian, $\hat{H}^{\text{eff}}_{\epsilon,T}$ Bastidas2018 , obtained as $\hat{H}^{\text{eff}}_{\epsilon,T}=i\hbar/T\log[\hat{\mathcal{F}}_{\epsilon}]$. By shifting to the graph theory framework, we can provide a visual and intuitive understanding of the system and use its well-established tools to understand physical phenomena. We discuss the application of percolation to unveil the complex dynamics of time-crystals, even though such ideas can be widely applicable to any periodically-driven (Floquet) system.

The effective Hamiltonian, $\hat{H}_{\epsilon,T}^{\text{eff}}$, is defined on a finite and discrete Hilbert space (dimension $N_{\mathcal{H}}$), and can be written as a tight-binding model using the configuration basis states $\{|i\rangle\}$, as

[TABLE]

where $\mathcal{E}_{i}$ is the energy of configuration $|i\rangle$, and $K_{ij}$ is the transition energy between configurations $|i\rangle$ and $|j\rangle$ (see Fig. 1.b). The Floquet states $|\Phi_{s}\rangle$ are linear combinations of the configurations $|i\rangle$. In general, the quasienergies $\lambda_{s}$ are complicated functions of the parameters $\mathcal{E}_{i}$ and $K_{ij}$. In this work, we consider a system with $n$ spin-$1/2$ particles and the configurations $\{|i\rangle\}$ with $i=0,1,\ldots,2^{n}-1$ can be chosen as the $N_{\mathcal{H}}=2^{n}$ elements of the computational basis, i.e. the product states of the eigenstates of $\sigma_{l}^{z}$. In general, any configuration can be written as $|i\rangle=(\sigma^{+}_{1})^{j_{1}}\cdots(\sigma^{+}_{n})^{j_{n}}|\downarrow\downarrow\cdots\downarrow\rangle$, where $i=(j_{1}j_{2}\ldots j_{n})_{2}$ is the binary decomposition of the integer $i$ and $\sigma^{+}_{l}|\downarrow\rangle_{l}=|\uparrow\rangle_{l}$. By using this notation, the states are denoted as $|0\rangle=|\downarrow\downarrow\cdots\downarrow\rangle$ and $|2^{n}-1\rangle=|\uparrow\uparrow\cdots\uparrow\rangle$. The basis configurations have been chosen to satisfy the relation $|2^{n}-1-i\rangle=\sigma^{x}_{1}\sigma^{x}_{2}\cdots\sigma^{x}_{n}|i\rangle$, such that the states $|i\rangle$ and $|2^{n}-1-i\rangle$ are always linked by the global $\pi$ rotation from $\hat{H}_{1}$ in Eq. 1, and will constitute the $2^{n}$ nodes of the graph. In addition, such nodes will be linked according to a percolation rule Roy2019 , where a non-weighted edge is considered to be active between nodes $i$ and $j$ only if the condition

[TABLE]

is satisfied. This rule eliminates the off-resonant transitions, allowing to effectively visualize the significant transitions in the Hamiltonian (see Sec..1 of the Methods for further detail).

In order to clarify the percolation on the graph, we define clusters as subsets of nodes which are connected with a path. The graph is percolated when all nodes belong to a single cluster Potter2015 . We want to stress that the nodes represent many-body states $|i\rangle$ and not physical spin sites $l$, so the graph spans through the full Hilbert –configuration– space (see Fig. 1.c). This percolation rule was proven to be useful to detect many-body localization to thermal phase transitions in undriven spin systems Roy2019 .

**Time crystal graphs

**Under the framework of Floquet graphs, the dynamics of the described $2T$-DTC can now be investigated. The effective Hamiltonian $\hat{H}^{\text{eff}}_{\epsilon,T}$ derived from Eq. (2) can be represented as a graph using the percolation rule from Eq. (4). The structure of the graph when $\epsilon=0$ presents decoupled dimers (two connected nodes), a signature of $2T$-DTCs crystals, as shown in the example presented in Fig. 1.c. We note that this visual characterization of the crystal time order is not unique of $2T$-DTC and can be extrapolated to crystals of different periodicity. For example, $3T$-DTC presents decoupled trimers; $4T$-DTC presents decoupled tetragons, and so on. Our results hence reveal that, in the case of $2T$-DTCs and in the absence of error, there are decoupled two-dimensional subspaces. For example if there is an active link between the configurations $|0\rangle$ and $|2^{n}-1\rangle$ defined above, there is an energy gap $\hbar\pi/T$ between the quasienergies $\lambda_{0}$ and $\lambda_{2^{n}-1}$. The corresponding Floquet states are given by $|\Phi_{0}\rangle\approx(|0\rangle+|2^{n}-1\rangle)/\sqrt{2}$ and $|\Phi_{2^{n}-1}\rangle\approx(|0\rangle-|2^{n}-1\rangle)/\sqrt{2}$, which are maximally entangled GHZ states. The same behaviour is observed for other configurations, $|i\rangle$ and $|2^{n}-1-i\rangle$, as all the paired Floquet states $|\Phi_{i}\rangle$, $|\Phi_{2^{n}-1-i}\rangle$ have an energy difference of $\hbar\pi/T$ Khemani2016 .

**Melting of the $2T$-DTC

**Consider now the case where our system is initialised to $|\psi(0)\rangle=|2^{n}-1\rangle=|\uparrow\uparrow\cdots\uparrow\rangle$ Zhang2017 . This state is a superposition of Floquet states, $|\psi(0)\rangle\approx(|\Phi_{0}\rangle-|\Phi_{2^{n}-1}\rangle)/\sqrt{2}$ and in the graph is represented as a dimer between the nodes $|0\rangle$ and $|2^{n}-1\rangle$. Such a dimer remains decoupled and isolated for modest levels of error (up to a ${\sim}3\%$ error in the rotation added as a non-zero value of $\epsilon$ in $\hat{H}_{1}$, see supplementary video). In Fig. 2 we present the effect of moderate ($\epsilon=0.012$) and large ($\epsilon=0.1$) levels of rotation errors $\epsilon$ for a $2T$-DTC graph with $n=8$ spins. When the error increases there are couplings between the different subspaces and larger clusters appear in the graph (Fig. 2.a). If we keep increasing the error, the crystal will melt and its associated graph will take the form of a single highly –percolated– connected cluster (Fig. 2.d). The presence/no presence of such a dimer indicates whether the system is still in a crystalline phase or, otherwise, has been melted. This demonstrates how our graph strategy, among other things, can track the robustness of the system in a very visual and intuitive manner. Importantly, we note that not all dimers show the same robustness, and some of them survive to higher error values. This, in fact, allows us to identify the most convenient state initialization: the initial state corresponding to one of the configurations of the most robust dimer will make the crystal phase last longer.

**Conserved quantities

**Let us now analyze what factors are determinant in the robustness of each pair of configurations or, alternatively, how they cluster as the crystal melts. With our system represented as a graph, we can now focus on how to interpret its structural properties in terms of conserved quantities. To do so, it is of fundamental importance to note that the system is switching stroboscopically between two configurations, $|i\rangle$ and $|2^{n}-1-i\rangle$. If we prepare the system in a symmetry broken state in terms of its total magnetization and observe the system every two periods ($2T$), the system will effectively remain in a manifold of states with a fixed quasienergy. For no error, the graph will be formed by decoupled single nodes. In that sense, if we sample the dynamics every two periods, the system is in an effective equilibrium state determined by conserved quantities that are defined stroboscopically. Now let us construct the conserved quantities in the absence of error and show how they are destroyed when clusters in the graph are formed (Fig. 3).

To calculate the conserved quantities, let us consider the square of the Floquet operator of Eq. (2), where $\hat{H}^{\text{eff}}_{\epsilon,2T}$ is the effective Hamiltonian now over two periods ($\hat{H}^{\text{eff}}_{\epsilon,2T}=i\hbar/2T\log[\hat{\mathcal{F}}^{2}_{\epsilon}]$). In that stroboscopic framework ($2T$), the Floquet operator for no error can be reduced to the following expression:

[TABLE]

In the absence of error, $\epsilon=0$, the disorder cancels out exactly and the effective Hamiltonian reads $\hat{H}_{\epsilon=0,2T}^{\text{eff}}=\hbar/2\sum_{lm}J_{lm}^{z}\sigma_{l}^{z}\sigma_{m}^{z}$. Despite its apparent simplicity, the effective Hamiltonian contains very relevant information. In particular, it preserves the local magnetization from which it follows the preservation of the total number of domain walls of the basis states $\hat{\mathcal{N}}=\sum_{l}(1-\sigma_{l}^{z}\sigma_{l+1}^{z})$ and parity $\hat{\Pi}=\prod_{l}\sigma^{x}_{l}$ such that $[\hat{H}_{\epsilon=0,2T}^{\text{eff}},\hat{\mathcal{N}}]=[\hat{H}_{\epsilon=0,2T}^{\text{eff}},\hat{\Pi}]=0$. Therefore our conserved quantities are classified according to parity and the total number of domain walls. If we prepare the system in the initial state $|\psi(0)\rangle=|2^{n}-1\rangle$, the system will remain stroboscopically within the symmetry multiplet with period-2 quasienergy $\lambda^{2T}_{1}=\hbar/2\sum_{lm}J_{lm}^{z}$ and zero domains walls $\langle\hat{\mathcal{N}}\rangle=0$. As shown in Fig. 1.c, in the absence of error, all the dimers conserve the number of domain walls (represented as a colormap) and their quasienergies. The aforementioned initial state, however, breaks the spatial parity symmetry $\hat{\Pi}$, which leads to a non-zero value of the correlation function $\langle\sigma^{z}_{l}(\tau)\rangle=\langle\psi(0)|\sigma^{z}_{l}(\tau)\sigma^{z}_{l}(0)|\psi(0)\rangle\neq 0$, as observed in the experiment Zhang2017 .

**Breaking symmetries

**For zero error, the $2^{n}$ quasienergies are local integrals of motion and the Hilbert space can be classified by using the symmetries that preserve the aforementioned conserved quantities. In turn, the effective Hamiltonian $\hat{H}^{\text{eff}}_{\epsilon,2T}$ is block diagonal (see Fig. 3.b) with $n$ sub-blocks given by the conservation of the number of domain walls. In the absence of error $(\epsilon=0)$, and after $2T$, one can find a classical quasienergy surface in phase space associated to the effective Hamiltonian $\hat{H}_{\epsilon=0,2T}^{\text{eff}}$. The idea is to represent the spins as vectors in a three-dimensional space. By using polar coordinates, one can derive the classical Hamiltonian in a $n$-dimensional phase space

[TABLE]

where $\theta_{l}$ is the polar angle. From this one can show that for zero error $\epsilon=0$ all configurations are stable fixed points. The presence of the error brings the system far from the classical limit and quantum effects become dominant. In such scenario the configurations are not fixed points of the dynamics anymore and they become unstable (see Sec. III and IV of the Supplementary Material for more detail). To investigate this, we use the Baker-Campbell-Hausdorff formula to obtain, up to a first order approximation, the effective Hamiltonian for a non-zero value of $\epsilon$:

[TABLE]

This allows us to understand how $\epsilon$ destroys the symmetries by coupling different symmetry multiplets. This can be observed from Fig. 3.b-e, as off-diagonal entries start populating the effective Hamiltonian matrix due to the error. This situation resembles the Kolmogorov-Arnold-Moser (KAM) theory in classical mechanics, which describes how invariant tori are broken in phase space under the effect of perturbations Lazutkin1993 , with our perturbation being $\epsilon$. First of all, the transverse field term $-\hbar g\epsilon/4\sum_{l}\sin{(B_{l}T)}\sigma_{l}^{y}$ breaks the parity symmetry. Most importantly, the transverse field term $-\hbar g\epsilon/4\sum_{l}(\cos{(B_{l}T)}+1)\sigma_{l}^{x}$, while preserving parity, melts the ferromagnetic order and breaks the conservation of the number of domain walls. One of the most interesting aspects of our work is the striking similarity with KAM theory in classical systems: a small perturbation can destroy integrals of motion Lazutkin1993 . The most fragile integrals of motion are periodic unstable orbits in phase space as discussed in the context of Eq. (6). This explains the robustness of the time crystal and how clusters form in the graph, as we increase the error. The time crystal is robust against small errors because the last integral of motion that is destroyed by the perturbation is the configuration with the lowest number of domain walls: this is, $|\Psi(0)\rangle=|2^{n}-1\rangle$ (see darkest blue dimer on the top-right corner of Fig. 2.a). From now on we will use this model as an example to illustrate the applicability of our method.

**Preferential attachment

**Let us now put the focus on the network topology for moderate levels of rotation error ($\epsilon\sim 0.012$). In the previous section we have seen that the quantum terms that appear in Eq. ( Simulating complex quantum networks with time crystals) break the conservation of number of domain walls. The error lifts the degeneracies present in the spectrum and new transitions between close states appear. In terms of the graph, the nodes with the same or similar number of domain walls connect following a preferential attachment mechanism. As most nodes in the network have $n/2$ or $n/2-1$ ($n$: even) number of domain walls, and they have very close quasienergies, these nodes easily connect to each other with a small value of error. This leads to the appearance of the large degree hub nodes as well as the heavy tailed degree distributions shown in Fig. 2.c. The tail of these distributions can be fit to a power law distribution, which is the characteristic of scale-free networks. Following Clauset et al.’s recommendations clauset2009power , we further test the goodness of fit of the power law against the lognormal distribution for each of these distributions, and observe that the power law favours over the lognormal distribution (see Sec. .2 of the Methods for further detail).

In addition, exploring the behaviour of the average degree of the graph obtained from $\hat{H}^{\text{eff}}_{\epsilon,T}$ further captures the existence of the preferential attachment mechanism and the non-trivial melting process of the $2T$-DTC. From Fig. 4, we can confirm that the nodes with 5 or 6 domain walls (where the $n=8$ case is considered here) tend to acquire neighbours preferentially and have more degree already from a small error $\epsilon<0.02$. This is the region where we can observe the scale-free behaviour emerging. When the error is increased, the preferential behaviour becomes weaker, and eventually most of the nodes tend to have large degree, closer to the structure of random networks, which indicates the crystal has melted.

This behaviour not only has consequences for the dynamics of the time crystals but can also be used to perform quantum simulation of networks with exotic properties such as scale-free-like ones Barabasi2002 . Note that the scale-free-like networks shown in Fig. 2 are only observed in the effective Hamiltonian obtained from single period $\hat{H}^{\text{eff}}_{\epsilon,T}$ and not in the one from two periods $\hat{H}^{\text{eff}}_{\epsilon,2T}$, which shows Poisson distributions, typical for random networks Migdal2013 ; Barabasi2002 . We observe that a small error has little effect to the effective Hamiltonian $\hat{H}^{\text{eff}}_{\epsilon,2T}$ (as shown in Fig. 3.b)). This can be understood from the fact that each $\pi$-rotation changes the direction of the rotation around the z-axis caused by the $B^{z}$ term. Therefore at time $2T$ the effect of the $B^{z}$ term is effectively cancelled, limiting the number of transitions happening. On the other hand, this is not the case at time $T$, and off-diagonal terms appear in the matrix form of $\hat{H}^{\text{eff}}_{\epsilon,T}$ as shown in Fig. 2.b). This indicates that transitions to a larger fraction of states can occur, which in turn yields a scale-free-like network for moderate values of $\epsilon$.

**CONCLUSION

**Our approach of representing discrete time crystals in terms of graphs is key to understand the structure of such an exotic phase of matter and to envisage prospective applications. We explicitly show that by translating a $2T$-DTC into a graph theory language we can study in detail how the crystal order disappears as an increasing error melts the time crystal. Among others, it allowed us to identify the crucial role symmetries and conserved quantities play in the resilience of the crystal, by observing the preferential attachment mechanism present in the formation of clusters for this specific model of time crystal. Crucially, the nature of the obtained graphs for moderate levels of error suggests that such systems could be used as scale-free-like network simulators, giving rise to a novel application of such devices in the field of complex networks Paparo2013 ; Faccin2014 . Importantly, because our networks span the configuration space, such simulation could be done with moderate numbers of qubits and thus available noisy intermediate-scale quantum (NISQ) platforms (ranging from ion traps to superconducting qubit chips) could be used for its implementation. Structural information about the network (i.e. the degree distribution) using DTCs in NISQ devices could be experimentally obtained by exploiting a quantum walk in the configuration space. One could track down the number of configurations visited dynamically in a given time by measuring a set of $l$-point spin correlation functions, a measure that is experimentally accessible with current technology. The number of configurations visited in a given time for different initial conditions is related to the degree of the different nodes of the graph, unveiling the complexity of the network (see Sec. V of the Supplementary Material for more detail on our proposed experimental protocol). Future work will include further investigation and characterization on the nature of these networks as well as the study of additional phenomena present in time crystals in terms of graphs. We believe that the use of this formalism will lead not only to a more complete and deep understanding of discrete time crystals and their related phenomena, but also be very advantageous in the study of periodically driven quantum systems.

Acknowledgements:— We thank A. Sakurai, M. Hanks, T. Haug, Y. Naka Renoust, and J. Schmiedmayer for valuable discussions. This work was supported in part from the Japanese program Q-LEAP, the MEXT KAKENHI Grant-in-Aid for Scientific Research on Innovative Areas Science of hybrid quantum systems Grant No.15H05870 and the JSPS KAKENHI Grant No. 19H00662. This project was also made possible through the support of a grant from the John Templeton Foundation (JTF 60478). The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the John Templeton Foundation. Author contributions: M.P.E, T.O. and V.M.B conceived the initial idea of the research. All authors contributed to develop the idea and analyse the system. M.P.E, T.O., V.M.B and B.R. performed numerical simulations. All authors contributed to discussions of the results and the development of the manuscript. K.S., W. J. M. and K.N. supervised the whole project. The authors declare that they have no competing interests. Data Availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Methods

.1 Floquet theory: Percolation rule and graph structure

One of the main tools used in our work is the percolation rule, which establishes when there is a link between two nodes of a graph representing the effective Hamiltonian. In this section we explain in detail the motivation of the percolation rule and its interpretation in terms of Floquet theory. In order to have an intuitive picture of the percolation rule, let us consider only two configurations $|i\rangle$ and $|j\rangle$ with energies $E_{i}$ and $E_{j}$ in the absence of drive, respectively. For simplicity, let us assume that the drive induces transitions only between the aforementioned configurations. Within the two-state approximation, the Hamiltonian from Eq.(1) reads:

[TABLE]

In general, the transitions occur when the driving is resonant with the energy level spacing between two configurations, i.e., when $E_{j}-E_{i}=\hbar k\omega$, where $k$ is an integer and $\omega=2\pi/T$ is the frequency of the drive. Therefore, when the drive is switched on, the quasienergies $\lambda_{s_{1}}$ and $\lambda_{s_{2}}$ corresponding to the energies $E_{i}$ and $E_{j}$, exhibit an anticrossing [22]. In the language of Floquet theory, the latter is referred to as a $m$-photon resonance because the quasienergies are defined modulo $\omega$.

If we restrict ourselves to the configurations $|i\rangle$ and $|j\rangle$, the effective Hamiltonian can be written as

[TABLE]

From this we can see that the energy gap for the quasienergies scales as $\delta\lambda=|\lambda_{s_{1}}-\lambda_{s_{2}}|=\sqrt{(\mathcal{E}_{i}-\mathcal{E}_{j})^{2}+(K_{ij})^{2}}$. In addition, we can derive the expressions for the eigenstates $|\Phi_{s_{1}}\rangle=\cos\theta_{ij}|i\rangle+\sin\theta_{ij}|j\rangle$ and $|\Phi_{s_{2}}\rangle=-\sin\theta_{ij}|i\rangle+\cos\theta_{ij}|j\rangle$, where $\cos\theta_{ij}=|\mathcal{E}_{i}-\mathcal{E}_{j}|/\delta\lambda$ and $\sin\theta_{ij}=|K_{ij}|/\delta\lambda$. The percolation rule $|\mathcal{E}_{j}-\mathcal{E}_{i}|<|K_{ij}|$ from Eq. (4) implies that when the link between the $i$-th and $j$-th nodes is active, the Floquet states $|\Phi_{s_{1}}\rangle$ and $|\Phi_{s_{2}}\rangle$ with quasienergies $\lambda_{s_{1}}$ and $\lambda_{s_{2}}$ are delocalized in the basis $|i\rangle$ and $|j\rangle$. On the contrary, if the percolation rule is violated, i.e., $|\mathcal{E}_{j}-\mathcal{E}_{i}|>|K_{ij}|$, the Floquet states are localized and no link between the $i$-th and $j$-th nodes is drawn in the graph.

Now let us consider a more detailed derivation of the percolation rule by using a locator-like expansion for driven systems. Here we explicitly derive the percolation rule and motivate its origin in terms of resonances. Let us first shortly discuss the most relevant aspects of the locator expansion. In the Anderson model, one can expand the resolvent in terms of the hopping strength. For strong disorder, this can be understood as a perturbative expansion of the localized wavefunctions in the hopping term of the Hamiltonian, which is referred to as the locator expansion. The latter usually diverges due to resonances [35].

In driven quantum systems, the Hamiltonian is periodic in time and all the important information is contained in the Floquet operator $\hat{\mathcal{F}}=\hat{U}(T)$, which is the evolution operator in a period $T$ of the drive. Let us begin be defining the discrete transformation

[TABLE]

The operator $\hat{G}_{\text{F}}(z)$ now plays the role of the resolvent of the Floquet operator $\hat{\mathcal{F}}$. Next we can express the Floquet operator in terms of the effective Hamiltonian, as $\hat{\mathcal{F}}=e^{-\mathrm{i}\hat{H}^{\text{eff}}_{\epsilon,T}T}$, where $\hat{H}^{\text{eff}}_{\epsilon,T}=\hat{H}_{0}+\hat{V}$ with $\hat{H}_{0}=\sum_{i}\mathcal{E}_{i}|i\rangle\langle i|$ and $\hat{V}=\sum_{i,j}K_{ij}|i\rangle\langle j|$.

Although Eq.(11) look different than the resolvent for undriven systems [35], they share the same structure if we consider the expansion

[TABLE]

in terms of the generating function for the Bernoulli numbers $B_{m}$ [36].

The next step is to use the perturbative expansion of Eq. (.1) in the hopping term $\hat{V}$

[TABLE]

This allows us to calculate the transition probability between two configurations $|i\rangle$ and $|j\rangle$, as follows

[TABLE]

where $\delta_{i,j}$ is the Kronecker delta. To investigate the physical meaning of the percolation rule, let us consider the role of the complex variable $z=\varepsilon+\mathrm{i}\delta$ when $\varepsilon=\mathcal{E}_{i}$. In this case, we can see that the relevant information of the series expansion is given by the ratio $|K_{ij}|/|\mathcal{E}_{i}-\mathcal{E}_{j}|$. This resembles the locator expansion for the undriven Anderson model. The main difference here is that that $|i\rangle$ and $|j\rangle$ are not lattice sites but configurations of spins. In any case, the series expansion (14) exhibit divergences whenever the condition $|\mathcal{E}_{i}-\mathcal{E}_{j}|<|K_{ij}|$ is satisfied. Keeping in mind the explanation of the percolation rule in terms of two level systems, here the divergence is a signature of ”hibridization” of the configurations $|i\rangle$ and $|j\rangle$. Thus, when the configurations are close to resonance, the states become delocalized under the action of the hopping term. On the contrary, when $|\mathcal{E}_{i}-\mathcal{E}_{j}|>|K_{ij}|$, the states remain localized.

.2 Power law fitting of the degree distributions

Let us first define the concept of degree within the context of graph theory. The degree of a vertex v belonging to a graph G is the number of edges connected to that vertex v. The degree distribution is therefore a probability distribution of vertex degrees over the whole graph G. Here we explain how the degree distributions of the graphs obtained from the effective Hamiltonian of the $2T$-DTC can be fit to a power law function $p(k)\propto k^{-\beta}$. We have followed the method by Clauset et al. [31]. First, the degree distributions $P(k)$ are obtained by measuring log-binned histograms of the degrees $k$ of the graph. We have obtained the distributions from 100 realizations of the graph with different disorder values for statistical convergence. Generally, the power law behaviour is observed in the large degree tail of the distribution rather than the entire domain. Therefore, we estimate the lower bound $k=k_{\text{min}}$ where the best power law fit can be obtained. This is done by estimating the power law exponent $\beta$ by applying maximum likelihood method on a certain portion $k\geq k_{\text{cutoff}}$ of the distribution. Once the exponent is estimated, the Kolmogorov-Smirnov distance between the distribution and the fit is computed. We compute this with various $k_{\text{cutoff}}>0$ and the $k_{\text{min}}$ is chosen at the $k_{\text{cutoff}}$ where the Kolmogorov-Smirnov distance is minimized. Once the $k_{\text{min}}$ and $\beta$ is estimated, we examine the goodness of fit of the power law using a comparative test. We consider an alternative probability distribution function that may be good fit for our data, and apply a likelihood ratio test between the power law and the alternative function. Here we chose the lognormal function as the alternative, which is another heavy-tailed function. From this comparative test, we conclude that the power law function well describes the degree distributions of our graphs.

**Supplementary materials for “Simulating complex quantum networks with time crystals”

** M. P. Estarellas,1,111These authors contributed equally. T. Osada,2,1,∗, V. M. Bastidas3,∗, B. Renoust4,1,5, K. Sanaka2, W. J. Munro3,1 and K. Nemoto1,5

1*National Institute of Informatics, 2-1-2 Hitotsubashi, Chiyoda-ku, Tokyo 101-8430, Japan

2Tokyo University of Science, 1-3 Kagurazaka, Shinjuku, Tokyo, 162-8601, Japan

3NTT Basic Research Laboratories & Research Center for Theoretical Quantum Physics,

3-1 Morinosato-Wakamiya, Atsugi, Kanagawa, 243-0198, Japan

4Osaka University, Institute for Datability Science, 2-8 Yamadaoka, Suita, Osaka Prefecture 565-0871, Japan

5Japanese-French Laboratory for Informatics, CNRS UMI 3527,

2-1-2 Hitotsubashi, Chiyoda-ku, Tokyo 101-8430, Japan

Contents

1. .1 Floquet theory: Percolation rule and graph structure
2. .2 Power law fitting of the degree distributions
3. I A roadmap through the supplementary materials
4. II Level statistics of the quasienergies and manybody localization
5. III Power spectrum of time series: Stability of the time-crystalline phase
6. IV Semiclassical description of a discrete time crystal: a stability analysis of the fixed points in the absence of error $\epsilon=0$
7. V Quantum walks in the configurations space: Experimental protocol to simulate complex quantum networks

I A roadmap through the supplementary materials

The purpose of this supplementary material is to provide the reader with the necessary tools to understand the main manuscript. In section II we focus on the relation between the robustness of time crystals and many-body localization (MBL). We explain in detail how to obtain the level statistics for the quasienergy gaps, which is a well-known diagnostic of MBL. In section III we show how to obtain signatures of the melting process by calculating the power spectrum of a time series of measured observables. Our approach is experimentally feasible and provides clear evidence that some configurations are more stable than others. In section IV we provide a semiclassical description of the stability of the different spin configurations. Finally, a very important aspect of our results is that there is an emergent scale-free behavior in the graphs associated to the DTC. With this purpose in mind, in section V we provide an experimental protocol to simulate complex quantum networks in near-term quantum devices. Our approach is based on a quantum walk in the configuration space, which can be realized by measuring $l$-point spin correlation functions.

II Level statistics of the quasienergies and manybody localization

In classical mechanics, integrable systems are characterized by the existence of integrals of motion restricting the motion to orbits that are constrained by conservation rules. The classical harmonic oscillator, for example, exhibits a two-dimensional phase space with coordinates $q$ and $p$ representing position and momentum, respectively. In the absence of driving, the energy of the oscillator is conserved. Thereby, given an initial condition $(q_{0},p_{0})$, the classical orbits are constrained to the energy surface $E(q,p)=p^{2}/2M+M\Omega q^{2}/2=p_{0}^{2}/2M+M\Omega q_{0}^{2}/2$, where $M$ is the mass and $\Omega$ is the oscillation frequency of the oscillator. For any value of the energy, the classical orbits are periodic with a period $T_{\text{orbit}}=2\pi/\Omega$. In the quantum world, quantum signatures of the periodic orbits appear in the spacings $\delta_{n}=E_{n+1}-E_{n}=\hbar\Omega$ between the eigenvalues $E_{n}=\hbar\Omega(n-1/2)$ of the Hamiltonian $\hat{H}=\hbar\Omega(\hat{a}^{\dagger}\hat{a}-1/2)$, with $n$ labelling each of the system’s eigenstates and $\hat{a}$ ($\hat{a}^{\dagger}$) being the anhiliation (creation) operator. Interestingly, if we calculate the probability distribution of level spacings $\delta_{n}$, this would be a Dirac delta distribution centered at $\Omega$. When we break conservation rules, however, the integrability of the system is compromised because the conserved quantities are destroyed and the level spacing distribution largely deviates from a simple Dirac-delta type distribution. This statistical analysis is one of the cornerstones of recent developments in the theory of thermalization, which is closely related to random matrix theory and quantum signatures of chaos [37, 38].

Recently, there has been an increasing interest in the study of level statistics as a diagnosis for thermalization and manybody localization. In the context of quantum manybody systems, the interactions distribute the energy between all the available states, because they usually break conserved quantities at the single-particle level, leading to thermalization and ergodic behavior. In systems where there is a well defined semiclassical limit, this behavior is associated with a fully chaotic phase space [37, 38]. As we discussed above, in quantum mechanics, level statistics is a powerful tool to determine the nature of the dynamics. In terms of level statistics, this means that the statistical behavior of an ergodic system will largely deviate from the quantum harmonic oscillator case described above. In particular, the destruction of conserved quantities induces strong correlations between the energy eigenvalues that are reflected in a strong level repulsion [38]. In the case of disordered manybody systems, there is a competition between disorder and interactions. Disorder can suppress ergodicity giving rise to the manybody localized phase, which is characterized by an extensive number of conserved quantities [39, 40, 41]. In the absence of drive, one is interested in the distribution of ratios $P(r)$ with $r_{s}=\min(\delta_{s},\delta_{s+1})/\max(\delta_{s},\delta_{s+1})\leq 1.0$ with $\delta_{s}=E_{s+1}-E_{s}$, where $E_{s}$ are the energy eigenvalues sorted in order of increasing value [38]. In the MBL phase, the levels are uncorrelated because the manybody states are highly localized in space and the level statistics follows a Poissonian distribution

[TABLE]

When the interactions overcome the effect of disorder, the aforementioned conserved quantities are destroyed and the system is able to explore all the configuration space with a constant energy. This is reflected in a strong level repulsion because some states are delocalized. As a consequence, in the case of real Hamiltonians, the system follows a universal level statistics

[TABLE]

associated to the Gaussian Orthogonal ensemble (GOE) of random matrices [38, 42].

In our manuscript, we are interested in the dynamics of DTCs under the effect of a rotation error as shown in the Hamiltonian from Eq. (1). Below we show that the time-crystalline order is directly related to the presence of an MBL phase. In analogy to the quantum harmonic oscillator, in the absence of error, the time crystal has several conserved quantities and the motion is periodic. However, discrete time crystals are quantum phases of matter that appear in periodically-driven quantum systems and energy is not conserved anymore. Therefore, we have to work with gaps $\delta_{s}=\lambda_{s+1}-\lambda_{s}$, where $-\hbar\omega/2<\lambda_{s}<\hbar\omega/2$ are the quasienergies and $\omega$ is the frequency of the drive. Under the effect of a small rotation error, there is a small coupling between different symmetry multiplets and the quasienergy level statistics follows a Poissonian behavior as in Eq. (S1). As a consequence, the system is in the MBL phase, which protects the system from heating up to infinite temperatures [43, 9]. When the error is increased, the conserved quantities are destroyed and the Floquet states become highly delocalized in the configuration space. The DTC melts and an MBL-to-ergodic phase transition takes place. The statistics of levels in this case is given by

[TABLE]

which is the same statistics as the circular orthogonal ensemble of random matrices [43, 44].

One of the most intriguing aspects of level statistics is its intimate relation to quantum signatures of chaos and random matrix theory. At this stage, it is important to clarify a technical point concerning the transition between regular and chaotic motion and its intimate relation to the MBL-to-ergodic transition in manybody systems. As we discussed above, the MBL phase is characterized by an extensive number of conserved quantities, which leads to a Poissonian level statistics. Under the effect of a perturbation, these conserved quantities are destroyed and the system has a transition to the ergodic phase. This process resembles the destruction of tori in phase space within the framework of KAM theorem. For small perturbations, some tori are more susceptible to be destroyed and as we increase the strength of the perturbation, the phase space becomes mixed: there are regular islands associated to the surviving tori and a chaotic sea. When all the tori are destroyed, the dynamics is fully chaotic. In quantum systems, signatures of this chaotic dynamics lead to a universal behavior of the level statistics giving rise to distributions $P_{\text{GOE}}(r)$ and $P_{\text{COE}}(r)$ for undriven and driven systems, respectively. Universality means that these distributions do not depend on microscopic details of the dynamics but are solely determined by symmetries of the system. However, quantum signatures of a mixed phase space are not universal. In the context of time crystals, when we increase the error, we see a crossover between $P_{\text{Poisson}}(r)$ and $P_{\text{COE}}(r)$, because the error destroys conserved quantities. As we increase the rotation error, we observe that the level statistics is closer to $P_{\text{COE}}(r)$. Figure S1 depicts the results for the level statistics in the case of a time crystal for errors $\epsilon=0.01$ and $\epsilon=0.1$. Our results therefore show that time crystalline order is protected by MBL.

III Power spectrum of time series: Stability of the time-crystalline phase

In this section our aim is to provide an additional experimentally-accessible way to investigate the stability of the different configurations. The latter will shed light on the melting mechanism of the time crystal and the emergent scale-free behavior. Our idea is based on a quantity that can be easily measured in several platforms. We just require to measure the expectation value of the total magnetization $M_{i,Z}(mT)=1/n\sum_{r=1}^{n}\langle\psi(mT)|\sigma_{r}^{z}|\psi(mT)\rangle$ at stroboscopic times, where $|\psi(0)\rangle=|i\rangle$. This quantity is related to the populations, as follows

[TABLE]

where $(j_{1},j_{2},\ldots,j_{n})_{2}$ is the binary decomposition of the integer number $j$ for all the configurations.

Next, for a given initial configuration $i$, one can record the measurements $M_{i,Z}(mT)$, thus obtaining a time series

[TABLE]

where $N$ is the number of periods that the system evolve. With the data for the time series at hand, we can calculate the discrete Fourier transformation

[TABLE]

where $\omega_{k}=2\pi k/NT$ and $k\in[0,N-1]$. For a given rotation error $\epsilon$, the associated power spectrum

[TABLE]

tell us how strong is the contribution of a the $k-$th harmonic to the time series. Note that for convenience, we have arranged the Fourier coefficients in the form of a vector $\bm{V}_{i,\epsilon}$. In the absence of error, all the configurations show a magnetization that varies as $M_{i,Z}(mT)=(-1)^{m}M_{i,Z}(0)=e^{\mathrm{i}n\pi}M_{i,Z}(0)$ and the power spectrum shows a single peak at a frequency $\omega_{N/2}$, which is half the frequency of the drive. The latter is a signature of the subharmonic response of the manybody system. As we discussed in Sec.II, the periodic motion is related to local integrals of motion. The DTC is protected agains heating due to MBL. However, as we increase the rotation error $\epsilon$, some configurations are more affected than others and they loose their periodic character. On the contrary, configurations like $|2^{n}-1\rangle=|\uparrow\uparrow\uparrow\ldots\uparrow\uparrow\rangle$ and $|0\rangle=|\downarrow\downarrow\downarrow\ldots\downarrow\downarrow\rangle$ are the most stable ones. To have a quantitative measurement of this statement that can be easily accessed in experiments, we propose to measure the fidelity $F_{i}^{\text{ps}}(\epsilon)$ of the power spectrum, defined as follows

[TABLE]

This very simple quantity tell us how far is the power spectrum from the one of a time crystal without error. In figure S2 we depict the fidelity as a function of the configuration number $i$ and the error $\epsilon$. This shows a clear evidence that not all the configurations are stable under the effect of a rotation error.

IV Semiclassical description of a discrete time crystal: a stability analysis of the fixed points in the absence of error $\epsilon=0$

In this section, we present numerical evidence of the stability of the different configurations from a semiclassical perspective. We start by considering the effective Hamiltonian for two periods of the drive

[TABLE]

which describes the dynamics of the system at stroboscopic times $t_{n}=2nT$.

After visualizing the classical orbits on the Bloch sphere, it is instructive to investigate a limiting case of the system. Here we focus on the system in the absence of error. In this case, there is an exact mapping of Eq.( Simulating complex quantum networks with time crystals) to classical Hamiltonian

[TABLE]

Configurations such as $\theta^{c}_{l}=0,\pi$ are fixed points of the stroboscopic dynamics. The latter can be derived from the first derivative of the energy function $\partial_{\theta_{j}}\mathcal{H}_{0,2T}(\theta_{1},\ldots,\theta_{n})=0$. For example, the quantum $|2^{n}-1\rangle=|\uparrow\uparrow\uparrow\ldots\uparrow\uparrow\rangle$ has a classical counterpart $(\theta^{c}_{1},\theta^{c}_{2},\ldots,\theta^{c}_{n})=(0,0,\ldots,0)$. Given a fixed point, we can explore its stability properties by investigating the eigenvalues of the Jacobian matrix with elements $\mathcal{J}_{i,j}=\partial_{\theta_{i}}\partial_{\theta_{j}}\mathcal{H}_{0,2T}(\theta_{1},\ldots,\theta_{n})|_{\theta=\theta^{c}}$, where the derivatives are evaluated at the critical points $\theta^{c}$ such that $\partial_{\theta_{j}}\mathcal{H}_{0,2T}(\theta_{1},\ldots,\theta_{n})=0$. For simplicity, let us focus on a particular example of a lattice with four sites and open boundary conditions. In the case of the configuration $(\theta^{c}_{1},\theta^{c}_{2},\theta^{c}_{3},\theta^{c}_{4})=(0,0,0,0)$ these matrix elements are given by

[TABLE]

and the equilibrium is stable [45]. For any configuration of the form $\theta^{c}_{l}=0,\pi$ the off-diagonal elements of the Jacobian matrix are zero. Let us consider the configuration $(\theta^{c}_{1},\theta^{c}_{2},\theta^{c}_{3},\theta^{c}_{4})=(0,0,\pi,\pi)$ with one domain wall, where $\mathcal{J}_{1,1}=-J^{z}_{1,2}+J^{z}_{1,3}+J^{z}_{1,4}$, $\mathcal{J}_{2,2}=-J^{z}_{1,2}+J^{z}_{2,3}+J^{z}_{2,4}$, $\mathcal{J}_{3,3}=-J^{z}_{1,3}+J^{z}_{2,3}-J^{z}_{3,4}$ and $\mathcal{J}_{4,4}=J^{z}_{1,4}+J^{z}_{2,4}-J^{z}_{3,4}$. Due to the power-law decay of the couplings, all the diagonal elements are negative and the fixed point is stable. Finally, for a configuration with three domain walls $(\theta^{c}_{1},\theta^{c}_{2},\theta^{c}_{3},\theta^{c}_{4})=(0,\pi,0,\pi)$, the matrix elements are $\mathcal{J}_{1,1}=J^{z}_{1,2}-J^{z}_{1,3}+J^{z}_{1,4}$, $\mathcal{J}_{2,2}=J^{z}_{1,2}+J^{z}_{2,3}-J^{z}_{2,4}$, $\mathcal{J}_{3,3}=-J^{z}_{1,3}+J^{z}_{2,3}+J^{z}_{3,4}$ and $\mathcal{J}_{4,4}=J^{z}_{1,4}-J^{z}_{2,4}+J^{z}_{3,4}$. In this case all the eigenvalues are positive. The elements of the Jacobian give us information about the local curvature of the energy function at the critical points. As long as they have the same sign, the fixed points are stable. For a non-zero value of the error, configurations such as $\theta^{c}_{l}=0,\pi$ are not fixed points of the dynamics anymore, but they belong to trajectories in phase space.

V Quantum walks in the configurations space: Experimental protocol to simulate complex quantum networks

So far, we have discussed spectral properties of driven systems and its consequences for the time crystal. The percolation rule discussed in the main text provides us with graphical way to represent localization properties of Floquet states. Therefore, the connectivity of the graph tell us information about how many configurations “participate” in a given Floquet state. When the latter are delocalized, the gaps between quasienergies become highly correlated and the system becomes ergodic as we discussed in section II. This establish an intriguing relation between the melting of a discrete time crystal and the MBL-to-ergodic transition as we increase the rotation error. During this process, the graph exhibits a scale-free behavior characteristic of complex quantum networks. Here, we describe in detail an experimental protocol to faithfully extract structural information about the simulated complex quantum networks in Noisy intermediate-scale quantum devices.

The cornerstone of our approach is the representation of the effective Hamiltonian

[TABLE]

in terms of the configuration basis states $\{|i\rangle\}$. Here $\mathcal{E}_{i}$ is the energy of configuration $|i\rangle$, and $K_{ij}$ is the transition energy between configurations $|i\rangle$ and $|j\rangle$. The quasienergies $\lambda_{s}$ have a complicated dependence on the parameters $\mathcal{E}_{i}$ and $K_{ij}$ that determine localization properties of the Floquet states $|\Phi_{s}\rangle$. The latter are linear combinations of configurations $|i\rangle$ with $i=0,1,2\ldots,2^{n}-1$, where $n$ is the number of sites of the spin chain. Interestingly, although the effective Hamiltonian describes a manybody system, in terms of configurations, the dynamics resemble a quantum walk in a quantum complex network.

Experimentally, one can initially prepare the state of the system in a given configuration such as $|\psi(0)\rangle=|i\rangle=|\uparrow\uparrow\uparrow\ldots\downarrow\downarrow\rangle$, which represents a node of the complex network. At stroboscopic times $t_{n}=nT$, the dynamics are governed by the effective Hamiltonian $\hat{H}^{\text{eff}}_{\epsilon,T}$ and the evolution of the state reads $|\psi(nT)\rangle=e^{-\mathrm{i}\hat{H}^{\text{eff}}_{\epsilon,T}nT/\hbar}|\psi(0)$. Inherently, the most relevant aspects of the dynamics can be recovered from this discrete evolution. Next, let us focus on the density matrix

[TABLE]

where $|\psi(nT)\rangle=\sum_{i}A_{i}(nT)\|i\rangle$. This expression give us already a recipe to measure the dynamics of the quantum walk in existing Noisy intermediate-scale quantum devices. This requires a modest set of measurements and it is straightforward with current experimental feasibilities. The idea is to measure the populations $\rho_{i,i}(nT)=|A_{i}(nT)|^{2}$ of the density matrix stroboscopically, to determine which configurations are occupied. This information can be accessed by measuring the correlation functions

[TABLE]

corresponding to all the configurations $i=0,1,2,\dots,2^{n}-1$. Here $(i_{1},i_{2},\ldots,i_{n})_{2}$ is the binary decomposition of the integer number $i$. For example, if we want to measure the population of the configuration $|2^{n}-1\rangle=|\uparrow\uparrow\uparrow\ldots\uparrow\uparrow\rangle$, we need to measure the correlation

[TABLE]

because $2^{n-1}=(1,1,\ldots,1)_{2}$. Similarly, to measure the population of $|0\rangle=|\downarrow\downarrow\downarrow\ldots\downarrow\downarrow\rangle$, we need to measure the correlation

[TABLE]

as $0=(0,0,\ldots,0)_{2}$. Interestingly, the aforementioned correlations can be reconstructed from measurements of all the possible combinations of single-point correlators $C_{Z}^{r_{1}}(nT)=\text{tr}\left[\hat{\rho}(nT)\sigma_{r_{1}}^{z}\right]$, two-point correlators $C_{Z}^{r_{1},r_{2}}(nT)=\text{tr}\left[\hat{\rho}(nT)\sigma_{r_{1}}^{z}\sigma_{r_{2}}^{z}\right]$ up to n-point correlators $C_{Z}^{r_{1},r_{2},\ldots,r_{n}}(nT)=\text{tr}\left[\hat{\rho}(nT)\sigma_{r_{1}}^{z}\sigma_{r_{2}}^{z}\cdots\sigma_{r_{n}}^{z}\right]$. In trapped ions, it is possible to measure these correlations by resorting to measurements on single ions. In superconducting qubits, it is even possible to perform full tomography for qubits arrays up to $N=12$ qubits and measurements of these correlations are within reach with current technology [46]. After recording the data of the populations $|A_{i}(nT)|^{2}$ for all the configurations $i$ at stroboscopic times $t_{n}=nT$, one can generate a space-time plot to depict the data $\left[nT,i,\rho_{i,i}(nT)\right]$. The latter is depicted in Fig. S3. The color scale represent the populations and one can visualize how the initial states travel through the configuration space.

The experimental protocol mentioned above provides us with a way of accessing very relevant information related to its degree distribution. To have a simple picture of this idea in mind, let us consider the case of a DTC in the absence of rotation error. In this case, if we prepare an initial configuration $|\psi(0)\rangle=|\uparrow\uparrow\uparrow\ldots\downarrow\downarrow\rangle$, this correspond to a Floquet cat state, because it is a quantum superposition of two Floquet states with quasienergies $\pm\hbar\pi/T$. Therefore, the expectation value of the total transverse magnetization $M_{i,Z}(nT)=1/n\sum_{r=1}^{n}\langle\psi(nT)|\sigma_{r}^{z}|\psi(nT)\rangle$ will show an oscillatory behavior with a period $2T$. In the configuration space, the quantum walk will show dynamics between just two configurations $|2^{n}-1\rangle=|\uparrow\uparrow\uparrow\ldots\uparrow\uparrow\rangle$ and $|0\rangle=|\downarrow\downarrow\downarrow\ldots\downarrow\downarrow\rangle$. If we now calculate how many configurations ”participate” in the quantum state $|\psi(nT)\rangle$, then we have a measure of how connected the network in the configuration space is. To be more precise, let us define the participation ratio

[TABLE]

which can be obtained from the measurement of the correlation function Eq. (V).

When approximating $\hat{H}^{\text{eff}}_{\epsilon,T}$ up to a first order in the error $\epsilon$, one can easily see that the local rotation that connects different resonant configurations (and therefore generates the clusters in the graph) has an approximate tunneling time of $\tau\sim\frac{T}{g\epsilon T_{1}}$. If we obtain $P_{i}(\tau)$ for a given initial $i$-th node or configuration, one can estimate how many configurations have been visited after the evolution time $\tau$, and thus estimate the degree of this node. If we now repeat the process for all the possible initial configurations $i=0,1,2,\ldots,2^{n}-1$ one can obtain a distribution for the participation ratio, $PR$. In Fig. S4 we present the distributions obtained from this experimental protocol for the same system that was depicted in Fig. 2 (network of $n=8$ sites and a rotation error of $\epsilon=0.012$ and $\epsilon=0.1$). We note that the shape of the distributions obtained from the calculation of $P_{i}(\tau)$, a) heavy-tailed with a power-law fit and b) normal distribution, matches the ones of Fig. 2 obtained directly from the degree distribution of the network. In a nutshell, by measuring the aforementioned correlation functions, we can simulate the network and obtain important information without calculating the unitary operator and without full quantum process tomography.