Simulating long-range hopping with periodically-driven superconducting qubits

Quantum computers are a leading platform for the simulation of many-body physics. This task has been recently facilitated by the possibility to program directly the time-dependent pulses sent to the computer. Here, we use this feature to simulate quantum lattice models with long-range hopping. Our approach is based on an exact mapping between periodically driven quantum systems and one-dimensional lattices in the synthetic Floquet direction. By engineering a periodic drive with a power-law spectrum, we simulate a lattice with long-range hopping, whose decay exponent is freely tunable. We propose and realize experimentally two protocols to probe the long tails of the Floquet eigenfunctions and to identify a scaling transition between weak and strong long-range couplings. Our work offers a useful benchmark of pulse engineering and opens the route towards quantum simulations of rich nonequilibrium effects.

A common assumption of many-body physics is that particles can interact only with their neighbors. Similarly, in canonical lattice models, particles are often assumed to hop only between neighboring sites. While generically correct in condensed matter physics, this assumption is not applicable to gravitational models, unscreened Coulomb interactions, and synthetic materials, where the couplings can decay as a power-law. Longrange couplings can lead to counter-intuitive effects, such as violations of thermodynamic identities [1][2][3][4][5] and anomalous mechanisms of information spreading [6][7][8].
For systems with power-law decaying couplings, the transition between long range and short range correlations occurs at α = d, where α is the decay power of the potential and d is the dimension [9][10][11]. This transition is further affected by disorder [12,13], thermal [14,15] and quantum [16] fluctuations, and competing orders [17,18], giving rise to a wide range of classical and quantum phases with peculiar scaling laws.
The experimental study of this transition requires a simulator where α is not set by the physical decay of elementary forces and can be tuned continuously. This requirement is partially fulfilled by trapped ions, where phonon-mediated interactions can be used to simulate long-range quantum spin models, and α can be tuned within a limited range [19,20]. Here, we propose and realize an alternative approach, based on periodically driven (Floquet) quantum models. In these systems, the frequency components of the wavefunctions give rise to an effective one-dimensional lattice and the harmonics of the drive correspond to arbitrary hopping terms [21]. Intuitively, this synthetic Floquet dimension describes the number of photons absorbed or emitted from the drive. In the theoretical limit of a classical drive, such dimension is infinite in both directions.
To realize our protocol, we use a feature of noisy intermediate-scale quantum (NISQ) computers that has been recently made available via the cloud, namely pulse engineering. This term refers to the possibility of feeding arbitrary time-dependent pulses to the device [22][23][24]. Pulse engineering is commonly used to characterize the qubits and to optimize the fidelity of a target gate, or of a set of gates [25][26][27][28][29][30][31][32][33][34][35][36][37]. Periodic drives have been recently used to create long-range couplings in the physical space [38] and to simulate a two dimensional topological insulator [39]. In this work, we use pulse engineering to simulate long-range couplings in the Floquet space.
Floquet theorem [40] guarantees that the time evolution of a system governed by a time-periodic Hamiltonian,Ĥ(t + τ ) =Ĥ(t), can be written as |ψ(t) = n c n e −iµnt |φ n (t) . Here, µ n are the quasienergies, and |φ n (t) are the Floquet functions. These functions are periodic in time, |φ n (t + τ ) = |φ n (t) , and can be expanded in a discrete Fourier series, where ω p = 2π/τ is the pump frequency. By integrating the Schrödinger equation over one period, one obtains a time independent linear system of equationsĤ F |φ n = µ n |φ n , whereĤ F is the Hamiltonian in an enlarged infinite dimensional Floquet space and |φ n is a vector in this space, with m|φ n = |φ n (m) . The central block of H F is given by whereĤ (m) = τ −1 τ 0Ĥ (t)e imωpt dt is the m th discrete Fourier component ofĤ(t) [41]. In the limit of large ω p , the off-diagonal terms ofĤ F can be treated perturbatively, giving rise to the well-known Magnus expansion. Note that Eq. (1) is invariant under the gauge transformation |φ n (m) → e iωpt |φ(m + 1) , and henceĤ F is invariant to discrete translations, m → m + 1,Ĥ F → H F + ω p . Textbook examples of Floquet systems often involve sinusoidal drives, whereĤ (m) is non zero only for m = ±1, andĤ F corresponds to a tilted lattice with nearestneighbor couplings. In this case, all the eigenstates ofĤ F are localized in Floquet space. In analogy to Bloch oscillations of electrons in an electric field, quantum particles oscillate in the Floquet dimension, leading to a periodic alternation of energy emission and absorption (Rabi oscillations) [42,43]. The situation does not change qualitatively when a finite number of harmonics is added. Hence, in most physical situations it is sufficient to focus only on a small segment of the Floquet space [44]. The opposite situation is encountered if the system is periodically driven by delta-functions in time, giving rise to celebrated kicked models [45], whereĤ (m) is independent on m. The HamiltonianĤ F corresponds to a model with all-to-all couplings and its eigenfunctions are completely delocalized in the Floquet space. In this case, a truncation of the Floquet space can lead to unphysical predictions [46].
Aiming at quantum simulations with superconducting qubits, we specifically consider the spin model where ω z and h x (t) are tunable fields, and σ α are Pauli matrices [47]. For h x (t) = h x (t + τ ) one obtains a Floquet Hamiltonian of the form of Eq.
This Hamiltonian is equivalent to a two-band Wannier-Stark ladder and is schematically shown in Fig. 1, where, for convenience, we have de-pictedĤ (m) for |m| ≤ 2 only. To study long-range hop-ping models, we introduce a time-dependent drive with a power-law frequency spectrum, This driving field interpolates between a kicked model at α = 0 and a weak sinusoidal drive for α 1. At intermediate α's, the Fourier transform of Eq. (4) corresponds to periodic kicks with a finite width. The resulting time evolution can be easily computed in two limiting cases: (i) For α = 0, the time evolution over one period is U (τ ) = U z U x , where U z = exp(−iω z τσ z /2) corresponds to the evolution between the kicks and U x = exp(−ih 0 τσ x /2) describes the kicks. Note that the same stroboscopic time evolution can be obtained by applying alternatively magnetic fields in the x and z direction for finite times. Such bang-bang [48,49] protocol and the kicked model are characterized by the same stroboscopic evolution (which is the evolution over an integer number of periods), but lead to a different micromotion (i.e. the intermediate evolution during the time periods). Due to this difference, these two protocols are mapped to different lattice models: the bang-bang protocol corresponds to a short-range hopping model, while the kicked protocol is mapped to a long-range model. (ii) For ω z = 0, the Hamiltonian corresponds to a time-dependent magnetic field in the x direction, whose time evolution is The driving field of Eq. (4) undergoes a scaling transition between short-range and long-range as a function of α. In general, a system is said to be short ranged (or long-ranged) if the integral over space of the coupling is finite (or infinite). In systems with short-range couplings the total energy is extensive (i.e. proportional to the size of the system), while for long-ranged couplings, the total energy grows faster than the system size. See, for example, Refs. [16,17] for the critical properties of this transition, in the context of particles with power-law interactions. In our mapping to a one-dimensional lattice, the Fourier component h x . In the case of the power-law spectrum introduced in Eq. (4), this quantity diverges for α < 1. Hence, the system under consideration undergoes a transition between short range and long range at α = 1.
In our mapping of periodically driven dynamics to long-range hopping, h sum corresponds to the effective drive strength at integer multiples of τ , h sum = h x (nτ ) with integer n, which is always finite in a real experiment. To address this problem, we introduce a physical cutoff M , by assuming that h  Our experimental setup is a single-qubit quantum computer, available on the cloud, the IBM Quantum processor [50] ibmq_armonk v2.4.0 which is one of the IBM Quantum Canary processors. The typical physical parameters that we use are ω z = 2π × 250 kHz, h sum = 2π ×1.2 MHz, and ω p = 2π ×300 kHz, or τ = 3.3 µs. This choice enables us to have δt τ min[T 1 , T 2 ], where δt = 0.000 22 µs is the smallest programmable time step and T 1 = 160 µs and T 2 = 280 µs are the decay and coherence times. In our experiment, we prepare the qubit in the |↓ state, set ω z to a fixed value, apply the drive h x (t) for time t, measure the expectation value s z (t) = σ z over 8192 shots, and repeat the procedure for 740 time steps between t = 0 and T = 5τ , see Fig. 2(b).
In the following, we propose and realize two complementary methods to experimentally probe the long-range nature of the model. The first method is based on the power-law scalings of the eigenstates of the enlarged Floquet HamiltonianĤ F . In the absence of a cutoff, the Floquet functions are composed of a central peak at m = 0 and long power-law tails, |φ n (m)| 2 ∼ |m| −x , see Fig. 3(a). The value of x can be analytically derived by the following scaling argument: If we consider the product of the Note that for all α ≥ 0, one has m |φ n (m)| 2 < ∞, ensuring that the Floquet functions |φ n (t) are normalizable. In the presence of a cutoff M , the Floquet functions decay exponentially for m > M , see Fig. 3(b). To preserve the normalization condition, the weight lost at m M appears as a pronounced peak at m ≈ M . This peak indicates that the Floquet functions are affected by the cutoff only at distances m ∼ M and appears sharp on a logarithmic scale. Its integrated weight is proportional to m>M m −2−2α ≈ 1/M 1+2α and rapidly decreases with α. This peak can be experimentally probed by studying the Fourier components The key feature that we observe is a pronounced peak at m = M . This peak is washed out as α increases and becomes invisible at α > 1, probing the scaling transition at α = 1.
Our second method to probe the scaling of the Floquet functions relies on the relation between the time derivatives of physical observables and the moments of |φ n (m) . According to Eq. (1), the q th time derivative of the expectation value of an operatorÔ in the state |φ n (t) , ( d/dt ) q Ô (t) n , equals to A rigorous bound to this quantity can be obtained by assuming thatÔ satisfies | φ n |Ô|φ m | ≤ Ô , where Ô is the operator norm [51]. This requirement is generically satisfied by operators that act on a finite Hilbert space, such as spin operators, and it does not apply to some operators of a infinite Hilbert space, such as the position or momentum operators of an oscillator, whose matrix elements can be arbitrary large. In this case, Using the scaling relation of Eq. (5), one can show that the right hand side of Eq. (6) scales as M q−α for large M , see appendix B. This expression is finite for q < α and diverges for q ≥ α. In particular, for the kicked Hamiltonian (α = 0) the q = 0 series diverges, highlighting the above-mentioned truncation problem. For all α ≤ 1, the q = 1 series diverges, leading to a divergence of the time derivative d/dt Ô (t) . For larger values of α, the q = 1 series is convergent, but higher-order series (and, hence, higher order derivatives) are infinite.
Hence, in our model the time derivatives of physical observables can serve as order parameters of the scaling transition, in analogy to the role of two-point correlation functions for symmetry breaking phase transitions. Figure 5 compares the experimental measurement and numerical calculation of the time derivative of s z (t). The scaling behavior changes drastically at the transition between long-range and short-range couplings (α = 1). For α ≤ 1, ds z /dt (t = τ ) diverges with increasing M , while for α > 1 it saturates to a finite value, in agreement with our scaling argument (see also appendix B).
In conclusion, we used Floquet engineering to simulate a one dimensional lattice with power-law hopping. We studied the scaling properties of the Floquet eigenstates and determined the effects of the long tails on the expectation values of physical observables and their time derivatives. By realizing this model on a quantum computer, we demonstrated the experimental capability of controlling and measuring a large number (M = 30) of harmonics. We were able to probe the long tails of the Floquet functions, both in real time and in Fourier space, and to study their dependence on the exponent α.
Our experimental observations probe the scaling transition between long-range and short-range couplings.
In recent years, superconducting circuits were used to simulate interesting nonequilibrium effects [52,53], including many-body localization of photons [54] and quench dynamics in one [55] and two [56] dimensions. Our work adds an important tool to these simulators, namely pulse engineering for the realization of synthetic Floquet dimensions. Pulse engineering allows us to both realize Floquet Hamiltonians with arbitrary couplings, and to measure the dynamics within each time period, the so called micromotion. By measuring the time harmonic decomposition of the Floquet function, and confirm their scaling properties. Our approach goes beyond previous studies based on stroboscopic measurements, which probe the Floquet quasienergies, but are insensitive to the associated eigenfunctions. By implementing our protocol on a multi-qubit quantum computer, it will be possible to explore equilibrium and nonequilibrium phase transitions with long-range couplings.

