Generic dynamical features of quenched interacting quantum systems: Survival probability, density imbalance and out-of-time-ordered correlator

TL;DR

This paper investigates the universal dynamical features of quenched many-body quantum systems, revealing that certain decay behaviors and correlation phenomena are common to both random matrix models and realistic disordered systems.

Contribution

It provides analytical expressions for key dynamical quantities and demonstrates their applicability to realistic disordered models, highlighting universal features of nonintegrable quantum systems.

Findings

Power-law decay at intermediate times

Presence of correlation holes at long times

Shared dynamical features between random matrices and disordered models

Abstract

We study numerically and analytically the quench dynamics of isolated many-body quantum systems. Using full random matrices from the Gaussian orthogonal ensemble, we obtain analytical expressions for the evolution of the survival probability, density imbalance, and out-of-time-ordered correlator. They are compared with numerical results for a one-dimensional disordered model with two-body interactions and shown to bound the decay rate of this realistic system. Power-law decays are seen at intermediate times and dips below the infinite time averages (correlation holes) occur at long times for all three quantities when the system exhibits level repulsion. The fact that these features are shared by both the random matrix and the realistic disordered model indicates that they are generic to nonintegrable interacting quantum systems out of equilibrium. Assisted by the random matrix analytical results, we propose expressions that describe extremely well the dynamics of the realistic chaotic system at different time scales.

Equations

$$W_{n_0}(t)={\left|⟨\Psi (0)|\Psi (t)⟩\right|}^2={\left|\sum _{\alpha }{\left|C_{n_0}^{(\alpha )}\right|}^2{e}^{-i{E}_{\alpha }t}\right|}^2,$$

$$I(t)=\frac{4}{L}\sum _{k=1}^L⟨\Psi (0)|S_k^z(0)S_k^z(t)|\Psi (0)⟩.$$

$$O_{toc}(t)=\frac{32(L-2)!}{L!\mathcal{N}}\sum _{n,k,k^{\prime }}⟨{\varphi }_n|S_{k^{\prime }}^z(t)S_k^z(0)S_{k^{\prime }}^z(t)S_k^z(0)|{\varphi }_n⟩,$$

$$\rho (E)=\frac{2\mathcal{N}}{\pi \epsilon }\sqrt{1-{\left(\frac{E}{\epsilon }\right)}^2},$$

$$W_{n_0}^{\text{FRM}}(t)=\frac{1-{\overline{W}}_{n_0}^{\text{FRM}}}{\mathcal{N}-1}\left[4\mathcal{N}\frac{{\mathcal{J}}_1^2(\epsilon t)}{(\epsilon t{)}^2}-b_2\left(\frac{\epsilon t}{4\mathcal{N}}\right)\right]+{\overline{W}}_{n_0}^{\text{FRM}},$$

$$⟨W_{n_0}(t)⟩=\frac{1-⟨{\overline{W}}_{n_0}⟩}{\mathcal{N}-1}\left[\mathcal{N}\frac{g(t)}{g(0)}-b_2\left(\frac{wt}{\mathcal{N}}\right)\right]+⟨{\overline{W}}_{n_0}⟩,$$

$$I^{\text{FRM}}(t)=\frac{I(0)-{\overline{I}}^{\text{FRM}}}{\mathcal{N}-1}\left[4\mathcal{N}\frac{{\mathcal{J}}_1^2(\epsilon t)}{(\epsilon t{)}^2}-b_2\left(\frac{\epsilon t}{4\mathcal{N}}\right)\right]+{\overline{I}}^{\text{FRM}},$$

Full text

Introduction. Nonequilibrium dynamics of isolated many-body quantum systems is a highly interdisciplinary subject covering a broad range of physics scales, from string theory and black holes to condensed matter and atomic physics. The connection between black hole physics and unitary quantum dynamics emerges from holographic dualities Maldacena (1998). On the experimental side, unitary quantum dynamics is investigated with cold atoms Kinoshita et al. (2006); Schreiber et al. (2015); Bordia et al. ; Kaufman et al. (2016), ion traps Jurcevic et al. (2014); Smith et al. (2016), and nuclear magnetic resonance platforms Gärttner et al. ; Wei et al. .

