Measuring out-of-time-order correlators on a nuclear magnetic resonance quantum simulator

The idea of the out-of-time-order correlator (OTOC) has recently emerged in the study of both condensed matter systems and gravitational systems. It not only plays a key role in investigating the holographic duality between a strongly interacting quantum system and a gravitational system, but also diagnoses the chaotic behavior of many-body quantum systems and characterizes the information scrambling. Based on the OTOCs, three different concepts -- quantum chaos, holographic duality, and information scrambling -- are found to be intimately related to each other. Despite of its theoretical importance, the experimental measurement of the OTOC is quite challenging and so far there is no experimental measurement of the OTOC for local operators. Here we report the measurement of OTOCs of local operators for an Ising spin chain on a nuclear magnetic resonance quantum simulator. We observe that the OTOC behaves differently in the integrable and non-integrable cases. Based on the recent discovered relationship between OTOCs and the growth of entanglement entropy in the many-body system, we extract the entanglement entropy from the measured OTOCs, which clearly shows that the information entropy oscillates in time for integrable models and scrambles for non-intgrable models. With the measured OTOCs, we also obtain the experimental result of the butterfly velocity, which measures the speed of correlation propagation. Our experiment paves a way for experimentally studying quantum chaos, holographic duality, and information scrambling in many-body quantum systems with quantum simulators.


I. INTRODUCTION
The out-of-time-order correlator (OTOC), given by is proposed as a quantum generalization of a classical measure of chaotic behaviors [1,2]. HereĤ is the system Hamiltonian andB(t) = e iĤtB e −iĤt , and ... β denotes averaging over a thermal ensemble at temperature 1/β = k B T . For a many-body system with local opera-torsÂ andB, the exponential deviation from unity of a normalized OTOC, i.e. F (t) ∼ 1 − #e λLt , gives rise to the Lyapunov exponent λ L . Quite remarkably, it is found recently that the OTOC also emerges in a different system that seems unrelated to chaos, that is, the scattering of shock waves nearby the horizon of a black hole and the information scrambling there [3][4][5]. A Lyapunov exponent of λ L = 2π/β is found there. Later it is also found that the quantum correction from the string theory always makes the Lyapunov exponent smaller [5]. Thus it leads to a conjecture that 2π/β is an upper bound of the Lyaponuv exponent, which is later proved for generic quantum systems [6]. This is a profound theoretical result. If a quantum system is exactly holographic dual to a black hole, its Lyapunov exponent will saturate the bound; and a more nontrivial speculation is that if the Lyapunov exponent of a quantum system saturates the bound, it will possess a holographic dual to a gravity model with a black hole. A concrete quantum mechanics model, now known as the Sachdev-Ye-Kitaev model, has been shown to fulfill this conjecture [2,7,8]. This establishes a profound connection between the existence of holographic duality and the chaotic behavior in many-body quantum systems [9]. perature limit (i.e. β = 0), connection between the OTOC and the growth of entanglement entropy in quantum many-body systems has also been discovered quite recently [10,11]. The OTOC can also characterize many-body localized phases, which are not even thermalized [10,[12][13][14][15].
Despite of the significance of the OTOC revealed by recent theories, experimental measurement of the OTOC remains challenging. First of all, unlike the normal timeordered correlators, the OTOC cannot be related to conventional spectroscopy measurements, such as ARPES, neutron scattering, through the linear response theory. Secondly, direct simulation of this correlator requires the backward evolution in time, that is, the ability of completely reverse the Hamiltonian, which is extremely challenging. One experimental approach closely related to time-reversal of quantum systems is the echo technique [16], and the echo has been studied extensively for both non-interacting particle systems and many-body systems to characterize the stability of quantum evolution in the presence of perturbations [17][18][19] and the physics is already quite close to OTOC. Recently it has been proposed that the OTOC can be measured using echo techniques [20]. In addition, there also exists several other theoretical proposals based on the interferometric approaches [21][22][23]. However, none of them has been experimentally implemented so far.
Here, we adopt a different approach to measure the OTOC. To make our approach work, some extent of "local control" is required. A universal quantum computer fulfills this need with having a "full local control" of the system -that is, a universal set of local evolutions can be realized, and this set of local evolutions can build up any unitary evolution of the many-body system, both forward and backward evolution in time. That is to say, we shall use a quantum computer to perform the measurement of the OTOC. In fact, historically, one of the key motivations to develop quantum computers is to simulate the dynamics of many-body quantum systems [24], and quantum simulation of many-body dynamics has been theoretically shown to be efficient with practical algorithms proposed [25]. Here the quantum computer we use is liquid-state NMR with molecules. In this work, we report measurements of OTOCs on a NMR quantum simulator. We should stress that, on one hand, our approach is universal and can be applied to any system that has full local quantum control, including superconducting qubit and trapped-ion; on the other hand, this experiment is currently limited to a small size not because of our scheme but because of the scalability issue of quantum computer.

