Timescales in the quench dynamics of many-body quantum systems: Participation ratio vs out-of-time ordered correlator

TL;DR

This paper investigates the relationship between participation ratio and out-of-time ordered correlators (OTOCs) in many-body quantum systems, revealing that while individual OTOCs do not grow exponentially, their sum decreases exponentially, providing insights into wave packet dynamics and quantum chaos.

Contribution

It demonstrates that the sum of all projection-OTOCs decreases exponentially, linking OTOC behavior to the exponential growth of the participation ratio in quantum chaos.

Findings

Sum of projection-OTOCs decreases exponentially over time

Participation ratio grows exponentially in quantum chaotic systems

Provides a new perspective on wave packet dynamics and information scrambling

Abstract

We study quench dynamics in the many-body Hilbert space using two isolated systems with a finite number of interacting particles: a paradigmatic model of randomly interacting bosons and a dynamical (clean) model of interacting spins-$1/2$. For both systems in the region of strong quantum chaos, the number of components of the evolving wave function, defined through the number of principal components $N_{pc}$ (or participation ratio), was recently found to increase exponentially fast in time [Phys. Rev. E 99, 010101R (2019)]. Here, we ask whether the out-of-time ordered correlator (OTOC), which is nowadays widely used to quantify instability in quantum systems, can manifest analogous time-dependence. We show that $N_{pc}$ can be formally expressed as the inverse of the sum of all OTOC's for projection operators. While none of the individual projection-OTOC's shows an exponential behavior, their sum decreases exponentially fast in time. The comparison between the behavior of the OTOC with that of the $N_{pc}$ helps us better understand wave packet dynamics in the many-body Hilbert space, in close connection with the problems of thermalization and information scrambling.

Figures

Equations

$$H=H_0+V,$$

$$H_0=\sum _k{E}_k^0|k⟩⟨k|$$

$$H=\sum _{\alpha }{E}^{\alpha }|\alpha ⟩⟨\alpha |,$$

$$H=\sum _{s=1}^M{ϵ}_s{a}_s^{†}{a}_s+\sum _{s_1,s_2,s_3,s_4=1}^M{V}_{s_1{s}_2{s}_3{s}_4}{a}_{s_1}^{†}{a}_{s_2}^{†}{a}_{s_3}{a}_{s_4},$$

$$\mathcal{D}=\frac{(N+M-1)!}{N!(M-1)!}.$$

$$\frac{(M-1)(M+2)}{2}\le \mathcal{N}\le N(M-1)\left[1+\frac{(N-1)(M-2)}{4}\right].$$

$$H$$

$$+$$

$$\mathcal{D}=\frac{L!}{N!(L-N)!}.$$

$$|\psi (0)⟩=|k_0⟩=\sum _{\alpha }{C}_{k_0}^{\alpha }|\alpha ⟩,$$

$$P_k(t)$$

$$=$$

$$P_{k_0}(t)$$

$${\rho }_{k_0}(E)\equiv \sum _{\alpha }|C_{k_0}^{\alpha }{|}^2\delta (E-E_{\alpha })$$

$$N_{pc}(t)=\frac{1}{{\sum }_k{P}_k(t{)}^2}=\frac{1}{{\sum }_k|⟨k|\psi (t)⟩{|}^4}.$$

$$F_{v,w}(t)=⟨{\hat{w}}^{†}(t)\hat{v}(0{)}^{†}\hat{w}(t)\hat{v}(0)⟩$$

$$|x(t)⟩=\hat{w}(t)\hat{v}(0)|k_0⟩$$

$$|y(t)⟩=\hat{v}(0)\hat{w}(t)|k_0⟩,$$

$$\begin{array}[]{lll}F_{k,k_{0}}(t)&=\left\langle k_{0}\,\right|e^{iHt}\left|k\,\right\rangle\left\langle k\,\right|e^{-iHt}\left|k^{\prime}\,\right\rangle\times\\ &\\ &\left\langle k^{\prime}\,\right|e^{iHt}\left|k\,\right\rangle\left\langle k\,\right|e^{-iHt}\left|k^{\prime}\,\right\rangle\left\langle k^{\prime}\,\right|k_{0}\rangle\\ &\\ &=\left\langle k_{0}\,\right|e^{iHt}\left|k\,\right\rangle\left\langle k\,\right|e^{-iHt}\left|k_{0}\,\right\rangle\times\\ &\\ &\hskip 10.52737pt\left\langle k_{0}\,\right|e^{iHt}\left|k\,\right\rangle\left\langle k\,\right|e^{-iHt}\left|k_{0}\,\right\rangle\\ &\\ &=|\left\langle k\,\right|e^{-iHt}\left|k_{0}\,\right\rangle|^{4}&\end{array}$$

$$[N_{pc}(t){]}^{-1}$$

$$O_{toc}(t)=\sum _{k\ne k_0}{F}_{k,k_0}(t).$$

$${\mathcal{M}}_1(k_0)=\left\{k\ne k_0,1\le k\le \mathcal{D},|⟨k|H|k_0⟩\ne 0\right\}.$$

$$W_1(t)=\sum _{k\in {\mathcal{M}}_1(k_0)}|⟨k|\psi (t)⟩{|}^2.$$

$${\left[\overline{N_{pc}^{\infty }}\right]}^{-1}=\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^{T}dt\sum _k|⟨k|e^{-iHt}|k_0⟩{|}^4$$

$$=2\sum _k{\left(\sum _{\alpha }|C_{k_0}^{\alpha }{|}^2|C_k^{\alpha }{|}^2\right)}^2-\sum _{\alpha }|C_{k_0}^{\alpha }{|}^4\sum _k|C_k^{\alpha }{|}^4.$$

