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

We investigate the entire evolution of the Dicke model, which is a two-degree-of-freedom interacting spin-boson model of great experimental interest. Our objects of study are the quantum and classical survival probabilities of initial coherent states and the corresponding classical evolution of the Wigner distribution in phase space. We show that major aspects of the system are uncovered by analyzing its long-time dynamics, such as whether the initial state is in a regular or chaotic region, in the vicinity of a separatrix, or yet close to an unstable periodic orbit. We demonstrate that a ratio of two between the quantum and classical asymptotic values of the survival probability is a clear indicator of maximal quantum ergodicity. In this case, the quantum survival probability develops a correlation hole, that is nonexistent in its classical version and results in a longer equilibration time for the quantum dynamics. These findings are corroborated by our analytical expressions for the survival probability and the equilibration time.

We investigate the entire evolution of the Dicke model, which is a two-degree-of-freedom interacting spin-boson model of great experimental interest. Our objects of study are the quantum and classical survival probabilities of initial coherent states and the corresponding classical evolution of the Wigner distribution in phase space. We show that major aspects of the system are uncovered by analyzing its long-time dynamics, such as whether the initial state is in a regular or chaotic region, in the vicinity of a separatrix, or yet close to an unstable periodic orbit. We demonstrate that a ratio of two between the quantum and classical asymptotic values of the survival probability is a clear indicator of maximal quantum ergodicity. In this case, the quantum survival probability develops a correlation hole, that is nonexistent in its classical version and results in a longer equilibration time for the quantum dynamics. These findings are corroborated by our analytical expressions for the survival probability and the equilibration time.
At long times, equilibration eventually takes place. In classical mechanics, the mixing properties of chaotic dynamics have provided a fundamental mechanism to explain the equilibration process and the ergodic properties of physical systems. The passage to the quantum domain entails new phenomena such as superpositions, quantum interferences, and the effects of the quantized spectra [16][17][18][19][20][21][22][23][24] and quantum scars [25][26][27][28][29][30][31][32][33][34][35]. Even though isolated quantum 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 the observable reaches a point, where it simply fluctuates around its asymptotic value, and these fluctuations decrease with the system size [36][37][38][39][40].
The quantity that we use for our studies is the probability to find the initial state at a later time, the so-called survival probability or return probability. We consider initial coherent states, which enables a direct comparison between the exact quantum evolution and its classical description obtained with the truncated Wigner approximation (TWA) [64][65][66]. This comparison allows for the identification of features in the equilibration process that are purely quantum. This includes quantum fluctuations and the effects of correlations between the eigenvalues. Both are associated with the discreteness of the spectrum and therefore manifest themselves at long times.
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 the different features of the quantum dynamics that emerge at different time scales.
We place initial coherent states in the regular and chaotic region and show how they can be distinguished by analyzing the behavior of the survival probability at long times. The quantum and classical relaxation times do not coincide in the chaotic regime, but they do in the regular case. We also demonstrate how the ratio between the asymptotic values of the quantum and classical survival probabilities tells us about the proximity of the initial state to the separatrix of the model or to an unstable periodic orbit. In addition, this ratio can be used to gauge the degree of scarring of the initial state. One of our strongest results is to show that a ratio equal to two indicates maximal quantum ergodicity. For this limit, we have an analytical expression for the entire evolution of the survival probability and for the equilibration time.
This paper is organized as follows. The core of the work is in Sec. IV, 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 to this analysis, Sec. II describes the Dicke Hamiltonian and its classical limit, as well as the initial coherent states studied. In Sec. III, 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 Sec. V.

II. DICKE MODEL
The Dicke model [41] represents a set of N two-level atoms with atomic transition frequency ω 0 interacting with a single mode of a radiation field with frequency ω. It is described by the following Hamiltonian,Ĥ where = 1,â (â † ) is the bosonic annihilation (creation) operator of the field mode,Ĵ x,y,z = 1 2 N k=1σ k x,y,z are collective pseudo-spin operators given by the sum of Pauli matricesσ x,y,z , and γ is the spin-boson interaction strength. When γ reaches a critical value γ c = √ ωω 0 /2, a secondorder quantum phase transition takes place [42,43]. There, the system 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Ĵ 2 =Ĵ 2 x +Ĵ 2 y +Ĵ 2 z determine the different invariant subspaces. We work with the maximum value j = N /2, which defines a symmetric atomic subspace that includes the ground state. The HamiltonianĤ D commutes also with the parity operatorΠ = e iπΛ , whereΛ =â †â +Ĵ z + j1 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 .

