Research ArticlePHYSICS

Simulating complex quantum networks with time crystals

See allHide authors and affiliations

Science Advances  16 Oct 2020:
Vol. 6, no. 42, eaay8892
DOI: 10.1126/sciadv.aay8892


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.


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 (14). 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 (79). 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).


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+T)=Ĥ(t), with T = 2π/ω being the period of that drive. Our approach builds from the calculation of the Floquet operator F̂=Û(T)=T̂ exp (i0TĤ(τ)dτ/), which is the evolution operator within one period from time t = 0 to time t = T (2123), with T̂ being the time ordering operator. In Floquet theory, one wishes to solve the eigenvalue problem F̂Φs=eiλsT/Φs. The eigenvectors ∣Φs⟩ are known as Floquet states and −ℏω/2 ≤ λs ≤ ℏω/2 are the quasienergies. At discrete times tn = nT, the evolution can be understood in terms of an effective “time-independent” Hamiltonian Ĥeff such that F̂=eiĤeffT/. 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 (22, 24, 25).

DTCs of period 2T

We now move to a well-known period-2 DTC (2T-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 = T1 + T2Ĥ(t)={Ĥ1g(1ϵ)Σlσlx0<t<T1Ĥ2ΣlmJlmzσlzσmz+ΣlBlzσlzT1<t<T(1)

Here, {σlx,σly,σlz} are the usual Pauli operators on the lth spin, JlmzJ0/lmα is the long-range interaction between spins l and m that takes the form of an approximate power law decay with a constant exponent α, and Blz[0,W] is a random longitudinal field. From the interplay between the driving, interactions, and disorder, an Many-Body Localization phase arises, which is a key ingredient that allows the existence of DTCs. 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 section SII for a more detailed description). The parameter g satisfies the condition 2gT1 = π such that Û1=exp (iĤ1T1/) becomes a global π pulse around the x axis, with a rotation error ϵ also introduced. The Floquet operator is then given byF̂ϵ=Ûϵ(T)=exp (iĤ2T2)exp (iĤ1T1)(2)

Crucially, the Floquet operator depends on the error ϵ, which determines how the time crystal will melt. In our work, we choose α = 1.51, J0T2 = 0.06 with a disorder strength WT2 = π, 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 Ĥ1.

Fig. 1 Obtaining the associated graph of a 2T-DTC.

(A) On the left, diagram of the 2T-DTC dynamics with no rotation error, ϵ. The initial state (∣Ψ(0)〉 = ∣ ↑ ↑ ↑ … ↑ 〉), represented with the green arrows pointing up, is recovered after two periods of the driving protocol. From the first period, we obtain the unitary, U(T), that will be used as the Floquet operator, F̂ϵ=0, to derive the effective Hamiltonian, Ĥϵ=0,Teff. On the right, fidelity of evolving state of an n = 8 sites 2T-DTC against its initial state, F(t) = ∣〈Ψ(0)∣ Ψ(t)〉∣2, showing the 2T periodicity of the dynamics. (B) The effective Hamiltonian, Ĥϵ=0,Teff, represented as a tight-binding matrix. ℰi and Kij are the energies of the ∣i〉 configuration (see the main text for full description) and transition energy between configurations ∣i〉 and ∣j〉, respectively. The right panel shows the entries of the effective Hamiltonian matrix, where only the diagonal and antidiagonal entries are nonzero. (C) After applying the percolation rule to the effective Hamiltonian, we obtain an adjacency matrix, which is, in turn, represented as a graph with the nodes being each of the 2n configuration basis set of the Hilbert space. In the right, for no rotation error, the crystal order of the 2T-DTC can be observed as 2n−1 decoupled dimers. Here, the elements of each dimer are the configurations ∣i〉 and ∣2n − 1 − i〉, which are related by a global π rotation along the x axis. All the nodes with the same color have the same number of domain walls with the color map gradient going from dark blue (0 domain walls) to dark red (n − 1 domain walls).

Floquet graphs

Let us now introduce the concept of a Floquet graph as the graphical representation of the effective Hamiltonian, Ĥϵ,Teff (26), obtained as Ĥϵ,Teff=i/T log[F̂ϵ]. 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, Ĥϵ,Teff, is defined on a finite and discrete Hilbert space (dimension N) and can be written as a tight-binding model using the configuration basis states { ∣ i⟩} asĤϵ,Teff=ΣiEiii+Σi,jKijij=ΣsλsΦsΦs(3)where ℰi is the energy of configuration ∣i〉 and Kij is the transition energy between configurations ∣i〉 and ∣j〉 (see Fig. 1B). The Floquet states ∣Φs〉 are linear combinations of the configurations ∣i〉. In general, the quasienergies λs are complicated functions of the parameters ℰi and Kij. In this work, we consider a system with n spin-1/2 particles, and the configurations {∣i〉} with i = 0,1, …, 2n − 1 can be chosen as the N = 2n elements of the computational basis, i.e., the product states of the eigenstates of σlz. In general, any configuration can be written as i=(σ1+)j1(σn+)jn, where i = (j1j2jn)2 is the binary decomposition of the integer i and σl+l=l. By using this notation, the states are denoted as ∣0〉 = ∣ ↓ ↓ ⋯ ↓ 〉 and ∣2n − 1〉 = ∣ ↑ ↑ ⋯ ↑ 〉. The basis configurations have been chosen to satisfy the relation 2n1i=σ1xσ2xσnxi, such that the states ∣i〉 and ∣2n − 1 − i〉 are always linked by the global π rotation from Ĥ1 in Eq. 1 and will constitute the 2n nodes of the graph. In addition, such nodes will be linked according to a percolation rule (27), where a nonweighted edge is considered to be active between nodes i and j only if the conditionEjEi<Kij(4)is satisfied. This rule eliminates the off-resonant transitions, allowing one to effectively visualize the relevant transitions in the Hamiltonian (see the “Floquet theory: Percolation rule and graph structure” section in Methods for further details).

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 Ĥϵ,Teff derived from Eq. 2 can be represented as a graph using the percolation rule from Eq. 4. The structure of the graph when ϵ = 0 presents decoupled dimers (two connected nodes), a signature of 2T-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 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〉 and ∣2n − 1〉 defined above, there is an energy gap ℏπ/T between the quasienergies λ0 and λ2n−1. The corresponding Floquet states are given by Φ0(0+2n1)/2 and Φ2n1(02n1)/2, which are maximally entangled Greenberger-Horne-Zeilinger (GHZ) states. The same behavior is observed for other configurations, ∣i〉 and ∣2n − 1 − i〉, as all the paired Floquet states ∣Φi〉, ∣Φ2n − 1 − i〉 have an energy difference of ℏπ/T (9).

Melting of the 2T-DTC

Consider now the case where our system is initialized to ∣ψ(0)〉 = ∣2n − 1〉 = ∣ ↑ ↑ ⋯ ↑ 〉 (10). This state is a superposition of Floquet states, ψ(0)(Φ0Φ2n1)/2, and, in the graph, is represented as a dimer between the nodes ∣0〉 and ∣2n − 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 Ĥ1, see movie S1). In Fig. 2, we present the effect of moderate (ϵ = 0.012) and large (ϵ = 0.1) levels of rotation errors ϵ 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. 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.