$$F_{k,k_0}(t)$$

$$F_{k,k_0}(t)\simeq t^4{H}_{k,k_0}^4+o(t^6)\mathrm{for}k\in {\mathcal{M}}_1(k_0).$$

$$⟨⟨F_{k,k_0}(t)⟩⟩$$

$$F_{k,k_0}(t)\simeq \frac{1}{16}{t}^8{\left[\sum _{k^{\prime }\in {\mathcal{M}}_1}{H}_{k,k^{\prime }}{H}_{k^{\prime },k_0}\right]}^4\mathrm{for}k\in {\mathcal{M}}_2(k_0).$$

Full text

I Introduction

There is currently great interest in the study of non-equilibrium quantum dynamics of isolated systems with many interacting particles. This is partially justified by significant experimental progress that makes possible the study of the coherent evolution of many-body quantum systems for long times Kaufman et al. (2016); Lewis-Swan et al. ; Wei et al. . Yet, despite important analytical and experimental advances, several questions remain open. A timely discussion refers to the conditions Borgonovi et al. (2016); D’Alessio et al. (2016) and timescales Borgonovi et al. (2019); Schiulaz et al. ; . for the onset of equilibration and thermalization that can emerge without the influence of an environment. When studying these topics, one should distinguish systems at the thermodynamic limit, addressed by mean-field theories Eriksson et al. (2018), from systems with a finite number of particles. The latter situation emerges in experiments with cold atoms and ion traps, where the number of particles can be small and controlled.

Analytical breakthroughs in the study of many-body quantum dynamics have been recently achieved in high energy physics Maldacena et al. (2017), where quantum systems without gravity are equated to classical gravitational systems in a higher spatial dimension. A quantity that became central in many of these studies is the out-of-time-order correlator (OTOC), first introduced in the semiclassical analysis of superconductivity in Ref. Lar . Existing analytical results for the evolution of the OTOC have been obtained by taking the average in the canonical ensemble Sekino and Susskind (2008); Kitaev (a); Maldacena and Stanford (2016); Maldacena et al. (2016), thus assuming implicitly the thermodynamic limit. The present work focuses on the dynamics of finite isolated systems with interacting Bose or Fermi particles and employs the OTOC to describe the gradual spreading of the initial wave packet in the many-body Hilbert space.

The OTOC can be measured experimentally with nuclear magnetic resonance platforms and ion traps Gärttner et al. (2017); Li et al. (2017); Niknam et al. . Among various applications, it has been used to quantify the spread of quantum information Swingle (2018) and the exponential instability of quantum systems that have a chaotic classical counterpart, as supported by semiclassical analysis Rammensee et al. (2018); Rodolfo A. Jalabert (2018). This has given birth to another method to detect chaos in quantum dynamics, a goal pursued by several earlier works Peres (1996); Levstein et al. (1998); Cucchietti et al. (2002); Gorin et al. (2006); Elsayed and Fine (2015).

The quantum-classical correspondence between the exponential growth rate of the OTOC and the classical Lyapunov exponent has being numerically corroborated for finite systems with few degrees of freedom, such as one-body chaotic systems Rozenbaum et al. (2017, ) and the Dicke model with two degrees of freedom Chávez-Carlos et al. (2019). However, little is known about this correspondence for finite quantum systems with many interacting particles. Studies of the OTOC have contributed to a significant renewed interest in the problem of the quantum-classical correspondence for chaotic systems, which is a study initiated about 40 years ago with the investigation of one-body chaos.

In the paradigmatic Kicked Rotator (KR) model, it was found numerically Casati et al. (1979) and explained analytically Chirikov et al. (1981); Shepelyansky (1983) that there are two timescales on which one can speak of the quantum-classical correspondence for the dynamics of wave packets. One is the timescale due to the Ehrenfest theorem according to which the center of the wave packet in phase space follows, for some time, the corresponding classical trajectories. In the case of strong chaos, the timescale $t_{E}$ for this correspondence was analytically studied in Refs. Berman and Zaslavsky (1978); Zaslavsky (1981) and shown to be proportional to $\ln(1/\hbar)$, where $\hbar$ stands for an effective dimensionless Planck constant. The other timescale, $t_{D}$, is due to the dynamical localization occurring in the momentum space of the KR Casati et al. (1979); Chirikov et al. (1981, 1988); Izrailev (1990). The second moment of the wave packet in momentum space nicely mimics classical diffusion on the timescale $t_{D}\propto 1/\hbar^{2}$, which is much longer than $t_{E}$. It was later argued that this localization may be compared with the Anderson localization in 1D disordered models with long-range hopping Fishman et al. (1982) and the localization in quasi-1D random models described by band random matrices Casati et al. (1990, 1991); Zyczkowski et al. (1992); Izrailev (1995); Izrailev et al. (1996).

The importance of these old results obtained for the KR is two-fold. First, they show that the classical diffusion coefficient is related to the localization length of the quasienergy eigenfunctions in momentum space Chirikov et al. (1981); Shepelyansky (1986), which is a pure quantum concept. Second, they demonstrate that the timescale for the quantum-classical correspondence can be very different for different observables. As mentioned above, global observables, such as the second moment of the probability distribution in momentum space, can coincide with their classical counterparts on a timescale much larger than that defined by the Ehrenfest theorem. This point is of special relevance for studies of the evolution of observables in many-body systems. A question of particular interest is how the number $N$ of quantum particles enters the characteristic timescales involved in the scrambling of information, equilibration, and thermalization Borgonovi et al. (2019); Schiulaz et al. ; . .

