Scrambling and thermalization in a diffusive quantum many-body system

Out-of-time ordered (OTO) correlation functions describe scrambling of information in correlated quantum matter. They are of particular interest in incoherent quantum systems lacking well defined quasi-particles. Thus far, it is largely elusive how OTO correlators spread in incoherent systems with diffusive transport governed by a few globally conserved quantities. Here, we study the dynamical response of such a system using high-performance matrix-product-operator techniques. Specifically, we consider the non-integrable, one-dimensional Bose-Hubbard model in the incoherent high-temperature regime. Our system exhibits diffusive dynamics in time-ordered correlators of globally conserved quantities, whereas OTO correlators display a ballistic, light-cone spreading of quantum information. The slowest process in the global thermalization of the system is thus diffusive, yet information spreading is not inhibited by such slow dynamics. We furthermore develop an experimentally feasible protocol to overcome some challenges faced by existing proposals and to probe time-ordered and OTO correlation functions. Our study opens new avenues for both the theoretical and experimental exploration of thermalization and information scrambling dynamics.

Dynamical correlations of many-body quantum systems provide direct information about many-body excitations [1], describe quantum phases and transitions [2], and characterize certain topological aspects [3,4]. The dynamical response of a many-body system to a local perturbation is obtained from a time ordered correlation function, Ŵ (t)V (0) , which describes the relaxation of the many-body system following the initial excitation by the operatorV that is then probed at later times byŴ . However, in general such time-ordered correlation functions cannot capture the spread of information across a quantum system, especially in a regime where quasiparticles are not well-defined.
Recently, it has been proposed that spreading or "scrambling" of quantum information across all the system's degrees of freedom can be characterized by out-of-time ordered (OTO) correlation functions: Ŵ † (t)V † (0)Ŵ (t)V (0) [5][6][7][8][9][10]. These correlation functions appear as the out-of-time ordered part of |[Ŵ (t),V (0)]| 2 and hence predict the growth of the squared commutator betweenŴ (t) andV (0). OTO correlators are thus capable of describing a quantum analogue of the butterfly effect in classical chaotic systems, which characterizes the spread of local excitations over the whole system. At short times, OTO correlators are expected to grow exponentially with a rate characterized by the Lyapunov exponent λ L . The Lyapunov exponent has been conjectured to be bounded by 0 ≤ λ L ≤ 2πT [9]. This bound is saturated in strongly coupled field theories with a gravity dual [6] and in disordered models describing a strange metal [7,11,12]. By contrast, λ L does not fully saturate the bound for a critical Fermi surface [13] and is parametrically smaller in Fermi liquids or other weakly coupled states [12,14,15].
Here, we study both time-ordered and OTO correlators in a diffusive many-body system by considering the concrete example of the non-integrable, one-dimensional Bose-Hubbard model. Thus far, it is a largely open question, how OTO cor-distance (i-j) correlators measure the scrambling of information across a quantum state. We compute OTO correlators Fij(t) = c † j (t)c † i cj(t)ci in the 1D Bose-Hubbard model at high temperature T = 4J for interactions U = J, chemical potential µ = 0, and system size L = 30. In the high temperature regime, well-defined quasiparticles cease to exist. However, the OTO correlator Fij exhibits a light-cone spreading of information. (b) The breakdown of well-defined quasiparticles is demonstrated by the one-particle Green's function Gij(t) = c † j (t)ci , which quickly decays to zero within τ J ∼ 0.6. The lifetime is thus shorter than the hopping rate, indicating a regime of incoherent transport. relators spread in diffusive systems with a few globally conserved quantities [10,13,15,16]. In our work, we study this question by performing matrix-product operator (MPO) based simulations of the Bose-Hubbard model at high temperatures, at which well defined quasi-particles cease to exist. We demonstrate that in this regime the time-ordered oneparticle correlation functions are strongly incoherent and feature rapidly decaying excitations, whereas the OTO correlators indeed describe the ballistic spreading of information across the quantum system (see Fig. 1). In contrast to the linear light-cone spreading of quantum information, the eventual global thermalization of the closed system takes parametrically longer, due to hydrodynamic power-laws resulting from globally conserved quantities. For example, we show that the local density correlation function decays as ∼ 1/ √ Dt, describing diffusion in one dimension with the corresponding diffusion constant D. Thus, the time scales associated with the spread of information and with global thermalization are different.
Despite their usefulness to characterize interacting manybody systems theoretically, it remains a challenge to experimentally measure such dynamical correlation functions in real space and time [17,18], as required to observe information spreading. Here, we propose generic experimental protocols to characterize both time-ordered and OTO correlators via local many-body interferometry. Our proposal to measure OTO correlators is unique because it overcomes some of the challenges that recently proposed protocols exhibit at finite temperatures and because it eliminates the scaling problems associated with global many-body interferometry [19,20]. Furthermore, our protocol does not require an ancillary atom to switch between different system Hamiltonians [21][22][23] and directly works with massive bosonic and fermionic particles (see also [24][25][26][27][28]). We show that our protocol not only enables the measurement of dynamical correlation functions but also rather generic static correlation functions (including off-diagonal ones), thus opening the way for a full state-tomography of many-body quantum states.