Driven by different purposes, studies of black hole information loss Roberts et al. (2015); Papadodimas and Raju (2015); Cotler et al. , quantum chaos Roberts and Stanford (2015); Maldacena et al. (2016), thermalization in isolated quantum systems Kinoshita et al. (2006); Kaufman et al. (2016); Borgonovi et al. (2016), many-body localization Schreiber et al. (2015); Fan et al. ; Wei et al. , quantum correlations Gärttner et al. , and quantum speed limits Mandelstam and Tamm (1945); Margolus and Levitin (1998); Muga et al. (2009) consider similar dynamical quantities. They include the survival probability, density imbalance, and out-of-time-ordered correlator (OTOC). Our goal is to characterize the evolution of these quantities at different time scales.

Given the complexity of out-of-equilibrium many-body quantum systems, we take the same approach as Wigner when studying heavy nuclei and use full random matrices (FRM) from the Gaussian orthogonal ensemble (GOE). These are matrices filled with random real numbers and constrained by time-reversal symmetry. The model is unrealistic, as it assumes simultaneous and infinite-range interactions among all particles. But it allows for the derivation of analytical expressions for the observables of interest.

The analysis of the FRM model assists in the identification of general features and bounds for the evolution of realistic systems. The analytical expressions obtained with FRM reveal different behaviors at different time scales. After determining the generic causes of these behaviors, one can propose expressions for the dynamics of realistic chaotic many-body quantum systems.

We compare the analytical expressions for FRMs with numerical results for the one-dimensional (1D) Heisenberg spin-1/2 model with on-site disorder. This system has been extensively studied in the context of many-body localization Santos et al. (2004); Dukesz et al. (2009); Nandkishore and Huse (2015). It shows a chaotic regime for small disorder Avishai et al. (2002); Santos (2004), which justifies the comparison with FRMs. The rate of the evolution is faster in the FRM case, but the overall dynamical behavior is similar for both models.

The basis of our analysis is the survival probability. It gives the probability of finding the initial state later in time and has been investigated since the early days of quantum mechanics Khalfin (1958). It is a main quantity in the studies of quantum speed limits Muga et al. (2009) and decay processes of unstable systems Fonda et al. (1978). More recently, it became central to the analysis of localization in noninteracting Ketzmerick et al. (1992); Huckestein and Schweitzer (1994) and interacting Torres-Herrera and Santos (2015, 2017) systems. The survival probability is also related del Campo et al. to the analytic continuation of the partition function used to study conformal field theories with holographic duals Dyer and Gur-Ari and to describe the time behavior of large anti-de Sitter black holes Maldacena (2003); Papadodimas and Raju (2015); Cotler et al. .

Our analytical expression for the survival probability for the FRM model covers the entire evolution at all different time scales. Following the same steps for its derivation, we find analytical expressions for the density imbalance and the OTOC. The density imbalance is measured in experiments with cold atoms Schreiber et al. (2015); Bordia et al. . The OTOC Maldacena et al. (2016) quantifies the degree of noncommutativity in time between two Hermitian operators that commute at time $t=0$ and has been studied experimentally Gärttner et al. . Guided by the derivations with FRM, we propose expressions that match very well the numerical evolution of the realistic spin model.

