Quantum vs classical dynamics in a spin-boson system: manifestations of spectral correlations and scarring

Experimental advances in the studies of isolated quantum systems have resulted in ever longer coherence times [1–4], prompting theoretical and experimental analysis of long-time quantum dynamics. Recent works on quantum chaos have shown that the short-time exponential growth of out-of-time-ordered correlators (OTOCs) [5–11] is not a universal signature of chaos, but can emerge also near critical points [12–14]. This has offered another motivation to switch the attention to long-time dynamics [13, 15].

At large times, equilibration eventually occurs, and a main question in studies of nonequilibrium quantum dynamics is how long it takes for this to happen. In classical mechanics, the mixing properties of chaotic dynamics have provided a fundamental mechanism to explain the equilibration process as well as the ergodic properties of physical systems. In quantum mechanics, even though isolated systems are described by linear equations, one can still talk about equilibration in the sense of saturation of the dynamics, that is, the evolution of observables reaches a point where it simply fluctuates around its asymptotic value, while these fluctuations decrease with the system size [16–20]. Since the passage from the classical to the quantum domain entails new phenomena, such as superpositions, quantum interferences, the effects of universal spectral correlations [21–29], and quantum scars [30–40], one may ask what differences and similarities one should find between the classical and quantum approaches to equilibrium. As we discuss here, the quantum–classical correspondence breaks down at long times, when effects that are purely quantum come to light.

The present work investigates the entire quantum and classical evolutions of the Dicke model, with the purpose of identifying general features. This two-degree-of-freedom interacting spin-boson model was introduced to explain the collective phenomenon of superradiance [41–46], a phenomenon that has been experimentally studied with cold atoms in optical cavities [47–52]. Depending on the parameters and excitation energy, the model presents regular and chaotic domains [45, 53–58]. It has been used in studies of nonequilibrium dynamics [27, 59–63] and as a paradigm of the ultra-strong coupling regime in several systems [64–66]. Experimentally, the model can be studied by means of cavity assisted Raman transitions [50, 67] and with trapped ions [68, 69].

By comparing the entire quantum and classical evolutions, we provide a broad picture of the dynamics and find when and why the quantum–classical correspondence no longer holds. We place initial coherent states in the regular and chaotic regions of the model and study the probability to find the initial state at a later time, the so-called survival probability or return probability. We choose coherent states, because they enable a direct comparison between the exact quantum evolution and its classical description obtained with the truncated Wigner approximation (TWA) [70–72]. An advantage of using the survival probability is that it recognizes differences in the quantum and classical dynamics that are not limited to the quantum fluctuations after saturation.

Results are presented for various initial states and four representative ones are studied at length. They are selected according to their level of delocalization in the energy eigenbasis. At short times, the quantum and classical evolutions coincide. At long times, four general cases are singled out, as listed below. They are distinguished according to the proximity to a separatrix in the regular regime and according to the manifestations of quantum chaos that are observed: scarring or ergodicity.

• (a)  
  In the regular regime, if the initial state is far from a separatrix, the quantum and classical evolutions coincide up to equilibration. Beyond this point, the classical survival probability reaches a constant value, while the quantum survival probability fluctuates around the classical asymptotic value.
• (b)  
  In the regular regime, if the initial state is close to a separatrix, small differences between the quantum and classical evolutions appear at long times due to a tunneling effect that exists only in the quantum domain.
• (c)  
  In the chaotic region, the quantum survival probability equilibrates faster than the classical one, when the quantum initial state has large components in scarred eigenstates. Quantum scars are signatures in the Hilbert space of the presence of classical unstable periodic orbits.
• (d)  
  In the chaotic region, if the initial state is highly delocalized in the energy eigenbasis (ergodic), the quantum survival probability takes longer to equilibrate than its classical counterpart. This is because the quantum survival probability equilibrates only after passing through the correlation hole, which is a dynamical manifestation of spectral correlations that is nonexistent in the classical limit.

We show that the ratio between the asymptotic values of the quantum and classical survival probabilities can be used as a metric to distinguish the four cases above. A ratio equal to two indicates maximal quantum ergodicity, as in (d). For this case, we derive an analytical expression for the evolution of the quantum survival probability and its equilibration time.

We also compare the results for the survival probability with the classical evolution of the Wigner distribution in phase space. With this parallel, we gain a deeper understanding of specific features of the quantum dynamics that emerge at different time scales.

This paper is organized as follows. The core of the work is in section 5, where a detailed comparative study of the classical and quantum evolution of the survival probability is presented for the regular and chaotic regimes. In preparation for this analysis, section 2 describes the Dicke Hamiltonian and its classical limit, and section 3 describes the initial coherent states studied. In section 4, we introduce the survival probability, its relation with the local density of states (LDoS), and its classical approximation using the TWA. Our conclusions are presented in section 6.

The Dicke model [41] represents a set of $\mathcal{N}$ two-level atoms with atomic transition frequency ω₀ interacting with a single mode of a radiation field with frequency ω. It is described by the following Hamiltonian,

$$\begin{equation}{\hat{H}}_{\mathrm{D}}=\omega {\hat{a}}^{{\dagger}}\hat{a}+{\omega }_{0}{\hat{J}}_{z}+\frac{2\gamma }{\sqrt{\mathcal{N}}}{\hat{J}}_{x}\left({\hat{a}}^{{\dagger}}+\hat{a}\right),\end{equation} \tag{ 1 }$$

where ℏ = 1, â (â^†) is the bosonic annihilation (creation) operator of the field mode, ${\hat{J}}_{x,y,z}=\frac{1}{2}{\sum }_{k=1}^{\mathcal{N}}{\hat{\sigma }}_{x,y,z}^{k}$ are collective pseudo-spin operators given by the sum of the Pauli matrices ${\hat{\sigma }}_{x,y,z}$, and γ is the spin-boson interaction strength. When γ reaches a critical value ${\gamma }_{\mathrm{c}}=\sqrt{\omega {\omega }_{0}}/2$, a second-order quantum phase transition takes place in the system [42–45]. It goes from the normal phase (γ < γ_c), where the ground state is characterized by all atoms in their ground state and no photons, to the the superradiant phase (γ > γ_c), where the ground state has a macroscopic population of photons and excited atoms.

The eigenvalues j(j + 1) of the total spin operator ${\hat{\mathbf{J}}}^{2}={\hat{J}}_{x}^{2}+{\hat{J}}_{y}^{2}+{\hat{J}}_{z}^{2}$ determine the different invariant subspaces. We work with the maximum value $j=\mathcal{N}/2$, which defines a symmetric atomic subspace that includes the ground state. The Hamiltonian Ĥ_D commutes also with the parity operator $\hat{{\Pi}}={\mathrm{e}}^{\mathrm{i}\pi \hat{{\Lambda}}}$, where $\hat{{\Lambda}}={\hat{a}}^{{\dagger}}\hat{a}+{\hat{J}}_{z}+j\hat{1}$ represents the total number of excitations with eigenvalues Λ = n + m + j. Here, n indicates the number of photons and m + j is the number of excited atoms, m being the eigenvalue of the operator Ĵ_z.

2.1. Classical limit

The corresponding classical Hamiltonian is obtained using Glauber coherent states for the bosonic sector [55, 56, 58, 73–75],

$$\begin{equation}\vert q,p\rangle ={\mathrm{e}}^{-\left(j/4\right)\left({q}^{2}+{p}^{2}\right)}\;{\mathrm{e}}^{\left[\sqrt{j/2}\left(q+\mathrm{i}p\right)\right]{\hat{a}}^{{\dagger}}}\vert 0\rangle ,\end{equation} \tag{ 2 }$$

and Bloch coherent states for the pseudo-spin sector,

$$\begin{equation}\vert Q,P\rangle ={\left(1-\frac{{Z}^{2}}{4}\right)}^{j}{\mathrm{e}}^{\left[\left(Q+\mathrm{i}P\right)/\sqrt{4-{Z}^{2}}\right]{\hat{J}}_{+}}\vert j,-j\rangle ,\end{equation} \tag{ 3 }$$

where Z² = Q² + P². The canonical j-independent variables (q, p) and (Q, P) are associated with the photonic and atomic degrees of freedom, respectively. The state |0⟩ denotes the photon vacuum and |j, −j⟩, the state with all atoms in their ground state. The rescaled classical Hamiltonian h_cl is given by (see appendix A for details),

$$\begin{align}\hfill {h}_{\mathrm{c}\mathrm{l}}& \equiv \frac{\langle q,p;Q,P\vert {\hat{H}}_{\mathrm{D}}\vert q,p;Q,P\rangle }{j}\hfill \\ \hfill & =\frac{\omega }{2}\left({q}^{2}+{p}^{2}\right)+\frac{{\omega }_{0}}{2}{Z}^{2}+2\gamma Qq\sqrt{1-\frac{{Z}^{2}}{4}}-{\omega }_{0}.\hfill \end{align} \tag{ 4 }$$

The rescaled classical Hamiltonian and its four-dimensional phase space $\mathcal{M}$ are independent of j. This is equivalent to working with an effective Planck constant ℏ_eff = 1/j [76].

2.2. Hamiltonian parameters

The Hamiltonian parameters ω, ω₀, and γ are chosen, so that the regular and chaotic regimes are clearly identified. The parameter γ controls the coupling between photons and atoms in the system and therefore the emergence of a chaotic region. We select the coupling strength in the superradiant phase, γ = 2γ_c, where the presence of chaotic behavior is guaranteed [58]. We choose ω = ω₀ for convenience, but other choices would not change the broad picture developed in this work.

For the system size j = 100 considered here, the normalized energy

$$\begin{equation}{\epsilon}=\frac{E}{{\omega }_{0}j}\end{equation} \tag{ 5 }$$

of the ground state is _GS = −2.125. The dynamics is regular from _GS up to ≈ −1.7 and chaotic above this point [58].

2.3. Numerical diagonalization