It was shown in Borgonovi et al. (2019) that when the eigenstates of a many-body quantum system are strongly chaotic, the number of principal components $N_{pc}$ (or participation ratio) involved in the dynamics of the wave function in the many-body Hilbert space increases exponentially fast in time. The growth rate was found to be $2\Gamma$, where $\Gamma$ is the energy width of the strength function. This function, introduced in nuclear physics and known in solid state physics as local density of states (LDOS), is defined by projecting an unperturbed many-body state onto the basis defined by the total Hamiltonian that includes the inter-particle interaction. Knowledge of the LDOS is very important in the analysis of quench dynamics, since its Fourier transform is the survival probability, which describes the decay of the initial state.

The exponential growth of $N_{pc}$ lasts for some time $t_{S}$ before the saturation of the dynamics, which happens due to the finite size of the many-body Hilbert space. It was found in Borgonovi et al. (2019) that, for a large number of particles, $N\gg 1$, the saturation time is approximately given by $t_{S}\propto N\hbar/\Gamma$. Since $t_{S}$ is proportional to the number of particles $N$, it can be much larger than the characteristic time for the depletion of the initial state given by $\hbar/\Gamma$. The timescale $t_{S}$ represents the time for thermalization Borgonovi et al. (2019), according to which an initial wave packet ergodically fills the energy shell Casati et al. (1993, 1996); Izrailev (2001). The spread of the initial state reflects the delocalization of the energy eigenstates, which is due to the strong inter-particle interactions Santos et al. (2012a, b). These states do not fill the whole Hilbert space, just the part defined by the inter-particle interaction.

In the present work, we explore the relationship between $N_{pc}$ and a particular kind of OTOC. The former quantifies the number of unperturbed many-body states that contribute to the evolution of the wave packet, while the OTOC measures the degree of non-commutativity in time between two different Hermitian operators. In the literature, these are usually taken as local operators in real space. Here, we use instead projection operators in the many-body Hilbert space, which are local in this space. We show that the inverse of the sum of all OTOC’s coincides with $N_{pc}$.

In our analysis, we distinguish between two categories of OTOC’s: the autocorrelator, where both projections are made on the initial state, and the case involving a projection onto a many-body state other than the initial state, referred to as projection-OTOC. While the autocorrelator decays exponentially as $e^{-2\Gamma t}$, we find that a single projection-OTOC does not exhibit exponential behavior. However, when we look at the sum of all projection-OTOC’s, we find a non-monotonic behavior in time, where an initial growth is followed by an exponential decay. This decay happens within the time interval of the exponential increase of $N_{pc}$.

We consider two models, the well-known two-body random ensemble (TBRE) with a finite number of bosons interacting randomly and a dynamical (deterministic) one-dimensional (1D) spin-$1/2$ model with nearest and next-nearest neighbor couplings only. The TBRE falls into the broader category of the so-called embedded ensembles, which have been thoroughly studied since the 1970’s in the context of nuclear physics and quantum chaos Brody et al. (1981); Kota (2001, 2014). The Sachdev-Ye-Kitaev (SYK) models Kitaev (b); Sachdev and Ye (1993), which have received increasing attention in high energy physics, are also examples of embedded random ensembles. For both models that we study, we choose parameters for which the eigenstates involved in the dynamics are composed by a very large number of unperturbed many-body states.

The paper is organized as follows. In Sec. II, we describe the two models considered. Section III presents the relationship between OTOC and $N_{pc}$. In Sec. IV, we show analytical as well as numerical results for both the TBRE and the spin model. In Sec. V, we summarize our results and discuss some possible future directions.

II Models and Quench Dynamics

We consider a bosonic TBRE and a 1D spin-1/2 system, both of them described by the Hamiltonian

[TABLE]

where

[TABLE]

stands for the unperturbed (integrable) part of the total Hamiltonian $H$, with

[TABLE]

and $V$ represents the two-body interactions. In what follows we set $\hbar=1$. We focus on the case where the perturbation $V$ is sufficiently strong, so that a large part of the energy spectrum of $H$ contains chaotic eigenstates.

Since our study concentrates on the dynamics occurring in the unperturbed many-body space of chaotic systems, a definition of what we mean by quantum chaos is in order. For one-body systems, it is common lore to associate quantum chaos with level statistics described by full random matrices. However, in realistic finite many-body models, not all eigenstates are random vectors, as in full random matrices, and not all of them are involved in the dynamics. Therefore, spectrum statistics obtained by taking into account all eigenvalues is not the best way to characterize the dynamics, which is only due to those eigenstates that are present in an initially excited wave packet. Our approach to quantum chaos is linked with the structure of the eigenstates. They are called chaotic when they are fully delocalized in the energy shell and are composed of many uncorrelated components (see, for example, Refs. Santos et al. (2012a, b)).

II.1 Two-body random ensemble

The TBRE describes $N$ identical bosons occupying $M$ single-particle levels; the latter are specified (and reordered) by random energies $\epsilon_{s}$. The mean spacing $\langle\epsilon_{s}-\epsilon_{s-1}\rangle\equiv\delta=1$ sets the energy scale defining the width of the unperturbed energy spectrum, $NM\delta$. The choice to have random single-particle energies is not a necessary condition for the results obtained below. It is used to remove the degeneracy in the unperturbed many-body spectrum.

The Hamiltonian of the TBRE is written as,

[TABLE]

where $a_{s}$ ($a_{s}^{\dagger}$) is the annihilation (creation) operator on the single-particle energy level $\epsilon_{s}$, so the number operator $n_{s}=a^{\dagger}_{s}a_{s}$ gives the probability for the occupation of the $s$-th single-particle energy level, $n_{s}/N$. The two-body matrix elements $V_{s_{1}s_{2}s_{3}s_{4}}$ are Gaussian random entries with zero mean and variance ${\cal V}^{2}$. The Hamiltonian conserves the total number of bosons, so the analysis is done for a single subspace of dimension