The short-time dynamics of the survival probability is controlled by the Fourier transform of the envelope of the energy distribution of the initial state, the so-called local density of states (LDOS). When the perturbation that takes the system out of equilibrium is strong, the LDOS is similar to the density of states (DOS). The DOS for the FRM has a semicircle shape, which leads to a decay $\propto{\cal J}_{1}^{2}(t)/t^{2}$, where ${\cal J}_{1}(t)$ is the Bessel function of first kind Torres-Herrera and Santos (2014a); Torres-Herrera et al. (2014); Torres-Herrera and Santos (2014b); Torres-Herrera et al. (2016). The initial decay of the density imbalance follows the same behavior, whereas the OTOC goes as ${\cal J}^{4}_{1}(t)/t^{4}$. For the spin system, where only two-body interactions exist, the decay is slower. In this case, maximally spread-out LDOS reach Gaussian shapes Torres-Herrera and Santos (2014a); Torres-Herrera et al. (2014); Torres-Herrera and Santos (2014b); Torres-Herrera et al. (2016); Torres-Herrera and Santos (2014c); Izrailev and Castañeda-Mendoza (2006), resulting in Gaussian decays.

The envelope of the oscillations of the term involving the Bessel function decays as $1/t^{3}$ for the survival probability Torres-Herrera et al. (2015); Távora et al. (2016, 2017) and imbalance and as $1/t^{6}$ for the OTOC. These behaviors emerge when the tails of the DOS fall with the square root of the energy Távora et al. (2016, 2017). In the spin model, the tails of the DOS decay slowly to its energy bounds, which yields smaller power-law exponents.

For long times, but still shorter than the inverse of the mean level spacing (Heisenberg time), the survival probability for both the FRM and the spin model shows a dip below its saturation value, known as correlation hole Leviandier et al. (1986); Guhr and Weidenmüller (1990); Wilkie and Brumer (1991); Alhassid and Levine (1992); Gorin and Seligman (2002). This is an explicit dynamical manifestation of level repulsion in systems with discrete spectra Torres-Herrera and Santos (2017, ). For yet longer times, the survival probability eventually saturates. Its increase from the bottom of the hole to saturation is nearly linear. We show that the correlation hole appears also for the imbalance and the OTOC.

Hamiltonians and dynamical quantities. We consider Hamiltonians $H=H_{0}+JV$ that have an unperturbed part $H_{0}$ and a perturbation $V$ of strength $J$. We set $J=1$ and $\hbar=1$.

For the 1D spin-1/2 model with onsite disorder, $L$ sites, and periodic boundary conditions, $H_{0}=\sum_{k=1}^{L}h_{k}S_{k}^{z}$ and $V=\sum_{k=1}^{L}\vec{S}_{k}\vec{S}_{k+1}$, where $\vec{S}_{k}$’s are the spin operators on site $k$. The amplitudes $h_{k}$ of the static magnetic fields are random numbers from a uniform distribution $[-h,h]$. The total spin in the $z$ direction ${\cal S}^{z}=\sum_{k}S_{k}^{z}$ is conserved. We study the largest subspace ${\cal S}^{z}=0$, which has dimension ${\cal N}=L!/(L/2)!^{2}$.

When $h=0$ or $h>h_{c}$, where $h_{c}$ is the critical point for spatial localization, the eigenvalues can cross, and the level spacing distribution is Poissonian as typical of integrable models. For $0<h<h_{c}$, the eigenvalues become correlated and repel each other. The level spacing distribution is intermediate between the Wigner-Dyson and the Poissonian distributions. The best agreement with the Wigner-Dyson distribution for ${\cal N}=12870$ occurs at $h\sim 0.5$ Torres-Herrera and Santos (2017).

In the FRM model, $H_{0}$ is the diagonal part of the matrix, and $V$ consists of the off-diagonal elements. In the FRM from the GOE, the matrix elements $H_{nm}$ are random numbers from a Gaussian distribution with mean zero. The variance of the elements of $V$ is $\sigma^{2}$, and for $H_{0}$, it is $2\sigma^{2}$. Due to the rotational symmetry, $H_{nm}=H_{mn}=H^{*}_{mn}$ Guhr et al. (1998). As in the spin model, ${\cal N}$ is the size of the matrix.

The system is initially in one of the eigenstates $|\phi_{n}\rangle$ of $H_{0}$. The dynamics starts by switching on the perturbation abruptly. The evolution of the initial state $|\Psi(0)\rangle=|\phi_{n_{0}}\rangle$ is dictated by $H$, $|\Psi(t)\rangle=e^{-iHt}|\Psi(0)\rangle$. The eigenvalues and eigenstates of $H$ are denoted by $E_{\alpha}$ and $|\psi_{\alpha}\rangle$. The dynamical quantities investigated are listed below.