The effective Planck's constant 1/j determines the resolution that a coherent initial state has in the classical phase space. This means that a large j is necessary for the quantum dynamics to reflect the classical effects and for the purely quantum properties to be identified. Large systems are needed to resolve the structure of the classical phase space [77]. We are able to consider a large system size, j = 100, because we employ an efficient basis that guarantees the convergence of the eigenvalues and eigenstates for a broad part of the bounded spectrum of the Dicke model [78, 79]. The dimension of the truncated Hilbert space in this efficient basis is given by dim = (2j + 1)(N_max + 1), where N_max is an upper bound to the modified bosonic subspace. For our case study, N_max = 300, so the dimension of the truncated Hilbert space is dim = 60 501. This ensures

$$\begin{equation}{N}_{\mathrm{c}}=30\;825\end{equation} \tag{ 6 }$$

converged eigenstates and eigenenergies, which range from the ground state energy _GS = −2.125 up to a truncation energy _T = 0.853.

The selected initial states in this work are the Glauber–Bloch coherent states |q, p; Q, P⟩. They allow for a direct connection between the quantum states and the coordinates (q, p, Q, P) in the classical phase space of h_cl, and also for a relatively simple calculation of the Wigner distribution needed for the evaluation of the classical dynamics. To perform the numerical quantum calculations, the coherent states |q, p; Q, P⟩ are expanded in the energy eigenstates through the efficient basis [57].

To choose the initial states, we restrict ourselves to the hyperplane p = 0 and solve the second-degree equation h_cl(q, p, Q, P) = ω₀ in q. This equation has two solutions, q₋ and q₊, with q₋ ⩽ q₊. The four representative initial states that we use are centered at (q, p, Q, P) = (q₊, 0, Q₀, P₀). Two of them have energy in the regular region, _R = −1.8, and the other two have the energy shell fully covered by chaotic trajectories, _C = −0.5. The criteria for our choices and the specific values of the coordinates (Q₀, P₀) are described below. In section 5, we provide a detailed analysis of the classical and quantum evolution of these four representative initial states. The results are then confirmed for several other initial states in the regular (see figure 6) and chaotic region (see figure 9).

3.1. Regular regime: _R = −1.8

Poincaré sections at p = 0 of the classical Hamiltonian at energy _R = −1.8 are shown in figure 1(a). The phase space splits in three regions of regular trajectories, whose borders are defined by the separatrix marked with large black dots. We choose state I at (Q₀, P₀) = (1, 0) and indicate it with a blue dot in figure 1(a). State II is close to the separatrix, (Q₀, P₀) = (1.2, 0), and is shown with a red dot in the same figure.

Zoom In Zoom Out Reset image size

The structure of the phase space reflects the quasi-conserved quantities of the Dicke Hamiltonian at low energies [80]. They can be identified by means of an adiabatic approximation, where the dynamics separates into two parts, one fast-evolving mode and a slow one, effectively decoupling the boson and pseudo-spin dynamics [81, 82]. For the parameters considered here, a quasi-constant of motion is given by the nutation angle of the pseudo-spin, which precesses fast around an axis whose direction oscillates slowly in a way dictated by the slow bosonic variables. State I is at the center of the slow-boson regular region located between 0.85 ⩽ Q ⩽ 1.16, where the boson modes are expressed by very small nutation angles and large amplitudes of the precession axis' oscillations. The rightmost region in figure 1(a) corresponds to the mode of the fast pseudo-spin degree of freedom, where the dynamics has maximal nutation angles and very small oscillations of the precession axis instead.

Quasi-integral behaviors are destroyed in non-linear systems due to the growth and consequent overlap of nonlinear resonances between the integrable modes. These non-linear resonances arise locally in the phase space [83] and are separated by a separatrix, around which a stochastic layer is formed [84]. Chaos emerges from it, destroying the remnant regular surfaces and opening the phase space to the diffusion of trajectories, which start wandering through the whole region of the resonance overlap. In the Poincaré sections of figure 1(a), one can see non-linear resonances between the adiabatic modes described above, specifically in the moon-shaped region of regular trajectories rotating around the point (Q, P) ≃ (1.3, 0). State II is in this region, but very close to the separatrix.

In figure 1(b), we plot the participation ratio P_R of the coherent states centered in each point of the Poincaré surface and projected into the energy eigenstates,

$$\begin{equation}{P}_{\mathrm{R}}=\frac{1}{{\sum }_{k}\vert {c}_{k}{\vert }^{4}},\end{equation} \tag{ 7 }$$

where c_k = ⟨ E_k|q₊, 0; Q₀, P₀⟩ and Ĥ_D|E_k⟩ = E_k|E_k⟩. This quantity measures the level of delocalization of the initial state in the energy eigenbasis. It varies from the minimum value, P_R = 1, when the initial state coincides with an eigenstate and is therefore maximally localized, to the largest value, P_R = N_c, when all eigenstates participate equally in the evolution of the initial state. In the case of a random vector, when the components c_k are uncorrelated random numbers from a Gaussian distribution, P_R = N_c/3.

The participation ratio is strongly correlated with the underlying classical dynamics [57, 85]. For instance, the stochastic layer appearing in the classical dynamics around the separatrix in figure 1(a) is associated with large values of P_R in figure 1(b). The two initial states that we chose in the regular regime have the lowest and largest P_R for that energy surface. These extremal cases allow us to identify the fundamental mechanisms underlying the evolution of initial coherent states.

3.2. Chaotic regime: _C = −0.5

The survival probability, S_P, is the probability of finding an evolved quantum state back in the initial state |Ψ(0)⟩ = ∑_kc_k|E_k⟩,

$$\begin{equation}{S}_{\mathrm{P}}\left(t\right)=\vert \langle {\Psi}\left(0\right)\vert {\Psi}\left(t\right)\rangle {\vert }^{2}={\left\vert \sum _{k}\vert {c}_{k}{\vert }^{2}\;{\mathrm{e}}^{-\mathrm{i}{E}_{k}t}\right\vert }^{2}.\end{equation} \tag{ 8 }$$

By introducing the local density of states (LDoS) or strength function, that is the energy distribution weighted by the components |c_k|² of the initial state,

$$\begin{equation}\mathcal{G}\left(E\right)=\sum _{k}\vert {c}_{k}{\vert }^{2}\delta \left(E-{E}_{k}\right),\end{equation} \tag{ 9 }$$

we can also write the survival probability as the squared norm of the Fourier transform of $\mathcal{G}\left(E\right)$

$$\begin{equation}{S}_{\mathrm{P}}\left(t\right)={\left\vert \int \mathrm{d}E\;\mathcal{G}\left(E\right){\mathrm{e}}^{-\mathrm{i}Et}\right\vert }^{2}.\end{equation} \tag{ 10 }$$

The evolution of the survival probability shows different behaviors at different time scales [25, 26, 86, 87]. By smoothing the LDoS, one gets insight on how its structure affects the dynamics at different times. The smoothing is done through a finite resolution function, given by

$$\begin{equation}{\rho }_{\mathrm{T}}\left(E\right)=\sum _{k}\vert {c}_{k}{\vert }^{2}{{\Upsilon}}_{\mathrm{T}}\left(E-{E}_{k}\right),\end{equation} \tag{ 11 }$$

where ϒ_T(E − E_k) = (T/π)sinc[(E − E_k)T], and sinc(x) = sin(x)/x. The time resolution T = π/Δ reflects aspects of the LDoS that are of order Δ in energy.

4.1. Time scales of the survival probability

Figures 2(a) and (b) and 3(a), (b), (d), (e), (g), (h) show the smoothed LDoS, ρ_T(E), of the initial coherent states described in section 3 for different time resolutions (equation (11)). In figures 2(c) and (d) and 3(j) and (k), we show the infinite-time resolution LDoS, $\mathcal{G}\left(E\right)$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1.1. Survival probability: initial decay

For all cases, up to a time resolution T* that depends on the state, we find a very good Gaussian distribution for the coherent states profiles,

$$\begin{equation}\rho \left(E\right)=\frac{1}{\sqrt{2\pi {\sigma }^{2}}}\;{\mathrm{e}}^{-{\left(E-{E}_{\mathrm{c}}\right)}^{2}/\left(2{\sigma }^{2}\right)}.\end{equation} \tag{ 12 }$$

The distributions are centered at the energy E_c of the initial state and have width given by the energy standard deviation σ, which can be calculated numerically or even analytically [88–90].

According to equation (10), the Gaussian envelope of the LDoS leads to an initial Gaussian decay of the survival probability,

$$\begin{equation}{S}_{\mathrm{P}}\left(t\right)={\mathrm{e}}^{-{\sigma }^{2}{t}^{2}},\end{equation} \tag{ 13 }$$

which is consistent with the universal quadratic behavior, S_P(t ≪ σ⁻¹) ≈ 1 − σ²t², for very short times.

In addition to the initial coherent states, we show in figures 3(c), (f), (i), and (l) the smoothed and infinite-time resolution LDoS of a random initial state. Motivated by the high level of delocalization of the initial coherent state IV in the chaotic regime, the random state is built with energy components generated randomly around the Gaussian distribution of state IV as

$$\begin{equation}\vert {c}_{k}^{\left(r\right)}{\vert }^{2}=\frac{{r}_{k}\rho \left({E}_{k}\right)}{\mathcal{A}\nu \left({E}_{k}\right)}.\end{equation} \tag{ 14 }$$

Above, ρ(E_k) is obtained from equation (12), ν(E_k) is given by the equation for the density of states (DoS) provided in the appendix A, r_k are random numbers from an exponential distribution P(r) = λ e^−λr, and $\mathcal{A}={\sum }_{q}{r}_{q}\rho \left({E}_{q}\right)/\nu \left({E}_{q}\right)$ is a normalization constant. The division by the DoS is done to compensate for the different energy densities and to guarantee a smooth enveloping distribution ρ(E). The exponential distribution for generating the random numbers is used because if we evaluate the numbers