[TABLE]

Throughout the paper, we fix the number of single-particles levels, $M=11$, and we vary the number of particles $N$ from 4 to 8. That corresponds to a size ${\cal D}$ of the many-body space ranging from 1001 up to 43758. The strength $V$ of the inter-particle interaction is chosen so that ${\cal V}=0.4$ to have a large energy region with strongly chaotic eigenstates Borgonovi and Izrailev (2017). The eigenstates $\left|k\,\right\rangle$ of $H_{0}$ constitute the unperturbed many-body basis (also called mean-field basis) in which we study the dynamics of the wave packets and in second quantized form they can be written as $\left|n_{1},...,n_{s},...,n_{M}\,\right\rangle$ where $n_{s}$ is the number of bosons in the $s$-th single-particle energy level.

The TBRE Hamiltonian matrix is very sparse, because only a fraction of the unperturbed many-body states of $H_{0}$ are directly connected by the two-body interaction $V$. The number of non-zero off-diagonal matrix elements ${\cal N}$ depends on the particularly chosen matrix line, but it is generally much smaller than the total matrix dimension $\cal{D}$. It is not possible to give a general analytical expression for ${\cal N}$, but upper and lower bounds as a function of $N,M$ have been estimated as follows Borgonovi and Izrailev (2017),

[TABLE]

In particular, the minimal number of directly coupled states, which is independent of $N$, is obtained when all $N$ particles occupy only one single-particle energy level. Another feature of the TBRE matrices is their band-like structure, which causes the eigenstates close to the ground state to be much less delocalized than the states closer to the center of the spectrum.

The TBRE was originally developed to explain the statistical properties of complex systems with interacting Fermi-particles, such as highly excited nuclei and molecules French and Wong (1970); Brody et al. (1981). It was later applied to systems of interacting bosons, to which, in the dilute limit, many aspects of energy spectra and eigenstates are similar to those of systems of random interacting fermions. To date, it has been extensively investigated for fermions Flambaum and Izrailev (1997); Altshuler et al. (1997) and for bosons Kota (2001); Kota and Sahu (2001); Kota (2014); Benet and Weidenmüller (2003). This model is a particular case of the embedded ensembles with $q$-body interactions. When $q=2$ we have the TBRE and when $q=N$, we recover the full random matrices.

In contrast to the standard ensembles of full random matrices, TBREs are much closer to realistic physical systems, since they take into account the two-body nature of the interactions, the type of interacting particles (fermions or bosons), the strength of the inter-particles interaction, and the properties of single-particle spectra.

II.2 Dynamical spin-1/2 model

The 1D spin-1/2 model that we study here is dynamical, that is it has no random elements. The Hamiltonian is given by

[TABLE]

The first part of this Hamiltonian contains only nearest-neighbor couplings and it is associated with the mean field $H_{0}$. The second part describes next-nearest-neighbor couplings and represents the perturbation $V$. Differently from the previous model, $V$ is a local interaction in space. The Pauli matrices $\sigma^{x,y,z}_{s}$ act on site $s$; $L$ is the number of sites which is chosen even; the coupling constant $J=1$ sets the energy scale; $\Delta$ stands for the anisotropy of the interaction, and $\lambda$ is the ratio between next-nearest-neighbor and nearest-neighbor couplings Santos (2009); Torres-Herrera et al. (2015).

The Hamiltonian conserves the total spin in the $z$-direction, ${\cal S}^{z}=\sum_{s=1}^{L}\sigma_{s}^{z}/2$. In what follows we consider the subspace ${\cal S}^{z}=-1$, which has $N=L/2-1$ excitations (up-spins) and dimension

[TABLE]

The unperturbed Hamiltonian $H_{0}$ is integrable, but as $\lambda$ increases, $H$ crosses over to the chaotic regime Santos et al. (2012a, b). For the parameters considered here, system size $L=16$, number of up-spins $N=7$, (so ${\cal D}=11440$), anisotropy $\Delta=0.48$, and $\lambda=1$, the model is strongly chaotic in a large region of the spectrum.

II.3 Quench Dynamics

To study the dynamics, we prepare the system in an unperturbed state $\left|k_{0}\,\right\rangle$,

[TABLE]

where $C_{k_{0}}^{\alpha}=\langle\alpha|k_{0}\rangle$ and $\left|\alpha\,\right\rangle$ are the exact energy eigenstates. The initial state $\left|\psi(0)\,\right\rangle$ evolves under the full Hamiltonian $H$ when the interaction $V$ is turned on. We consider initial states that have energy $E_{k_{0}}=\left\langle k_{0}\,\right|H\left|k_{0}\,\right\rangle$ away from the edges of the spectrum of $H$.

We notice that the initial state for the spin model is not a site-basis vector (computational basis vector) for which the spin on each site either points up or down in the $z$-direction, but it is instead an eigenstate of $H_{0}$. In analogy with the TBRE, we refer to these states as the unperturbed many-body basis.

The probability to find the evolved state in a basis state $\left|k\,\right\rangle$ at the time $t$ is given by

[TABLE]

The particular case where $k=k_{0}$ corresponds to the survival probability (also known as return probability), which can be written as

[TABLE]

where

[TABLE]