A. Classical Limit of Dicke Hamiltonian
The corresponding classical Hamiltonian is obtained using Glauber coherent states for the bosonic sectors [52,53,67,68], and Bloch coherent states for the pseudo-spin sectors, where Z 2 = Q 2 +P 2 . The canonical j-independent variables (q, p) and (Q, P ) are associated with the photonic and atomic degrees of freedom, respectively. The ket |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), h cl ≡ q, p; Q, P |Ĥ D |q, p; Q, P j (4) The rescaled classical Hamiltonian and its four-dimensional phase space M are independent of j. This is equivalent to working with an effective Planck constant eff = 1/j [69]. Depending on the parameters and energies, the model displays chaotic or regular dynamics. As a case study, we choose ω = ω 0 , j = 100, and a coupling strength in the superradiant phase γ = 2γ c . For these parameters, the normalized energy of the ground state is GS = −2.125. The dynamics is regular from GS up to ≈ −1.7 and chaotic above this point [53].

B. Initial States
Our initial states are the Glauber-Bloch coherent states |q, p; Q, P . They allow for a direct connection between the quantum states and the points 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 select the initial states, we restrict ourselves to the hyperplane p = 0 and solve the seconddegree equation h cl (q, p, Q, P ) = ω 0 in q. This equation has two solutions, q − and q + , with q − ≤ q + . Our initial states are centered at (q, p, Q, P ) = (q + , 0, Q 0 , P 0 ). Two of them have energy in the regular region, R = −1.8, and two of them have the energy shell fully covered by chaotic trajectories, C = −0.5. The criteria for our choices and the specific values of Q 0 and P 0 are described below.
Poincaré sections at p = 0 of the classical Hamiltonian at energy R = −1.8 are shown in Fig. 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 0 , P 0 ) = (1, 0) and indicate it with a blue dot in Fig. 1 (a). State II is close to the separatrix, (Q 0 , P 0 ) = (1.2, 0), and is shown with a red dot in the same figure.
The structure of the phase space reflects the quasi-conserved quantities of the Dicke Hamiltonian at low energies [70]. 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 [71]. For the parameters considered here, a quasiconstant 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 Fig. 1 (a) corresponds to modes 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 [72] and are separated by a separatrix, around which a stochastic layer is formed [73]. 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 Fig. 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 Fig. 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, where c k = E k |q + , 0; Q 0 , P 0 andĤ|E k = E k |E k . The participation ratio is strongly correlated with the underlying classical dynamics [53,74]. For instance, the stochastic layer appearing in the classical dynamics around the separatrix is associated with larger values of P R , which indicate states that are more delocalized in the energy eigenbasis. Poincaré sections for the selected energy C = −0.5 are shown in Fig. 1 (c). They reveal a region of hard chaos, where chaotic trajectories, all with the same positive Lyapunov exponent, densely fill the whole phase space. From these Poincaré sections, no particular region can be identified, but the participation-ratio map in Fig. 1 (d) provides a richer picture, with coherent states showing different levels of spreading in the energy eigenbasis. To analyze how the structure of the initial states affects the dynamics and equilibration, we select two initial states. State III, indicated with the cyan point in Fig. 1 (d), is located in the region with small values of P R , at (Q 0 , P 0 ) = (1.75, 0). State IV, marked with a red point in Fig. 1 (d), is in the region of large values of P R , at (Q 0 , P 0 ) = (−1.25, 0.75). These two states are representative of the typical cases used in the studies of the dynamics of chaotic regions in [28].

