Scrambling and thermalization in a diffusive quantum many-body system

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 $\hat{W}$ and $\hat{V}$ in time

$$\begin{eqnarray}&&C(t)=-\langle | [\hat{W}(t),\hat{V}(0)]{| }^{2}\rangle .\end{eqnarray} \tag{ 2 }$$

In a semiclassical picture, the commutator in equation (2) can be replaced by Poisson brackets. Then, for the choice of $\hat{W}={p}_{j}$ and $\hat{V}={p}_{i}$, this quantity reduces to $C{(t)\sim \langle (\partial {p}_{j}(t)/\partial {q}_{i}(0))}^{2}\rangle$. Therefore, the correlation function C(t) describes the sensitivity of the time evolution and is expected to grow exponentially at short times $\sim \,\exp [{\lambda }_{{\rm{L}}}t]$, with a rate ${\lambda }_{{\rm{L}}}$ that resembles the Lyapunov exponent in classical chaotic dynamics. Rewriting these momenta and coordinates as combinations of creation and annihilation operators, equation (2) generically consists of OTO correlators of the form

$$\begin{eqnarray}&&{{ \mathcal F }}_{{ij}}(t)=\langle {c}_{j}^{\dagger }(t){c}_{i}^{\dagger }{c}_{j}(t){c}_{i}\rangle .\end{eqnarray} \tag{ 3 }$$