is the LDOS, that is the energy distribution weighted by the components $|C^{\alpha}_{k_{0}}|^{2}$ of the initial state. The subscript $k_{0}$ in Eq. (10) stresses the important point that the LDOS depends on the initial state $\left|k_{0}\,\right\rangle$. As evident from Eq. (9), the survival probability is the Fourier transform of the LDOS. The inverse of the width $\Gamma$ of the LDOS gives the characteristic decay time of $P_{k_{0}}(t)$.

The maximal size of the LDOS, obtained when $H_{0}$ is negligible and $H\sim V$, defines the energy shell, which is only a part of the total energy spectrum. The shape of the energy shell depends on the density of states, which in systems with few-body interactions typically has a Gaussian form Brody et al. (1981). The eigenstates of $H$ written in the unperturbed basis are chaotic when they fill the energy shell completely and the components $C_{k}^{\alpha}$ are random numbers following the Gaussian envelope of the energy shell Santos et al. (2012a, b).

To quantify how the initial state spreads in time, in the many-body Hilbert space, we compute the number of principal components (also known as participation ratio),

[TABLE]

For the TBRE, we use the notation $\langle\langle N_{pc}(t)\rangle\rangle$ to indicate average over the random configurations of the two-body interaction.

III OTOC for projection operators and number of principal components

The OTOC for two Hermitian operators ${\hat{w}}$ and ${\hat{v}}$ is defined as,

[TABLE]

where ${\hat{w}}(t)=e^{iHt}{\hat{w}}(0)e^{-iHt}$ is the operator in the Heisenberg representation. In the literature, $\langle.\rangle$ originally referred to the average over the canonical ensemble, but later, averages over all states of an unperturbed Hamiltonian or over one particular initial state $\left|k_{0}\,\right\rangle$, as we do here, have also been considered.

Written in terms of the initial state, the OTOC has a clear physical meaning, which can be explained as follows. Let us define the two states,

[TABLE]

and

[TABLE]

which represents the action of the two operators taken in the reversed order. The state $\left|x(t)\,\right\rangle$ is obtained by first applying ${\hat{v}}$, then evolving forward with the full Hamiltonian for time $t$, applying ${\hat{w}}$, and finally evolving backward for the same time $t$. For $\left|y(t)\,\right\rangle$, the order is exchanged: first the evolution is forward, then ${\hat{w}}$ is applied, followed by the backward evolution, and finally the application of ${\hat{v}}$. Thus, $F_{v,w}(t)$ quantifies the decay of the overlap between these two states, $\langle y(t)|x(t)\rangle$, caused by the exchanged action of the two operators ${\hat{v}}(0)$ and ${\hat{w}}(t)$. It probes the way ${\hat{v}}$ and ${\hat{w}}$ inhibit the cancellation between forward and backward evolution. Equivalently, $F_{v,w}(t)$ measures the degree of non-commutativity between the two operators.

The OTOC is related to the $N_{pc}$ when in Eq. (12) we use projection operators in the unperturbed many-body states, ${\hat{w}}(0)=\left|k\,\right\rangle\left\langle k\,\right|$, ${\hat{v}}(0)=\left|k^{\prime}\,\right\rangle\left\langle k^{\prime}\,\right|$, and compute the expectation value in the initial state $\left|k_{0}\,\right\rangle$. This gives,

[TABLE]

Since ${\hat{v}}(0)\left|k_{0}\,\right\rangle=\delta_{k^{\prime},k_{0}}\left|k^{\prime}\,\right\rangle$, it is clear that to have a non-zero correlation function one needs to choose ${\hat{v}}(0)=\left|k_{0}\,\right\rangle\left\langle k_{0}\,\right|$. Comparing the equation above with Eq. (11), one sees that

[TABLE]

In the above, we separate $k=k_{0}$ from $k\neq k_{0}$. We refer to $F_{k,k_{0}}(t)$ for $k\neq k_{0}$ as projection-OTOC’s, while the autocorrelation function $F_{k_{0},k_{0}}(t)=\left\langle k_{0}\,\right|e^{-iHt}\left|k_{0}\,\right\rangle|^{4}=P_{k_{0}}(t)^{2}$ is simply the squared survival probability. We denote by $O_{toc}(t)$ the extensive sum over all projection-OTOC’s,

[TABLE]

The inverse of the $N_{pc}$ is therefore $O_{toc}(t)$ plus the squared survival probability.

IV Analytical estimates and numerical results

We now have the tools to compare the results for the $N_{pc}$ and the OTOC for the TBRE and the dynamical spin-1/2 model in the strongly chaotic regime. As mentioned above, the initial states have energy $E_{k_{0}}=\left\langle k_{0}\,\right|H\left|k_{0}\,\right\rangle$ far from the edges of the spectrum.

IV.1 TBRE: Number of principal components and OTOC

For the TBRE, we focus on initial states, where all particles are on a single level, which we choose to be the fifth level, such as in $\left|0,0,0,0,N,0,0,0,0,0\,\right\rangle$. States of this kind have $E_{k_{0}}$ close to the center of the band. This choice of initial state is made, because the number of directly coupled matrix elements is minimal and independent of $N$. The number of states directly coupled with the initial state together with the strength of the perturbation determine the width of the LDOS and thus the decay rate of the survival probability.

In Fig. 1, we confirm that for the chosen perturbation and initial states, the survival probability decays exponentially and the decay rate is approximately independent of the number of particles. Needless to say, for very short time, $t\ll\Gamma^{-1}$, the survival probability decays quadratically in time, as given by perturbation theory. This behavior is subsequently followed by a region of exponential decay with rate $\Gamma$, as seen in Fig. 1. This rate defines the timescale $t_{\Gamma}=1/\Gamma$ for the depletion of the initial state Borgonovi et al. (2019). At this point, the probability to be in the initial state is reduced by a factor $1/e$.

