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

## INTRODUCTION

Symmetries are of utmost importance in condensed matter physics and statistical mechanics owing to the strong relation between quantum phases of matter and symmetry breaking (*1*–*4*). Among the zoo of symmetry-broken phases, discrete time crystals (DTCs) play a fundamental role on their own due to the type of symmetry involved (*5*). 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 (*6*), and a discrete version of time crystals can be realized by periodically driven quantum systems (*7*–*9*). 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 (*5*).

The exploration of time crystals is a very active field of research and several experimental realizations with trapped ions (*10*), dipolar spin impurities in diamond (*11*), ordered dipolar many-body systems (*12*), ultracold atoms (*13*), and nuclear spin-1/2 moments in molecules (*14*) have been achieved. Yet, an intuitive and complete insight into the nature of time crystals and its characterization, as well as a set of proposed applications, is lacking. We here provide 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 (*15*, *16*). 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 (*17*, *18*). This type of networks, far from being regular or random, contains nontrivial topological structures present in many biological, social, and technological systems. Small-world and scale-free networks are two of the most popular examples, the latter 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 (*19*), the development of new algorithms in deep learning (*14*), or the analysis of genetic and neural structures in biological systems (*20*).

## 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 *T* = 2π/ω being the period of that drive. Our approach builds from the calculation of the Floquet operator *t* = 0 to time *t* = *T* (*21*–*23*), with * _{s}*⟩ are known as Floquet states and −ℏω/2 ≤ λ

*≤ ℏω/2 are the quasienergies. At discrete times*

_{s}*t*=

_{n}*nT*, the evolution can be understood in terms of an effective “time-independent” Hamiltonian

*22*,

*24*,

*25*).

### DTCs of period 2*T*

We now move to a well-known period-2 DTC (2*T*-DTC) model (*10*, *15*) 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}

Here, *l*th spin, *l* and *m* that takes the form of an approximate power law decay with a constant exponent α, and *g* satisfies the condition 2*gT*_{1} = π such that *x* axis, with a rotation error ϵ also introduced. The Floquet operator is then given by

Crucially, the Floquet operator depends on the error ϵ, which determines how the time crystal will melt. In our work, we choose α = 1.51, *J*_{0}*T*_{2} = 0.06 with a disorder strength *WT*_{2} = π, which are similar values used in the recent experiment (*10*). The key feature of this unitary evolution is that at each period *T*, as shown in Fig. 1A, the system evolves between the states ∣↑↑⋯↑〉 and ∣ ↓↓⋯ ↓ 〉, when the state is initialized at one of these two. This periodic evolution is robust against moderate errors ϵ in

### Floquet graphs

Let us now introduce the concept of a Floquet graph as the graphical representation of the effective Hamiltonian, *26*), obtained as

The effective Hamiltonian, *N _{ℋ}*) and can be written as a tight-binding model using the configuration basis states { ∣

*i*⟩} as

*is the energy of configuration ∣*

_{i}*i*〉 and

*K*is the transition energy between configurations ∣

_{ij}*i*〉 and ∣

*j*〉 (see Fig. 1B). The Floquet states ∣Φ

*〉 are linear combinations of the configurations ∣*

_{s}*i*〉. In general, the quasienergies λ

*are complicated functions of the parameters ℰ*

_{s}_{i}and

*K*. In this work, we consider a system with

_{ij}*n*spin-1/2 particles, and the configurations {∣

*i*〉} with

*i*= 0,1, …, 2

*− 1 can be chosen as the*

^{n}*N*

_{ℋ}= 2

*elements of the computational basis, i.e., the product states of the eigenstates of*

^{n}*i*= (

*j*

_{1}

*j*

_{2}…

*j*)

_{n}_{2}is the binary decomposition of the integer

*i*and

*− 1〉 = ∣ ↑ ↑ ⋯ ↑ 〉. The basis configurations have been chosen to satisfy the relation*

^{n}*i*〉 and ∣2

*− 1 −*

^{n}*i*〉 are always linked by the global π rotation from

*nodes of the graph. In addition, such nodes will be linked according to a percolation rule (*

^{n}*27*), where a nonweighted edge is considered to be active between nodes

*i*and

*j*only if the condition

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 (*28*). We want to stress that the nodes represent many-body states ∣*i*〉 and not physical spin sites *l*, so the graph spans through the full Hilbert—configuration—space (see Fig. 1C). This percolation rule was proven to be useful to detect many-body localization to thermal phase transitions in undriven spin systems (*27*).