$$\begin{equation}{r}_{k}^{\text{CS}}=\frac{\mathcal{A}\nu \left({E}_{k}\right)\vert {c}_{k}{\vert }^{2}}{\rho \left({E}_{k}\right)},\end{equation} \tag{ 15 }$$

for the components c_k of the coherent state IV in the energy eigenbasis, we obtain the histogram in figure 3(m), which is very well fitted with an exponential distribution.

By comparing the smoothed LDoS of the state III in figure 3(d), the state IV in figure 3(e), and the random state in figure 3(f), all of them with time resolution ω₀T = 2, one notices that the smoothed LDoS of the random state already deviates from a Gaussian, which must affect its short-time dynamics. While the initial coherent states III and IV are expected to remain on the Gaussian decay at this time scale, the random state should already diverge from it, as we indeed confirm in section 5. For an even higher time resolution, as ω₀T = 3 in figures 3(g)–(i), the smoothed LDoS of state III also deviates from the Gaussian, and so does its Gaussian decay (see section 5).

4.1.2. Survival probability: intermediate times

The behavior of the survival probability at intermediate times for the different initial coherent states can be anticipated from the infinite-time resolution LDoS in figures 2(c) and (d) and 3(j) and (k). For the regular coherent state I (figure 2(c)), very few components have no negligible values and the LDoS is well described by a Gaussian distribution. For the regular coherent state II located close to the separatrix (figure 2(d)), the number of participating energy levels is larger, but they are still organized according to a set of different Gaussians with different amplitudes and centers [90]. In contrast, the components of the chaotic coherent states in figure 3 have a very different structure. For the coherent state III, with small P_R (figure 3(j)), the components are bunched around some specific energy levels, while the coherent state IV, with a large P_R (figure 3(k)), counts with the participation of most components. The consequences of these distributions to the evolution of the survival probability are discussed in section 5.

4.1.3. Survival probability: asymptotic values

The asymptotic value of the survival probability,

$$\begin{equation}{S}_{\mathrm{P}}^{\infty }=\langle {S}_{\mathrm{P}}{\left(t\right)\rangle }_{t\to \infty }=\underset{t\to \infty }{\mathrm{lim}}\frac{1}{t}{\int }_{0}^{t}\mathrm{d}{t}^{\prime }\;{S}_{\mathrm{P}}\left({t}^{\prime }\right),\end{equation} \tag{ 16 }$$

can be derived from

$$\begin{equation}{S}_{\mathrm{P}}\left(t\right)=\sum _{k\ne l}\vert {c}_{l}{\vert }^{2}\vert {c}_{k}{\vert }^{2}{\mathrm{e}}^{-\mathrm{i}\left({E}_{k}-{E}_{l}\right)t}+\sum _{k}\vert {c}_{k}{\vert }^{4}.\end{equation} \tag{ 17 }$$

In the absence of energy degeneracies, the first term on the right-hand side of the equation cancels out on average, so

$$\begin{equation}{S}_{\mathrm{P}}^{\infty }=\sum _{k}\vert {c}_{k}{\vert }^{4},\end{equation} \tag{ 18 }$$

which is the inverse of the participation ratio shown in equation (7). However, when the energy levels have degeneracies of degree d_k (that is, |E_k, m⟩ with m = 1, ‥, d_k), the first term in equation (17) contributes with additional terms to the asymptotic value, which is now given by

$$\begin{equation}{S}_{\mathrm{P}}^{\infty }=\sum _{{E}_{k}}{\left(\sum _{m=1}^{{d}_{k}}\vert {c}_{k,m}{\vert }^{2}\right)}^{2},\end{equation} \tag{ 19 }$$

where c_k,m are the components in the degenerate space of E_k. In the superradiant phase of the Dicke model, for the energy range going from the ground state to = −1, the spontaneous breaking of the parity symmetry [91] produces an energy spectrum with two-fold degeneracies⁶. For the initial coherent states in the regular region, these degeneracies affect their asymptotic values, while for the chaotic initial states, where the energy _C is well above = −1, the influence of the degeneracies is rather marginal.

4.2. Classical limit of the survival probability

To find the classical limit of the survival probability, we use the Wigner formalism [92]. Details are given in appendix B and appendix C. The basic idea is to use the overlap property between two arbitrary quantum states |A⟩ and |B⟩ to write the survival probability in terms of the Wigner function W [93],

$$\begin{equation}{\left\vert \left\langle A\vert B\right\rangle \right\vert }^{2}={\left(2\pi \hslash \right)}^{d}\int \mathrm{d}\boldsymbol{u}\;{W}_{A}\left(\boldsymbol{u}\right){W}_{B}\left(\boldsymbol{u}\right),\end{equation} \tag{ 20 }$$

where d represents the degrees of freedom of the system. We scale the Wigner function of our Glauber–Bloch initial coherent states and work with $w\left(\boldsymbol{u},t\right)=\left(1/j\right)W\left(\sqrt{j}q,\sqrt{j}p,Q,P,t\right)$, where u = (q, p, Q, P) is a point in the j-scaled phase space $\mathcal{M}$. This gives

$$\begin{equation}{S}_{\mathrm{P}}\left(t\right)={\left(\frac{2\pi }{j}\right)}^{2}{\int }_{\mathcal{M}}\mathrm{d}\boldsymbol{u}\;w\left(\boldsymbol{u},0\right)w\left(\boldsymbol{u},t\right).\end{equation} \tag{ 21 }$$

To finally obtain the classical limit, we use the TWA. As shown in appendix C, the short-time dependence of the Wigner function can be written in terms of the Hamiltonian flow ${\boldsymbol{\varphi }}^{t}:\mathcal{M}\to \mathcal{M}$. One finds that w(u, t) = w(φ^−t(u), 0) for short times, so the classical survival probability can be defined as

$$\begin{equation}{\mathfrak{S}}_{\mathrm{P}}\left(t\right)={\left(\frac{2\pi }{j}\right)}^{2}{\int }_{\mathcal{M}}\mathrm{d}\boldsymbol{u}\;w\left(\boldsymbol{u}\right)w\left({\boldsymbol{\varphi }}^{-t}\left(\boldsymbol{u}\right)\right),\end{equation} \tag{ 22 }$$

where w(u) = w(u, t = 0). This quantity is numerically constructed through a Monte Carlo method (see appendix C). Its asymptotic value, ${\mathfrak{S}}_{\mathrm{P}}^{\infty }={\left\langle {\mathfrak{S}}_{\mathrm{P}}\left(t\right)\right\rangle }_{t\to \infty }$, is obtained as in equation (16).

We compare the entire evolution of the quantum (equation (8)) and classical (equation (22)) survival probabilities for the four initial coherent states selected in section 3. This analysis allows us to identify features of the quantum evolution that reflects the classical dynamics and properties that are purely quantum. We find that:

• The S_P(t) for coherent initial states reaches values close to zero at short times. This behavior can be understood from the short-time evolution of the Wigner distribution in phase space, and is thus a classical effect.
• The long-time behavior of initial states with the same energy in the regular or chaotic region may differ according to the states' level of delocalization in the energy eigenbasis.
• Purely quantum properties of the survival probability include not only the quantum fluctuations of S_P(t) after the saturation of the dynamics, but also tunneling effects, the dynamical consequences of the phenomenon of quantum scarring, and the dynamical manifestation of spectral correlations in the form of the correlation hole.
• The saturation times of the quantum and classical dynamics in the chaotic regime do not coincide. Quantum equilibration is shorter, when the quantum initial state has large projections in scarred eigenstates, reflecting the proximity to unstable periodic orbits in phase space. Quantum equilibration is longer, when the initial state is highly delocalized and S_P(t) develops the correlation hole before saturation.
• The value of the ratio between the asymptotic values of the quantum and classical survival probabilities indicates whether the initial state is close to a separatrix and whether it is close to an unstable periodic orbit. A ratio equal to two indicates maximal quantum ergodicity.

5.1. Regular region

We analyze first the quantum S_P(t) and the classical ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$ for the initial coherent states I and II (figure 1(a)), which are in the regular regime (_R = −1.8).

5.1.1. Initial state I: center of the regular mode

An analytical equation is available for the survival probability of the initial coherent state I, located at the center of the slow-boson regular region [90]. It is given by

$$\begin{equation}{S}_{\mathrm{P}}\left(t\right)\approx \frac{{\omega }_{1}}{\sigma \sqrt{\pi }}\sum _{n=1}{\mathrm{e}}^{-{n}^{2}\left(\frac{{\omega }_{1}^{2}}{4{\sigma }^{2}}+\frac{{t}^{2}}{{t}_{\mathrm{D}}^{2}}\right)}\mathrm{cos}\left(n{\omega }_{1}t\right)+{S}_{\mathrm{P}}^{\infty },\end{equation} \tag{ 23 }$$

where the index n denotes the distance between the eigenenergies of the levels with no negligible components, ω₁ is related with the spacing between neighboring levels and is equal to the average of the classical frequency over the initial Wigner distribution, σ is the width of the LDoS, t_D = ω₁/(σ|e₂|) is the decay time, with e₂ being the anharmonicity of the participating spectrum, and

$$\begin{equation}{S}_{\mathrm{P}}^{\infty }=\frac{{\omega }_{1}}{2\sigma \sqrt{\pi }}.\end{equation} \tag{ 24 }$$

The time t_D comes from the Gaussian decay of the slowest component n = 1.

In the bottom panel of figure 4, we compare the full quantum evolution of S_P(t) with the analytical expression in equation (23) and the classical ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$. The latter two show perfect agreement and they are indistinguishable from the quantum result up to the decay time t_D. Soon after this point, S_P(t) shows quantum fluctuations around ${S}_{\mathrm{P}}^{\infty }$ that emerge due to the discreteness of the spectrum and which the analytical expression in equation (23) and ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$ are unable to reproduce.

Zoom In Zoom Out Reset image size