IV.1.1 Number of principal components

The parameter $\Gamma$ is at the basis of a phenomenological cascade model Borgonovi et al. (2019), that describes in a coarse-grained way the spreading of the initial many-body state in the many-body Hilbert space. The basic idea is to analyze the dynamics at different time steps, each being associated with the probability to find the system in a specific subset of unperturbed many-body states, referred to as a “class”. The class that contains only the initial state is the ${\cal M}_{0}(k_{0})$ class and the probability to be in this class is just the survival probability $P_{k_{0}}(t)$. ${\cal M}_{1}(k_{0})$ is the set of all unperturbed states directly coupled to the initial state,

[TABLE]

The probability to be in this class is defined as

[TABLE]

The subset with states coupled to $\left|k_{0}\,\right\rangle$ in second order of perturbation theory is ${\cal M}_{2}(k_{0})$, and so on. This description of the dynamics in terms of the spread of the wave packet in the many-body Hilbert space was also explored in Altshuler et al. (1997); Flambaum and Izrailev (2001). With this picture, we obtained in Borgonovi et al. (2019) approximate rate equations for the probability to find the system in each class. The sum of the square of these probabilities gives the inverse of the number of principal components $N_{pc}$. Our analysis predicted an exponential growth for $N_{pc}$ with exponent $2\Gamma$, which was verified numerically. This is shown in Fig. 2(a) for different initial states with increasing number of particles.

It is important to remark that the exponential increase of the number of principal components continues beyond $t_{\Gamma}$. At long times, since the many-body Hilbert space is finite, $N_{pc}(t)$ finally saturates to an equilibrium value, which is obtained by taking the infinite time average,

[TABLE]

An estimate of the saturation time $t_{S}$ can be obtained by equating $e^{2\Gamma t_{S}}\simeq\overline{N_{pc}^{\infty}}$. We showed in Ref. Borgonovi et al. (2019) that for $M,N\gg 1$, this estimate is given by $t_{S}\sim Nt_{\Gamma}$. This result is seen clearly in Fig. 2 (b), together with a linear fit. The values for $t_{S}$ are obtained from the intersections in Fig. 2 (a) between the exponential curve and the horizontal lines, which indicate the saturation values from Eq. (17). We note that the saturation time $t_{S}$ was shown to coincide with the time necessary for the onset of the Bose-Einstein distribution for single-particle occupation numbers (for details see Borgonovi and Izrailev (2019)). One can therefore identify $t_{S}$ with the thermalization time.

IV.1.2 Out-of-time ordered correlator

We now proceed with the analysis of the OTOC and comparison with $N_{pc}$. The OTOC behavior at short time can be obtained with the expansion,

[TABLE]

where $H_{k,k_{0}}=\langle k|H|k_{0}\rangle$. For $k\neq k_{0}$, there are different behaviors, as listed below.

(i) The first one corresponds to $k\in{\cal M}_{1}(k_{0})$, for which one gets,

[TABLE]

Taking the average over disorder realizations in the TBRE, we come to the following estimate,

[TABLE]

To obtain the last line above, we took into account that $H_{k,k^{\prime}}$ are Gaussian variables with zero mean and variance ${\cal V}^{2}$.

(ii) For the case $k\in{\cal M}_{2}(k_{0})$, one has a $t^{8}$ behavior,

[TABLE]

(iii) For the projection-OTOC’s of higher-order classes, where $\langle k|H|k_{0}\rangle=\langle k|H^{2}|k_{0}\rangle=0$, the initial numerical power-law growth gives a $t^{12}$ behavior.

The behaviors $t^{4}$, $t^{8}$ and $t^{12}$ for the various projection-OTOC’s are shown in Fig. 3, respectively as dashed, full and dot-dashed lines. Perturbation theory is approximately valid for $t<t_{\Gamma}$. In the region marked by the exponential growth of the $N_{pc}$, that is $t_{\Gamma}<t<t_{S}$, the OTOC’s have a non-generic and non-monotonous behavior. For $t>t_{S}$, the OTOC’s just show fluctuations around some equilibrium value.

In Fig. 4, we examine the behavior of the sum of all projection-OTOC’s [Eq. (15)]. Our figure shows the time dependence of the $\langle\langle O_{toc}(t)\rangle\rangle$ for different numbers of particles. We can see that it reaches a maximum approximately at $t_{\Gamma}$ (vertical orange line), when the probability to be in the initial state is reduced by a factor $1/e$. After this point, $\langle\langle O_{toc}(t)\rangle\rangle$ decays exponentially, with an exponent between $\Gamma$ and $2\Gamma$ (actually $1.2\Gamma$ for this set of initial states). This exponent comes out from the sum of many different contributions from states belonging to different classes, and it cannot be obtained by taking into account the first-class states only. We note that extensive sums of local operators were also used in the analysis of the OTOC in Ref. Kukuljan et al. (2017), where it is argued that only the sum, and not a single local observable, can exhibit indefinite exponential growth in the thermodynamic limit.

We do not have yet a theory to extract the exponential decay rate for $O_{toc}(t)$. It should be possible to associate the characteristic decay time for the sum $\sum_{k\in{\cal M}}F_{k,{k_{0}}}(t)$ of projection-OTOC’s that belong to a specific class ${\cal M}$ to the scrambling time of the correlations during the flow from one class to the other. The timescale $t_{S}$ would emerge as a result of the summation of all different timescales associated to all classes. We leave this study to a future work. We note that the exponential decay of the out-of-time order correlators was recently obtained analytically for chaotic quantum maps Garcia-Mata et al. . In that work, the approach to the stationary value was found to occur with a rate determined by the Ruelle-Pollicot resonances.