I. RESULTS
We study dynamical correlation functions of the onedimensional Bose-Hubbard model focusing mainly on the incoherent intermediate to high temperature regime. The Hamiltonian of the system is given bŷ where J is the tunneling matrix element, U the interaction strength, and µ the chemical potential. The bosonic creation (annihilation) operator on site i is denoted as c † i (c i ) and the local particle number operator isn i = c † i c i . At zero temperature and commensurate filling, the Bose-Hubbard model exhibits a quantum phase transition from a gapped Mott insulating phase with short range correlations at strong interactions to a compressible superfluid phase with power-law correlations at weak interactions [29]. At finite temperatures the system is a correlated, normal fluid. We compute the dynamical correlation functions at finitetemperature for systems up to L = 50 sites using MPO techniques. The presented results are evaluated for virtual bond dimension 200 to 400 and the local bosonic Hilbert space is truncated to three states, which is sufficient to render the system nonintegrable. The presented results are checked for convergence with respect to the MPO bond dimension and system size; see Methods III A for details on the numerical simulations.

A. Spread of Quantum Information
Recently, OTO correlation functions have been proposed as a useful diagnostic tool to quantify the dynamical spreading of quantum entanglement and quantum chaos in many-body systems. OTO correlators describe the growth of the commutator between two local operatorsŴ andV in time In a semiclassical picture, the commutator in Eq. (2) can be replaced by Poisson brackets. Then, for the choice of W = p j andV = p i , this quantity reduces to C(t) ∼ (∂p j (t)/∂q i (0)) 2 . Therefore, the correlation function C(t) describes the sensitivity of the time evolution and is expected to grow exponentially at short times ∼ exp[λ L t], with a rate λ L that resembles the Lyapunov exponent in classical chaotic dynamics. Rewriting these momenta and coordinates as combinations of creation and annihilation operators, Eq. (2) generically consists of OTO correlators of the form Below we mainly consider the quantum statistical average . . . = tr[ρ . . .] over an initial thermal state with weights distributed according to the Gibbs ensembleρ = e −Ĥ/T /Z where Z is the partition function and we set the Boltzmann constant k B to one. Alternatively, the average can also be performed with respect to an arbitrary initial state, for example a pure stateρ = |ψ 0 ψ 0 |. For thermalizing systems, it is then expected that an effective temperature is approached at late times which depends on the energy density imprinted on the system by the initial state [30][31][32].
OTO correlators F ij evaluated at comparatively high temperatures T = 4J, interactions U = J, and chemical potential µ = 0 are shown in Fig. 1 (a) as a function of time t and distance (i−j). Despite the high temperature, the OTO correlator F ij unveils a pronounced light-cone spreading of the information across the quantum state for |i − j| 7. For larger distances the light cone seems to grow super-ballistically, which we, however, attribute to the finite MPO bond dimension considered in the numerical simulations; see Methods III A. OTO correlators are in that respect challenging to simulate with MPO techniques, because they directly reflect the fast spreading of entanglement.
The OTO correlator F ij (t) should be contrasted to the timeordered single-particle Green's function which is shown in Fig. 1 (b). In the incoherent transport regime, where well-defined quasiparticles do not exist, the Green's function G ij (t) rapidly decays in time. Therefore, it is not capable of characterizing the spread of quantum information or entanglement across the state which is generically not linked to the transport of quasi-particles [33]. For the chosen parameters (U = J, µ = 0, T = 4J), we find that the quasiparticle lifetime is approximately τ J ∼ 0.6 and hence shorter than the microscopic hopping rate, which indicates incoherent transport. By contrast, the out-of-time ordered structure of F ij (t) reveals a well defined linear spread of quantum information despite the high temperature. We now characterize the OTO correlators F ij (t) in detail. To this end, we subtract n inj from the F ij (t) and consider its relative change: F r ij (t) = |F ij (t) − n inj |/ n inj . Examples for the reduced OTO correlator F r ij (t) are shown in Fig. 2 for interaction U = J and different temperatures T . The reduced OTO correlator F r ij (t) starts off at zero, forms the light-cone plateau, and approaches the steady-state value as an exponential.
From the light-cone spread of the OTO correlator, we extract two velocities [ Fig. 3 (a)]: (1) The light-cone velocity v lc , which we define by the space-time region where the reduced OTO correlators F r ij (t) surpasses a small threshold of 0.05% of its final value. (2) The butterfly velocity v b , which we define by the space-time region where the OTO correlator attains a large fraction (20%) of its final value. We find that v b does not significantly depend on this cutoff, as long as it is chosen to a sizeable fraction; see Methods The reduced OTO correlator F r ij (t) is expected to grow exponentially on a timescale set by the butterfly velocity vb with a rate that defines the Lyapunov exponent λL. In our system, the regime of exponential growth is restricted to a rather small time range, see also Fig. 8. Our data suggests that the Lyapunov exponent λL is parametrically smaller than the conjectured upper bound 2πT and increases slowly as the temperature T is lowered. The data is shown for interaction strength U = {1, 3, 9}J.
from below by the zero temperature Luttinger liquid velocity; see Fig. 3 (b). The butterfly velocity v b is systematically lower than v lc and is almost independent of temperature. The butterfly velocity determines the time scale t scr for scrambling information across the many-body quantum state which is linear in system size t scr ∼ L/v b . Based on results from holography, it has been argued in Ref. [16] that the light-cone and the butterfly velocity should be quite generally the same. This should be contrasted to our results for the Bose-Hubbard model and to a study of non-relativistic non-Fermi liquids [13]. In both cases the butterfly velocity has been found to be smaller than the light-cone velocity. For details on the analysis of the OTO correlators; see Methods III B.
We characterize the initial dynamics of the OTO correlator by its growth rate λ L around the space-time cone set by the butterfly velocity: and refer to this rate as Lyapunov exponent. However, we note that in contrast with predictions from strongly coupled field theories [9] or disordered SYK models [7], there appears to be no parametrically large regime of exponential growth in our model, see Fig. 8, which hints at more complex initial dynamics. Finding the analytic form for the initial growth over an extended space-time region remains an outstanding challenge. Alternatively, one can extract the growth rate Λ lc of the reduced OTO correlator around a space-time region set by the light-cone velocity v lc which we find to be larger than λ L ; see Methods III B. We show the Lyapunov exponent λ L as a function of temperature for different values of the interaction strength U and chemical potential µ = 0 in Fig. 4. It has been conjectured that the Lyapunov exponent is bounded by 2πT , which is the value it assumes in a strongly coupled field theory with a gravity dual [9]. In our system, λ L is parametrically lower than this bound and increases slowly when lowering the temperature. Moreover, we find that the depen- Conserved quantities restrict the approach of a closed quantum system to global equilibrium, thus, rendering global thermalization a slow process. In the Bose-Hubbard model the total particle number is conserved leading to diffusive power-law tails in the connected density correlator Dt, where D is the diffusion constant. By contrast, at low temperature T = J the diffusive regime has not yet been reached within the numerically accessible times and the correlations rather decay ballistically Cn(0, t) ∼ 1/t. The slow relaxation of the hydrodynamic modes leads to the global thermalization time scale tth ∼ L 2 /D that is parametrically larger than the scrambling time scale tscr ∼ L/vb of quantum information.
dence of the Lyapunov exponent on the interaction strength U is small with slightly larger values of λ L for intermediate interaction strength, U = 3J, which is in the vicinity of the quantum critical point.