Below we mainly consider the quantum statistical average $\langle \ldots \rangle ={\rm{tr}}[\hat{\rho }\ldots ]$ over an initial thermal state with weights distributed according to the Gibbs ensemble $\hat{\rho }={{\rm{e}}}^{-\hat{H}/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 $\hat{\rho }=| {\psi }_{0}\rangle \langle {\psi }_{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–32].

OTO correlators ${{ \mathcal F }}_{{ij}}$ evaluated at comparatively high temperatures $T=4J$, interactions U = J, and chemical potential $\mu =0$ are shown in figure 1(a) as a function of time t and distance $(i-j);$ for OTO correlators in the coherent Luttinger Liquid regime see [33]. Despite the high temperature, the OTO correlator ${{ \mathcal F }}_{{ij}}$ unveils a pronounced light-cone spreading of the information across the quantum state for $| i-j| \lesssim 7$. For larger distances the light cone seems to exhibit a plateau, which we, however, attribute to the finite MPO bond dimension considered in the numerical simulations; see methods 3.1. OTO correlators are in that respect challenging to simulate with MPO techniques, because they directly reflect the fast spreading of entanglement.

The OTO correlator ${{ \mathcal F }}_{{ij}}(t)$ should be contrasted to the time-ordered single-particle Green's function

$$\begin{eqnarray}&&{{ \mathcal G }}_{{ij}}(t)=\langle {c}_{j}^{\dagger }(t){c}_{i}\rangle ,\end{eqnarray} \tag{ 4 }$$

which is shown in figure 1(b). In the incoherent transport regime, where well-defined quasiparticles do not exist, the Green's function ${{ \mathcal 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 [34]. For the chosen parameters ($U=J,\mu =0,T=4J$), we find that the quasiparticle lifetime is approximately $\tau J\sim 0.6$ and hence shorter than the microscopic hopping rate, which indicates incoherent transport. By contrast, the OTO structure of ${{ \mathcal F }}_{{ij}}(t)$ reveals a well defined linear spread of quantum information despite the high temperature.

We now characterize the OTO correlators ${{ \mathcal F }}_{{ij}}(t)$ in detail. To this end, we subtract $\langle {\hat{n}}_{i}{\hat{n}}_{j}\rangle$ from the ${{ \mathcal F }}_{{ij}}(t)$ and consider its relative change: ${{ \mathcal F }}_{{ij}}^{r}(t)=| {{ \mathcal F }}_{{ij}}(t)-\langle {\hat{n}}_{i}{\hat{n}}_{j}\rangle | /\langle {\hat{n}}_{i}{\hat{n}}_{j}\rangle$. Examples for the reduced OTO correlator ${{ \mathcal F }}_{{ij}}^{r}(t)$ are shown in figure 2 for interaction U = J and different temperatures T. The reduced OTO correlator ${{ \mathcal F }}_{{ij}}^{r}(t)$ starts off at zero, forms the light-cone plateau, and approaches the steady-state value as an exponential.

Zoom In Zoom Out Reset image size

From the light-cone spread of the OTO correlator, we extract two velocities (figure 3(a)): (1) the light-cone velocity ${v}_{\mathrm{lc}}$, which we define by the space–time region where the reduced OTO correlators ${{ \mathcal F }}_{{ij}}^{r}(t)$ surpasses a small threshold of 0.05% of its final value. (2) The butterfly velocity ${v}_{{\rm{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}_{{\rm{b}}}$ does not significantly depend on this cutoff, as long as it is chosen to a sizeable fraction; see methods figure 7. The light cone velocity ${v}_{\mathrm{lc}}$ increases with temperature T and is bounded from below by the zero temperature Luttinger liquid velocity; see figure 3(b). The butterfly velocity ${v}_{{\rm{b}}}$ is systematically lower than ${v}_{\mathrm{lc}}$ and is almost independent of temperature. The butterfly velocity determines the time scale ${t}_{\mathrm{scr}}$ for scrambling information across the many-body quantum state which is linear in system size ${t}_{\mathrm{scr}}\sim L/{v}_{{\rm{b}}}$. Based on results from holography, it has been argued in [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 3.2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Previous studies of strongly coupled field theories [9] or disordered SYK models [7] predict an exponential growth in the initial dynamics of the reduced OTO correlator of the form ${{ \mathcal F }}_{x}^{r}(t)\sim \exp [{\lambda }_{{\rm{L}}}(t-x/{v}_{{\rm{b}}})]$ with butterfly velocity ${v}_{{\rm{b}}}$ and a growth rate ${\lambda }_{{\rm{L}}}$, which is refered to as Lyapunov exponent. Both systems display a clear separation between the collision time τ and the scrambling time ${t}_{\mathrm{scr}}$, at which the OTO correlator assumes an appreciably large value. Hence, in these systems the exponential growth occurs up to parametrically late times determined by a large parameter N, which controls the approach to a semiclassical limit. We tested the results of our numerical simulations against this prediction and find that they are incompatible with an exponential growth in time inasmuch as there appears to be no parametrically large regime of exponential growth in our model, see figure 3. This is a consequence of the average particle density, which is controlled by the chemical potential, being typically on the order of one for the chosen parameters. This effectively restricts the accessed local Hilbert space dimension and hence there is no large separation between collision and scrambling time [35]. Therefore, we cannot unambiguously define a Lyapunov exponent. Finding the analytic form for the initial growth in our system remains an outstanding challenge.

Nonetheless, one can estimate a growth rate ${\lambda }_{{\rm{L}}}$ by linearizing the OTO correlator around the space–time region set by the butterfly velocity, which is shown as a function of temperature for different values of the interaction strength U and chemical potential $\mu =0$ in figure 4. We note however that the concrete values for the Lyapunov exponent depend on the space–time region around which we linearize, see figure 8. It has been conjectured that the Lyapunov exponent is bounded by $2\pi T$, which is the value it assumes in a strongly coupled field theory with a gravity dual [9]. In our system, ${\lambda }_{{\rm{L}}}$ is parametrically lower than this bound and increases slowly when lowering the temperature. Moreover, we find that the dependence of the Lyapunov exponent on the interaction strength U is small with slightly larger values of ${\lambda }_{{\rm{L}}}$ for intermediate interaction strength, $U=3J$, which is in the vicinity of the quantum critical point.

Zoom In Zoom Out Reset image size

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 [36, 37]

$$\begin{eqnarray}&&{\partial }_{t}n-D{{\rm{\nabla }}}^{2}n=0,\end{eqnarray} \tag{ 5 }$$

where D is the diffusion constant. The connected density correlation function ${C}_{n}(x-x^{\prime} ,t-t^{\prime} )={\rm{Re}}[\langle n(x,t)n(x^{\prime} ,t^{\prime} )\rangle -\langle n(x)\rangle \langle n(x^{\prime} )\rangle ]$ relates the density at space–time (x, t) to the density at $(x^{\prime} ,t^{\prime} )$ via $n(x,t)\sim \int {\rm{d}}{t}^{\prime} {\rm{d}}{x}^{\prime} {C}_{n}(x-x^{\prime} ,t-t^{\prime} )n(x^{\prime} ,t^{\prime} )$ and in a hydrodynamic regime is expected to be of the form

$$\begin{eqnarray}&&{C}_{n}(x,t)\cong \displaystyle \frac{\tilde{C}}{\sqrt{4\pi D| t| }}{{\rm{e}}}^{-\tfrac{{x}^{2}}{4D| t| }},\end{eqnarray} \tag{ 6 }$$

with $\tilde{C}=\int {\rm{d}}{x}\,{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), figure 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^{\prime} =0$) density correlation function. At high temperatures $T\gtrsim 4J$ the local correlator exhibits a diffusive power-law decay ${C}_{n}(x=0,t)\sim 1/\sqrt{{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 table 1. 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/\sqrt{{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 [37].

Table 1.  Diffusion constant D and the ratio $D{\lambda }_{{\rm{L}}}/{v}_{{\rm{B}}}^{2}$ for different temperatures T. The errors as indicated in the parentheses are errors from the fits.

T/J	$D/({a}^{2}J)$	$D{\lambda }_{{\rm{L}}}/{v}_{{\rm{b}}}^{2}$
4	14.29(27)	7.2(6)
6	11.69(10)	6.0(4)
8	10.42(04)	5.4(3)
10	9.79(01)	5.1(3)

It has been proposed that the diffusion constant is related to the butterfly velocity ${v}_{{\rm{b}}}$ and the Lyapunov exponent ${\lambda }_{{\rm{L}}}$ via $D\sim {v}_{{\rm{b}}}^{2}/{\lambda }_{{\rm{L}}}$ [38–41], where $1/{\lambda }_{{\rm{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{\lambda }_{{\rm{L}}}/{v}_{{\rm{b}}}^{2}\sim 5.5$ for temperatures $T\gtrsim 6J;$ see table 1, 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}_{\mathrm{th}}\sim {L}^{2}/D$. Experimentally measuring OTO correlators (section 1.3) and density correlators (section 3.3.3) will make it possible to further check these holographic predictions.

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 [42–44] 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. In the following, we sketch the main features of the protocols and refer to the methods section 3.3 for technical details on the implementation.

Global interferometry precisely yields the square modulus of the single-particle Green's functions ${{ \mathcal G }}_{{ij}}^{\mathrm{gl}}(t)=| {{ \mathcal G }}_{{ij}}(t){| }^{2}$ and the OTO correlators ${{ \mathcal F }}_{{ij}}^{\mathrm{gl}}(t)=| {{ \mathcal F }}_{{ij}}(t){| }^{2}$ for pure initial states (see section 3.3.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 [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 section 3.3.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 ${{ \mathcal G }}_{{ij}}^{\mathrm{loc}}(t)\sim {\rm{Im}}[\langle {a}_{j}^{\dagger }(t){a}_{i}\rangle \langle {a}_{j}(t){a}_{i}^{\dagger }\rangle ]$ and OTO correlator ${{ \mathcal F }}_{{ij}}^{\mathrm{loc}}(t)\sim {\rm{Im}}[\langle {a}_{j}^{\dagger }(t){a}_{i}^{\dagger }{a}_{j}(t){a}_{i}\rangle \langle {a}_{j}^{\dagger }(t){a}_{i}{a}_{j}(t){a}_{i}^{\dagger }\rangle ]$. However, we demonstrate in section 3.3.2 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 straightforward.

Our numerical simulations are based on finite-temperature, time-dependent MPO [45–50].

For the density correlations, we evaluate [49]

$$\begin{eqnarray}&&{\langle {\hat{n}}_{{\ell }}(t){\hat{n}}_{j}(0)\rangle }_{\beta }\mathop{=}\limits^{\mathrm{def}}\displaystyle \frac{1}{Z}\,{\rm{tr}}[{{\rm{e}}}^{-\beta \hat{H}}({{\rm{e}}}^{{\rm{i}}{t}\hat{H}}{\hat{n}}_{{\ell }}{{\rm{e}}}^{-{\rm{i}}{t}\hat{H}}){\hat{n}}_{j}]=\displaystyle \frac{1}{Z}\,{\rm{tr}}\left[\left({{\rm{e}}}^{{\rm{i}}\tfrac{t}{2}\hat{H}}{{\rm{e}}}^{-\beta \hat{H}/2}{\hat{n}}_{{\ell }}{{\rm{e}}}^{-{\rm{i}}\tfrac{t}{2}\hat{H}}\right)\left({{\rm{e}}}^{-{\rm{i}}\tfrac{t}{2}\hat{H}}{\hat{n}}_{j}{{\rm{e}}}^{-\beta \hat{H}/2}{{\rm{e}}}^{{\rm{i}}\tfrac{t}{2}\hat{H}}\right)\right],\end{eqnarray} \tag{ 7 }$$

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 ${{\rm{e}}}^{-\beta \hat{H}/2}{\hat{n}}_{{\ell }}$ and ${\hat{n}}_{j}{{\rm{e}}}^{-\beta \hat{H}/2}$, respectively, and then performing a real-time evolution up to $\tfrac{t}{2}$ and $-\tfrac{t}{2}$. By exploiting the time translation invariance, $\langle {\hat{n}}_{{\ell }}(t){\hat{n}}_{j}(0){\rangle }_{\beta }=\langle {\hat{n}}_{{\ell }}(t/2){\hat{n}}_{j}(-t/2){\rangle }_{\beta }$, 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 ${{\rm{e}}}^{-\beta \hat{H}/2}$, we employ a second-order Suzuki–Trotter decomposition with imaginary time step ${\rm{\Delta }}\tau$ (typically ${\rm{\Delta }}\tau J=0.025$) after splitting the Hamiltonian into even and odd bonds, as described in [45]. The real-time evolution proceeds by Liouville steps $\hat{A}(t+{\rm{\Delta }}t)={{\rm{e}}}^{{\rm{i}}{\rm{\Delta }}t\hat{H}}\hat{A}(t){{\rm{e}}}^{-{\rm{i}}{\rm{\Delta }}t\hat{H}}$. For each of the steps we combine a fourth-order partitioned Runge–Kutta method [51] with even–odd bond splitting of the Hamiltonian. As noted in [49], 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 ${\rm{\Delta }}{tJ}=0.125$ or 0.25).

For the OTO correlators $\langle {c}_{j}^{\dagger }(t){c}_{{\ell }}^{\dagger }{c}_{j}(t){c}_{{\ell }}{\rangle }_{\beta }$, a regrouping analogous to equation (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

$$\begin{eqnarray}&&{\langle {c}_{j}^{\dagger }(t){c}_{{\ell }}^{\dagger }{c}_{j}(t){c}_{{\ell }}\rangle }_{\beta }=\displaystyle \frac{1}{Z}\,{\rm{tr}}[({{\rm{e}}}^{{\rm{i}}{t}\hat{H}}{{\rm{e}}}^{-\beta \hat{H}}{c}_{j}^{\dagger }{{\rm{e}}}^{-{\rm{i}}{t}\hat{H}}{c}_{{\ell }}^{\dagger })({{\rm{e}}}^{{\rm{i}}{t}\hat{H}}{c}_{j}{{\rm{e}}}^{-{\rm{i}}{t}\hat{H}}{c}_{{\ell }})]\end{eqnarray} \tag{ 8 }$$

and time-evolve both ${{\rm{e}}}^{-\beta \hat{H}}{c}_{j}^{\dagger }$ and c_j up to time t. Subsequent application of the site-local operators ${c}_{{\ell }}^{\dagger }$ 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 mainly 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 figure 6 we compare the data obtained for the same simulation parameters but different maximal bond dimensions. The MPO bond dimension of 20 apparently leads to a spread of quantum information from distances $| i-j| \approx 6$ to $| i-j| \approx 10$ within a duration $\delta t\sim 0$, which manifests itself in a plateau in the light-cone, see figure 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.

We describe in detail how we determine the light-cone velocity ${v}_{\mathrm{lc}}$, the butterfly velocity ${v}_{{\rm{b}}}$, and the Lyapunov exponent. The light-cone velocity is defined as the ratio of the distance $| i-j|$ and the time at which the reduced OTO correlator ${{ \mathcal F }}_{{ij}}^{r}(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 ${{ \mathcal F }}_{{ij}}^{r}(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 ${{ \mathcal F }}^{* }$ of the reduced OTO correlator ${{ \mathcal F }}_{{ij}}^{r}(t)$ in figure 7. For large values of ${{ \mathcal F }}^{* }$, the velocity converges toward a constant. Hence, the butterfly velocity will be largely insensitive to the precise choice of ${{ \mathcal F }}^{* }$ as long as it is large enough. For the definition of ${v}_{{\rm{b}}}$, we consider the specific value of ${{ \mathcal F }}_{{\rm{b}}}^{* }=0.2$.

In the limit ${{ \mathcal F }}^{* }\to 0$, there is a strong dependence of v on the choice of the threshold. The light-cone velocity ${v}_{\mathrm{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 ${{ \mathcal F }}_{\mathrm{lc}}^{* }=0.0005;$ see inset in figure 7.

As described in section 1.1, from the SYK model and strongly coupled field theories one could expect that the OTO correlator grows exponentially on a timescale set approximately by the butterfly velocity. We thus fit the exponential function

$$\begin{eqnarray}&&{{ \mathcal F }}_{x}^{r}(t)=a\cdot {{\rm{e}}}^{{\lambda }_{{\rm{L}}}\left(t-\tfrac{x}{{v}_{{\rm{b}}}}\right)}\end{eqnarray} \tag{ 9 }$$

to the numerical data simultaneously for distances $1\leqslant | i-j| \leqslant 5$ within the range $-2.5\leqslant \mathrm{log}{{ \mathcal F }}_{{ij}}^{r}(t)\leqslant -1$. The butterfly velocity ${v}_{{\rm{b}}}$ is determined as described above with the threshold ${{ \mathcal F }}^{* }=0.2$, which lies well within the interval of the considered data points; see figure 8(a) for an exemplary plot. We note, however, that our data does not support an exponential growth of the OTO correlator over parametrically large times, since the collision time τ is similar to the scrambiling time ${t}_{\mathrm{scr}}$ in our system. This is demonstrated by extracting the growth rate ${{\rm{\Lambda }}}_{\mathrm{lc}}$ from linearizing the OTO correlator around the light-cone velocity within the range $-14\leqslant \mathrm{log}{{ \mathcal F }}_{{ij}}^{r}(t)\leqslant -4.5$, figure 8(b), which yields larger rates. In particular, for the parameters shown in figure 8, we obtain ${\lambda }_{{\rm{L}}}=2.9(1)$ and ${{\rm{\Lambda }}}_{\mathrm{lc}}=10.3(5)$, respectively.

We elaborate on the different protocols outlined in section 1.3, which can be used to experimentally access the theoretical findings presented in this work.

3.3.1. 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 $\hat{H}$. We first focus on the real-time and spatially resolved single-particle Green's functions ${{ \mathcal G }}_{{ij}}^{\mathrm{gl}}(t)$, which can be measured by the following protocol, figure 9(a): (1) initially, prepare two identical copies of a pure state $| \psi \rangle \otimes | \psi \rangle$. 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}| \psi \rangle \otimes | \psi \rangle$. (2) The system evolves in time for a period $t,\exp [-{\rm{i}}\hat{H}t]{c}_{i}| \psi \rangle \otimes \exp [-{\rm{i}}\hat{H}t]| \psi \rangle$. (3) Create a hole on site j of the right system

$$\begin{eqnarray}&&| {\psi }_{{\rm{l}}}(t)\rangle \otimes | {\psi }_{{\rm{r}}}(t)\rangle \equiv {{\rm{e}}}^{-{\rm{i}}\hat{H}t}{c}_{i}| \psi \rangle \otimes {c}_{j}{{\rm{e}}}^{-{\rm{i}}\hat{H}t}| \psi \rangle .\end{eqnarray} \tag{ 10 }$$

We abbreviate the sequence of the operations (1)–(3) as $\hat{O}(i,j)$ and illustrate the corresponding quantum circuit in the bottom of the blue box in figure 9(a). (4) Finally, measure the swap operator $\hat{{ \mathcal V }}$, which interchanges the particles between the left and the right subsystem

$$\begin{eqnarray}&&\langle \hat{{ \mathcal V }}\rangle ={\rm{tr}}[| {\psi }_{{\rm{r}}}(t)\rangle \langle {\psi }_{{\rm{l}}}(t)| \otimes | {\psi }_{{\rm{l}}}(t)\rangle \langle {\psi }_{{\rm{r}}}(t)| ]=| \langle {\psi }_{{\rm{r}}}(t)| {\psi }_{{\rm{l}}}(t)\rangle {| }^{2}=| \langle {c}_{j}^{\dagger }(t){c}_{i}\rangle {| }^{2}={{ \mathcal G }}_{{ij}}^{\mathrm{gl}}(t).\end{eqnarray} \tag{ 11 }$$

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, 42, 43].

Zoom In Zoom Out Reset image size

OTO correlation functions are measured in a similar fashion, figure 9(b). To begin with, we recycle the first three steps of the Green's function protocol, compiled in $\hat{O}(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 [52]. Furthermore, by appropriately tuning the drive frequency of a modulated optical lattice, the sign of the hopping matrix element can be flipped [53]. Combining these already established experimental techniques, the global sign of the Hamiltonian is inverted. As a next step, the operations $\hat{O}(j,i)$ are applied again, leading to the time evolved state

$$\begin{eqnarray}&&| {\psi }_{{\rm{l}}}(t)\rangle \otimes | {\psi }_{{\rm{r}}}(t)\rangle \equiv {{\rm{e}}}^{{\rm{i}}\hat{H}t}{c}_{j}{{\rm{e}}}^{-{\rm{i}}\hat{H}t}{c}_{i}| \psi \rangle \otimes {c}_{i}{{\rm{e}}}^{{\rm{i}}\hat{H}t}{c}_{j}{{\rm{e}}}^{-{\rm{i}}\hat{H}t}| \psi \rangle .\end{eqnarray} \tag{ 12 }$$

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 $| \psi \rangle$ 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–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 figure 9, the system is prepared in the state ${\rho }_{{\rm{l}}}(t)\otimes {\rho }_{{\rm{r}}}(t)$, where ${\rho }_{\alpha }(t)$ is a generic density matrix. The measurement of the swap operator $\hat{{ \mathcal V }}$ yields [43]

$$\begin{eqnarray}\langle \hat{{ \mathcal V }}\rangle & = & {\rm{tr}}\hat{{ \mathcal V }}{\rho }_{{\rm{l}}}(t)\otimes {\rho }_{{\rm{r}}}(t)={\rm{tr}}\hat{{ \mathcal V }}\displaystyle \sum _{\mu \nu }{\rho }_{{\rm{l}},\mu }(t){\rho }_{{\rm{r}},\nu }(t)| \mu \rangle \langle \mu | \otimes | \nu \rangle \langle \nu | ={\rm{tr}}\displaystyle \sum _{\mu \nu }{\rho }_{{\rm{l}},\mu }(t){\rho }_{{\rm{r}},\nu }(t)| \nu \rangle \langle \mu | \otimes | \mu \rangle \langle \nu | \\ & = & \displaystyle \sum _{\mu }{\rho }_{{\rm{l}},\mu }(t){\rho }_{{\rm{r}},\mu }(t)={\rm{tr}}{\rho }_{{\rm{l}}}(t){\rho }_{{\rm{r}}}(t).\end{eqnarray} \tag{ 13 }$$

For pure states, ${\rho }_{{\rm{l}},{\rm{r}}}(t)=| {\psi }_{{\rm{l}},{\rm{r}}}(t)\rangle \langle {\psi }_{{\rm{l}},{\rm{r}}}(t)|$, we directly obtain equation (11). 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

$$\begin{eqnarray}&&\langle \hat{{ \mathcal V }}\rangle =\displaystyle \sum _{\mu \nu }\,{\rho }_{\mu }{\rho }_{\nu }\langle \mu | {c}_{i}^{\dagger }{c}_{j}(t)| \nu \rangle \langle \nu | {c}_{j}^{\dagger }(t){c}_{i}| \mu \rangle .\end{eqnarray} \tag{ 14 }$$

By contrast, the desired modulus square of the thermal Green's function would be

$$\begin{eqnarray}&&| \langle {c}_{i}{(t)}^{\dagger }{c}_{j}\rangle {| }^{2}=\displaystyle \sum _{\mu \nu }{\rho }_{\mu }{\rho }_{\nu }\langle \mu | {c}_{i}^{\dagger }{c}_{j}(t)| \mu \rangle \langle \nu | {c}_{j}^{\dagger }(t){c}_{i}| \nu \rangle .\end{eqnarray} \tag{ 15 }$$

Hence, at high temperatures, equation (14) 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.

3.3.2. 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 tunneling Hamiltonian

$$\begin{eqnarray}&&{\hat{H}}_{{\rm{l}}}^{{\rm{BS}}}=-{J}_{{\rm{BS}}}({a}_{{\rm{l}}}^{\dagger }{b}_{{\rm{l}}}+{b}_{{\rm{l}}}^{\dagger }{a}_{{\rm{l}}}),\end{eqnarray} \tag{ 16 }$$

where ${a}_{{\rm{l}}}^{\dagger }$ (${b}_{{\rm{l}}}^{\dagger }$) creates a particle in the left (right) system. The unitary evolution under equation (16), ${{\rm{BS}}}_{{\rm{l}}}(\tau )=\exp [-{\rm{i}}{\hat{H}}_{{\rm{l}}}^{{\rm{BS}}}\tau ]$, for time ${\tau }_{{\rm{BS}}}=\pi /4{J}_{{\rm{BS}}}$ defines a 50%-50% beam splitter operation

$$\begin{eqnarray}&&\left(\begin{array}{c}{\bar{a}}_{{\rm{l}}}\\ {\bar{b}}_{{\rm{l}}}\end{array}\right)=\mathop{\underbrace{\displaystyle \frac{1}{\sqrt{2}}\left(\begin{array}{cc}1 & -{\rm{i}}\\ -{\rm{i}} & 1\end{array}\right)}}\limits_{={\mathrm{BS}}_{{\rm{l}}}}\left(\begin{array}{c}{a}_{{\rm{l}}}\\ {b}_{{\rm{l}}}\end{array}\right).\end{eqnarray} \tag{ 17 }$$

Furthermore, the phase of the beam splitter can be adjusted by applying a field gradient between the left and the right system ${\hat{H}}_{{\rm{l}}}^{{\rm{F}}}=\tfrac{h}{2}({b}_{{\rm{l}}}^{\dagger }{b}_{{\rm{l}}}-{a}_{{\rm{l}}}^{\dagger }{a}_{{\rm{l}}})$ for a duration ${\tau }_{f},R({\tau }_{f})=\exp [-{\rm{i}}{\hat{H}}_{{\rm{l}}}^{{\rm{F}}}{\tau }_{f}]$:

$$\begin{eqnarray}{\overline{\mathrm{BS}}}_{{\rm{l}}}={R}^{\dagger }(\phi ){\mathrm{BS}}_{{\rm{l}}}\,R(\phi )=\displaystyle \frac{1}{\sqrt{2}}\left(\begin{array}{cc}1 & -{{\rm{ie}}}^{-{\rm{i}}\phi }\\ -{{\rm{ie}}}^{{\rm{i}}\phi } & 1\end{array}\right),\end{eqnarray} \tag{ 18 }$$

where $\phi =h{\tau }_{f}$.

The time ordered Green's function for a system prepared in an arbitrary density matrix can be measured by the following sequence (figure 10): (1) apply a beam splitter operation on site i for a short duration $\tau {J}_{\mathrm{BS}}\ll 1$. In that limit, the unitary evolution can be linearized ${{\rm{BS}}}_{{\rm{l}}}(\tau )=1-{\rm{i}}{\hat{H}}_{{\rm{l}}}^{{\rm{BS}}}\tau +O({J}_{{\rm{BS}}}^{2}{\tau }^{2})$. (2) Let the two copies evolve for the physical time t. (3) Apply a 50%-50% beam splitter operation on site j with a phase that is detuned from the first one by $\phi =\pi /2$. (4) Finally, the density difference $\delta {\bar{n}}_{j}$ between the right and the left subsystem is measured. This leads to the following measurement outcome

$$\begin{eqnarray}&&{{ \mathcal G }}_{{ij}}^{\mathrm{loc}}(t)=\langle {{\rm{BS}}}_{i}^{\dagger }(\tau ){{\rm{e}}}^{{\rm{i}}\hat{H}t}{\overline{{\rm{BS}}}}_{j}^{\dagger }\delta {\bar{n}}_{j}{\overline{{\rm{BS}}}}_{j}{{\rm{e}}}^{-{\rm{i}}\hat{H}t}{{\rm{BS}}}_{i}(\tau )\rangle .\end{eqnarray} \tag{ 19 }$$

We first calculate the densities after the beam splitter operation ${\overline{\mathrm{BS}}}_{j}$, which gives

$$\begin{eqnarray}&&{\bar{n}}_{{\rm{l}},j}={\bar{a}}_{j}^{\dagger }{\bar{a}}_{j}=\displaystyle \frac{1}{2}({a}_{j}^{\dagger }-{b}_{j}^{\dagger })({a}_{j}-{b}_{j}),\end{eqnarray} \tag{ 20a }$$

$$\begin{eqnarray}&&{\bar{n}}_{{\rm{r}},j}={\bar{b}}_{j}^{\dagger }{\bar{b}}_{j}=\displaystyle \frac{1}{2}({a}_{j}^{\dagger }+{b}_{j}^{\dagger })({a}_{j}+{b}_{j}).\end{eqnarray} \tag{ 20b }$$

Computing the density difference between the right and the left system, we find

$$\begin{eqnarray}&&\delta {\bar{n}}_{j}={\bar{n}}_{{\rm{l}},j}-{\bar{n}}_{{\rm{r}},j}=-({b}_{j}^{\dagger }{a}_{j}+{a}_{j}^{\dagger }{b}_{j}).\end{eqnarray} \tag{ 21 }$$

Considering now that the duration of the first beam splitter operation on site i is short and using the particle number conservation, we obtain

$$\begin{eqnarray}&&{{ \mathcal G }}_{{ij}}^{\mathrm{loc}}(t)=4{J}_{{\rm{BS}}}\tau {\rm{Im}}\{\langle {a}_{j}^{\dagger }(t){a}_{i}\rangle \langle {a}_{j}(t){a}_{i}^{\dagger }\rangle \}+O({J}_{{\rm{BS}}}^{3}{\tau }^{3}).\end{eqnarray} \tag{ 22 }$$

The conventional time ordered one-body correlation function is defined as ${{ \mathcal G }}_{{ij}}(t)=\langle {a}_{j}^{\dagger }(t){a}_{i}\rangle$. In our protocol, the imaginary part of the product of a particle and a hole correlation function is measured. However, we argue below that this observable carries related information as the time-ordered correlation function ${{ \mathcal G }}_{{ij}}(t);$ see also figure 11(b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

OTO correlators are measured by a straight forward extension; figure 10(b): (1) apply a beam splitter operation for a short duration $\tau {J}_{\mathrm{BS}}\ll 1$ at site i. (2) Let the system evolve for a physical time t. (3) Use single-site addressing to remove a particle on site j in both copies. (4) Flip the sign of the Hamiltonian $\hat{H}\to -\hat{H}$, as suggested in the previous section. Let the system evolve in time for the duration t. (5) Apply the 50%-50% beam splitter operation ${\overline{\mathrm{BS}}}_{i}$ on site i. Evaluating these steps, we find

$$\begin{eqnarray}{{ \mathcal F }}_{{ij}}^{\mathrm{loc}}(t) & = & \langle {\mathrm{BS}}_{i}^{\dagger }(\tau ){{\rm{e}}}^{{\rm{i}}\hat{H}t}{a}_{j}^{\dagger }{b}_{j}^{\dagger }{{\rm{e}}}^{-{\rm{i}}\hat{H}t}{\overline{\mathrm{BS}}}_{i}^{\dagger }\delta {\bar{n}}_{i}\times \,{\overline{\mathrm{BS}}}_{i}{{\rm{e}}}^{{\rm{i}}\hat{H}t}{a}_{j}{b}_{j}{{\rm{e}}}^{-{\rm{i}}\hat{H}t}{\mathrm{BS}}_{i}(\tau )\rangle \\ & = & 4{J}_{\mathrm{BS}}\tau {\rm{Im}}\{\langle {a}_{j}^{\dagger }(t){a}_{i}^{\dagger }{a}_{j}(t){a}_{i}\rangle \langle {a}_{j}^{\dagger }(t){a}_{i}{a}_{j}(t){a}_{i}^{\dagger }\rangle \}+O({J}_{\mathrm{BS}}^{3}{\tau }^{3}).\end{eqnarray} \tag{ 23 }$$

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 ${{ \mathcal F }}_{{ij}}^{\mathrm{loc}}(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 ${{ \mathcal F }}_{{ij}}^{\mathrm{loc}}(t)$ (see figure 11(a)) from which light-cone and butterfly velocities can be extracted. In figure 11(a), we once again attribute the plateau in the light-cone, which starts at $| i-j| \gtrsim 7$, to the finite MPO bond dimension of 400. Similarly, ${{ \mathcal G }}_{{ij}}^{\mathrm{loc}}(t)$ starts off at zero but then develops a peak that quickly decays; figure 11(b). From that we determine the quasiparticle lifetime $\tau J\sim 0.32$ which corresponds roughly to half the lifetime obtained for the Green's function ${{ \mathcal 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 $\delta {\bar{n}}_{i}$ determine one-body correlation functions of the original many-body state:

$$\begin{eqnarray}&&\langle \delta {\bar{n}}_{i}\delta {\bar{n}}_{j}\rangle =2\langle {a}_{i}^{\dagger }{a}_{j}\rangle \langle {a}_{i}{a}_{j}^{\dagger }\rangle .\end{eqnarray} \tag{ 24 }$$

Here we used that the left and the right initial states are identical. Higher order static correlation functions in the creation ${a}_{i}^{\dagger }$ and annihilation operators a_i are straightforwardly obtained by measuring higher order correlators in $\delta {\bar{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 as

$$\begin{eqnarray}&&\hat{\rho }={ \mathcal N }\displaystyle \sum _{{i}_{1},\ldots {i}_{N}}{r}_{{i}_{1},\ldots {i}_{N}}{\hat{\sigma }}_{{i}_{1}}\otimes ...\,\otimes \,{\hat{\sigma }}_{{i}_{N}},\end{eqnarray} \tag{ 25 }$$

where ${ \mathcal N }$ is a normalization constant and the ${\hat{\sigma }}_{{i}_{j}}$ constitute a suitable basis [54]. 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}_{{i}_{1},\ldots {i}_{N}}$ and thereby to reconstruct the density matrix, which paves the way for the full state tomography of quantum states with massive particles.

3.3.3. 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,\omega )$, which is the spatial and temporal Fourier transform of the density correlator ${C}_{n}(x,t)$, can be measured with Bragg spectroscopy [55, 56]. 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,\omega )$. Diffusion manifests itself in the wavevector and frequency resolved structure factor $S(k,\omega )$ as Lorentzian peaks with half-width-half-maximum that scales as Dq².

It is furthermore possible to measure the dissipative response $\langle [{n}_{i}(t),{n}_{j}]\rangle$ to a local perturbation of the system in a quantum gas microscope. To this end, a local potential $\delta H={n}_{j}\delta \mu$ is created at site j by applying a laser for a short time τ, yielding the time evolution $\exp [-{\rm{i}}\delta H\tau ]\sim 1-{\rm{i}}\delta H\tau +O(\delta {\mu }^{2}{\tau }^{2})$. Measuring the density at site i after the unitary time evolution for duration t we obtain

$$\begin{eqnarray}&&{\chi }_{{ij}}^{\mathrm{loc}}(t)=\langle {n}_{i}(t)\rangle +{\rm{i}}\delta \mu \tau \langle [{n}_{i}(t),{n}_{j}]\rangle +O(\delta {\mu }^{2}{\tau }^{2}).\end{eqnarray} \tag{ 26 }$$

In equilibrium, the fluctuation-dissipation theorem provides an exact relation between $\langle [{n}_{i}(t),{n}_{j}]\rangle$ and $\langle {n}_{i}(t){n}_{j}\rangle$. The accurate measurement of the former therefore enables the observation of diffusive response in the dynamical density correlator.