The exponential decay of $\langle\langle O_{toc}(t)\rangle\rangle$ for $t_{\Gamma}<t<t_{S}$ indicates that the extensive sum of OTOC’s plays an important role in the exponential growth of the number of principal components beyond $t_{\Gamma}$. In Fig. 5, we compare the two terms appearing in the denominator of $N_{pc}$, that is $\langle\langle O_{toc}(t)\rangle\rangle$ and $\langle\langle P_{k_{0}}(t)^{2}\rangle\rangle$, for the case with $N=8$ particles. Initially $\langle\langle N_{pc}(t)^{-1}\rangle\rangle$ is entirely dominated by the squared survival probability. Later, due to the different decay rates for $\langle\langle O_{toc}(t)\rangle\rangle$ and $\langle\langle P_{k_{0}}(t)^{2}\rangle\rangle$, these two contributions become of the same order of magnitude and they eventually cross.

As seen in Fig. 5, for the system size and set of initial states considered, the crossing between $\langle\langle O_{toc}(t)\rangle\rangle$ and $\langle\langle P_{k_{0}}(t)^{2}\rangle\rangle$ occurs after the saturation time $t_{S}$. As a result, the relaxation of $\langle\langle N_{pc}(t)\rangle\rangle^{-1}$ to its infinite time-average value is entirely due to the saturation of $\langle\langle O_{toc}(t)\rangle\rangle$. The two saturate roughly at the same time. In contrast, the squared survival probability reaches its stationary value at a timescale much larger than $t_{S}$.

Figure 6 illustrates the timescale for the relaxation of the survival probability. By comparing this time with the saturation time $t_{S}$ for $N_{pc}$ shown in Fig. 2, we can see that the former is more that two orders of magnitude larger. This is due to the presence of the so-called correlation hole (see Torres-Herrera and Santos (2017a, b); Torres-Herrera et al. (2018) and references therein), which is a dip below the saturation value. This hole is clearly visible for the survival probability, but it is not so evident for $\langle\langle N_{pc}(t)\rangle\rangle$ (for a comparison see Ref. Schiulaz et al. ). The correlation hole ends at the Heisenberg time, beyond which there are only fluctuations around the infinite-time average, given by $\sum_{\alpha}|C_{k_{0}}^{\alpha}|^{4}$.

IV.2 Spin-1/2 model: Number of principal components and OTOC

For the spin model, we do not perform any average, since the Hamiltonian has no random elements and a single initial state with energy $E_{k_{0}}\approx-0.5$ is considered. The results are very similar to those presented in Fig. 2, Fig. 3, and Fig. 4.

Figure 7 (a) shows the number of principal components, which grows as $e^{2\Gamma t}$ in the time interval $t_{\Gamma}<t<t_{S}$. In Fig. 7 (b), we depict the behavior of some projection-OTOC’s. They show power-law growths proportional to $t^{4}$ and $t^{8}$ for $t<t_{\Gamma}$, as seen also in Fig. 3. The behaviors become non-monotonic for $t_{\Gamma}<t<t_{S}$. From the figure, it is clear that states belonging to the first class (those having a $t^{4}$ initial growth) reach their maximal value before the states in the second class (those with a $t^{8}$ behavior). Since they reach the maximum at different times, they start to decay at different times, so we might expect a complicated behavior in the time region $t_{\Gamma}<t<t_{S}$. However, as clear from Fig. 7 (c), in this time interval, the extensive sum of all projection-OTOC’s actually decays exponentially before saturation, with $\alpha=0.66$ in $e^{-\alpha\Gamma t}$. The result is similar to the one observed in Fig. 3 for the TBRE.

The results for the spin model corroborate that for $t_{\Gamma}<t<t_{S}$, the sum given by $O_{toc}(t)$ contributes to the exponential behavior of $N_{pc}$, despite the fact that individually, the projection-OTOC’s do not show any sign of exponential behavior in this time interval. We find a different decay exponent $\alpha$ from TBRE case. It is not clear at this point what this exponent might depend on, such as number of particles, energy of the initial state, and connectivity of the model. We leave this point for future investigations.

We notice that even though $H_{0}$ for the spin model can be solved with the Bethe ansatz, this is not at all trivial. Thus, we obtain numerically the eigenstates $\left|k\,\right\rangle$, used as the basis to write $H$. As a result, all matrix elements of $H$ become non-zero. To identify which elements correspond to effective couplings between the unperturbed states, we use a threshold $\xi=0.1$, that is, we assume that $\left|k\,\right\rangle$ is directly coupled with $\left|k^{\prime}\,\right\rangle$ only if $H_{k,k^{\prime}}>\xi|H_{k,k}-H_{k^{\prime},k^{\prime}}|$.

V Discussion

We studied the relationship between the out-of-time ordered correlator (OTOC) and the number of principal components $N_{pc}$ (or participation ratio), and their relevance to the relaxation process of many-body quantum systems. Two chaotic models were considered: One model belongs to the two-body random ensemble (TBRE), where randomness is introduced ad hoc as random couplings between many-body unperturbed states, and the other is a clean system of spin-$1/2$ particles on a linear chain with non-random two-body interactions.

In a recent work Borgonovi et al. (2019), we had shown that, starting with a single many-body state of the unperturbed Hamiltonian $H_{0}$, the effective number of unperturbed many-body states participating in the dynamics, dictated by the perturbed Hamiltonian $H=H_{0}+V$, increases exponentially in time. This happens when the inter-particle interactions are sufficiently strong and the many-body eigenstates are superpositions of many effectively pseudo-random components, which is a main feature of strong quantum chaos. The quantity employed to characterize the spread of the initial wave packet in the Hilbert space was the number of principal components $N_{pc}$.