B. Thermalization
Closed quantum systems approach their global equilibrium only very slowly, due to the slow evolution of observables that overlap with conserved quantities. In the Bose-Hubbard model (1), energy, lattice momentum, and total particle number are conserved. From hydrodynamics we infer that, for example, the conserved particle number leads to a diffusion equation of the density [34,35] ∂ t n − D∇ 2 n = 0, where D is the diffusion constant. The connected density cor- and in a hydrodynamic regime is expected to be of the form withC = dx C n (x, 0). Whereas local equilibrium is approached after a few scattering events, attaining global equilibrium is restricted due to the relaxation of such conserved quantities, which have to be transported over long distances. At comparatively low temperatures (T = J), Fig. 5 (a), the ballistic spread of sound modes dominates the dynamics of the connected density correlator in the numerically accessible time regime. However, at high temperatures (T = 10J), (b), the density correlator approaches diffusive transport after a few hopping scales and attains a finite value in the region between the sound modes. To be more quantitative, we study the local (x − x = 0) density correlation function. At high temperatures T 4J the local correlator exhibits a diffusive power-law decay C n (x = 0, t) ∼ 1/ √ Dt, (c). For this parameter set we extract the diffusion constant D = 9.79(1)a 2 J for T = 10J and D = 14.29(27)a 2 J for T = 4J, where a is the lattice spacing; see Tab. I. The decrease of the diffusion constant with increasing temperature is somewhat counterintuitive. We attribute this behavior to the fact that the calculations are performed in the grand-canonical ensemble. Hence the particle density depends on the temperature and, in particular, increases with temperature in the chosen parameter regime. We note that the connected density correlator does not exhibit pronounced hydrodynamic long time tails, which could result from higher order gradient corrections to the diffusion equation and mask the 1/