II. NMR QUANTUM SIMULATION OF THE OTOC
The system we will simulate is an Ising spin chain model, whose Hamiltonian is written aŝ whereσ x,y,z i are Pauli matrices on the i-site. The parameter values g = 1, h = 0 correspond to the traverse field Ising model, where the system is integrable. The system is non-integrable whenever both g and h are nonzero. We simulate the dynamics governed by the system HamiltonianĤ, and measure the OTOCs of operators that are initially acting on different local sites. The time dynamics of the OTOCs are observed, from which entanglement entropy of the system and butterfly velocities of the chaotic systems are extracted.

A. Physical System
The physical system to perform the quantum simulation is the ensemble of nuclear spins provided by Iodotrifluroethylene (C 2 F 3 I) which is dissolved in d-chloroform, see Fig. 1(a) for the sample's molecular structure. For this molecule, the 13 C nucleus and the three 19 F nuclei ( 19 F 1 , 19 F 2 and 19 F 3 ) constitute a four-qubit quantum simulator. Each nucleus corresponds to a spin site of the Ising chain, as shown in Fig. 1(b). In experiment, the sample is placed in a static magnetic field alongẑ direction, resulting in the following form of system Hamilto- HereR = 1,Rx(−π/2),Ry(π/2) forÂ =σ z 1 ,σ y 1 ,σ x 1 , respectively.
where ω 0i /2π is the Larmor frequency of spin i, J ij is the coupling strength between spins i and j. The values of these system parameters are given in Appendix A. The system is controlled by radio-frequency (r.f.) pulses, and the corresponding control Hamiltonian goeŝ where ω 1 (t) and φ(t) denote the amplitude and the emission phase of the r.f. field respectively. The control pulse shape can be elaborately monitored to realize desired dynamic evolution. Actually, such a system has been demonstrated complete controllability [26], which guarantes that arbitrary system evolution can be implemented on it. Our experiments were carried out on a Bruker AV-400MHz spectrometer (9.4 T) at temperature T = 305 K.