For strong perturbation, namely $H_{0}\sim V$, we found that $N_{pc}(t)$ increases as $e^{2\Gamma t}$, where $\Gamma$ is the width of the LDOS. Our numerical data, as well as the analytical estimates, showed that this exponential behavior holds up to the saturation time $t_{S}\sim Nt_{\Gamma}$, where $N$ is the number of particles for the TBRE and number of excitations for the spin model. This timescale is larger than the time $t_{\Gamma}\sim 1/{\Gamma}$ for the effective decrease of the survival probability.

In the present paper, we showed that $N_{pc}$ is the square of the survival probability plus the sum of all projection-OTOC’s. For the latter, the operators ${\hat{w}}$ and ${\hat{v}}$ in Eq. (12) are projection operators in the many-body Hilbert space, ${\hat{w}}$ being the projection on a state other than the initial state.

Our semi-analytical description of $N_{pc}(t)$ was based on the spread of the initial wave packet into different classes of unperturbed many-body states. At the shortest timescale, only the many-body states of $H_{0}$ directly coupled to the initial state by the two-body interactions get excited. Later in time, the wave packet propagates to those states which are coupled to the initial state in the second order of perturbation theory, and even later, higher orders are reached successively. This dynamics may be compared with the spread (mixing) of packets of classical trajectories in phase space: initially the whole phase space is scarcely occupied, but as time grows it gets more densely occupied. Within this picture, the projection-OTOC’s describe the flow of the wave packet probability between specific classes. At short time, each one increases as $t^{4},t^{8},t^{12}$, depending on the class the OTOC is associated with, and in accordance with our analytical estimates. After reaching a maximal value, the projection-OTOC’s decay to a stationary value given by the infinite-time average value. In the course of this process, none of the individual projection-OTOC’s shows an exponential behavior. It is only the sum of the projection-OTOC’s over all classes that decays exponentially for $t>t_{\Gamma}$. This non-monotonic behavior contrasts with that for the autocorrelation function (squared survival probability), which decays as $e^{-2\Gamma t}$ already at short times.

It should be possible to associate to the sum of the projection-OTOC’s belonging to a specific class ${\cal M}$, a characteristic decay time that represents the scrambling time for that class. The saturation of the entire dynamics at $t_{S}$ happens after the saturation of the projection-OTOC’s for all classes. After the time $t_{S}\sim Nt_{\Gamma}$, the system is fully equilibrated (thermalized) in a finite but very large domain of the unperturbed basis.

We finish this conclusion with a discussion about the quantum-classical correspondence for chaotic many-body systems. For this, we recall that the LDOS, which has width $\Gamma$, has a well defined classical limit with width $\Gamma_{cl}$ F. Borgonovi and Izrailev (1998); G.A.Luna-Acosta et al. (2000); Izrailev (2001); Luna-Acosta et al. (2001, 2002). Our results show that for $N_{pc}(t)$, which is a global observable, the timescale $t_{S}$ over which one can speak of exponential instability diverges in the thermodynamic limit, provided the semiclassical limit $\Gamma\to\Gamma_{cl}$ is done before $N\rightarrow\infty$. This suggests that there may be global observables for which the quantum-classical correspondence remains indefinitely in the thermodynamic limit.

The divergence of $t_{S}$ does not contradict the conventional picture of the Ehrenfest theorem, according to which the timescale of the quantum-classical correspondence for one-body chaotic systems is very small, $t_{E}\sim\ln(1/\hbar)$. As shown for the KR model, this is the timescale for a local observable, but there is another timescale, $t_{D}\sim 1/\hbar$, corresponding to the dynamical localization in the momentum space, which is related to a global observable. Therefore, the timescales for the quantum-classical correspondence depend on the choice of the observable and can vary significantly from one observable to another. Our study for many-body models focused on the global observable $N_{pc}$, rather than on local observables. There is not yet any direct comparison between $N_{pc}$ and a classical analog. We suggested in Ref. Borgonovi et al. (2019) that such comparison will have to be done with the use of the Kolmogorov-Sinai entropy, which is the main characteristic of the dynamics for classical many-body systems, whose dynamics occurs in a $2N$ dimensional phase space.

One should mention that the quantum diffusion in the KR is not a “true” diffusion as that occurring in classical systems. As shown in Shepelyansky (1983), the quantum diffusion is completely reversible, despite the presence of small, but finite errors associated with any numerical calculation. This is at variance with classical diffusion, which is non-reversible due to the exponential sensitivity with respect to unavoidable computation errors. This is a distinctive property of the observed quantum-classical correspondence for the wave packet width in the momentum space. One can conjecture that a similar picture should arise for many-body chaos. Even though the quantum-classical correspondence may look very good for global observables (for the number of principal components in our case), quantum properties such as local quantum correlations and entanglement may still be present during the relaxation process and even at thermalization. In fact, it was recently shown numerically and semi-analytically in Ref. Borgonovi and Izrailev (2019) that the Bose-Einstein distribution for occupation numbers emerges on the same timescale as the thermalization time $t_{S}$. This implies the coexistence of classical and quantum features in the dynamics on a very large timescale $t>t_{S}\sim N/\Gamma$. The quantum-classical correspondence for many-body systems is a challenging problem that requires further studies.

Acknowledgements.

F.B. acknowledges support by the I.S. INFNDynSysMath. F.M.I. acknowledges financial support from CONACyT (Grant No. 286633). L.F.S. is supported by the U.S. National Science Foundation (NSF) Grant No. DMR-1603418.