III. SURVIVAL PROBABILITY
The survival probability, S P , is the probability of finding an evolved quantum state back in the initial state |Ψ(0) = k c k |E k , By introducing the local density of states (LDoS) or strength function, that is the energy distribution weighted by the components |c k | 2 of the initial state, we can also write the survival probability as the squared norm of the Fourier transform of G(E) as The evolution of the survival probability shows different behaviors at different time scales [20,21,75,76]. 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 where , and sinc(x) = sin(x)/x. The time resolution T = π/∆ reflects aspects of the LDoS that are of order ∆ in energy.

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, 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 [77][78][79]. According to Eq. (9), the Gaussian envelope of the LDoS leads to the initial Gaussian decay of the survival probability, which is consistent with the universal quadratic behavior, S P (t σ −1 ) ≈ 1 − σ 2 t 2 , for very short times.
In addition to the initial coherent states, we show in Figs. 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 |c (r) Above, ρ(E k ) is obtained from Eq. (11), ν(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 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 for the components c k of the coherent state IV in the energy eigenbasis, we obtain the histogram in Fig. 3 (m), which is very well fitted with an exponential distribution. By comparing the smoothed LDoS of the state III in Fig. 3 (d), the state IV in Fig. 3 (e), and the random state in Fig. 3 (f), all of them with time resolution ω 0 T = 2, one notices that the smoothed LDoS of the random state already deviates from a Gaussian, which must affect its shorttime 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 Sec. IV. For an even higher time resolution, as ω 0 T = 3 in Figs. 3 (g)-(i), the smoothed LDoS of state III also deviates from the Gaussian, and so does its Gaussian decay (see Sec. IV).  Fig. 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 [79]. In contrast, the components of the chaotic coherent states in

Survival Probability: Asymptotic Values
The asymptotic value of the survival probability, can be derived from In the absence of energy degeneracies, the first term on the right-hand-side of the equation cancels out on average, so which is the inverse of the participation ratio shown in Eq. (6). However, when the the energy levels have degeneracies of degree d k (that is, |E k , m with m = 1, .., d k ), the first term in Eq. (16) contributes with additional terms to the asymptotic value, which is now given by 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 [80] produces an energy spectrum with two-fold degeneracies 1 . 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.

B. Classical Limit of the Survival Probability
To find the classical limit of the survival probability, we use the Wigner formalism [81]. 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 [82], 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(u, 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 ϕ t : M → M. One finds that w(u, t) = w(ϕ −t (u), 0) for short times, so the classical survival probability can be defined as where w(u) = w(u, t = 0). This quantity is numerically constructed through a Monte Carlo method. Its asymptotic value, S ∞ P = S P (t) t→∞ , is obtained with Eq. (15).

IV. CLASSICAL AND QUANTUM DYNAMICS
By comparing the behaviors of the quantum and classical survival probabilities for the four initial coherent states selected in Sec. II B, we are able: -To understand why S P (t) for coherent initial states reaches values close to zero at short times.
-To track down the origins of the different long-time behaviors of initial states that have the same energy.
-To identify properties of the survival probability that are purely quantum and which are seen only at long times.
-To explain why the saturation times of the quantum and classical dynamics in the chaotic regime do not coincide.

A. Regular Region
We analyze first the quantum S P (t) and the classical S P (t) for the initial coherent states I and II [ Fig. 1 (a)], which are in the regular regime ( R = −1.8).

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 [79]. It is given by where the index n denotes the distance between the eigenenergies of the levels with no negligible components, ω 1 is related with the spacing between neighboring levels, σ is the width of the LDoS, t D = ω 1 /(σ|e 2 |) is the decay time, with e 2 being the anharmonicity of the spectrum, and The time t D comes from the Gaussian decay of the slowest component n = 1.
In the bottom panel of Fig. 4, we compare the full quantum evolution of S P (t) with the analytical expression in Eq. (22) and the classical S P (t). 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 ∞ P that emerge due to the discreteness of the spectrum and which the analytical expression in Eq. (22) and S P (t) are unable to reproduce.
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 Fig. 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, 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) [83]): -At ω 0 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 ω 0 t = 2 (second column), the distribution is effectively outside the initial region and the survival probability becomes close to zero 2 .
-After a full period, the distribution comes back to the starting region at ω 0 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 ω 0 t = 100 (fourth column).
-Classically, at ω 0 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 S P (t) reaches the constant value S ∞ P , while S P (t) fluctuates around the asymptotic value S ∞ P ≈ S ∞ P .

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 [ Fig. 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 [ Fig. 1 (a)]. An analytical expression can also be obtained for this case by generalizing Eq. (22), where in addition to contributions from Gaussians, we also need to take into account interferences between them [79].
In the bottom panel of Fig. 5, we show the quantum and classical survival probability for the state II. The agreement between the two is good, but contrary to Fig. 4, the values of S P (t) are slightly smaller than those for S P (t). This is better seen with the temporal averages of the curves, shown in Fig. 5 in blue for S P (t) and in red for S P (t).
FIG. 5. Top panels: Classical evolution of the Wigner distribution projected onto the Q-P plane (top row) and q-p plane (bottom row) for the initial coherent state II from Fig. 1 (a). The initial distribution (small black circle) is inside the available phase space (big black circle). Each column represents a time, as indicated. (Animations can be found in the SM [83].) Bottom panel: Quantum survival probability (light gray solid line), its temporal average (blue solid line), classical survival probability (dark gray solid line), and its temporal average (red solid line) for the initial coherent state II. The temporal averages are computed for temporal windows of constant size in the logarithmic scale. The horizontal blue dotted line is the quantum asymptotic value S ∞ P and the red dashed line corresponds to the classical asymptotic value S ∞ P .
To better understand the small differences between S P (t) and S P (t), 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 Fig. 1 (a). We plot the values of S ∞ P and S ∞ P in Fig. 6. One sees that S ∞ P is smaller than S ∞ P 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 S P (t) to a dynamic tunneling effect that takes place in the quantum regime [84], but is absent in the classical limit. 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 Fig. 4 and explains specific features of S P (t), such as the revival at ω 0 t = 6.5. At very long times, such as ω 0 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 Fig. 5 is one order of magnitude smaller than in Fig. 4.

B. Chaotic Region
We now study the quantum S P (t) and the classical S P (t) for the initial coherent states III and IV [ Fig. 1 (d)], which are in the chaotic regime ( R = −0.5). We start the analysis with state III, which has a smaller number of contributing energy eigenbasis (smaller P R ) than state IV.

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 suggested by the low value of P R in Fig. 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 Fig. 7. There is excellent agreement from ω 0 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 Fig. 7. Analogously to the discussions about Fig. 4 and Fig. 5, as the originally localized Wigner distribution at ω 0 t = 0 (first column) FIG. 7. Top panels: Classical evolution of the Wigner distribution projected onto the Q-P plane (top row) and q-p plane (bottom row) for the initial coherent state III from Fig. 1 (d). The initial distribution (small black circle) is inside the available phase space. Each column represents a time, as indicated. (Animations can be found in the SM [83].) The black arrows indicate an enhancement of the distribution around an unstable periodic orbit (see text). Bottom panel: Quantum survival probability (light gray solid line), its temporal average (blue solid line), classical survival probability (dark gray solid line), and its temporal average (red solid line) for the initial coherent state III. The temporal averages are computed for temporal windows of constant size in the logarithmic scale. The horizontal blue dotted line is the quantum asymptotic value S ∞ P and the horizontal red dashed line corresponds to the classical asymptotic value S ∞ P . The inset shows the classical approximation (TWA) of the survival probability in a semi-log scale, which makes clear that the oscillations are periodic. moves outside the starting region at ω 0 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 ω 0 t = 5.6. (third column). These recurrences are periodic, as confirmed with the inset in the bottom panel of Fig. 7.
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 S P (t) 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), The asymptotic result S ∞ P of the quantum survival probability is more than twice this value. Large S ∞ P , 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 ω 0 t = 400 (fifth column in the top panel of Fig. 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 S P (t).
For longer times, the periodic orbits produce opposing effects in the quantum and classical regimes. Classically, we observe a slow decay of S P (t) towards 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 [28,29]. In the quantum evolution, the scarring decreases the equilibration time. As discussed in Ref. [27] and as can be seen in Fig. 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 separate 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 [28,29]. 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 S P (t).

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 Sec. III A, 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 [22] and Appendix F for details) S (r) where η = 2 √ πσν c and ν c = ν(E c ) is the DoS [see Eq. (A4) in Appendix A] evaluated in the center of the energy profile E c , σ is the width of the LDoS, the b 2 function is the Gaussian orthogonal ensemble (GOE) two-level form factor studied in random matrix theory [85], and D = 2/ν c is the mean level spacing of the correlated eigenvalues. By comparing the quantum survival probability with S  times, S P (t) and S P (t) overlap. There is also perfect agreement between the Gaussian decay of the temporal average of S P (t), S P (t) and S (r) P (t), as expected from the great similarity between the smoothed LDoS at low time resolution (ω 0 T = 0.2) of the coherent and the random state [ Fig. 3 (b) and (c)]. However, the curves for S P (t) and S (r) P (t) diverge after ω 0 t ∼ 0.2. One sees that S P (t) remains on the Gaussian decay, while S The quantum and classical survival probabilities continue their Gaussian decay for ω 0 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 Fig. 8. As the initially localized Wigner distribution at ω 0 t = 0 (first column) moves out of the starting region at ω 0 t = 2 and ω 0 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 ω 0 t ≈ 30 (marked with a vertical line in the large bottom panel of Fig. 8), the classical survival probability attains its asymptotic equilibration value S ∞ P and the Wigner distribution covers the entire energy shell, as seen in the top panels of Fig. 8 for ω 0 t = 31.4 and ω 0 t = 400 (fourth and fifth columns). In contrast, the quantum S P (t) continues to raise, leading to a longer equilibration time. For ω 0 t > 30, the temporal averaged S P (t) coincides again extremely well with the analytical expression for S (r) P (t). During this ramp toward equilibration, the dynamics is controlled by the b 2 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 [43,52,68,86], which justifies the use in Eq. (25) of the same form of the b 2 function from random matrix theory. At this time scale the dynamics becomes universal. The region of the ramp is referred to in Fig. 8 as "GOE correlations" and is commonly known as correlation hole [16-22, 87, 88]. This is a quantum effect associated with the discreteness of the spectrum, which does not appear for the classical S P (t).
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 (r) P (t) (see Ref. [22] and Appendix F), where δ S P is a small parameter indicating that the values of S P (t > t r ) are already within the quantum fluctuations around S ∞ P . This time is marked with a vertical line in Fig. 8. The asymptotic value of the quantum survival probability agrees with the saturation point of S where r and r 2 are the first and second moments of the exponential distribution used to generate the random numbers in Eq. (13), and the ratio 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. One sees that due to this ratio, the value of S ∞ P is twice as large as the asymptotic value of the classical survival probability, which is given by Eq. (24). If the fluctuations of the energy components of the LDoS of the random state were absent and r 2 / r 2 was one, S ∞ P 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 [27,29,89]), which effectively limits the ergodicity of the quantum evolution [90] when compared to the classical dynamics.
We close this discussion with an interesting observation. The minimum value of S (r) P (t), reached right after the Gaussian decay, coincides exactly with S ∞ P . That is, before the ramp towards equilibration, S (r) P (t) stabilizes at the value associated with the ergodicity of the classical evolution. Why the random state is able to reach this greatest level of spreading to later contract to the value S ∞ P of maximal quantum ergodicity is an open question to us.

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 Fig. 8.
In the central panel of Fig. 9, we show the ratio for a mesh of initial coherent states distributed over the same Poincaré surface as in Fig. 1 (c). In the equation above, S (r),∞ P is given by the analytical expression in Eq. (27), 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). There are few ratios that are smaller than one (dark color) and therefore affected by unstable periodic orbits. The ratio R is twice the factor F introduced in [28] to gauge the fraction of phase space explored by the quantum state.
Around the central panel of Fig. 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 Eq. (25) for the random initial state (orange solid line). Panels (a), (b), (d)-(g) make it clear that a ratio R ≈ 1 leads to a generic behavior of the quantum S P (t), well described by S (r) P (t). In all these panels, there appears a ramp towards equilibration associated with the presence of correlated eigenvalues. These random-like coherent states contrast with those where ratio R < 1, which are affected by quantum scarring, as in Fig. 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.