(i) The survival probability is given by

[TABLE]

where $C^{(\alpha)}_{n_{0}}=\langle\psi_{\alpha}|\Psi(0)\rangle$.

(ii) The imbalance of the spin density for all sites is computed as in Luitz et al. (2016); Lee et al. ,

[TABLE]

(iii) In terms of spin operators, the OTOC that we calculate is similar to the one in Fan et al. ,

[TABLE]

where we average over all pairs of sites $k^{\prime}>k$. In the thermal ensemble average, all states $|\phi_{n}\rangle$ of the subspace ${\cal N}$ are assumed to contribute equally.

Survival probability.– We can write Eq. (1) in terms of the Fourier transform of the spectral autocorrelation function as $W_{n_{0}}(t)=\int G(E)e^{-iEt}dE+\overline{W}_{n_{0}}$, where $G(E)=\sum_{\alpha_{1}\neq\alpha_{2}}|C^{(\alpha_{1})}_{n_{0}}|^{2}|C^{(\alpha_{2})}_{n_{0}}|^{2}\delta(E-E_{\alpha_{1}}+E_{\alpha_{2}})$ and $\overline{W}_{n_{0}}=\sum_{\alpha}|C_{n_{0}}^{(\alpha)}|^{4}$ is the infinite time average.

In the GOE FRM model, the eigenstates are random vectors, so $\langle\overline{W}_{n_{0}}\rangle_{\text{FRM}}=\overline{W}_{n_{0}}^{\text{FRM}}=3/({\cal N}+2)$, where $\langle.\rangle_{\text{FRM}}$ represents the ensemble average. Since the eigenvalues and eigenstates are statistically independent, $G(E)$ is separated into $\langle\sum_{\alpha_{1}\neq\alpha_{2}}|C^{(\alpha_{1})}_{n_{0}}|^{2}|C^{(\alpha_{2})}_{n_{0}}|^{2}\rangle_{\text{FRM}}=1-\overline{W}_{n_{0}}^{\text{FRM}}$ and $\langle\delta(E-E_{\alpha_{1}}+E_{\alpha_{2}})\rangle_{\text{FRM}}=\int\delta(E-E_{\alpha_{1}}+E_{\alpha_{2}})R_{2}(E_{\alpha_{1}},E_{\alpha_{2}})dE_{\alpha_{1}}dE_{\alpha_{2}}/[{\cal N}({\cal N}-1)]$, where $R_{2}(E_{\alpha_{1}},E_{\alpha_{2}})$ is the two-point correlation function. $R_{2}$ splits in the one-point correlation function, which is simply the DOS, and the two-level cluster function Mehta (1991). As ${\cal N}\rightarrow\infty$, the DOS converges to the Wigner semicircle law

[TABLE]

where $2\varepsilon$ is the length of the spectrum.

The Fourier transform of the semicircle leads to a term $\propto{\cal J}_{1}(\varepsilon t)/t$ Torres-Herrera and Santos (2014a). The Fourier transform of the two-level cluster function gives the two-level form factor $b_{2}(Dt/2\pi)$, where $D$ is the mean level spacing Mehta (1991); Guhr et al. (1998). In the large ${\cal N}$ limit, $D\approx 1/\rho(0)$. Therefore,

[TABLE]

where $b_{2}(\overline{t})=[1-2\overline{t}+\overline{t}\ln(1+2\overline{t})]\Theta(1-\overline{t})+\{-1+\overline{t}\ln[(2\overline{t}+1)/(2\overline{t}-1)]\}\Theta(\overline{t}-1)$ and $\Theta$ is the Heaviside step function.

In Fig. 1 (a), we compare Eq. (4) with the numerical results for the GOE FRM. The agreement is excellent; the two curves can hardly be distinguished.