√
Dt decay. This seems to be a particular property of the density correlator, as we find at high-temperatures pronounced t −3/4 corrections in the energy-density correlation function (not shown), in agreement with Ref. [35].
It has been proposed that the diffusion constant is related to the butterfly velocity v b and the Lyapunov exponent λ L via [36][37][38][39], where 1/λ L is a bound for the local thermalization time in which the system is able to attain local equilibrium characterized by a local temperature and local chemical potential that varies between different regions in space. From our simulations, we obtain coefficients of the order Dλ L /v 2 b ∼ 5.5 for temperatures T 6J; see Tab. I, which seems to suggest a connection between the spread of information and local thermalization, as suggested by calculations for holographic matter. However, clearly global thermalization is a parametrically slower process than information scrambling and takes for systems of size L times of the order t th ∼ L 2 /D. Experimentally measuring OTO correlators (Sec. I C) and density correlators (App. A 3) will make it possible to further check these holographic predictions.

C. Measuring Dynamical Time-Ordered and Out-of-Time Ordered Correlators
We develop two generic interferometric protocols that measure time-ordered as well as OTO correlation functions for systems of bosons or fermions in an optical lattice. The first is based on globally interfering two many-body states and the second on local interference. The quantum interference of two copies of the many-body state is realized by local beam splitter operations. Variants of this approach have been proposed to study Rényi entropies [40][41][42] and have been demonstrated experimentally using a quantum gas microscope [19,20]. Both protocols that we propose consist only of elements which have already been used in experiments. The two protocols are complementary and each of them has its own advantages.
Global interferometry precisely yields the square modulus of the single-particle Green's functions G gl ij (t) = |G ij (t)| 2 and the OTO correlators F gl ij (t) = |F ij (t)| 2 for pure initial states (see Supplemental Information Sec. A 1). In a thermalizing system, effective finite temperatures can be obtained with the help of quenches from pure initial states. However, generic high temperature initial states are not accessible with this protocol. The measurement schemes proposed e.g. in Ref. [23] face similar challenges. The global interferometry protocol is furthermore limited to rather small system sizes, since the many-body wave function overlap has to be measured, which requires an extensive number of beam splitter operations.
These limitations are overcome by the second proposed protocol which uses local interferometry (see Supplemental Information Sec. A 2 for details on the implementation). In this protocol, only two local beam-splitter operations are required irrespective of the system size, and only local density differences between the two copies have to be measured. Furthermore, initial thermal density matrices can be studied as well.
This local approach yields a slightly amended two-point correlation function G loc However, we demonstrate in the Supplemental Material that these correlators carry much of the same information as the ones we discussed previously. Static correlation functions are accessible with the same techniques and an extension to higher order correlators is straight forward.