The remarkable agreement between the quantum and classical results allows for a more intuitive classical interpretation of the behavior of the survival probability. In the top panels of figure 4, we show the classical evolution of the Wigner distribution in the phase space. The top row gives the distribution projected onto the Q–P plane and the bottom row, the distribution on the q–p plane. Each column corresponds to a specific instant of time, as indicated. The classical survival probability measures the percentage of the Wigner distribution (in color in the figure) that is inside of the starting region occupied at Ω₀ t = 0 (see appendix D for a detailed explanation). This starting region is marked with a black outline which is inside a bigger one marking the available phase space. The parallel between the evolution of the survival probability and of the Wigner distribution in the phase space goes as follows (animations are available in the supplemental material (SM) (https://stacks.iop.org/NJP/22/063036/mmedia)):

• At ω₀t = 0 (first column), the distribution is entirely within the starting region. Concurrently, the value of the survival probability in the bottom panel is one.
• As the distribution moves out of the starting region, the survival probability follows a Gaussian decay. At ω₀t = 2 (second column), the distribution is effectively outside the initial region and the survival probability becomes close to zero⁷.
• After a full period given by t = 2π/ω₁, the distribution comes back to the starting region at ω₀t = 6.5 (third column) and we have the first revival of S_P(t).
• After several periods, the distribution further spreads over the classical orbit and the amplitude of the revivals of the survival probability decreases, as for example at ω₀t = 100 (fourth column).
• Classically, at ω₀t = 5000 (fifth column), the Wigner distribution ends up filling homogeneously the region covered by the trajectories of the phase-space points of the initial distribution. At this point the classical ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$ reaches the constant value ${\mathfrak{S}}_{\mathrm{P}}^{\infty }$, while S_P(t) fluctuates around the asymptotic value ${S}_{\mathrm{P}}^{\infty }\approx {\mathfrak{S}}_{\mathrm{P}}^{\infty }$.

5.1.2. Initial state II: close to the separatrix

The LDoS of the initial coherent state II presents different Gaussian distributions with different amplitudes, centers and widths (figure 2(d)). This indicates that this state activates the two adiabatic modes (both the slow boson mode and the fast pseudo-spin mode), as well as the non-linear resonances between them. This is indeed expected, since state II is very close to the separatrix (figure 1(a)). An analytical expression can also be obtained for this case by generalizing equation (23), where in addition to contributions from Gaussians, we also need to take into account interferences between them [90].

In the bottom panel of figure 5, we show the quantum and classical survival probability for the state II. The agreement between the two is good, but contrary to figure 4, the values of S_P(t) are slightly smaller than those for ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$. This is better seen with the temporal averages of the curves, shown in figure 5 in blue for S_P(t) and in red for ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$.

Zoom In Zoom Out Reset image size

To better understand the small differences between S_P(t) and ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$, we calculate their asymptotic values for initial coherent states (q, p, Q, P) = (q₊, 0, Q_i, 0), where Q_i covers all possible points of the selected Poincaré surface in figure 1(a). We plot the values of ${S}_{\mathrm{P}}^{\infty }$ and ${\mathfrak{S}}_{\mathrm{P}}^{\infty }$ in figure 6. One sees that ${S}_{\mathrm{P}}^{\infty }$ is smaller than ${\mathfrak{S}}_{\mathrm{P}}^{\infty }$ near the separatrix regions (Q = 0.845, 1.182, and 1.428), while the two values get closer near the centers of the adiabatic-modes regions (Q = 1.0, 1.5) and of the center of the non-linear resonances (Q = 1.31). Since a coherent state located in the separatrix has a Wigner distribution defined in regions classically not connected by the trajectories, we attribute the small difference between S_P(t) and ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$ to a dynamic tunneling effect that takes place in the quantum regime [94], but is absent in the classical limit.

Zoom In Zoom Out Reset image size

Even though the quantum–classical agreement for state II is not exact, the comparison between the results for the survival probability in the bottom panel and the evolution of the Wigner distribution in the top panels is similar to that presented in figure 4 and explains specific features of S_P(t), such as the revival at ω₀t = 6.5. At very long times, such as ω₀t = 5000 in the figure, the Wigner distribution fills homogeneously the classical region that is covered by the trajectories of the initial distribution. Since this region is larger than in the case of state I, the asymptotic value in figure 5 is one order of magnitude smaller than in figure 4.

5.2. Chaotic region

We now study the quantum S_P(t) and the classical ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$ for the initial coherent states III and IV (figure 1(d)), which are in the chaotic regime (_C = −0.5). We start the analysis with state III, which has a smaller number of contributing energy eigenstates (smaller P_R) than state IV.

5.2.1. Initial state III: low P_R

Despite being in the chaotic region, initial coherent states such as state III count with a relatively small number of contributing energy eigenstates, as indicated by the low value of P_R in figure 1(d). As we show below, these states lead to large recurrences at intermediate times and, contrary to the regular case, the equilibration time for the classical and quantum dynamics no longer coincide. The equilibration time for the classical dynamics is now longer than the equilibration time for the quantum S_P(t).

The quantum and classical survival probabilities of the initial coherent state III are shown in the bottom panel of figure 7. There is excellent agreement from ω₀t = 0 up to times beyond the decaying recurrences. Initially both curves decay on a Gaussian to values close to zero and then revivals appear. These two features can be well understood by studying the classical evolution of the Wigner distribution shown in the top panels of figure 7. Analogously to the discussions about figures 4 and 5, as the originally localized Wigner distribution at ω₀t = 0 (first column) moves outside the starting region at ω₀t = 2 (second column), the survival probability becomes effectively zero. The recurrences are connected with the return of the distribution to the starting region, such as at ω₀t = 5.6 (third column). These recurrences are periodic, as confirmed with the inset in the bottom panel of figure 7.

Zoom In Zoom Out Reset image size

Soon after the last revival, the quantum–classical correspondence breaks down. At this point the quantum S_P(t) reaches an asymptotic value, around which one finds large quantum fluctuations, while the classical ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$ continues decreasing and equilibrates at a longer time. The asymptotic value of the classical survival probability is obtained using the TWA and the ergodic hypothesis (see appendix E),

$$\begin{equation}{\mathfrak{S}}_{\mathrm{P}}^{\infty }=\frac{1}{2\sqrt{\pi }\sigma {\nu }_{\mathrm{c}}}.\end{equation} \tag{ 25 }$$

The asymptotic result ${S}_{\mathrm{P}}^{\infty }$ of the quantum survival probability is more than twice this value.

Large ${S}_{\mathrm{P}}^{\infty }$, or equivalently small P_R, in the chaotic region is usually associated with the phenomenon of scarring, which is indeed the case here. The classical Wigner distribution at ω₀t = 400 (fifth column in the top panel of figure 7) is enhanced in a small closed region indicated with black arrows in the figure. This reveals the presence of unstable periodic orbits of relatively short period. These orbits are responsible for the short-time periodic revivals of the quantum S_P(t) and classical ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$.

For longer times, the periodic orbits produce opposing effects in the quantum and classical regimes. Classically, we observe a slow decay of ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$ toward its asymptotic value and thus a long classical equilibration time if compared to the quantum case. Recurrences imply that the dynamics revisits part of the phase space that was initially covered instead of exploring new regions, which slows down the full spread over the phase space [33, 34]. In the quantum evolution, the scarring decreases the equilibration time. As discussed in reference [32] and as can be seen in figure 3(g), the proximity of an initial state to unstable periodic orbits of short period produces a particular structure in the distribution of the energy eigenbasis components. Many of these components are very small, while the large ones are organized in bunches separated in energy by ΔE_s ≈ 2π/τ, where τ is the period of the classical unstable periodic orbit. These highly populated eigenstates are scarred, meaning that they are concentrated in the phase space around the unstable periodic orbits [33, 34]. Because of this concentration, the phase space available for the quantum evolution of the initial coherent state is effectively shrunk, resulting in an smaller equilibration time for S_P(t) than for ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$ and an asymptotic quantum value larger than the classical one, ${S}_{\mathrm{P}}^{\infty }{ >}2{\mathfrak{S}}_{\mathrm{P}}^{\infty }$.

5.2.2. Initial state IV: high P_R

The behavior of the survival probability presented in the previous subsection is not general. Most initial coherent states in the chaotic region are similar to state IV, being highly delocalized in the energy eigenbasis. In fact, as discussed in section 4.1, the energy distribution of state IV is comparable to that of a random state. This latter state is very useful, because one can derive an analytical expression for its survival probability. The expression for the average over an ensemble of initial random states is given by (see [27] and appendix F for details)

$$\begin{equation}{S}_{\mathrm{P}}^{\left(r\right)}\left(t\right)=\frac{1-{S}_{\mathrm{P}}^{\left(r\right),\infty }}{\eta -1}\left[\eta \;{\mathrm{e}}^{-{\sigma }^{2}{t}^{2}}-{b}_{2}\left(\frac{Dt}{2\pi }\right)\right]+{S}_{\mathrm{P}}^{\left(r\right),\infty },\end{equation} \tag{ 26 }$$

where $\eta =2\sqrt{\pi }\sigma {\nu }_{\mathrm{c}}$ and ν_c = ν(E_c) is the DoS (see equation (A.4) in appendix A) evaluated in the center of the energy profile E_c, σ is the width of the LDoS, the b₂ function is the Gaussian orthogonal ensemble (GOE) two-level form factor studied in random matrix theory [95], D = 2/ν_c is the mean level spacing of the correlated eigenvalues, and ${S}_{\mathrm{P}}^{\left(r\right),\infty }=2/\eta$. By comparing the quantum survival probability with ${S}_{\mathrm{P}}^{\left(r\right)}\left(t\right)$ and the classical ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$, we can explain the different behaviors of S_P(t) found at different time scales.

In the bottom panel of figure 8, we show the quantum S_P(t) of state IV (light gray solid line), its time average (blue solid line), the classical ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$ (dark gray solid line), its time average (red solid line), and the time average of the analytical expression ${S}_{\mathrm{P}}^{\left(r\right)}\left(t\right)$ (orange solid line). At short times, S_P(t) and ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$ overlap. There is also perfect agreement between the Gaussian decay of the temporal average of ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$, S_P(t) and ${S}_{\mathrm{P}}^{\left(r\right)}\left(t\right)$, as expected from the great similarity between the smoothed LDoS at low time resolution (ω₀T = 0.2) of the coherent and the random state (figures 3(b) and (c)). However, the curves for S_P(t) and ${S}_{\mathrm{P}}^{\left(r\right)}\left(t\right)$ diverge after ω₀t ∼ 0.2. One sees that S_P(t) remains on the Gaussian decay, while ${S}_{\mathrm{P}}^{\left(r\right)}\left(t\right)$ reaches a plateau. This divergence can also be understood from the smoothed LDoS, but now at higher time resolution (ω₀T = 2 and 3) (see figures 3(e), (f), (h) and (i)): while the smoothed LDoS for the state IV is still a good Gaussian, that of the random state shows deviations. This discrepancy implies that the components of the coherent state IV are not exactly random.

Zoom In Zoom Out Reset image size

The quantum and classical survival probabilities continue their Gaussian decay for ω₀t > 0.2 and reach values close to zero, which can once again be understood from the analysis of the classical evolution of the Wigner distribution shown in the top panels of figure 8. As the initially localized Wigner distribution at ω₀t = 0 (first column) moves out of the starting region at ω₀t = 2 and ω₀t = 8.3 (second and third columns), the survival probability becomes effectively zero. Beyond these times and contrary to what one sees for the regular regime and for state III, the quantum–classical correspondence breaks down before the quantum equilibration, as we explain next.

At ω₀t ≈ 30 (marked with a vertical line in the large bottom panel of figure 8), the classical survival probability attains its asymptotic equilibration value ${\mathfrak{S}}_{\mathrm{P}}^{\infty }$ and the Wigner distribution covers the entire energy shell, as seen in the top panels of figure 8 for ω₀t = 31.4 and ω₀t = 400 (fourth and fifth columns). In contrast, the quantum S_P(t) continues to raise, leading to a longer equilibration time. For ω₀t > 30, the temporal averaged S_P(t) coincides again extremely well with the analytical expression for ${S}_{\mathrm{P}}^{\left(r\right)}\left(t\right)$. During this ramp toward equilibration, the dynamics is controlled by the b₂ function, which reflects the correlations between the energy levels in the chaotic regime. For the chaotic Dicke model, these correlations are equivalent to those in the spectrum of GOE random matrices [45, 54, 75, 96], which justifies the use in equation (26) of the same form of the b₂ function from random matrix theory. At this time scale the dynamics becomes universal. The region of the ramp is referred to in figure 8 as 'GOE correlations' and is commonly known as correlation hole [21–27, 97, 98]. This is a quantum effect associated with the discreteness of the spectrum, which does not appear for the classical ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$.

In contrast to the scarred state III, the quantum survival probability of the state IV reaches its asymptotic value at a time t_r that is longer than the classical relaxation time, and is given by the same expression obtained with the analytical equation for ${S}_{\mathrm{P}}^{\left(r\right)}\left(t\right)$ (see reference [27] and appendix F),

$$\begin{equation}{t}_{\mathrm{r}}=\frac{\pi {\nu }_{\mathrm{c}}}{2\sqrt{6\;{\delta }_{{S}_{\mathrm{P}}}}},\end{equation} \tag{ 27 }$$

where ${\delta }_{{S}_{\mathrm{P}}}$ is a small parameter indicating that the values of S_P(t > t_r) are already within the quantum fluctuations around ${S}_{\mathrm{P}}^{\infty }$. This time is marked with a vertical line in figure 8.

The asymptotic value of the quantum survival probability agrees with the saturation point of ${S}_{\mathrm{P}}^{\left(r\right)}\left(t\right)$ (see appendix F),

$$\begin{equation}{S}_{\mathrm{P}}^{\infty }={S}_{\mathrm{P}}^{\left(r\right),\infty }=\frac{\langle {r}^{2}\rangle }{\langle {r\rangle }^{2}}\frac{1}{2\sqrt{\pi }\sigma {\nu }_{\mathrm{c}}}=\frac{1}{\sqrt{\pi }\sigma {\nu }_{\mathrm{c}}},\end{equation} \tag{ 28 }$$

where ⟨r⟩ and ⟨r²⟩ are the first and second moments of the exponential distribution used to generate the random numbers in equation (14), and their ratio

$$\begin{equation}\frac{\langle {r}^{2}\rangle }{\langle {r\rangle }^{2}}=2\end{equation} \tag{ 29 }$$

is determined by the fluctuations of the energy components of the LDoS of the random state with respect to the exact Gaussian envelope for state IV. Interestingly, one sees that due to this ratio, the value of ${S}_{\mathrm{P}}^{\infty }$ is twice as large as the asymptotic value of the classical survival probability, which is given by equation (25). If the fluctuations of the energy components of the LDoS of the random state were absent and ⟨r²⟩/⟨r⟩² was one, ${S}_{\mathrm{P}}^{\infty }$ would coincide with the classical result. The origin of these fluctuations is rooted in a remaining structure present even in strongly chaotic eigenstates (the so-called nodal structure [32, 34, 99]), which effectively limits the ergodicity of the quantum evolution [100] when compared to the classical dynamics.

We close this discussion with an interesting observation. The minimum value of ${S}_{\mathrm{P}}^{\left(r\right)}\left(t\right)$, reached right after the Gaussian decay, coincides exactly with ${\mathfrak{S}}_{\mathrm{P}}^{\infty }$. That is, before the ramp toward equilibration, ${S}_{\mathrm{P}}^{\left(r\right)}\left(t\right)$ stabilizes at the value associated with the ergodicity of the classical evolution. Why the random states are able to reach this greatest level of spreading to later contract to the value ${S}_{\mathrm{P}}^{\infty }$ of maximal quantum ergodicity is an open question to us.

5.2.3. Quantum asymptotic values and unstable periodic orbits

The purpose of this subsection is to show numerically that most initial coherent states in the chaotic region are indeed marginally affected by unstable periodic orbits, being well described by the behavior of the survival probability reported in figure 8.

In the central panel of figure 9, we show the ratio

$$\begin{equation}R=\frac{{S}_{\mathrm{P}}^{\left(r\right),\infty }}{{S}_{\mathrm{P}}^{\infty }}=\frac{2{\mathfrak{S}}_{\mathrm{P}}^{\infty }}{{S}_{\mathrm{P}}^{\infty }}\end{equation} \tag{ 30 }$$

for a mesh of several initial coherent states distributed over the same Poincaré surface as in figure 1(c). In the equation above, ${S}_{\mathrm{P}}^{\left(r\right),\infty }$ is given by the analytical expression in equation (28), which is the asymptotic value of an initial coherent state without the influence of unstable periodic orbits. We see that for most of the initial coherent states the ratio R is very close to one (light color), which implies that $\frac{{S}_{\mathrm{P}}^{\infty }}{{\mathfrak{S}}_{\mathrm{P}}^{\infty }}\approx 2$. There are few values of R that are smaller than one (dark color), $\frac{{S}_{\mathrm{P}}^{\infty }}{{\mathfrak{S}}_{\mathrm{P}}^{\infty }}{ >}2$, which indicates initial states affected by unstable periodic orbits. The ratio R is twice the factor $\mathcal{F}$ introduced in [33] to gauge the fraction of phase space explored by the quantum state.

Zoom In Zoom Out Reset image size

Around the central panel of figure 9, several plots of the quantum (blue solid line) and classical (red solid line) survival probabilities averaged over temporal windows are shown together with the analytical expression in equation (26) for random initial states (orange solid line). Panels (a), (b), (d)–(g) make it clear that a ratio $R\approx 1\;\left(\frac{{S}_{\mathrm{P}}^{\infty }}{{\mathfrak{S}}_{\mathrm{P}}^{\infty }}\approx 2\right)$ leads to a generic behavior of the quantum S_P(t), well described by ${S}_{\mathrm{P}}^{\left(r\right)}\left(t\right)$. In all these panels, there appears a ramp toward equilibration associated with the presence of correlated eigenvalues. These random-like coherent states contrast with those where ratio $R{< }1\;\left(\frac{{S}_{\mathrm{P}}^{\infty }}{{\mathfrak{S}}_{\mathrm{P}}^{\infty }}{ >}2\right)$, which are affected by quantum scarring, as in figures 9(c) and (h). These latter cases lead to recurrences in the evolution of the survival probability and a strong influence of the unstable periodic orbits in the quantum dynamics.

This work provides a comprehensive study of the dynamics and equilibration process of the Dicke model in the regular and chaotic regimes. The classical dynamics obtained via the truncated Wigner approximation (TWA) explains some features of the quantum evolution. We show that the quantum and classical survival probabilities agree extremely well up to the point where one of the two saturates. After equilibration, the classical ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$ reaches a constant value, while the quantum S_P(t) shows fluctuations due to the discreteness of the energy spectrum.

By comparing the entire quantum and classical evolutions, we identify properties other than the quantum fluctuations that are also purely quantum and analyze how they affect the equilibration times. We build a broad picture of the dynamics that counts with two general cases in the regular region and two in the chaotic region. These four scenarios are distinguished via the ratio between the asymptotic values of the quantum and classical survival probabilities, $\frac{{S}_{\mathrm{P}}^{\infty }}{{\mathfrak{S}}_{\mathrm{P}}^{\infty }}$. We itemize below our main findings.

• At short times, the survival probability of coherent states reaches values close to zero. This reflects the localization of the Wigner distribution in phase space and the fact that it moves out from its starting region. Later revivals of S_P(t) are associated with the return of the Wigner distribution to its starting classical region.
• The long-time dynamics in the regular regime depends on the proximity to a separatrix. For an initial state close to a separatrix, dynamical tunneling between classically disconnected phase-space regions leads to an asymptotic value of the quantum survival probability slightly lower than the classical one, $\frac{{S}_{\mathrm{P}}^{\infty }}{{\mathfrak{S}}_{\mathrm{P}}^{\infty }}< sim 1$, while away from the separatrix the asymptotic values coincide, $\frac{{S}_{\mathrm{P}}^{\infty }}{{\mathfrak{S}}_{\mathrm{P}}^{\infty }}=1$. This difference holds even though both states have the same energy.
• The equilibration times of the classical and quantum dynamics coincide in the regular regime, but differ in the chaotic region. In the latter case, dynamical manifestations of quantum chaos associated with quantum scars and spectral correlations separate the curves of the classical and quantum survival probability, one saturating earlier than the other.
• For initial coherent states in the chaotic region and affected by quantum scarring, the equilibration of the quantum dynamics happens before the classical one. In this case, the number of energy eigenstates that participate in the quantum evolution is small. This effectively shrinks the phase space available for the quantum evolution as compared to what is available for the classical evolution, so quantum equilibration occurs sooner. This is reflected in the value of the ratio $\frac{{S}_{\mathrm{P}}^{\infty }}{{\mathfrak{S}}_{\mathrm{P}}^{\infty }}{ >}2$. The larger this ratio is, the larger the degree of scarring.
• For initial states highly delocalized in the chaotic regime (ergodic states), the equilibration of the quantum dynamics happens after the classical one. This is due to the emergence of the correlation hole, which is caused by the correlations between the eigenvalues and is nonexistent in the classical limit. In this case, the quantum evolution of the survival probability at long times coincides with that for random initial states. This analogy allows us to derive analytically the asymptotic value ${S}_{\mathrm{P}}^{\infty }$ and the relaxation time. We find that $\frac{{S}_{\mathrm{P}}^{\infty }}{{\mathfrak{S}}_{\mathrm{P}}^{\infty }}=2$. This value signals the onset of maximal quantum ergodicity. The fact that this ratio is not 1 indicates that the quantum evolution is more restricted than the classical one due to remaining structures of the quantum states. We emphasize that the correlation hole is a universal behavior that emerges in any chaotic quantum model, provided the initial state is non-scarred.

Tunneling, quantum scars, and the correlation hole are purely quantum properties that are not exclusive to the Dicke model. We therefore expect our results to be applicable to other models with regular and chaotic regimes and with a quantum–classical correspondence.

We acknowledge the support of the Computation Center—ICN, in particular to Enrique Palacios, Luciano Díaz, and Eduardo Murrieta, and valuable conversations with Jonathan Torres and Jorge Chávez-Carlos. We acknowledge financial support from Mexican CONACyT project CB2015-01/255702, DGAPA-UNAM projects IN109417 and IN104020. LFS is supported by the NSF Grant No. DMR-1936006. LFS and JGH thank the hospitality of the Aspen Center for Physics and the Simons Center for Geometry and Physics at Stony Brook University, where some of the research for this paper was performed.

To obtain the classical Hamiltonian in equation (4), we use the Glauber (equation (2)) and Bloch (equation (3)) coherent states expressed in terms of the general parameters $\alpha ,z\in \mathbb{C}$ as a tensor product of the form,

$$\begin{align}\hfill \vert \alpha ,z\rangle =& \vert \alpha \rangle \otimes \vert z\rangle \hfill \\ \hfill =& \frac{{\mathrm{e}}^{-\vert \alpha {\vert }^{2}/2}}{{\left(1+\vert z{\vert }^{2}\right)}^{j}}\sum _{n=0}^{\infty }\sum _{m=-j}^{j}\sqrt{\left(\begin{pmatrix}2j\\ \hfill j+m\hfill \end{pmatrix}\right)}\frac{{\alpha }^{n}{z}^{j+m}}{\sqrt{n!}}\vert n\rangle \otimes \vert j,m\rangle ,\hfill \end{align} \tag{ A.1 }$$

and take the expectation value of the Dicke Hamiltonian Ĥ_D [73, 74],

$$\begin{align}\hfill {H}_{\mathrm{c}\mathrm{l}}=& \langle \alpha ,z\vert {\hat{H}}_{\mathrm{D}}\vert \alpha ,z\rangle \hfill \\ \hfill =& \omega \vert \alpha {\vert }^{2}-j{\omega }_{0}\frac{1-\vert z{\vert }^{2}}{1+\vert z{\vert }^{2}}+\gamma \sqrt{2j}\frac{z+{z}^{{\ast}}}{1+\vert z{\vert }^{2}}\left(\alpha +{\alpha }^{{\ast}}\right).\hfill \end{align} \tag{ A.2 }$$

Considering the harmonic oscillator $\alpha =\sqrt{\frac{j}{2}}\left(q+\mathrm{i}p\right)$ and the Bloch sphere $z=\sqrt{\frac{1+{j}_{z}}{1-{j}_{z}}}\;{\mathrm{e}}^{-\mathrm{i}\phi }$ parameters in terms of canonical variables (q, p) and (ϕ, j_z), we obtain

$$\begin{equation}{H}_{\mathrm{c}\mathrm{l}}=j\left[\frac{\omega }{2}\left({q}^{2}+{p}^{2}\right)+{\omega }_{0}{j}_{z}+2\gamma \sqrt{1-{j}_{z}^{2}}q\;\mathrm{cos}\left(\phi \right)\right].\end{equation} \tag{ A.3 }$$

To finally get equation (4), we perform a canonical transformation to the atomic variables (ϕ, j_z) → (Q, P), where $Q=\sqrt{2\left(1+{j}_{z}\right)}\mathrm{cos}\left(\phi \right)$ and $P=-\sqrt{2\left(1+{j}_{z}\right)}\mathrm{sin}\left(\phi \right)$ satisfy the Poisson bracket {Q, P} = 1, and rescale the overall classical Hamiltonian H_cl by j.

To obtain the representation of coherent states given by equations (2) and (3), we express the coherent state parameters α and z in terms of the canonical variables (q, p, Q, P), such that $\alpha =\sqrt{\frac{j}{2}}\left(q+\mathrm{i}p\right)$ and $z=\frac{1}{\sqrt{4-{Z}^{2}}}\left(Q+\mathrm{i}P\right)$ with Z² = Q² + P².

With the classical Hamiltonian H_cl, a semi-classical approximation to the DoS can be found by integrating the phase space volume available for a given energy. It reads [75, 101]

$$\begin{equation}\nu \left({\epsilon}\right)=\frac{2j}{\omega }\begin{cases}\frac{1}{\pi }{\int }_{{y}_{-}}^{{y}_{+}}\mathrm{d}y\;{\mathrm{cos}}^{-1}\left(\sqrt{\frac{2\left(y-{\epsilon}\right)}{{\overline{\gamma }}^{2}\left(1-{y}^{2}\right)}}\right),\quad {{\epsilon}}_{\mathrm{G}\mathrm{S}}{\leqslant}{\epsilon}{< }-1\quad \hfill \\ \frac{1+{\epsilon}}{2}+\frac{1}{\pi }{\int }_{{\epsilon}}^{{y}_{+}}\mathrm{d}y\;{\mathrm{cos}}^{-1}\left(\sqrt{\frac{2\left(y-{\epsilon}\right)}{{\overline{\gamma }}^{2}\left(1-{y}^{2}\right)}}\right),\quad \vert {\epsilon}\vert {\leqslant}1\quad \hfill \\ 1,\quad {\epsilon}{ >}1\quad \hfill \end{cases},\end{equation} \tag{ A.4 }$$

where ${y}_{{\pm}}=-{\overline{\gamma }}^{-1}\left({\overline{\gamma }}^{-1}\mp \sqrt{2\left({\epsilon}-{{\epsilon}}_{0}\right)}\right)$ and $\overline{\gamma }=\gamma /{\gamma }_{\mathrm{c}}$. The ground state energy is given by _GS = −1 for the normal phase, and by ${{\epsilon}}_{\mathrm{G}\mathrm{S}}={{\epsilon}}_{0}=-\frac{1}{2}\left({\overline{\gamma }}^{2}+{\overline{\gamma }}^{-2}\right)$ for the superradiant phase. In figure A1, we compare the DoS obtained numerically with equation (A.4), showing excellent agreement. The figure also shows the energies selected for our studies throughout this paper.

Zoom In Zoom Out Reset image size

As we see from appendix A, the coherent states for the Dicke model are the product of Glauber and Bloch coherent states (equation (A.1)). Thus, the associated Wigner function is the product of the Wigner functions associated to the Glauber and Bloch states. For a Glauber state |α⟩ with $\alpha \left(q,p\right)=\sqrt{\frac{j}{2}}\left(q+\mathrm{i}p\right)$, the Wigner function is given by a normal distribution

$$\begin{equation}{w}_{{q}_{0},{p}_{0}}\left(q,p\right)=\frac{j}{\pi }\;{\mathrm{e}}^{-j{d}^{2}},\end{equation} \tag{ B.1 }$$

with $d=\sqrt{{\left(q-{q}_{0}\right)}^{2}+{\left(p-{p}_{0}\right)}^{2}}$. For a Bloch coherent state, the Wigner function in variables (θ, ϕ) may be written as a sum of Legendre polynomials P_k(x) [102]

$$\begin{equation}{w}_{{\theta }_{0},{\phi }_{0}}\left(\theta ,\phi \right)=\frac{\left(2j\right)!}{4\pi }\sum _{k=0}^{2j}\sqrt{\frac{\left(2k+1\right)}{\left(2j-k\right)!\left(2j+k+1\right)!}}{P}_{k}\left(\mathrm{cos}\;{\Theta}\right),\end{equation} \tag{ B.2 }$$

where Θ is the angle between (θ, ϕ) and (θ₀, ϕ₀) obtained from

$$\begin{equation}\mathrm{cos}\;{\Theta}=\mathrm{cos} \theta \;\mathrm{cos}\;{\theta }_{0}+\mathrm{sin} \theta \;\mathrm{sin}\;{\theta }_{0}\;\mathrm{cos}\left(\phi -{\phi }_{0}\right).\end{equation} \tag{ B.3 }$$

For large j values this is very well approximated by a normal distribution on the Bloch sphere

$$\begin{equation}{w}_{{\theta }_{0},{\phi }_{0}}\left(\theta ,\phi \right)\approx \frac{j}{\pi }\;{\mathrm{e}}^{-j{{\Theta}}^{2}}.\end{equation} \tag{ B.4 }$$

This approximation is very useful, because sampling from normal distributions is easy and computationally cheap. The complete coherent Wigner function is just the product of the Wigner functions given by equations (B.1) and (B.4),

$$\begin{equation}{w}_{{\boldsymbol{u}}_{\mathbf{0}}}\left(\boldsymbol{u}\right)={\left(\frac{j}{\pi }\right)}^{2}{\mathrm{e}}^{-j\left({d}^{2}+{{\Theta}}^{2}\right)}.\end{equation} \tag{ B.5 }$$

The temporal evolution of the Wigner function is governed by the so-called Moyal equation⁸ [103]

$$\begin{equation}\frac{\partial w}{\partial t}\left(\boldsymbol{u},t\right)={\left\{w\left(\boldsymbol{u},t\right),{h}_{\mathrm{c}\mathrm{l}}\left(\boldsymbol{u}\right)\right\}}_{M}.\end{equation} \tag{ C.1 }$$

Here, ${\left\{A,B\right\}}_{M}=\frac{2}{\hslash }A\;\mathrm{sin}\left[\frac{\hslash }{2}\left(\overleftarrow{{\partial }_{\boldsymbol{q}}} \overrightarrow {{\partial }_{\boldsymbol{p}}}-\overleftarrow{{\partial }_{\boldsymbol{p}}} \overrightarrow {{\partial }_{\boldsymbol{q}}}\right)\right]B$ represents the Moyal bracket. Taylor expanding the sine, one may write ${\left\{A,B\right\}}_{M}=\left\{A,B\right\}+\mathcal{O}\left({\hslash }^{2}\right)$ where $\left\{\cdot ,\cdot \right\}$ is the Poisson bracket. So

$$\begin{equation}\frac{\partial w}{\partial t}\left(\boldsymbol{u},t\right)=\left\{w\left(\boldsymbol{u},t\right),{h}_{\mathrm{c}\mathrm{l}}\left(\boldsymbol{u}\right)\right\}+\mathcal{O}\left({j}^{-2}\right).\end{equation} \tag{ C.2 }$$

If we ignore the j⁻² order terms in this equation, we are left with the classical Liouville equation. This is known as the truncated Wigner approximation (TWA) and yields the correct quantum evolution for small times.

Within this approximation, the Wigner function remains constant along classical trajectories in phase space, so the time dependence of w(u, t) may be written in terms of the Hamiltonian flow ${\boldsymbol{\varphi }}^{t}:\mathcal{M}\to \mathcal{M}$. This function describes the time evolution of an initial condition ${\boldsymbol{u}}_{\boldsymbol{0}}\in \mathcal{M}$ by

$$\begin{equation}\boldsymbol{u}\left(t\right)={\boldsymbol{\varphi }}^{t}\left({\boldsymbol{u}}_{\boldsymbol{0}}\right),\end{equation} \tag{ C.3 }$$

and it satisfies the one-parameter group identities φ⁰ = Id, ${\boldsymbol{\varphi }}^{-t}={\left({\boldsymbol{\varphi }}^{t}\right)}^{-1}$, and ${\boldsymbol{\varphi }}^{{t}_{1}+{t}_{2}}={\boldsymbol{\varphi }}^{{t}_{2}}{\circ}{\boldsymbol{\varphi }}^{{t}_{1}}$. Staying constant along classical trajectories means that for any pair of times t₁ and t₂

$$\begin{equation}w\left({\boldsymbol{\varphi }}^{{t}_{1}}\left(\boldsymbol{u}\right),{t}_{1}\right)=w\left({\boldsymbol{\varphi }}^{{t}_{2}}\left(\boldsymbol{u}\right),{t}_{2}\right).\end{equation} \tag{ C.4 }$$

In particular, taking t₂ = 0 and performing an adequate change of variables,

$$\begin{equation}w\left(\boldsymbol{u},t\right)=w\left({\boldsymbol{\varphi }}^{-t}\left(\boldsymbol{u}\right),0\right).\end{equation} \tag{ C.5 }$$

Inserting equation (C.5) into (21), we find equation (22).

Note that equation (22) may be written as

$$\begin{equation}{\mathfrak{S}}_{\mathrm{P}}\left(t\right)={\left\langle w\left({\boldsymbol{\varphi }}^{-t}\left(\boldsymbol{u}\right)\right)\right\rangle }_{w}.\end{equation} \tag{ C.6 }$$

Provided that w is everywhere positive, as it is for coherent states (see equation (B.5)), this expected value may be efficiently approximated using a Monte Carlo method,

$$\begin{equation}{\mathfrak{S}}_{\mathrm{P}}\left(t\right)\approx \frac{1}{M}\sum _{i=1}^{M}w\left({\boldsymbol{\varphi }}^{-t}\left({\boldsymbol{u}}_{\boldsymbol{i}}\right)\right),\end{equation} \tag{ C.7 }$$

where the points ${\boldsymbol{u}}_{\boldsymbol{i}}\in \mathcal{M}$ are randomly sampled from the initial distribution w, and M is a sufficiently large, albeit computationally accessible, integer.

For the plots of ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$ shown in figures 4, 5, 7, and 8, and for the animations available in the supplementary data, we use a value of M = 10⁷ to significantly reduce numerical noise. For figures 6 and 9, because temporal averages reduce the numerical error, we only need to use M = 3 × 10⁵.

An intuitive understanding of the classical survival probability may be obtained by considering a simple measurable set $\mathcal{S}\subseteq \mathcal{M}$ (a sphere, for example). Let ${\boldsymbol{1}}_{\mathcal{S}}$ be its indicator function⁹. Consider the normalized distribution ${\rho }_{\mathcal{S}}={\mathbf{1}}_{\mathcal{S}}/{V}_{\mathcal{S}}$ in lieu of the Wigner function. Here ${V}_{\mathcal{S}}=\int \mathrm{d}\boldsymbol{u}\;{\mathbf{1}}_{\mathcal{S}}\left(\boldsymbol{u}\right)$ is the volume of $\mathcal{S}$. In this case, to retain normalization, our prefactor becomes $1/{V}_{\mathcal{S}}$ instead of (2π/j)², and then, from equation (22),

$$\begin{equation}{\mathfrak{S}}_{\mathrm{P}}\left(t\right)=\frac{1}{{V}_{\mathcal{S}}}{\int }_{\mathcal{M}}\mathrm{d}\boldsymbol{u}\;{\boldsymbol{1}}_{\mathcal{S}}\left(\boldsymbol{u}\right){\boldsymbol{1}}_{\mathcal{S}}\left({\boldsymbol{\varphi }}^{-t}\left(\boldsymbol{u}\right)\right).\end{equation} \tag{ D.1 }$$

Since ${\mathbf{1}}_{\mathcal{S}}\left({\boldsymbol{\varphi }}^{-t}\left(\boldsymbol{u}\right)\right)={\mathbf{1}}_{{\boldsymbol{\varphi }}^{t}\left(\mathcal{S}\right)}\left(\boldsymbol{u}\right)$ and using that the product of the indicator functions of two sets is the indicator function of their intersection, we get

$$\begin{equation}{\mathfrak{S}}_{\mathrm{P}}\left(t\right)=\frac{1}{{V}_{\mathcal{S}}}{\int }_{\mathcal{M}}\mathrm{d}\boldsymbol{u}\;{\boldsymbol{1}}_{\mathcal{S}\cap {\boldsymbol{\varphi }}^{t}\left(\mathcal{S}\right)}\left(\boldsymbol{u}\right)=\frac{{V}_{\mathcal{S}\cap {\boldsymbol{\varphi }}^{t}\left(\mathcal{S}\right)}}{{V}_{\mathcal{S}}}.\end{equation} \tag{ D.2 }$$

This result is easy to interpret: ${\mathfrak{S}}_{\mathrm{P}}\left(t\right)$ is the percentage of the volume of $\mathcal{S}$ that is found back inside $\mathcal{S}$ after time t. In other words, it is the probability that a point sampled from $\mathcal{S}$ is back in $\mathcal{S}$ after time t [30, 104].

To obtain an expression for ${\mathfrak{S}}_{\mathrm{P}}^{\infty }$, we assume ergodicity. This is a reasonable assumption for chaotic behaviors. If the Hamiltonian flow φ^t is ergodic over the energy shells in phase space, then, for any real function f(u) of phase space, temporal averages in composition with the flow are equal to space averages over the corresponding energy shell, that is, for any fixed point $\boldsymbol{u}\in \mathcal{M}$ with energy E = H_cl(u),

$$\begin{align}\hfill {\left\langle f\left({\boldsymbol{\varphi }}^{t}\left(\boldsymbol{u}\right)\right)\right\rangle }_{t\to \infty }=& {\left\langle f\right\rangle }_{E}\hfill \\ \hfill =& \frac{{j}^{2}}{{\left(2\pi \right)}^{2}\nu \left(E\right)}{\int }_{\mathcal{M}}\mathrm{d}\boldsymbol{v}\;\delta \left({H}_{\mathrm{c}\mathrm{l}}\left(\boldsymbol{v}\right)-E\right)f\left(\boldsymbol{v}\right),\hfill \end{align} \tag{ E.1 }$$

where j⁻²(2π)²ν(E) is the volume of the energy shell for E/j in $\mathcal{M}$, and ν(E) is given by equation (A.4).

It is then straightforward to calculate

$$\begin{align}\hfill {\left\langle {\mathfrak{S}}_{\mathrm{P}}\left(t\right)\right\rangle }_{t\to \infty }=& {\left(\frac{2\pi }{j}\right)}^{2}{\int }_{\mathcal{M}}\mathrm{d}\boldsymbol{u}\;w\left(\boldsymbol{u}\right){\left\langle w\left({\boldsymbol{\varphi }}^{-t}\left(\boldsymbol{u}\right)\right)\right\rangle }_{t\to \infty }\hfill \\ \hfill =& {\left(\frac{2\pi }{j}\right)}^{2}{\int }_{\mathcal{M}}\mathrm{d}\boldsymbol{u}\;w\left(\boldsymbol{u}\right){\left\langle w\right\rangle }_{{H}_{\mathrm{c}\mathrm{l}}\left(\boldsymbol{u}\right)}\hfill \\ \hfill =& {\left(\frac{2\pi }{j}\right)}^{4}{\int }_{{E}_{GS}}^{\infty }\mathrm{d}E\;{\left\langle w\right\rangle }_{E}^{2}\nu \left(E\right).\hfill \end{align} \tag{ E.2 }$$

The last equality is obtained by using $1{=\int }_{{E}_{\mathrm{G}\mathrm{S}}}^{\infty }\mathrm{d}E\;\delta \left(E-{H}_{\mathrm{c}\mathrm{l}}\left(\boldsymbol{u}\right)\right)$ inside the integral, a change in the integration order, and a substitution of the value of ${\left\langle w\right\rangle }_{E}$. In the above, E_GS represents the ground state energy. The classical energy distribution for the state associated to w is

$$\begin{align}\hfill {\rho }_{\mathrm{c}\mathrm{l}}\left(E\right)=& {\int }_{\mathcal{M}}\mathrm{d}\boldsymbol{u}\;\delta \left(E-{H}_{\mathrm{c}\mathrm{l}}\left(\boldsymbol{u}\right)\right)w\left(\boldsymbol{u}\right)\hfill \\ \hfill =& {\left(2\pi \right)}^{2}{j}^{-2}\nu \left(E\right){\left\langle w\right\rangle }_{E},\hfill \end{align} \tag{ E.3 }$$

therefore

$$\begin{equation}{\left\langle {\mathfrak{S}}_{\mathrm{P}}\left(t\right)\right\rangle }_{t\to \infty }={\int }_{{E}_{\mathrm{G}\mathrm{S}}}^{\infty }\mathrm{d}E\;\frac{{\rho }_{\mathrm{c}\mathrm{l}}^{2}\left(E\right)}{\nu \left(E\right)}.\end{equation} \tag{ E.4 }$$

This result is exact, but we may approximate ν(E) (see equation (A.4)) by the value at the center of the classical energy distribution (ν(E) ≈ ν(E_c) = ν_c) to get

$$\begin{equation}{\left\langle {\mathfrak{S}}_{\mathrm{P}}\left(t\right)\right\rangle }_{t\to \infty }=\frac{1}{{\nu }_{\mathrm{c}}}{\int }_{{E}_{\mathrm{G}\mathrm{S}}}^{\infty }\mathrm{d}E\;{\rho }_{\mathrm{c}\mathrm{l}}^{2}\left(E\right).\end{equation} \tag{ E.5 }$$

By further approximating ρ_cl with the Gaussian distribution given by equation (12), we obtain equation (25).

An analytical expression for the survival probability averaged over an ensemble of random initial states was derived in [27]. The properties of the ensemble are as follows. Its members are constrained to have a smooth LDoS ρ(E) and their energy components are given by

$$\begin{equation}\vert {c}_{k}^{\left(r\right)}{\vert }^{2}=\frac{{r}_{k}\;\rho \left({E}_{k}\right)}{\mathcal{A}\;\nu \left({E}_{k}\right)},\end{equation} \tag{ F.1 }$$

where ν(E) is the DoS, the numbers r_k are randomly generated from a given probability distribution P(r), and $\mathcal{A}={\sum }_{q}{r}_{q}\rho \left({E}_{q}\right)/\nu \left({E}_{q}\right)$ is a normalization constant. This gives,

$$\begin{equation}{S}_{\mathrm{P}}^{\left(r\right)}\left(t\right)=\frac{1-{S}_{\mathrm{P}}^{\left(r\right),\infty }}{\eta -1}\left[\eta {S}_{\mathrm{P}}^{\text{bc}}\left(t\right)-{b}_{2}\left(\frac{Dt}{2\pi }\right)\right]+{S}_{\mathrm{P}}^{\left(r\right),\infty },\end{equation} \tag{ F.2 }$$

where the initial decay of the survival probability is dictated by

$$\begin{equation}{S}_{\mathrm{P}}^{\text{bc}}\left(t\right)={\left\vert \int \mathrm{d}E\rho \left(E\right){\mathrm{e}}^{-\mathrm{i}Et}\right\vert }^{2},\end{equation} \tag{ F.3 }$$

the effective dimension of the ensemble

$$\begin{equation}\eta =\frac{1}{\int \mathrm{d}E{\rho }^{2}\left(E\right)/\nu \left(E\right)}\approx \frac{{\nu }_{\mathrm{c}}}{\int \mathrm{d}E\;{\rho }^{2}\left(E\right)},\end{equation} \tag{ F.4 }$$

where ν_c = ν(E_c) is the density of states evaluated in the center of the energy profile, and the asymptotic value of the survival probability corresponds to

$$\begin{equation}{S}_{\mathrm{P}}^{\left(r\right),\infty }=\frac{\langle {r}^{2}\rangle }{\langle {r\rangle }^{2}}\frac{1}{\eta },\end{equation} \tag{ F.5 }$$

with ⟨rⁿ⟩ the nth moments of the distribution P(r).

The function b₂ is the two-level form factor of the GOE [95]

$$\begin{equation}{b}_{2}\left(\overline{t}\right)=\left[1-2\overline{t}+\overline{t}\;\mathrm{ln}\left(2\overline{t}+1\right)\right]{\Theta}\left(1-\overline{t}\right)+\left[\overline{t}\;\mathrm{ln}\left(\frac{2\overline{t}+1}{2\overline{t}-1}\right)-1\right]{\Theta}\left(\overline{t}-1\right),\end{equation} \tag{ F.6 }$$

where Θ is the Heaviside step function. The factor D in the argument of b₂ is the mean level spacing of the correlated eigenvalues.

For the random ensemble that we consider in this paper, ρ(E) has a Gaussian energy profile (see equation (12)) and the random numbers r_k are generated from the exponential distribution P(r) = λ e^−λr (see figure 3(m)) with ⟨rⁿ⟩ = n!/λⁿ, which implies that

$$\begin{equation}{S}_{\mathrm{P}}^{\text{bc}}\left(t\right)={\mathrm{e}}^{-{\sigma }^{2}{t}^{2}},\end{equation} \tag{ F.7 }$$

$$\begin{equation}\eta =2\sqrt{\pi }\sigma {\nu }_{\mathrm{c}},\end{equation} \tag{ F.8 }$$

$$\begin{equation}{S}_{\mathrm{P}}^{\left(r\right),\infty }=\frac{2}{\eta }=\frac{1}{\sqrt{\pi }\sigma {\nu }_{\mathrm{c}}}.\end{equation} \tag{ F.9 }$$

Because the correlations in the spectrum appears only for energy levels in the same parity sector, the mean level spacing of correlated eigenvalues is D = (D₊ + D₋)/2, where the mean level spacing for each parity sector, D_± = 1/ν_±, is given by the respective density of states. These densities are, in turn, given by ν_± = ν_c/2, with ν_c the density of states of the whole spectrum, yielding D = 2/ν_c. From the previous results for ${S}_{\mathrm{P}}^{\left(r\right),\infty }$, η, ${S}_{\mathrm{P}}^{\text{bc}}$, and D, we get equation (26).

To derive the relaxation time t_r of the ensemble-averaged survival probability, we consider the asymptotic form of b₂, which grows toward saturation following a power-law behavior

$$\begin{equation}{b}_{2}\left(\frac{t}{\pi {\nu }_{\mathrm{c}}}\right)\to \frac{{\pi }^{2}{\nu }_{\mathrm{c}}^{2}}{12{t}^{2}} \text{for}\;\frac{t}{\pi {\nu }_{\mathrm{c}}}\gg 1.\end{equation} \tag{ F.10 }$$

At this temporal scale, the contribution of the initial decay ${S}_{\mathrm{P}}^{\text{bc}}$ is negligible and the asymptotic form of (F.2) is given by

$$\begin{equation}{S}_{\mathrm{P}}^{\left(r\right)}\left(t\right)\to -\frac{1}{\eta }{b}_{2}\left(\frac{t}{\pi {\nu }_{\mathrm{c}}}\right)+{S}_{\mathrm{P}}^{\left(r\right),\infty }={S}_{\mathrm{P}}^{\left(r\right),\infty }\left(1-\frac{{\pi }^{2}{\nu }_{\mathrm{c}}^{2}}{24{t}^{2}}\right),\end{equation} \tag{ F.11 }$$

where in the last step we have used $\eta =2/{S}_{\mathrm{P}}^{\left(r\right),\infty }$. We define the relaxation time according to

$$\begin{equation}{S}_{\mathrm{P}}^{\left(r\right)}\left({t}_{\mathrm{r}}\right)=\left(1-{\delta }_{{S}_{\mathrm{P}}}\right){S}_{\mathrm{P}}^{\left(r\right),\infty },\end{equation} \tag{ F.12 }$$

where ${\delta }_{{S}_{\mathrm{P}}}$ is a small parameter determining the point where ${S}_{\mathrm{P}}^{\left(r\right)}\left(t\right)$ is already within the fluctuations around the asymptotic value. By substituting equation (F.11) in (F.12), we obtain equation (27).