V. CONCLUSIONS
Our comprehensive study of the dynamics and equilibration process of the Dicke model shows that the quantum survival probability of coherent states agrees extremely well with its classical limit obtained with the truncated Wigner approximation up to the point where either the classical or the quantum survival probability equilibrates. It is at this point, at very long-times, that one can identify major aspects of the model, such as those itemized below.
-Features of the equilibration process that are purely quantum and caused by the discreteness of the energy spectrum. One example is the quantum fluctuations around the asymptotic value of the quantum survival probability, which are present in both the regular and chaotic regime. The other is the onset of the correlation hole caused by the correlations between the eigenvalues, which is a universal behavior emerging in the chaotic regime for non-scarred initial states. -Whether the initial state is located in the regular or chaotic regime. The equilibration time of the classical and quantum dynamics coincides in the regular regime, but differ in the chaotic region. In the latter case, correlations between the eigenvalues or the proximity to a periodic orbit separates the curves of the classical and quantum survival probability, one saturating earlier than the other.
-Whether the initial state in the regular regime is close to the separatrix. For initial states located in the border of disconnected classical phase-space regions, the dynamical tunneling between these regions leads to an asymptotic value of the quantum survival probability slightly lower than the classical value, while away from the separatrix the asymptotic values coincide.
-The degree of scarring of the initial coherent states in the chaotic regime. For initial coherent states affected by quantum scarring, the ratio between the asymptotic value of the quantum survival probability and the asymptotic value of the classical survival probability is larger than 2 and the equilibration of the quantum dynamics happens before the classical one. The larger the ratio is, the larger the degree of scarring.
-The onset of maximal quantum ergodicity. In the chaotic regime, for initial coherent states unaffected by quantum scarring, the quantum evolution of the survival probability at long times coincides with that for random states. This analogy allows us to derive analytical results for S ∞ P and for the relaxation time. We find that the ratio between the asymptotic values of the quantum and classical survival probabilities equals 2 and classical equilibrium takes place before the saturation of the quantum dynamics. 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.
This work makes evident the great advantage of having access to the TWA in the phase space, which can be used to explain various features of the quantum evolution of the survival probability. This includes time intervals where the values of S P (t) are close to zero, revivals, and the presence of periodic orbits.
To obtain the representation of coherent states given by Eq. (2) and Eq. (3), we express the coherent state parameters α and z in terms of the canonical variables (q, p, Q, P ), such that α = j 2 (q + ip) and z = 1 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 [68,91] where y ± = −γ −1 γ −1 ∓ 2( − 0 ) andγ = γ/γ c . The ground state energy is given by gs = −1 for the normal phase, and by gs = 0 = − 1 2 γ 2 +γ −2 for the superradiant phase. In Fig. 10, we compare the DoS obtained numerically with Eq. (A4), showing excellent agreement. The figure also shows the energies selected for our studies throughout this paper. Hamiltonian parameters: ω = ω 0 = 1,γ = 2, and j = 100. Vertical black dotted and dashed lines indicate respectively the energies R = −1.8 and C = −0.5 selected for our studies. A truncated Hilbert space was employed using the basis of Refs. [92,93], ensuring 30825 converged eigenenergies, which range from the ground state energy gs = −2.125 to the truncation at T = 0.853.