II. DISCUSSION
We studied time-ordered as well as out-of-time ordered correlation functions in the one-dimensional Bose-Hubbard model and suggest different protocols to experimentally access them. At high temperatures, well-defined quasi-particles cease to exist and the time-ordered Green's function decays within short times. However, the spread of information is not necessarily linked to the transport of quasi-particles. Our numerical results for the out-of-time ordered correlators clearly indicate the ballistic spread of information even at high temperatures where transport is incoherent. In our onedimensional system, this linear spread sets the timescale for scrambling information across the quantum state to be proportional to the system size. Moreover, the existence of conserved quantities in the Bose-Hubbard model leads to diffusive behavior of the corresponding time-ordered correlation functions. Global thermalization therefore scales with the square of the system size and takes parametrically longer than scrambling quantum information.
For future work, it would be on the one hand interesting to develop analytical predictions for the growth of OTO correlators, which in our numerics deviates significantly from the simple exponential growth obtained in strongly coupled field theories, or for the bounds that characterize the information propagation and Lyapunov exponents. On the other hand, the numerical study of out-of-time ordered correlators in other interacting many-body systems, including Fermi-Hubbard models, spin models, or continuum Lieb-Liniger models, could be beneficial. Taking such routes could help to advance our fundamental understanding of information scrambling, transport, and thermalization.
Experimental measurements of both time-ordered and outof-time ordered correlators will be eminent for the investigation of the dynamical properties of many-body systems. We proposed two different protocols to measure correlation functions, which can be either static, time ordered, or outof-time ordered. The schemes are respectively based on the global and local interference of two copies of the many-body state of interest. The required techniques have already been demonstrated in experiments with synthetic quantum matter. An extension of the described experimental schemes to twodimensional systems is conceivable as well and could provide a valuable perspective on the dynamics of many-body quantum systems.