### Time crystal graphs

Under the framework of Floquet graphs, the dynamics of the described *2T*-DTC can now be investigated. The effective Hamiltonian *T*-DTC crystals, as shown in the example presented in Fig. 1C. We note that this visual characterization of the crystal time order is not unique to 2*T*-DTC and can be extrapolated to crystals of different periodicity. For example, 3*T*-DTC presents decoupled trimers, 4*T*-DTC presents decoupled tetragons, and so on. Our results hence reveal that, in the case of 2*T*-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〉 and ∣2* ^{n}* − 1〉 defined above, there is an energy gap ℏπ/

*T*between the quasienergies λ

_{0}and λ

_{2n−1}. The corresponding Floquet states are given by

*i*〉 and ∣2

*− 1 −*

^{n}*i*〉, as all the paired Floquet states ∣Φ

*〉, ∣Φ*

_{i}_{2n − 1 − i}〉 have an energy difference of ℏπ/

*T*(

*9*).

### Melting of the 2*T*-DTC

Consider now the case where our system is initialized to ∣ψ(0)〉 = ∣2* ^{n}* − 1〉 = ∣ ↑ ↑ ⋯ ↑ 〉 (

*10*). This state is a superposition of Floquet states,

*− 1〉. Such a dimer remains decoupled and isolated for modest levels of error (up to a ~3% error in the rotation added as a nonzero value of ϵ in*

^{n}*T*-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. 2A). 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. 2D). The presence/absence 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. We note that not all dimers show the same robustness, and some of them survive to higher error values. This 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*〉 and ∣2* ^{n}* − 1 −

*i*〉. 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 *T*), the Floquet operator for no error can be reduced to the following expression

In the absence of error, ϵ = 0, the disorder cancels out exactly and the effective Hamiltonian reads * ^{n}* − 1〉, the system will remain stroboscopically within the symmetry multiplet with 2

*T*quasienergy

*10*).

### 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

*n*subblocks given by the conservation of the number of domain walls. In the absence of error (ϵ = 0), and after 2

*T*, one can find a classical quasienergy surface in phase space associated to the effective Hamiltonian

*n*-dimensional phase space

*is the polar angle. From this, one can show that for zero error ϵ = 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 sections SIII and SIV for more details). To investigate this, we use the Baker-Campbell-Hausdorff formula to obtain, up to a first-order approximation, the effective Hamiltonian for a nonzero value of ϵ*

_{l}This allows us to understand how ϵ destroys the symmetries by coupling different symmetry multiplets. This can be observed from Fig. 3 (B to 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 (*29*), with our perturbation being ϵ. First of all, the transverse field term *29*). 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; that is, ∣ψ(0)〉 = ∣2* ^{n}* − 1〉 (see the darkest blue dimer on the top right corner of Fig. 2A). 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 (ϵ ∼ 0.012). In the previous section, we have seen that the quantum terms that appear in Eq. 7 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. 2C. The tail of these distributions can be fit to a power law distribution, which is a characteristic of scale-free networks. Following the recommendations of Clauset *et al.* (*30*), 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 favors over the lognormal distribution (see the “Power law fitting of the degree distributions” section in Methods for further details).

In addition, exploring the behavior of the average degree of the graph obtained from *T*-DTC. From Fig. 4, we can confirm that the nodes with five or six domain walls (where the *n* = 8 case is considered here) tend to acquire neighbors preferentially and have more degree already from a small error ϵ < 0.02. This is the region where we can observe the scale-free behavior emerging. When the error is increased, the preferential behavior becomes weaker, and eventually most of the nodes tend to have a large degree, closer to the structure of random networks, which indicates that the crystal has melted.

This behavior not only has consequences for the dynamics of the time crystals but also can be used to perform quantum simulation of networks with exotic properties such as scale-free–like ones (*18*). Note that the scale-free–like networks shown in Fig. 2 are only observed in the effective Hamiltonian obtained from a single period *18*, *31*). We observe that a small error has little effect to the effective Hamiltonian *z* axis caused by the *B ^{z}* term. Therefore, at time 2

*T*, the effect of the

*B*term is effectively cancelled, limiting the number of transitions happening. On the other hand, this is not the case at time

^{z}*T*, and off-diagonal terms appear in the matrix form of

## DISCUSSION

Our approach of representing DTCs 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 2*T*-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 previously unknown application of such devices in the field of complex networks (*32*, *33*). 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 section SV for more details 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 not only lead to a more complete and deep understanding of DTCs and their related phenomena but also be very advantageous in the study of periodically driven quantum systems.