B. Experimental Procedure
As schematically illustrated in Fig. 1(c), here we focus on β = 0 case and measuring OTOC mainly consists of the following parts.
1. Initial state preparation. This step aims at preparing an initial state with density matrixρ i ∝Â =σ α 1 , α = x, y or z.
1.1. The natural system is originally in the thermal equilibrium stateρ eq populated according to the Boltzmann distribution. In high-temperature approximation, ρ eq ≈ 1/2 4 (1 + 4 i=1 iσ z i ), where 1 is the identity and i ∼ 10 −5 denotes the equilibrium polarization of spin i. Because there is no observable and unitary dynamical effect on 1, effectively we writeρ eq = 4 i=1 iσ z i . 1.2. We engineer the system fromρ eq intoρ 0 = σ z 1 . This is accomplished in two steps: first to remove the polarizations of the spins except for that of F 2 by using selective saturation pulses, and then to transfer the polarization from F 2 to 13 C. Details of the method are described in Appendix B.
1.3. For initial stateρ 0 with α = x, y, we need to further rotate spin at site-1 by π/2 pulse around y or −x axes, respectively.
2. Implementing unitary evolution ofÛ (t) = e iĤtB e −iĤt . The key point is that according to the Trotter formula [25], the time evolution e −iĤt of the Ising spin chain of Eq. (2) can be approximately simulated through the decomposition e −iĤmτ ≈ e −iĤxτ /2 e −iĤzτ /2 e −iĤzzτ e −iĤzτ /2 e −iĤxτ /2 m for small enough τ . Here the dynamics is divided into m pieces with t = mτ , and Each propagator inside the bracket of Eq. (5) corresponds to either single-spin operation or coupled two-spin operation, and can be implemented through manipulatinĝ H NMR with r.f. controlĤ rf : single-spin operation terms are global rotations around x or z axis, which can be easily done through hard pulses; two-spin operation term e −iĤzzτ can be generated through some suitably designed pulse sequence based on the NMR refocusing techniques [27]. More details of the method are described in Appendix B. The reversal of Ising dynamics e iĤt can be done in the similar manner. Note in the case considered here,B is a local unitary operator on the site-N spin and B =σ γ N with γ = x, y, z that can be implemented by a selective r.f. π pulse on the site-N spin. Hence, for any given t, the total unitary evolution e iĤtB e −iĤt can be simulated. 3.
Readout. The OTOC is obtained by measuring the expectation value of the observableÔ = e iĤtB e −iĤtÂ e iĤtB e −iĤtÂ . For the infinite temperature β = 0, the equilibrium state of the many-body systemĤ is the maximally mixed state 1/2 4 . Since whenB is unitary,Û (t)ρ 0Û † (t) is a density matrix ρ(t) evolved from ρ 0 byÛ (t), as simulated in step 2. Finally Ô β=0 becomes measuring the expectation value ofÂ under ρ(t). Because that NMR detection is performed on a bulk ensemble of molecules, readout is an ensemble-averaged macrosopic measurement. When the system is prepared at state ρ(t), the expectation value of A can then be directly obtained from the spectrum. See Appendix B for details.

C. Results of OTOC
Two sets of typical experimental results of the OTOC at β = 0 are shown in Fig. 2. Here we normalize the OTOC by B † (0)B(0) Â † (0)Â(0) , and becauseÂ and B † commute at t = 0, the initial value of this normalized OTOC is unity. The experimental data (red points) agree very well with the theoretical results (blue curves). The sources of experimental errors include imperfections in state preparation, control inaccuracy, and decoherence. See Appendix C for more details. We also measure OTOC for other operators (Â =σ α 1 ,B =σ γ 4 with α, γ = x, y, z) and they all behave similarly. The experimental results are put in Appendix B. In both the integrable case (the first column in Fig. 2) and the non-integrable cases (the second and the third columns in Fig. 2), the early time behaviors look similar. That is, the OTOC starts to deviate from unity after a certain time (for the unit of time t, See Appendix D for details.). However, the long time behaviors are very different between the integrable and non-integrable cases. In the integrable case, after the decreasing period, the OTOC revives and recovers unity. This reflects that the system has well-defined quasi-particle. And there exists extensive number of integral of motions, which is related to the fact that an integrable system does not thermalize. While in the non-integrable case, the OTOC decreases to a small value and oscillates, which will not revive back to unity in a practical time scale. This relates to the fact that the information does scramble in a non-integrable system [11].