Fig. 2 Melting of a 2T-DTC using Ĥϵ,Teff.

(A) Graph representation obtained from Ĥϵ,Teff with n = 8 sites and ϵ = 0.012. Nodes are color-coded according to the number of domain walls of the corresponding configuration (see color map). For moderate levels of error, the nodes attach to each other according to their number of domain walls. This can also be observed in the effective Hamiltonian matrix of (B) with the basis ordered in increasing number of domain walls and delimited by the colored squares. As the nodes start to cluster due to the presence of error, some nonzero off-diagonal terms appear in the center of the matrix. For this level of error, some dimers survive [see the top right corner of (A)], serving as good indication of the robustness of the system and meaning that the crystal order is still present. In (C), the degree distributions of the corresponding graph of (A) with different system sizes (n = 8–12) averaged over 100 realizations of disorder are shown in both linear and logarithmic scale (see the inset), which display heavy-tailed distributions. The distributions are fitted with a power law curve (solid lines in the inset), indicating the presence of large degree hub nodes. (D) Graph representation obtained from Ĥϵ,Teff with n = 8 sites and ϵ = 0.1 and (E) its associated effective Hamiltonian matrix. As we increase the error, the system forms a single cluster. This can be seen from the appearance of many new off-diagonal entries in the Hamiltonian matrix as well as the presence of a giant component in the graph. In this scenario, the time crystal has melted completely and no crystal order is left (no presence of isolated dimers). The degree distributions, shown in (F), also indicate that the heavy-tailed nature is destroyed and approximate to a normal distribution.

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 ∣2n − 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).