## METHODS

### 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. To have an intuitive picture of the percolation rule, let us consider only two configurations ∣*i*〉 and ∣*j*〉 with energies *E _{i}* and

*E*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

_{j}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}*k*ω, where

*k*is an integer and ω = 2π/

*T*is the frequency of the drive. Therefore, when the drive is switched on, the quasienergies λ

_{s1}and λ

_{s2}, corresponding to the energies

*E*and

_{i}*E*, respectively, exhibit an anticrossing (

_{j}*21*). In the language of Floquet theory, the latter is referred to as an

*m*-photon resonance because the quasienergies are defined modulo ω.

If we restrict ourselves to the configurations ∣*i*〉 and ∣*j*〉, the effective Hamiltonian can be written as

From this, we can see that the energy gap for the quasienergies scales as _{s1}〉 = cos θ* _{ij}*∣

*i*〉 + sin θ

*∣*

_{ij}*j*〉 and ∣Φ

_{s2}〉 = − sin

*θ*∣

_{ij}*i*〉 + cos

*θ*∣

_{ij}*j*〉, where cos θ

*= ∣ ℰ*

_{ij}*− ℰ*

_{i}*∣ /δλ and sin θ*

_{j}*= ∣*

_{ij}*K*∣ /δλ. The percolation rule ∣ℰ

_{ij}*− ℰ*

_{j}*∣ < ∣*

_{i}*K*∣ from Eq. 4 implies that when the link between the

_{ij}*i*th and

*j*th nodes is active, the Floquet states ∣Φ

_{s1}〉 and ∣Φ

_{s2}〉 with quasienergies λ

_{s1}and λ

_{s2}are delocalized in the basis ∣

*i*〉 and ∣

*j*〉. In contrast, if the percolation rule is violated, i.e., ∣ℰ

*− ℰ*

_{j}*∣ > ∣*

_{i}*K*∣, the Floquet states are localized and no link between the

_{ij}*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 wave functions in the hopping term of the Hamiltonian, which is referred to as the locator expansion. The latter usually diverges due to resonances (*34*).

In driven quantum systems, the Hamiltonian is periodic in time and all the important information is contained in the Floquet operator *T* of the drive. Let us begin by defining the discrete transformation

The operator

Although Eq. 11 looks different from the resolvent for undriven systems (*34*), they share the same structure if we consider the expansion*B _{m}* (

*35*).

The next step is to use the perturbative expansion of Eq. 12 in the hopping term

This allows us to calculate the transition probability between two configurations ∣*i*〉 and ∣*j*〉, as follows_{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* = ε + *i*δ when ε = ℰ* _{i}*. In this case, we can see that the relevant information of the series expansion is given by the ratio ∣

*K*∣ / ∣ℰ

_{ij}*− ℰ*

_{i}*∣. This resembles the locator expansion for the undriven Anderson model. The main difference here is that that ∣*

_{j}*i*〉 and ∣

*j*〉 are not lattice sites but configurations of spins. In any case, the series expansion (Eq. 14) exhibits divergences whenever the condition ∣ℰ

*− ℰ*

_{i}*∣ < ∣*

_{j}*K*∣ is satisfied. Keeping in mind the explanation of the percolation rule in terms of two-level systems, here, the divergence is a signature of “hybridization” of the configurations ∣

_{ij}*i*〉 and ∣

*j*〉. Thus, when the configurations are close to resonance, the states become delocalized under the action of the hopping term. In contrast, when ∣ℰ

*− ℰ*

_{i}*∣ > ∣*

_{j}*K*∣, the states remain localized.

_{ij}### 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 2*T*-DTC can be fit to a power law function *p*(*k*) ∝ *k*^{−β}. We have followed the method by Clauset *et al.* (*30*). 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 behavior is observed in the large degree tail of the distribution rather than the entire domain. Therefore, we estimate the lower bound *k* = *k*_{min} where the best power law fit can be obtained. This is done by estimating the power law exponent β by applying the maximum likelihood method on a certain portion *k* ≥ *k*_{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*_{cutoff} > 0, and the *k*_{min} is chosen at the *k*_{cutoff} where the Kolmogorov-Smirnov distance is minimized. Once the *k*_{min} and β are 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 a 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

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/42/eaay8892/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank A. Sakurai, M. Hanks, T. Haug, Y. Naka Renoust, and J. Schmiedmayer for valuable discussions.

**Funding:**This work was supported, in part, by 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 analyze 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.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials 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.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).