ACKNOWLEDGMENTS
Data availability All MATLAB, Julia, Python, and QISKIT codes used to generate the experimental and theoretical data presented in this article are made freely available online in [57].
Acknowledgments We thank Eli Arbel, Angelo Russomanno, and Alessandro Silva for useful discussions. We thank the support received through the IBM's QISKIT Slack channels. MMR and EGDT are supported by the Israel Science Foundation grants No. 151/19 and 154/19. We acknowledge use of the IBM Q for this work. We acknowledge the access to advanced services provided by the IBM Quantum Researchers Program. The views expressed are those of the authors, and do not reflect the official policy or position of IBM or the IBM Quantum team.

Appendix A: Quasienergies crossings
In this section we provide a brief analysis of the first Floquet-Brillouin zone of the periodically driven system with long-range interactions with the Floquet Hamiltonian (see Eqs. (2) and (4) in the main text),Ĥ F , and look for any spectral signatures related to the scaling transition. Because we are dealing with a two-level system, the folded spectrum has just two quasienergies. In the extended representation, this pair of quasienergies is repeated every integer multiple of the pump frequency, ω p . Hence, we cannot observe a clustering of states, which is associated with a phase transition, but at most, a level crossing.
We numerically calculated the quasienergies in the first Floquet-Brillouin zone. Typical results of the quasienergies in the folded zone are shown in Fig. 7(a)-(c), for dif-ferent values of h 0 . We, indeed, find that the quasienergies can cross leading to exceptional points of the system. These crossings generically do not occur at α = 1 and are not related to the transition between long-range and short range couplings discussed in this work. This section provides a deeper analysis of theoretical predictions related to Fig. 5 in the main text. In Eq. (6), we have obtained an upper bound for any derivative of generic observables [58],Ô(t). We will now show how to use the scaling law obtained in Eq. (5) in Eq. (6). We begin by remembering that the eigenfunctions scales as We now separate the sum over m and m in Eq. (6) to two sums, one in which m or m equals zero and another where m = 0 and m = 0 using the scaling law of Eq. (B1) we can rewrite Eq. (6) as we will examine the scaling of the integral, so we will ignore any numerical constants for the sake of brevity At this point we need to preform the integral carefully. If q = 2α the integral has a logarithmic behavior. Otherwise, if q > 2α, the integral scales as M q−2α or saturates to a finite values of order ine if q < 2α. Hence, we find . (B7) Now we will focus on the sum in which either m or m equals zero. As before, we convert the sum into an integral Just like before, we need to preform the integral carefully. If q = α the integral has a logarithmic behavior. Otherwise, if q > α the integral scales as M q−α , or scales as 1 if q < α. x | . This quantity converges for all α > 1, and it can diverge only for α ≤ 1.
Due to physical limitations, in the experiment we were only able to realize our model only up to M = 30 (with accurate results). To accurately describe the scaling transition, we plot the theoretical value of | ds z /dt |, as a function of M in Fig. 7. As expected, we observe that this quantity diverges as M 1−α for α < 1 and remains finite for α > 1. At the scaling transition, for α = 1, the time derivative diverges logarithmically.