Simulating complex quantum networks with time crystals

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 (−i∫0TĤ(τ)dτ/ℏ), which is the evolution operator within one period from time t = 0 to time t = T (21–23), with T̂ being the time ordering operator. In Floquet theory, one wishes to solve the eigenvalue problem F̂∣Φs〉=e−iλsT/ℏ∣Φs〉. The eigenvectors ∣Φ_s⟩ are known as Floquet states and −ℏω/2 ≤ λ_s ≤ ℏω/2 are the quasienergies. At discrete times t_n = nT, the evolution can be understood in terms of an effective “time-independent” Hamiltonian Ĥeff such that F̂=e−iĤ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 = T₁ + T₂Ĥ(t)={Ĥ1≡ℏg(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, Jlmz≡J0/∣l−m∣α 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 2gT₁ = π 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, J₀T₂ = 0.06 with a disorder strength WT₂ = π, 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.

• Download high-res image
• Open in new tab
• Download Powerpoint

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=ΣiEi∣i〉〈i∣+Σi,jKij∣i〉〈j∣=Σsλs∣Φs〉〈Φs∣(3)where ℰ_i is the energy of configuration ∣i〉 and K_ij 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 K_ij. In this work, we consider a system with n spin-1/2 particles, and the configurations {∣i〉} with i = 0,1, …, 2ⁿ − 1 can be chosen as the N_ℋ = 2ⁿ 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 = (j₁j₂…j_n)₂ is the binary decomposition of the integer i and σl+∣↓〉l=∣↑〉l. By using this notation, the states are denoted as ∣0〉 = ∣ ↓ ↓ ⋯ ↓ 〉 and ∣2ⁿ − 1〉 = ∣ ↑ ↑ ⋯ ↑ 〉. The basis configurations have been chosen to satisfy the relation ∣2n−1−i〉=σ1xσ2x⋯σnx∣i〉, such that the states ∣i〉 and ∣2ⁿ − 1 − i〉 are always linked by the global π rotation from Ĥ1 in Eq. 1 and will constitute the 2ⁿ 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 condition∣Ej−Ei∣<∣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 ∣2ⁿ − 1〉 defined above, there is an energy gap ℏπ/T between the quasienergies λ₀ and λ_2ⁿ−1. The corresponding Floquet states are given by ∣Φ0〉≈(∣0〉+∣2n−1〉)/2 and ∣Φ2n−1〉≈(∣0〉−∣2n−1〉)/2, which are maximally entangled Greenberger-Horne-Zeilinger (GHZ) states. The same behavior is observed for other configurations, ∣i〉 and ∣2ⁿ − 1 − i〉, as all the paired Floquet states ∣Φ_i〉, ∣Φ_{2ⁿ − 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)〉 = ∣2ⁿ − 1〉 = ∣ ↑ ↑ ⋯ ↑ 〉 (10). This state is a superposition of Floquet states, ∣ψ(0)〉≈(∣Φ0〉−∣Φ2n−1〉)/2, and, in the graph, is represented as a dimer between the nodes ∣0〉 and ∣2ⁿ − 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.

• Download high-res image
• Open in new tab
• Download Powerpoint

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ⁿ − 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).

• Download high-res image
• Open in new tab
• Download Powerpoint

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)〉 = ∣2ⁿ − 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 2ⁿ 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σmz−ℏgϵ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)〉 = ∣2ⁿ − 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.

• Download high-res image
• Open in new tab
• Download Powerpoint

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 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 Ĥϵ,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 ϵ.