III. METHODS
In this section, we discuss the numerical method with which the calculations have been performed as well as the procedure to obtain the velocities v lc , v b and the Lyapunov exponent λ L .
For the density correlations, we evaluate [47] n (t)n j (0) β where β is the inverse temperature and Z the partition function. We construct the MPO approximation of the two terms in the parentheses by first computing e −βĤ/2n andn j e −βĤ/2 , respectively, and then performing a real-time evolution up to t 2 and − t 2 . By exploiting the time translation invariance, n (t)n j (0) β = n (t/2)n j (−t/2) β , the maximum simulated time has effectively been reduced by a factor two, which in turn reduces the required virtual bond dimension of the MPO.
To evaluate e −βĤ/2 , we employ a second-order Suzuki-Trotter decomposition with imaginary time step ∆τ (typically ∆τ J = 0.025) after splitting the Hamiltonian into even and odd bonds, as described in Ref. 43. The real-time evolution proceeds by Liouville stepsÂ(t + ∆t) = e i∆tĤÂ (t)e −i∆tĤ . For each of the steps we combine a fourth-order partitioned Runge-Kutta method [49] with even-odd bond splitting of the Hamiltonian. As noted in Ref. 47, the Liouville time evolution has the advantage that the virtual bond dimension does not increase outside the space-time cone set by Lieb-Robinson-type bounds. The high order decomposition also allows for relatively large time steps (in our case ∆tJ = 0.125 or 0.25).
For the OTO correlators c † j (t)c † c j (t)c β , a regrouping analogous to Eq. (7) would lead to four terms inside the trace, such that a straightforward contraction to evaluate the trace becomes computationally very expensive. Instead, we evaluate and time-evolve both e −βĤ c † j and c j up to time t. Subsequent application of the site-local operators c † and c does not affect the virtual bond dimension in the MPO representation.
In our simulations, we restrict the local Hilbert space to three states due to computational limitations. Since the average particle number per site is approximately one, this restriction should not qualitatively affect the simulation results. Moreover, truncating the local Hilbert space to three states is sufficient to render the system non-integrable, which is crucial to observe the thermalization behavior studied in this work.
Since OTO correlators are closely linked to the spreading of entanglement, it is challenging to simulate them using MPO techniques. In Fig. 6 we compare the data obtained for the same simulation parameters but different maximal bond dimensions. The MPO bond dimension of 20 leads to an apparently super-ballistic growth of the light-cone around |i − j| ≈ 5, see Fig. 6 (a). Increasing the bond dimension shifts this numerical artifact to larger distances. It is however exponentially costly to reach full convergence of the OTO correlator. In the analysis of the numerical data we therefore only considered small distances, where we checked that increasing the bond dimension does not alter the correlators.

B. Data Analysis
We describe in detail, how we determine the light-cone velocity v lc , the butterfly velocity v b , and the Lyapunov exponent λ L . The light-cone velocity is defined as the ratio of the distance |i − j| and the time at which the reduced OTO correlator F r ij (t) reaches a small threshold. The butterfly velocity, however, sets a scale for the time it takes to scramble information over the system and is therefore defined via the time at which F r ij (t) attains a large value of order one. The specific threshold one chooses to determine the butterfly velocity is thus somewhat arbitrary. We illustrate the dependence of the velocity v on the chosen threshold F * of the reduced OTO correlator F r ij (t) in Fig. 7. For large values of F * , the velocity converges toward a constant. Hence, the butterfly velocity will be largely insensitive to the precise choice of F * as long as it is large enough. For the definition of v b , we consider the specific value of F * b = 0.2. In the limit F * → 0, there is a strong dependence of v on the choice of the threshold. The light-cone velocity v lc is defined by the fastest spread of information through the system and is determined by the reduced OTO correlator attaining a small value. To fulfill this definition, we fix F * lc = 0.0005; see inset in Fig. 7.
As described in Sec. I A, the OTO correlator is expected to grow exponentially on a timescale set approximately by the butterfly velocity. We thus fit the exponential function to the numerical data simultaneously for distances 1 ≤ |i − j| ≤ 5 within the range −2.5 ≤ log F r ij (t) ≤ −1. The butterfly velocity v b is determined as described above with the threshold F * = 0.2, which lies well within the interval of the considered data points; see Fig. 8 (a) for an exemplary plot.
We note that the exponential growth of the reduced OTO correlator is in our model limited to a rather small dynamical range. Extracting by contrast the growth rate Λ lc of the OTO correlator from the timescale set by the light-cone velocity within the range −14 ≤ log F r ij (t) ≤ −4.5, Fig. 8 (b), yields larger rates. In particular, for the parameters shown in Fig. 8, we obtain λ L = 2.9(1) and Λ lc = 10.3(5), respectively.  (5). The dashed gray line denotes the threshold value F * used to determine the velocities vb and vlc, respectively. We obtain the exponents by fitting our data in a restricted regime around the threshold value F * to the predicted exponential growth, see text for details. However, we note that in our data the exponential growth is limited to a rather small time range. The errorbars shown in Fig. 4 correspond to errors obtained from such fits.