Appendix B: Wigner Function of Glauber and Bloch Coherent States
As we see from Appendix A, the coherent states for the Dicke model are the product of Glauber and Bloch coherent states [Eq. (A1)]. Thus, the associated Wigner function is the product of the Wigner functions associated to the Glauber and Bloch states. For a Glauber state |α with α(q, p) = j 2 (q + ip), the Wigner function is given by a normal distribution with d = (q − q 0 ) 2 + (p − p 0 ) 2 . For a Bloch coherent state, the Wigner function in variables (θ, φ) may be written as a sum of Legendre polynomials P k (x) [94] where Θ is the angle between (θ, φ) and (θ 0 , φ 0 ) obtained from cos Θ = cos θ cos θ 0 + sin θ sin θ 0 cos(φ − φ 0 ).
For large j values this is very well approximated by a normal distribution on the Bloch sphere 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 Eq. (B1) and Eq. (B4), and it satisfies the one-parameter group identities ϕ 0 = Id, ϕ −t = (ϕ t ) −1 , and ϕ t 1 +t 2 = ϕ t 2 • ϕ t 1 .
Staying constant along classical trajectories means that for any pair of times t 1 and t 2 w(ϕ t 1 (u), t 1 ) = w(ϕ t 2 (u), t 2 ). (C4) In particular, taking t 2 = 0 and performing an adequate change of variables, Inserting Eq. (C5) into Eq. (20), we find Eq. (21). Note that Eq. (21) may be written as Provided that w is everywhere positive, as it is for coherent states [see Eq. (B5)], this expected value may be efficiently approximated using a Monte Carlo method, where the points u i ∈ M are randomly sampled from the initial distribution w, and M is a sufficiently large, albeit computationally accessible, integer.