The initial evolution of $W_{n_{0}}^{\text{FRM}}(t)$ is controlled by the term with the Bessel function, which leads to oscillations that decay as $1/t^{3}$, as indicated by the dashed line in Fig. 1 (a). The correlation hole, corresponding to the full time interval where $W_{n_{0}}^{\text{FRM}}(t)$ is below $\overline{W}_{n_{0}}^{\text{FRM}}$, is caused by $b_{2}(\overline{t})$. As we approach the Heisenberg time, the hole fades away, and the dynamics eventually saturates at $\overline{W}_{n_{0}}^{\text{FRM}}$.

The correlation hole is a direct probe of short- and long-range correlations between the eigenvalues. For level statistics given by the Poissonian distribution, $b_{2}(\overline{t})=0$, and the hole is nonexistent.

In Fig. 1 (b), we compare Eq. (4) (the dotted line) with numerical results for the analytic continuation of the partition function, $|Z(\beta+it)|^{2}=\sum_{\alpha}\exp[-(\beta+it)E_{\alpha}]/Z(\beta)$ (the solid lines). As discussed in Ref. del Campo et al. , $|Z(\beta+it)|^{2}$ is analogous to the survival probability if one considers as initial state, a thermofield-double state, that is $|\Psi(0)\rangle=\sum_{\alpha}\exp(-\beta E_{\alpha}/2)|\psi_{\alpha}\rangle/\sqrt{Z(\beta)}$. As illustrated in Fig. 1 (b), the results for $|Z(\beta+it)|^{2}$ for GOE FRM show qualitative agreement with $W_{n_{0}}^{\text{FRM}}(t)$. The survival probability and $|Z(\beta+it)|^{2}$ for $\beta=0$ decay initially as ${\cal J}_{1}^{2}(\varepsilon t)/(\varepsilon t)^{2}$, and all curves in Fig. 1 (b) show correlation holes. However, this comparison has limitations, since in quench dynamics $C^{(\alpha)}_{n_{0}}$ cannot be chosen independently of $H_{0}$ and $H$ as performed for the thermofield state. Contrary to $|Z(\beta+it)|^{2}$, $W_{n_{0}}(t)$ depends on the quench protocol.

Figure 1 (c) depicts the survival probability for the spin model with different disorder strengths. The curves are averages over disorder realizations and $0.1{\cal N}$ initial states with energy in the middle of the spectrum. Even deep in the chaotic regime ($h=0.5$), the decay of $\langle W_{n_{0}}(t)\rangle$ is slower than that for the FRM model, being bounded by Eq. (4). This is caused by two related factors typical of realistic systems with two-body interactions: the Gaussian shape of the DOS Brody et al. (1981) and the lack of full ergodicity of the eigenstates.

Using as a reference the steps for the analytical derivation of $G(E)$ for FRM, namely that the $R_{2}$ function splits into the DOS and the two-level cluster function, we look for an expression that can reproduce the evolution of the chaotic spin model. We take into account the following features of the realistic system: (i) The Fourier transform of a Gaussian LDOS gives a Gaussian decay at short times $e^{-w^{2}t^{2}}$, where $w$ is the width of the energy distribution Torres-Herrera and Santos (2014a); Torres-Herrera et al. (2014); Torres-Herrera and Santos (2014b, c); Izrailev and Castañeda-Mendoza (2006), (ii) this distribution is bounded in energy Távora et al. (2016, 2017) and nearly constant at the edges, which causes a power-law behavior $\propto 1/t^{2}$, and (iii) the presence of level repulsion induces the correlation hole at long times. These aspects, together with the saturation of $\langle W_{n_{0}}(t)\rangle$, motivate the expression

[TABLE]

where $g(t)=e^{-w^{2}t^{2}}+A(1-e^{-w^{2}t^{2}})/(w^{2}t^{2})$ and $A$ is a fitting constant. Apart from the first term, which depends on the shape and tails of the energy distribution, Eq.(5) is equal to Eq.(4). It is impressive that, with a single fitting constant, our expression captures so well the entire evolution of $\langle W_{n_{0}}(t)\rangle$ for $h=0.5$ as seen in Fig. 1 (c).