Supplemental Information Appendix A: Details on the Experimental Protocols
We elaborate on the different protocols outlined in Sec. I C, which can be used to experimentally access the theoretical findings presented in this work.

Global Many-Body Interferometry
We consider a system of bosons or fermions in an optical lattice. At this point, we do not make any assumptions about the specific form of the HamiltonianĤ. We first focus on the real-time and spatially resolved single-particle Green's functions G gl ij (t), which can be measured by the following protocol, Fig. 9 (a): (1) Initially, prepare two identical copies of a pure state |ψ ⊗ |ψ . Remove a particle on site i in the left system by locally transferring the atom to a hyperfine state that is decoupled from the rest of the system or by transferring it to a higher band of the optical lattice, yielding c i |ψ ⊗ |ψ . (2) The system evolves in time for a period t, We abbreviate the sequence of the operations (1-3) asÔ(i, j) and illustrate the corresponding quantum circuit in the bottom of the blue box in Fig. 9 (a). (4) Finally, measure the swap operatorV, which interchanges the particles between the left and the right subsystem The expectation value of the swap operator is experimentally determined by a global 50%-50% beam splitter operation, which is realized by tunnel-coupling the left and the right system, followed by a measurement of the parity-projected particle number [19,40,41]. OTO correlation functions are measured in a similar fashion, Fig. 9b. To begin with, we recycle the first three steps of the Green's function protocol, compiled inÔ(i, j). As a second step, the sign of the Hamiltonian needs to be inverted globally. The sign of the interaction can be flipped by ramping the magnetic field across a Feshbach resonance, as demonstrated experimentally, for instance, in the realization of negative temperature states [50]. Furthermore, by appropriately tuning the drive frequency of a modulated optical lattice, the sign of the hopping matrix element can be flipped [51]. Combining these already established experimental techniques, the global sign of the Hamiltonian is inverted. As a next step, the operationsÔ(j, i) are applied again, leading to the time evolved state |ψ l (t) ⊗ |ψ r (t) ≡ e iĤt c j e −iĤt c i |ψ ⊗ c i e iĤt c j e −iĤt |ψ .
(A3) The square modulus of the OTO correlators is then obtained by measuring the wavefunction overlap of the left and the right system using beam splitters as discussed before.
For the measurement of both the Green's function and the OTO correlators, the initial state |ψ can be an arbitrary pure state, such as the ground state, or a simple product state. An effective finite temperature state can be obtained for quenches from initial pure states to some final Hamiltonian. In a thermalizing system [30][31][32], the effective temperature is then determined by the energy-density produced by the quantum quench. In the case of a thermal initial state, after the first three steps of our protocol, blue box in Fig. 9, the system is prepared in the state ρ l (t) ⊗ ρ r (t), where ρ α (t) is a generic density matrix. The measurement of the swap operatorV yields [41] For pure states, ρ l,r (t) = |ψ l,r (t) ψ l,r (t)|, we directly obtain Eq. (A2). However, at finite temperature, the measurement does not directly yield the square of the correlation function. In particular, we obtain for the Green's function protocol By contrast, the desired modulus square of the thermal Green's function would be (A6) Hence, at high temperatures, Eq. (A5) is suppressed by a factor 1/Z, where Z is the partition sum, and thus vanishes in the thermodynamic limit. A similar reasoning applies in the case of OTO correlators.