III. ENTROPY DYNAMICS
To better illustrate the different behaviors of the information dynamics in the two cases of integrable and non-integrable systems, we reconstruct the entanglement entropy of a subsystem from the measured OTOCs. Entanglement entropy has become an important quantity not only for quantum information processing, but also for describing a quantum many-body system, such as quantum phase transition, topological order and thermalization. However, measuring entanglement entropy is always challenging [28].
OTOC opens a new door for entropy measurement. An equivalence relationship between OTOCs at equilibrium and the growth of the 2nd Rényi entropy after a quench has recently been established [10], which states that In the left-hand side of Eq. (8), S A is the 2nd Rényi entropy of the subsystem A, after the system is quenched by an operatorÔ at time t = 0. That is, S In our experiment, we choose the quench operatorÔ ∝ (1+σ x 1 ) at the first site, and we take the first three sites as the subsystem A and the fourth site as the subsystem B, as marked in Fig. 1(b). In this setting, S A measures how much the quench operation induces additional correlation between the subsystems A and B.
The results of 2nd Rényi entropy S A are shown in Fig. 3. At short time, all three curves start to grow sig- A after a quench. A quench operator (1 +σ x 1 ) (up to a normalization factor) is applied to the system at t = 0, and the entropy is measured by tracing out the fourth site as the subsystem B. Different colors correspond to different parameters of g and h in the Ising spin model. The points are experimental data, the curves are theoretical calculations. nificantly after certain time. This demonstrates that it takes certain time for the perturbation applied at the first site to propagate to the subsystem B at the fourth site (see the discussion of butterfly velocity below). Then, for all three cases, S A s grow roughly linearly in time. This indicates that the extra information caused by the initial quench starts to scramble between subsystems A and B. The differences lie in the long-time regime. For the integrable model, the S A oscillates back to around its initial value after some time, which means that this extra information moves back to the subsystem A around that time window. As a comparison, such a large amplitude oscillation does not occur for the two non-integrable cases and the S A s saturate after growing. This supports the physical picture that the local information moves around in the integrable model, while it scrambles in the nonintegrable models [11].