Appendix D: Intuitive Interpretation of the Classical Survival Probability
An intuitive understanding of the classical survival probability may be obtained by considering a simple measurable set S ⊆ M (a sphere, for example). Let 1 S be its indicator function 4 . Consider the normalized distribution ρ S = 1 S /V S in lieu of the Wigner function. Here V S = du 1 S (u) is the volume of S. In this case, to retain normalization, our prefactor becomes 1/V S instead of (2π/j) 2 , and then, from Eq. (21), Since 1 S (ϕ −t (u)) = 1 ϕ t (S) (u) and using that the product of the indicator functions of two sets is the indicator function of their intersection, we get This result is easy to interpret: S P (t) is the percentage of the volume of S that is found back inside S after time t. In other words, it is the probability that a point sampled from S is back in S after time t [25,96].

Appendix E: Asymptotic Value of the Classical Survival Probability
To obtain an expression for S ∞ P , 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 u ∈ M with energy E = H cl (u), where j −2 (2π) 2 ν(E) is the volume of the energy shell for E/j in M, and ν(E) is given by Eq. (A4). It is then straightforward to calculate The last equality is obtained by using 1 = ∞ Egs dE δ(E − H cl (u)) inside the integral, a change in the integration order, and a substitution of the value of w E . In the above, E gs represents the ground state energy. The classical energy distribution for the state associated to w is By further approximating ρ cl with the Gaussian distribution given by Eq. (11), we obtain Eq. (24).
Appendix F: Analytical Expression of the Survival Probability for a Random Ensemble with a Gaussian Energy Profile An analytical expression for the survival probability averaged over an ensemble of random initial states was derived in [22]. 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 |c (r) where ν(E) is the DoS, the numbers r k are randomly generated from a given probability distribution P (r), and A = q r q ρ(E q )/ν(E q ) is a normalization constant. This gives, S (r) where the initial decay of the survival probability is dictated by the effective dimension of the ensemble where ν c = ν(E c ) is the density of states evaluated in the center of the energy profile, and r n are the n-th moments of the distribution P (r), and the asymptotic value of the survival probability corresponds to S (r),∞ P = r 2 r 2 1 η .
The function b 2 is the two-level form factor of the GOE [85] b 2 (t) =[1 − 2t +t ln(2t + 1)]Θ(1 −t)+ where Θ is the Heaviside step function. The factor D in the argument of b 2 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 Eq. (11)] and the random numbers r k are generated from the exponential distribution P (r) = λe −λr with r n = n!/λ n in Fig. 3, which implies that S bc P (t) = e −σ 2 t 2 , 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 (r),∞ P , η, S bc P , and D, we get Eq. (25).
To derive the relaxation time t r of the ensemble-averaged survival probability, we consider the asymptotic form of b 2 , which grows toward saturation following a power-law behavior b 2 t πν c → π 2 ν 2 c 12t 2 for t πν c 1.

(F10)
At this temporal scale, the contribution of the initial decay S bc P is negligible and the asymptotic form of (F2) is given by where in the last step we have used η = 2/S (r),∞ P . We define the relaxation time according to S (r) where δ S P is a small parameter determining the point where S