Fig. 3 Melting of a 2T-DTC using Ĥϵ,2Teff.

(A) Graph representation obtained from Ĥϵ,2Teff with n = 8 sites and a rotation error of ϵ = 0.012 with the nodes corresponding to the configuration basis set of the Hilbert space color-coded according to the number of domain walls of the corresponding state (see color map). For no or moderate levels of error, the nodes are decoupled, as shown in the bottom part of the graph, or forming clusters of small size. In (B), we present the effective Hamiltonian Ĥϵ,2Teff that presents mainly diagonal terms. For moderate levels of error, the nodes form small clusters. Again, they attach to each other according to their number of domain walls. In (C), the degree distributions of the corresponding graph of (A) with different system sizes (n = 8–12) averaged over 100 realizations of disorder are shown in both linear and logarithmic scale (see the inset), fitted with a Poisson distribution. This indicates that the corresponding graph is a random graph with low connectivity. (D) Graph representation obtained from Ĥϵ,2Teff with n = 8 and a rotation error of ϵ = 0.1. As we increase the error, the system forms a larger, highly connected cluster. (E) This can be seen from the increasing magnitude of the off-diagonal entries in the effective Hamiltonian matrix as well as the presence of a giant component in the graph. In this scenario, the time crystal has melted completely. (F) The degree distributions are well approximated to a Poisson distribution in the lower degree region, which is the characteristic of highly connected random networks.

To calculate the conserved quantities, let us consider the square of the Floquet operator of Eq. 2, where Ĥϵ,2Teff is the effective Hamiltonian now over two periods (Ĥϵ,2Teff=i/2T log [F̂ϵ2]). In that stroboscopic framework (2T), the Floquet operator for no error can be reduced to the following expressionF̂ϵ=02=exp (2iT2ΣlmJlmzσlzσmz)(5)

In the absence of error, ϵ = 0, the disorder cancels out exactly and the effective Hamiltonian reads Ĥϵ=0,2Teff=/2ΣlmJlmzσlzσmz. 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 N̂=Σl(1σlzσl+1z) and parity Π̂=lσlx such that [Ĥϵ=0,2Teff,N̂]=[Ĥϵ=0,2Teff,Π̂]=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 ∣ψ(0)〉 = ∣2n − 1〉, the system will remain stroboscopically within the symmetry multiplet with 2T quasienergy λ12T=/2ΣlmJlmz and zero domain walls Nˆ=0. As shown in Fig. 1C, in the absence of error, all the dimers conserve the number of domain walls (represented as a color map) and their quasienergies. The aforementioned initial state, however, breaks the spatial parity symmetry Π̂, which leads to a nonzero value of the correlation function σlz(τ)=ψ(0)σlz(τ)σlz(0)ψ(0)0, as observed in the experiment (10).

Breaking symmetries

For zero error, the 2n 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 Ĥϵ,2Teff is block diagonal (see Fig. 3B) with n subblocks given by the conservation of the number of domain walls. In the absence of error (ϵ = 0), and after 2T, one can find a classical quasienergy surface in phase space associated to the effective Hamiltonian Ĥϵ=0,2Teff. 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 an n-dimensional phase spaceHϵ=0,2Teff(θ1,θ2,,θn)=2ΣlmJlmzcos θl cos θm(6)where θl 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 ϵĤϵ,2Teff=T2TΣlmJlmzσlzσmzgϵT12TΣl[(cos (Bl2T2)+1)σlx+sin (Bl2T2)σly](7)

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 gϵ/4Σlsin (BlT)σly breaks the parity symmetry. The transverse field term gϵ/4Σl(cos (BlT)+1)σlx, 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 notable similarity with KAM theory in classical systems: A small perturbation can destroy integrals of motion (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)〉 = ∣2n − 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 Ĥϵ,Teff further captures the existence of the preferential attachment mechanism and the nontrivial melting process of the 2T-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.