IV. THE BUTTERFLY VELOCITY
The OTOC also provides a tool to determine the speed for correlation propagating. At t = 0,Â andB commute with each other since they are operators at different sites. As time grows, the higher order terms in the Baker-Campbell-Hausdorff formulâ becomes more and more important and some terms fail to commute withÂ, at which the normalized OTOC starts to drop. Thus, the larger the distance between sites for A andB, the later time the OTOC starts deviating from unity. In general, the OTOC behaves as where a and b are two non-universal constants, |x| denotes the distance between two operators. Here v B defines the butterfly velocity [5,11,[29][30][31]. It quantifies the speed of a local operator growth in time and defines a light cone for chaos, which is also related to the Lieb-Robinson bound [31,32]. In our experiment, we fixÂ at the first site, and movê B from the fourth site to the third site, and to the second site. From the experimental data, we can phenomenologically determine a characteristic time t d for the onset of chaos in each OTOC, i.e. the time that the OTOC starts departing from unity. By comparing the three different OTOCs in Fig. 4, it is clear that the closer the distance betweenÂ andB, the smaller t d . In the insets of Fig.  4(a) and (b), we plot t d as a function of the distance, and extract the butterfly velocity from the slope. We find that, for OTOC withÂ =σ z 1 andB =σ x i , v B = 2.10; and for OTOC withÂ =σ y 1 andB =σ z i , v B = 2.22. The butterfly velocity is nearly independent of the choice of local operators, which is kind of manifestation of the chaotic behaviour of the system.

V. OUTLOOK
OTOC provides a faithful reflection of the information scrambling and chaotic behaviour of quantum many-body systems. It goes beyond the normal order correlators studied in linear response theory, which only capture the thermalization behaviour of the system. Measuring the OTOC functions can reveal how quantum entanglement and information scrambles across all of the degrees of freedom in a system. In the future it will be possible to simulate more sophisticated systems that may possess holographic duality, with larger size and different β, to extract the corresponding Lyapunov exponents such that one can experimentally verify the connection between the upper bound of the Lyapunov exponent and the holographic duality.
We have used liquid-state NMR as a quantum simulator for the demonstration of OTOC measurement. NMR provides an excellent platform to benchmark the measurement ideas and techniques. Our work here represents a first and encouraging step towards further experimentally observing OTOCs on large-sized quantum systems. The present method can be readily translated to other controllable systems. For instance, in trapped-ion systems there have been realized high-fidelity execution of arbitrary control with up to five atomic ions [33]. Superconducting quantum circuits also allow for engineering on local qubits with errors at or below the threshold [34,35], hence offering another very promising experimental approach. Recent years' progress in these two quantum hardware platforms has been fast and astounding, particularly in the pursuit of fabrication of quantum computing architecture at large scale. It is reckoned that quantum simulators consisting of tens of or even hundreds of qubits are within reach in the near future [36][37][38][39]. Experimentalists will see the great opportunity of applying these technologies for studying quantum chaotic behaviors for much more complicated quantum many-body systems. Note added.-After finishing this work, we notice a related work [40], where OTOCs are measured in a trapped ion quantum magnet. 13  We use Iodotrifluroethylene dissolved in dchloroform [41].
The system Hamiltonian is given byĤ where ω 0i /2π is the Larmor frequency of spin i, J ij are the coupling strength between spins i and j. The values of parameters ω 0i and J ij are given in Fig. 5.
Appendix B: Experimental Procedure

Initialization
The system is required to be initialized intoρ 0 ∝σ z 1 from the equilibrium stateρ eq . We first exploit the steady state effect when a relaxing nuclear spin system is subjected to multiple-pulse irradiation [42]. To implement this, we apply the periodic sequence [π 1,2,4 − d] to the system, where π 1,2,4 means simultaneous π rotations on the spins 13 C, F 1 and F 3 , and d is a time delay parameter to be adjusted, see the first part of the circuit shown in Fig. 6(a). To do π 1,2,4 , we use a pulse which is composed of three frequency components, each Hermite-180 shaped in 500 segments, with a duration of 1 ms. With increasing the number of applied cycles, under the joint effects of relaxation and π reversions, the equilibrium Zeeman magnetizations σ z 1,2,4 gradually decay to zero. Only the magnetizationσ z 3 is retained at last as it is the fixed point to the periodic driving. We adjust the time interval d between the π pulses to achieve the bestquality steady state. In experiment, we set d = 25 ms and after more than 500 cycles we found that the system was effectively steered into a steady stateρ ss ∝σ z 3 (in this sample, we did not see observable Overhauser enhancement). Next, with a SWAP operation we transfer the polarization from the high-sensitivity F 2 nucleus to the low-sensitivity 13 C nucleus. Using the method, we finally get an initial stateρ 0 ∝σ z 1 . The resulting experimental spectrum is shown in Fig. 6(c).  6. (a) Quantum circuit that measures the OTOCs. The first part aims at reset an arbitrary state to the desired initial state. Here the time interval between the π pulses is 25 ms, the number of cycles is l = 500, and Gz denotes z axis gradient pulse. (b) Sequences for implementing the dynamics of e −iĤzz τ (left) and e iĤzz τ (right). The refocusing circuits are designed to generate the right amount of coupled evolution. (c) 13 C experimental spectrum for equilibrium state (blue) and stateρ0 (red) after a readout pulseR 1 y (π/2). They are shown at the same scale for comparison.

Simulating time evolution of Ising spin chain
According to Eq. (5) of the main text, the key ingredient in simulating the evolution of Ising HamiltonianĤ is to implement Here, except for e −iĤzzτ , all other four terms are global rotation around x (and z) axis, which can be easily done through hard pulses. e −iĤzzτ can be generated by manipulating the natural physical HamiltonianĤ NMR with a suitable refocusing scheme [43]. The basic idea is to evolve the system with the J-term inĤ NMR and then to use spin echoes to engineer the evolution. That is to say, for instance, for theσ z iσ z j term, when a transverse π pulses is applied to reverse the polarization of one of the two spins, the evolution is also reserved. Hence by designing a suitable refocusing scheme, the dynamics of H zz and −Ĥ zz can be efficiently simulated.
Although general and efficient refocusing scheme exists for anyσ zσz -coupled evolution [27], for the present task it is possible to find a much simplified circuit construction. Fig. 6(b) shows our ideal circuits. Let O 1 and O 2 define the reference frequency for carbon and fluorine channel respectively. Consider the refocusing circuit ( Fig. 6(b), left) for implementing e −iĤzzτ , it automatically refocuses the fluorine spins and decouples the terms J 31 , J 41 and J 43 , and the evolution of other terms should fulfil the following requirements to yield the right amount of evolution: The solution to the above system of equations is given by O 1 = ω 01 /2π = 15480.0 Hz, t 1 = 0.004935τ , t 2 = 0.009870τ and τ 1 = 0.000534τ . As to the refocusing circuit for implementing e iĤzzτ , we found that it suffices to just make slight changes to the circuit for −e iĤzzτ , as shown in the figure, and one can then reverse the dynamics of all terms. Now, the whole network for implementing Ising dynamics is expressed in terms of single-spin rotations and evolution of J-terms inĤ NMR . In practices, each single spin rotation is realized through a selective r.f. pulse of Gaussian shape, with a duration of 0.5 to 1 ms. We then conduct a compilation procedure to the sequence of selective pulses to eliminate the control imperfections caused by off-resonance and coupling effects up to the first order [44,45]. To further improve the control performance, we employ the gradient ascent pulse engineering (GRAPE) technique [46] on the complied sequences. Because that compilation procedure has the capability of directly providing a good initial start for subsequent gradient iteration, the GRAPE searching quickly finds out high performance pulse controls for the desired propagators. The obtained shaped pulses for different set of Hamiltonian parameters (g, h) all have the numerical fidelities above 0.999, and have been optimized with practical control field inhomogeneity taken into consideration.
The Ising dynamics to be simulated is discretized into 20 steps, with each time step of duration τ = 0.35 ms. Choosing different operators forÂ andB, we have experimentally measured the corresponding OTOC. All the experimental results are given in Fig. 7. The theoretical trajectories are plotted for comparison. Although some discrepancies between the data and the simulations remain, the experimental results reflect very well how OTOCs behave differently in the integrable and chaotic cases.

Readout
All the observations are made on the probe spin 13 C. Because we use an unlabeled sample in real experiment, the molecules with a 13 C nucleus are present at a concentration of about 1%. The NMR signal in high field is obtained from the precessing transverse magnetization of the ensemble of molecules in the sample: As the precession frequencies of different spins are distinguishable, they can be individually detected, e.g., we obtained the measurements of σ x 1 and σ y 1 at the 13 C Larmor frequency. To measure σ z 1 , we need to apply a π/2 rotation alongŷ. By fitting the 13 C spectrum, the real part and imaginary parts of the peaks are extracted, which corresponds to σ x 1 and σ y 1 , respectively.

Appendix C: Experimental Error Analysis
The sources of experimental errors include imperfections in initial state preparation, infidelities of the GRAPE pulses, r.f. inhomogeneity, and decoherences. We make analysis to the data set of the caseÂ =σ x 1 ,B = σ y 4 to get an understanding on the role of each type of error sources. We calculated the standard deviations σ exp := 20 i=1 ( Â i exp − Â i th ) 2 /20 for the experimental data, which are presented in Table I.
We have run the initialization process for 50 times and found that the fluctuation of the initial state polarization ofρ 0 is around 3.40%. The fluctuation is due to (i) error in state preparation; (ii) error in spectrum fitting. The latter can be inferred from the signal-to-noise ratio of the spectrum, which is estimated to be ≈ 2.13%.
All the GRAPE pulses for implementing e −iĤτ and e iĤτ , are of fidelities above 0.999. On such precision level, if we assume no other sources of error and assume that the pulse generator ideally generates these pulses, then the experimental results should match the theoretical predictions almost perfectly.  and 19 F channel. To understand to what extent the r.f. field inhomogeneity may affect the experimental results, we calculate the deviation of the dynamics based on a simple inhomogeneity model. The model assumes that the output power discrepancy of the r.f. fields is uniformly distributed between ±3%. The simulated results σ err inhomo are shown in Table I. Another major source of error comes from decoherence effects. We compare the experimental data to a simple phenomenological error model, i.e., the system undergoes uncorrelated dephasing channel, parameterized with a set of phase flip error probabilities {p i } i=1,2,3,4 per evolution time step t 0 . The density matrixρ is then, at each evolution step, subjected to the composition of the error channels E i for each qubit [47] ρ where with p i = (1 − e −t0/T2,i )/2 (see Fig. 5 for the values of T 2,i ). The results are presented in Table I. The results indicate that, with decoherence effects taken into account, the discrepancy between theoretical and experimental data for g = 1, h = 0 is expected to be larger than that of the other two cases, consistent with the experiment data.
In summary, we conclude that r.f. inhomogeneity and decoherence effects are two major sources of errors.