Local Many-Body Interferometry
Interfering two system copies globally requires beam splitter operations with high fidelity, as in each measurement for systems of size L the same number of beam splitter operations have to be applied. To overcome this challenge, we introduce an alternative protocol that is scalable since it only requires two beam splitter operations irrespective of the system size.
A local beam splitter operation on site l is realized by coupling the left and the right copy of the quantum system by a A detailed description of the protocl is givne in the text. HamiltonianĤ → −Ĥ, as suggested in the previous section. Let the system evolve in time for the duration t. (5) Apply the 50%-50% beam splitter operation BS i on site i. Evaluating these steps, we find This expression corresponds to the product of two OTO correlation function. Four point correlators in spin systems can be obtained with related protocols [18].
The OTO correlator F loc ij (t) obtained from local interference contains at the considered temperatures essentially the same information as the one we originally introduced. As the protocol measures the imaginary part of a product of two OTO correlators, it starts out at zero. The scrambling across the quantum state manifests itself in the linear propagation of a wave-packet in F loc ij (t) [see Fig. 11 (a)] from which light-cone and butterfly velocities can be extracted. In Fig. 11 (a), we once again attribute the superballistic spread of information, which kicks in at |i − j| 7, to the finite MPO bond dimension of 400. Similarly, G loc ij (t) starts off at zero but then develops a peak that quickly decays; Fig. 11 (b). From that we determine the quasiparticle lifetime τ J ∼ 0.32 which corresponds roughly to half the lifetime obtained for the Green's function G ij (t). This factor can be attributed to the fact that here the product of two correlation functions is measured.
With the protocols discussed so far, static one-body correlation functions can be measured by setting the physical time t = 0. Moreover, a generalization of the local protocol makes it possible to measure static correlations functions of arbitrary order. Specifically, correlators of δn i determine one-body correlation functions of the original many-body state: Here we used that the left and the right initial states are identical. Higher order static correlation functions in the creation a † i and annihilation operators a i are straightforwardly obtained by measuring higher order correlators in δn i . We emphasize that this protocol scales favorable with system size, and that correlators between arbitrary sites and of arbitrary order can be taken in a single shot by performing the beam splitter operations on the full system.
The density matrix describing the quantum state of a system can be expressed aŝ where N is a normalization constant and theσ ij constitute a suitable basis [52]. In the case of fermions or hard-core bosons, one possible choice for the basis are the Pauli matrices. The knowledge of correlators up to sufficient order makes it possible to determine the so-called Stokes parameters r i1,...i N and thereby to reconstruct the density matrix, which paves the way for the full state tomography of quantum states with massive particles.

Measuring Dynamical Density Correlators
In this section we discuss two different possibilities to measure dynamical density correlation functions and thereby observe their diffusive behavior.
The dynamic structure factor S(k, ω), which is the spatial and temporal Fourier transform of the density correlator C n (x, t), can be measured with Bragg spectroscopy [53,54]. In Bragg spectroscopy, the detuning of the two laser beams sets the frequency ω and the angle between the beams the transferred momentum k. A measurement of the absorption of the system as a function of k and ω directly maps out the dynamic structure factor S(k, ω). Diffusion manifests itself in the wavevector and frequency resolved structure factor S(k, ω) as Lorentzian peaks with half-width-half-maximum that scales as Dq 2 .
It is furthermore possible to measure the dissipative response [n i (t), n j ] to a local perturbation of the system in a quantum gas microscope. To this end, a local potential δH = n j δµ is created at site j by applying a laser for a short time τ , yielding the time evolution exp[−iδHτ ] ∼ 1 − iδHτ + O(δµ 2 τ 2 ). Measuring the density at site i after the unitary time evolution for duration t we obtain χ loc ij (t) = n i (t) + iδµτ [n i (t), n j ] + O(δµ 2 τ 2 ). (A17) In equilibrium, the fluctuation-dissipation theorem provides an exact relation between [n i (t), n j ] and n i (t)n j . The accurate measurement of the former therefore enables the observation of diffusive response in the dynamical density correlator.