The inset of Fig. 1 (c) confirms that $b_{2}$ is the appropriate function to describe the correlation hole also for the chaotic spin system. The $h=0.5$ curve follows closely the FRM analytical expression. This indicates that the long-time behavior of realistic chaotic many-body systems (before saturation) depends only on the correlations in the eigenvalues, not on details of the model, such as the shape of the DOS and structure of the eigenstates.

The origin of the $1/t^{3}$ decay for the FRM model is the square-root edge of the DOS. This power-law exponent is observed also for the Sachdev-Ye-Kitaev model Cotler et al. ; Bagrets et al. where the DOS is also a semicircle at the edges García-García and Verbaarschot (2016, 2017) and for ($1+1$)-dimensional conformal field theories with a gravity dual Dyer and Gur-Ari . Since field theories with holographic duals set bounds to certain dynamical coefficients Kovtun et al. (2005), one may speculate whether the $1/t^{3}$ behavior is a general bound to the decay of the survival probability and related quantities of generic lattice many-body quantum systems. If we replace the Gaussian distribution of the random entries of the FRM by distributions involving higher even powers, it is possible to achieve DOS whose tails go as $|E-E_{0}|^{\xi}$ where $\xi=3/2,5/2,\ldots$ and $E_{0}$ is the edge of the spectrum Brézin et al. (1978), which would lead to decays faster than $1/t^{3}$. Whether there may be realistic systems with such DOS is an open question.

Density imbalance. Level repulsion manifests itself not only as the correlation hole of the survival probability. It is revealed also in the long-time evolution of experimental observables such as the spin density imbalance.

The curves for the density imbalance for the FRM model and for the disordered spin system with different values of $h$ show a dip below the saturation value as illustrated in Fig. 2 (a). As $h$ increases above $0.5$ and the realistic system moves away from the chaotic region, the hole in Fig. 2 (a) shows the same features of the hole in Fig. 1 (c). It gets less deep, its time interval shrinks, and the moment when it first appears gets deferred to longer times. This is consistent with the fact that the long-range correlations in the eigenvalues diminish as the realistic system moves towards a localized phase. The depth of the correlation hole has been used to signal the metal-insulator transition in Refs. Torres-Herrera and Santos (2017, ).

To obtain an analytical expression for the density imbalance, we refer to the equation $O(t)=\int K(E)e^{-iEt}dE+\overline{O}$ for a general observable $O$, where $K(E)=\sum_{\alpha_{1}\neq\alpha_{2}}C^{(\alpha_{1})}_{n_{0}}C^{(\alpha_{2})}_{n_{0}}O_{\alpha_{1}\alpha_{2}}\delta(E-E_{\alpha_{1}}+E_{\alpha_{2}})$ with $O_{\alpha_{1}\alpha_{2}}=\langle\psi_{\alpha_{1}}|O|\psi_{\alpha_{2}}\rangle$ and $\overline{O}=\sum_{\alpha}|C^{(\alpha)}_{n_{0}}|^{2}O_{\alpha\alpha}$ is the infinite time average. In the FRM model, where the eigenvalues, eigenstates, and $O_{\alpha_{1}\alpha_{2}}$ are statistically independent, we can separate $K(E)$ into $\langle\sum_{\alpha_{1}\neq\alpha_{2}}C^{(\alpha_{1})}_{n_{0}}C^{(\alpha_{2})}_{n_{0}}O_{\alpha_{1}\alpha_{2}}\rangle_{\text{FRM}}=O(0)-\overline{O}^{\text{FRM}}$ and $\langle\delta(E-E_{\alpha_{1}}+E_{\alpha_{2}})\rangle_{\text{FRM}}$, already computed for Eq. (4).

Using the reasoning above, we obtain the following expression for the density imbalance:

[TABLE]

where $\overline{I}^{\text{FRM}}=2I(0)/({\cal N}+2)$. The result is very similar to that for the survival probability, leading also to the $1/t^{3}$ decay of the oscillations as seen in Fig. 2 (a).

The decay of the density imbalance for the spin model is bounded by Eq. (6). It shows a power-law behavior also in the chaotic domain, which indicates that algebraic decays are not exclusive to systems in the vicinity of a localized phase.

The relaxation of $I(t)$ for the disordered spin model was investigated in Ref. Luitz et al. (2016). There, a fitting function with nine free parameters was proposed for the intermediate times where the power-law behavior is observed. We add to this picture the description of the short- and long-time dynamics.

The imbalance for the spin system follows closely what happens for the survival probability. The initial decay, up to $wt\sim 2$, is Gaussian as shown in Fig. 2 (b).

The correlation hole emerges at long times and is shown in Fig. 2 (c). The numerical curve for $h=0.5$ is fitted with the function $At^{-B}-Cb_{2}\left(\frac{wt}{{\cal N}}\right)$, where $A$, $B$, and $C$ are fitting constants. We use the same $b_{2}(\overline{t})$ used for the survival probability in Fig. 1 (c). The agreement is extremely good, covering a large time interval all the way to saturation.

Out-of-time-ordered correlator. Analogous to what happens for the density imbalance, the evolution of the OTOC for the FRM model is initially very fast and later shows oscillations that decay as $1/t^{6}$. The OTOC involves the four-point correlation function $R_{4}(E_{\alpha_{1}},E_{\alpha_{2}},E_{\alpha_{3}},E_{\alpha_{4}})$ derived from the ensemble average $\langle\delta(E-E_{\alpha_{1}}+E_{\alpha_{2}}-E_{\alpha_{3}}+E_{\alpha_{4}})\rangle_{\text{FRM}}$. $R_{4}$ can be expressed as the determinant of a single spectral kernel which is known explicitly Mehta (1991). For short and intermediate times, the leading contribution to the Fourier transform of $R_{4}$ is proportional to ${\cal J}_{1}^{4}(\varepsilon t)/(\varepsilon t)^{4}$, which causes the $1/t^{6}$ decay. At long times, $b_{2}^{2}(Dt/2\pi)$ becomes dominant and causes the correlation hole.

The $1/t^{6}$ behavior of the OTOC is shown in Fig. 3 (a). The agreement between the numerical data and the analytical prediction from the FRM is very good. In Fig. 3 (b), the analytical curve for the FRM model is compared with the decay for the disordered spin system with $h=0.5$. The decay of the latter is slower and exhibits a Gaussian behavior for short times.

The survival probability, and therefore $I(t)$ and $O_{toc}(t)$, are not self-averaging Prange (1997). The size of the ensemble of random matrices needed to reasonably expose the correlation hole for the density imbalance and the OTOC is significantly larger than for $\langle W_{n_{0}}(t)\rangle$.

Conclusion. We have found analytical expressions for the evolution of the survival probability, density imbalance, and OTOC for a FRM model. These observables are central to theoretical and experimental studies of quantum systems out of equilibrium. The analytical findings were compared with numerical results for a 1D-disordered spin-1/2 system. The power-law decays, for intermediate times, and dips below the saturation values, for longer times, revealed by the FRM model appeared also for the chaotic spin model. The identification of these generic properties helped us finding and justifying functions that describe very well the numerical evolution of the spin model at different time scales. This approach can be used also for describing equivalent realistic lattice many-body quantum systems with level repulsion.

Acknowledgments. E.J.T.-H acknowledges funding from CONACyT and VIEP-BUAP, Mexico. He is also grateful to LNS-BUAP for allowing the use of their supercomputing facility. A.M.G.-G acknowledges partial financial support from a QuantEmX grant from ICAM and the Gordon and Betty Moore Foundation through Grant No. GBMF5305. L.F.S. was supported by the NSF Grant No. DMR-1603418.