Fig. 4 The average degree of the graph obtained from Ĥϵ,Teff with n = 8 sites, plotted against the error ϵ.

The nodes are distinguished by their number of domain walls and averaged individually for each value of domain walls, as well as averaging over 100 realizations of disorder. Notice that some of the curves exhibit a local maximum approximately between 0.01 < ϵ < 0.03. This is exactly the range in which we observe the scale-free behavior emerging in our system. Comparing each curve also confirms that the states with a lower number of domain walls are more robust against the error in a sense that they have lower degrees.

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 Ĥϵ,Teff and not in the one from two periods Ĥϵ,2Teff, which shows Poisson distributions, typical for random networks (18, 31). We observe that a small error has little effect to the effective Hamiltonian Ĥϵ,2Teff (as shown in Fig. 3B). This can be understood from the fact that each π rotation changes the direction of the rotation around the z axis caused by the Bz term. Therefore, at time 2T, the effect of the Bz 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 Ĥϵ,Teff as shown in Fig. 2B. 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 ϵ.


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


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 Ei and Ej 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Ĥ(t)={Ĥ1g(1ϵ)(ij+ji)0<t<T1Ĥ2Eiii+EjjjT1<t<T2(8)

In general, the transitions occur when the driving is resonant with the energy level spacing between two configurations, i.e., when EjEi = ℏ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 Ei and Ej, respectively, exhibit an anticrossing (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Ĥϵ,Teff=Eiii+Ejjj+Kijij+Kjiji=(9)λs1Φs1Φs1+λs2Φs2Φs2(10)

From this, we can see that the energy gap for the quasienergies scales as δλ=λs1λs2=(EiEj)2+(Kij)2. In addition, we can derive the expressions for the eigenstates ∣Φs1〉 = cos θiji〉 + sin θijj〉 and ∣Φs2〉 = − sin θiji〉 + cos θijj〉, where cos θij = ∣ ℰi − ℰj ∣ /δλ and sin θij = ∣ Kij ∣ /δλ. The percolation rule ∣ℰj − ℰi∣ < ∣ Kij∣ from Eq. 4 implies that when the link between the ith and jth 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 ∣ > ∣ Kij∣, the Floquet states are localized and no link between the ith and jth 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 F̂=Û(T), which is the evolution operator in a period T of the drive. Let us begin by defining the discrete transformationĜF(z)=Σl=0F̂neinz=11F̂eiz(11)

The operator ĜF(z) now plays the role of the resolvent of the Floquet operator F̂. Next, we can express the Floquet operator in terms of the effective Hamiltonian, as F̂=eiĤϵ,TeffT, where Ĥϵ,Teff=Ĥ0+V̂ with Ĥ0=ΣiEiii and V̂=Σi,jKijij.

Although Eq. 11 looks different from the resolvent for undriven systems (34), they share the same structure if we consider the expansionĜF(z)=Σm=0Bm(iT)m1m!(zĤϵ,Teff)m1iT(zĤϵ,Teff)+B1+(12)in terms of the generating function for the Bernoulli numbers Bm (35).

The next step is to use the perturbative expansion of Eq. 12 in the hopping term V̂1zĤϵ,Teff=1zĤ0+(1zĤ0)V̂(1zĤ0)+(13)

This allows us to calculate the transition probability between two configurations ∣i〉 and ∣j〉, as followsiĜF(z)jiδi,jT(zEi)+iKijT(zEi)(zEj)+B1δi,j+(14)where δ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 ∣Kij ∣ / ∣ℰi − ℰj∣. This resembles the locator expansion for the undriven Anderson model. The main difference here is that that ∣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 ∣ < ∣ Kij∣ 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 ∣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 ∣ > ∣ Kij∣, the states remain localized.

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) ∝ 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 = kmin 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 kkcutoff 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 kcutoff > 0, and the kmin is chosen at the kcutoff where the Kolmogorov-Smirnov distance is minimized. Once the kmin 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 material for this article is available at

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.


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.

Stay Connected to Science Advances

Navigate This Article