Chaos and Thermalization in the Spin-Boson Dicke Model

We present a detailed analysis of the connection between chaos and the onset of thermalization in the spin-boson Dicke model. This system has a well-defined classical limit with two degrees of freedom, and it presents both regular and chaotic regions. Our studies of the eigenstate expectation values and the distributions of the off-diagonal elements of the number of photons and the number of excited atoms validate the diagonal and off-diagonal eigenstate thermalization hypothesis (ETH) in the chaotic region, thus ensuring thermalization. The validity of the ETH reflects the chaotic structure of the eigenstates, which we corroborate using the von Neumann entanglement entropy and the Shannon entropy. Our results for the Shannon entropy also make evident the advantages of the so-called “efficient basis” over the widespread employed Fock basis when investigating the unbounded spectrum of the Dicke model. The efficient basis gives us access to a larger number of converged states than what can be reached with the Fock basis.


Introduction
The onset of thermalization in isolated quantum systems, which evolve unitarily and are described by pure states, was discussed in von Neumann's 1929 paper on the quantum ergodic theorem [1][2][3][4]. In Pechukas's 1984 paper "Remarks on Quantum Chaos" [5], one finds a brief summary of von Neumann's work and the subsequent developments against it, particularly by Bocchieri and Loinger [6], which, Pechukas claims, actually make von Neumann's results sharper [7]. In Pechukas's words, if one selects a state "at random" from an energy shell and determines, as a function of time, the probability that it lies in a "typical" subspace of the shell, the time average of this probability is liable to be much closer to the statistical expectation than is its instantaneous value, the more so the less degenerate the spectrum of the Hamiltonian.
The ideas in this quotation are connected with the notion of "typicality" [3,4], which has been one of the directions in studies of thermalization.
In retrospect, one can identify in the 1985 work by Jensen and Shakar [8] the seeds for what later became known as the eigenstate thermalization hypothesis (ETH) [9,10]. There, the authors study the infinite-time average and the eigenstate expectation value of the magnetization of a finite spin-1/2 chain for different initial states and different energy eigenstates, and show that both exhibit very good agreement with the microcanonical average when the system is chaotic. They state that in the chaotic system, "the magnetization is a fairly smooth and monotonic function of energy", and in this case, even an initial energy eigenstate will exhibit a constant value for the observable which is very close to that predicted by the statistical theory, while for the integrable Hamiltonian, the "deviations from equilibrium tend to be large". This means that an observable O should thermalize when its expectation value E k |Ô|E k is a smooth function of energy, since in this situation, its infinite-time average, O, will be very close to the microcanonical average, O mic ; that is [11] where |E k are the eigenstates of the system's Hamiltonian, O k,k = E k |Ô|E k are the diagonal elements of the operatorÔ, c k = E k |Ψ(0) are the coefficients of the initial state |Ψ(0) = ∑ k c k |E k , and W E,∆E is the number of eigenstates |E k contained in the energy window E k ∈ [E − ∆E, E + ∆E] with |E − E k | < ∆E. Under the assumption of the smooth behavior in energy, even a single eigenstate inside the microcanonical window should give O k,k very close to the microcanonical average; therefore, the term "eigenstate thermalization hypothesis" (ETH) was coined by Srednicki in his 1994 paper [12]. Notions of quantum chaos [8], random matrices [13], and Berry's conjecture [12] have been invoked to justify the validity of Equation (1). Starting with Rigol et al's 2008 paper [11], several numerical studies have confirmed Equation (1) for chaotic systems. Since then, various studies further elaborated the framework of the ETH to take into account the convergence of O and O mic as the system size increases [9,14], the analysis of the off-diagonal elements O k,k = E k |Ô|E k , known as off-diagonal ETH [9,15], the dependence on the energy of the initial state [16,17], and the structure of the eigenstates in realistic many-body quantum systems [18][19][20]. The analysis in Ref. [18] was inspired by the 1990s papers from previous members of the Novosibirsk school, including Casati [21], Flambaum [22,23], Izrailev [24,25], Shepelyansky [26], and Zelevinsky [27,28], who focused on the structure of the eigenstates to define quantum chaos and explain the onset of thermalization [29]. In realistic interacting many-body quantum systems, the energy eigenstates can be rather complicated, but they are not random vectors, as the eigenstates of Gaussian random matrices or the random superpositions of plane waves with random phases and Gaussian random amplitude stated by the Berry's conjecture [30]. For such random vectors, Equation (1) is trivially satisfied, but they are associated with unphysical models. In the above mentioned 1990s papers, chaotic eigenstates are those that, when written in the mean field basis, exhibit coefficients that are random variables following a Gaussian distribution around the envelope defined by the energy shell. They emerge away from the edges of the spectrum of quantum systems with many strongly interacting particles and ensure the validity of the Fermi-Dirac [24] and Bose-Einstein [31] distributions. These ideas were employed in Ref. [18] to corroborate the validity of the ETH for few-body observables when the eigenstates are chaotic.
The present article is dedicated to Professor Giulio Casati on the occasion of his 80th birthday in 2022 and to Professor Felix Izrailev's 80th birthday in 2021. Not only their past, but also their recent works [31][32][33][34][35][36][37][38] continue to have an enormous impact in the developments of quantum chaos and its connections with thermalization, quantum statistical mechanics, and the quantum-classical correspondence.
Here, we investigate the validity of the diagonal and off-diagonal ETH, in connection with the structure of the eigenstates, for the spin-boson Dicke model [39][40][41]. This system has a well-defined classical limit and two-degrees of freedom, so it does not quite fall within the requirement for thermalization of a large number of coupled degrees of freedom [13].
We analyze the diagonal and off-diagonal ETH for the number of photons and the number of excited atoms in the regular and chaotic regions of the Dicke model, and compare the results with the structure of the eigenstates analyzed with the entanglement entropy and the Shannon entropy. The latter is a delocalization measure that depends on the basis in which the eigenstates are written. The Hilbert space of the Dicke model is infinite, because its number of bosons is unbounded. Our results show that the values of the Shannon entropy grows rapidly with the eigenvalues when the eigenstates are written in the Fock basis, but are more restricted when the "efficient basis" is used, which makes evident the advantages of the latter when one wants to study large systems and high energies.
The paper is organized as follows. We introduce the Dicke model in Section 2 and briefly review the onset of classical and quantum chaos in Section 3. In Section 4, we study the diagonal and off-diagonal ETH and analyze the eigenstates in Section 5. Our conclusions are summarized in Section 6.

Dicke Model
The Dicke model [39] describes a system of N two-level atoms coupled with a single mode of a quantized radiation field. Settingh = 1, the Hamiltonian is given bŷ whereĤ F defines the field's energy,Ĥ A the energy of the two-level atoms, andĤ I the atomfield interaction energy. In the equations above,â † (â) is the bosonic creation (annihilation) operator of the single field mode andĴ + (Ĵ − ) is the raising (lowering) collective pseudospin operator, whereĴ ± =Ĵ x ± iĴ y andĴ x,y,z = (1/2) ∑ N k=1σ k x,y,z for the Pauli matriceŝ σ k x,y,z acting on the kth two-level atom. Since the squared total pseudo-spin operator given by j(j + 1), label the invariant subspaces in the Hilbert space. We work within the totally symmetric subspace, which includes the collective ground state and is defined by the maximum pseudo-spin value j = N /2. The Dicke HamiltonianĤ D commutes with the parity operator whereΛ =â †â +Ĵ z + j1 =n +n ex , the number of photons isn =â †â , and the number of excited atoms isn ex =Ĵ z + j1. Because of this symmetry, the eigenstates have one of two different parities,Π|E k = ±|E k . The three parameters of the Dicke Hamiltonian, ω, ω 0 , and γ, determine the energy scales of the system and the onset of chaos. The radiation frequency of the single-mode electromagnetic field (the boson) is given by ω, the energy splitting of each two-level atom is ω 0 , and γ is the coupling strength modulating the atom-field interaction. In the thermodynamic limit, N → ∞, the light-matter coupling divides the parameter space in the normal phase, when the strength γ in smaller than the critical value γ c = √ ωω 0 /2, and the superradiant phase, when γ > γ c [42,[72][73][74]. The Dicke model also presents regular and chaotic regions depending on the parameters and excitation energies [47]. We fix the coupling strength in the superradiant phase, γ = 2γ c = 1, and the resonant frequencies at ω = ω 0 = 1, so that we ensure that both the classical dynamics and the quantum eigenspectrum are chaotic at excitation energies ≥ −0.8 [47], where = E/j denotes the energy scaled to the system size j.
The Hilbert space of the Dicke model is infinite. In Appendix A, we provide details on how to diagonalize the system's Hamiltonian on an efficient basis. We consider three system sizes, j = 30, 60, 100, whose associated Hilbert space dimensions in the efficient basis are given by d EB D = 8 601, 30 371, 80 601, ensuring the convergence of the set of eigenstates from the ground-state energy GS = −2.125 until a truncated value T .

Chaos in the Dicke Model
The classical Hamiltonian of the Dicke model is presented in Appendix B. The generation of Poincaré sections and the computation of Lyapunov exponents for trajectories in phase space [47,75] reveal the onset of chaos for high energies and strong couplings. In Figure 1a, we show a classical map of the percentage of chaos as a function of excitation energies and coupling strengths. The percentage of chaos is defined as the ratio of the number of chaotic initial conditions over the total number of initial conditions used in the sample and is illustrated with a color gradient, where dark indicates that the majority of the initial conditions are regular and light indicates that most initial conditions are chaotic. Signatures of the onset of chaos in classical systems are found in the quantum domain, as suggested by Casati et al's 1980 paper [76] and Bohigas et al's 1984 work [77], and became known as "quantum chaos". The eigenvalues of a quantum system become correlated when its classical counterpart is chaotic. The degree of correlations between neighboring levels can be detected with the level spacing distribution [78] or the ratio of consecutive energy levels [79,80], where s k = E k+1 − E k is the nearest-neighbor spacing between energy levels and r k = s k /s k−1 ; while both short-and long-range correlations are measured with quantities such as the level number variance [81] or the spectral form factor [78]. In the case of the Dicke model, level spacing distributions [75,82,83] and spectral form factors [53,54] in agreement with random matrix theory (RMT) were verified for the energies and parameters associated with the onset of chaos in the classical limit.
Other tests of quantum chaos include Peres lattices [84] and the evolution of OTOCs [85][86][87]. The Peres lattice is a plot of the eigenstate expectation values of an observable as a function of the eigenvalues. It has been widely employed (not always under this name) in studies of ETH and provides visual evidence of the loss of integrability [75,88]. As for the OTOCs, their exponential growth rates are seen as quantum analogs of the classical Lyapunov exponents. This was indeed confirmed for the Dicke model in the chaotic region [55], although its exponential growth happens also in regular regions due to instability [57].
In Figure 1, we compare the onset of classical chaos, investigated in Figure 1a, with the degree of correlations between neighboring energy levels, as captured by the ratio of consecutive energies in Figure 1b. To smooth spectral fluctuations, an average is performed and denoted by r . In the regular regime, where the eigenvalues are uncorrelated and the level spacing distribution is Poissonian, r P ≈ 0.39, while in the chaotic region, where the eigenvalues are correlated as in RMT and the level spacing distribution follows the Wigner-Dyson distribution, r WD ≈ 0.53. Our results in Figure 1 exhibit an evident classical-quantum correspondence. There is a visible relationship between the classical map of Figure

Thermalization in the Dicke Model
For an isolated quantum system initially in a pure state, the evolution is governed by the unitary time operator as whereĤ|E k = E k |E k and c k = E k |Ψ(0) . Thus, the expectation value of a given operator O under states evolved in time can be calculated as follows where O k,k = E k |Ô|E k are the matrix elements of the operator expressed in the energy eigenbasis. The onset of thermalization according to the ETH happens for local observables if two assumptions are satisfied: (i) The infinite-time average, coincides with the microcanonical average, or more precisely, O approaches the thermodynamic average as the system size grows. This is sometimes referred to as "diagonal ETH". In the equation above, W E,∆E is the number of eigenstates |E k contained in the energy window As explained above, the ETH should be valid when the eigenstates are chaotic (ergodic), filling the energy shell. In this case, one obtains Gaussian distributions for few-body observables in many-body quantum systems, which can be understood as follows. Assume that the few-body observableÔ has N nonzero eigenvalues O n inÔ|O n = O n |O n , where N is large. We can project the energy eigenstates in the basis |O n as |E k = ∑ n C k O n |O n , and write If the eigenstates are fully chaotic (ergodic), then C k O n 's are independent Gaussian random numbers. According to the central limit theorem, a large sum of independent random numbers O n,n |C k O n | 2 will also follow a Gaussian distribution. Therefore, in the region where the eigenstates are chaotic, the distribution of O k,k should be Gaussian.
(ii) The fluctuations around O, which are determined by the phases e i(E k −E k )t , the coefficients c k , and the off-diagonal elements O k,k in Equation (10), decrease with system size and cancel out on average. This is sometimes referred to as "off-diagonal ETH". A Gaussian distribution of O k,k indicates that the eigenstates are strongly chaotic (ergodic), so condition (ii) should be satisfied if the energy of the initial state (of a system perturbed far from equilibrium) is in the chaotic region.
The relationship between the Gaussian distribution of the off-diagonal elements of a few-body observable and ergodicity can be understood as above, usingÔ|O n = O n |O n and |E k = ∑ n C k O n |O n in For chaotic eigenstates, C k O n 's are independent Gaussian random numbers, and the product of two independent Gaussian random numbers, (C k O n ) * C k O n , which appears in each term of Equation (14), is also an independent random number. According to the central limit theorem, a large sum of independent random numbers follows a Gaussian distribution, so in the region where the eigenstates are chaotic, the distribution of O k,k should be Gaussian. This has been confirmed for different chaotic systems with many degrees of freedom [15,89,90], but not in chaotic systems with one [91] or few particles [92], few degrees of freedom [93] or in many-body systems whenÔ is not few-body [94].

Diagonal ETH
To test the validity of the diagonal ETH, we compute the deviation of the eigenstate expectation values with respect to the microcanonical value [19,20] and also use a stronger test that takes into account the normalized extremal fluctuations and is given by [19] ∆ mic where max(O) and min(O) are taken from the same energy window E k ∈ [E − ∆E, E + ∆E] used in Equation (12). We consider as observables the number of excited atoms,n ex = J z + j1, and the number of photons,n =â †â , of the Dicke model. In Figure 2, we show the Peres lattices for the expectation values of the number of photons, n k,k = E k |n|E k (Figure 2a), and of the number of excited atoms, (n ex ) k,k = E k |n ex |E k (Figure 2b), for all eigenstates of the Dicke model with positive parity that range from the ground state energy, GS = −2.125, up to a maximal converged eigenstate with eigenenergy T = 1.755. In the regular region of low energies, the Peres lattices present a clear pattern related with the quasi-conserved quantities. Above ≈ −0.8, where the system becomes chaotic, the lattices become smoother in energy. It is visible that the spread of the expectation values in the chaotic region decreases as the system size increases; that is, for high excitation energies, the fluctuations are larger for j = 30 than for j = 100. It is also noticeable that the number of excited atoms fluctuate less than the number of photons. The reduction of the fluctuations with the increase of system size, which is a necessary condition for the validity of the ETH, is better quantified in Figure 2c,d, where we show ∆ mic [n] and ∆ mic [n ex ], respectively. To compute these quantities, we used moving windows of eigenenergies for each system size j = 30, 60, 100, which contain 100, 350, 900 energy levels. In the chaotic region, both deviations clearly decrease as j increases, while at low energies ( < −0.8), where chaos and regularity coexist, the fluctuations decrease very slowly or do not decrease at all.
The insets of Figure 2c,d contain the results for the normalized extremal fluctuation in Equation (16). Only the chaotic energy interval defined by ∈ [−0.5, 1.755] is shown. The reduction of the extremal fluctuations as j increases is also perceptible in these insets, thus confirming the validity of the diagonal ETH in the chaotic region of the Dicke model.
In Figure 3, we consider the energy interval ∈ [0.5, 1], where hard chaos manifests itself classically, and use only the positive parity sector of eigenstates. In Figure 3a,b, we show, respectively, the distribution of the diagonal elements of the number of photons and the number of excited atom, which are presented in a semi-logarithmic scale. The figures confirm that the distribution shape is Gaussian, as explained in Equation (13).

Off-Diagonal ETH
We now analyze the off-diagonal elements of the number of photons in Figure 3c and the number of excited atoms in Figure 3d. As discussed in Equation (14), a Gaussian distribution validates the off-diagonal ETH. This is indeed the case of Figure 3c,d.

Entropies of the Eigenstates of the Dicke Model
In the last section, we verified numerically the validity of the ETH in the chaotic region of the Dicke model using the matrix elements of the number of photons and the number of excited atoms. In this section, we use the von Neumann entanglement entropy and the Shannon entropy to analyze the structure of the eigenstates of the model. The von Neumann entanglement entropy can be regarded as the limit of higher-order observables in a replicated Hilbert space, and one may generalize the ETH to include higher-order statistical moments [95] arising from these replicated spaces, which allows one to explain the so-called Page correction [96] to the volume law. The von Neumann entanglement entropy has been linked to the onset of chaos in quantum systems [97][98][99][100][101]. The Shannon entropy is a basis-dependent quantity. We study this entropy with respect to both the Fock and an efficient basis.
The Dicke model is a bipartite system, whose Hilbert space is a tensor product of the atomic H A and bosonic H B sectors, H D = H A ⊗ H B . The bosonic sector is infinitedimensional, while the atomic one has dimension d A = 2j + 1. For a pure state expanded in the Fock basis |n; j, m z , the density matrix is the projector operatorρ = |Ψ Ψ|, and the reduced density matrix in the atomic sector is calculated as, c n,m z c * n,m z |j, m z j, m z |.
The von Neumann entanglement entropy is given by For numerical convenience, we trace out the infinite bosonic sector first, but from the Schmidt decomposition [102], the result for S En is the same if we would instead trace out the atomic sector first. (See Appendix C for a generalization of quantum entanglement to multipartite systems.) The standard basis to compute the entanglement entropy is the Fock basis, which is defined as a decoupled basis between photons and atoms |n; j, m z = |n ⊗ |j, m z (see Appendix A.1). Nevertheless, alternative bases can be used to calculate this entropy by selecting the correct partition of the system. In this work, in addition to the Fock basis, we use an efficient basis (see Appendix A.2 for a full description) that allows us to reach larger system sizes (j > 30) than we can reach with the Fock basis [103,104]. The efficient basis is defined as the tensor product |N; j, m x = |N m x ⊗ |j, m x , where |j, m x are the atomic pseudo-spin states rotated by an angle −π/2 around the y axis and |N m x are displaced Fock states, whose displacing depends explicitly on the atomic pseudo-spin eigenvalue m x . In the following general notation for a pure state, (x, y) = (n, m z ) stands for the Fock basis and (x, y) = (N, m x ) for the efficient basis.
To select the correct atomic sector starting with the efficient basis, we need to perform an adequate trace over the modified bosonic sector related with this basis, such that it is equivalent to the trace over the bosonic sector of the Fock basis. This can be accomplished by mapping one basis into the other (see Appendix A.3). As the Hilbert space associated with the atomic pseudo-spin states |j, m z is the same as the Hilbert space of the rotated ones |j, m x , the atomic sector for both subspaces is the same and the entanglement entropy can be computed properly. The Shannon entropy of the pure state in Equation (20) is given by

Results for the Entanglement Entropy
We start by comparing the von Neumann entanglement entropy for the Dicke model, which is nonintegrable, with the results for one of the integrable limits of the model, namely the Tavis-Cummings model (see Appendix D for the derivation of this model). This comparison was also made in Ref. [105].
In Figure 4, we plot the Peres lattice of the exponential of the von Neumann entanglement entropy in Equation (19) for the integrable Tavis-Cummings model (Figure 4a) and the nonintegrable Dicke model (Figure 4b). Results for all states from both parity sectors, from the ground-state of the Tavis-Cummings model up to a maximal converged eigenstate with eigenenergies k ∈ [ GS = −2.136, T = 2.5], are shown. As evident in Figure 4a, the values of the entropy for the integrable case show large fluctuations and patterns associated with regularity are visible. The fluctuations indicate that even states very close in energy may have very different structures, so ETH should not be satisfied. In contrast, S En becomes a smoother function of energy in the chaotic region of the Dicke model ( > −0.8), as seen in Figure 4b, which indicates the states close in energy are very similar. In the case of the Dicke model, regular patterns are restricted to the low energies, where the model is not chaotic.
We stress, however, that some isolated eigenstates located in the chaotic regime of the Dicke model present low entanglement. They should be related with strongly scarred states, which are known to exist in this model [61,62,106].  (19)) scaled to the atomic Hilbert-space dimension d A = 2j + 1 for eigenstates |E k from both parity sectors of the integrable Tavis-Cummings model (a) (see Equation (A23)) and the nonintegrable Dicke model (b) (see Equation (2)). The system size in both panels (a) and (b) is j = 30.
In Figure 5a, we analyze the von Neumann entanglement entropy for the eigenstates of the Dicke model, ranging from the ground state until a maximal converged eigenstate with eigenenergies k ∈ [ GS = −2.125, T = 1.841], for three values of the system size j = 30, 60, 100. As discussed in Section 4 and in Figure 2, thermalization requires the convergence of the infinite-time average of a few-body observable towards the microcanonical average in the thermodynamic limit, so scaling analysis needs to be performed. As seen in Figure 5a, the patterns at the low energies of the regular region of the Dicke model do not disappear as the system size increases, but in the chaotic region, the fluctuations clearly shrink as the system size increases. This behavior is quantified in Figure 5b, where we present the deviation of the entanglement entropy from the microcanonical average, as computed in Equation (15) for the three system sizes j = 30, 60, 100. The fluctuations decay to values close to zero for energies above ≈ −1.2, indicating the transition to a region where the majority of the eigenstates are chaotic. In the inset of Figure 5b, we present the extremal fluctuations calculated with Equation (16) for the chaotic energy interval only, ∈ [−0.5, 1.841]. We avoid the regular region, where the changes in the extremal values are abrupt. The inset confirms the decrease in fluctuations with the increase in system size.
By comparing Figure 5b with Figure 3c,d, we observe that the deviation from the microcanonical average of the von Neumann entanglement entropy decreases to zero more abruptly as compared to the the number of photons and excited atoms. This different behavior between entropies and expectation values of observables is not well understood yet and motivates future work.

Results for the Shannon Entropy
We proceed with the analysis of the structure of the eigenstates of the Dicke model making use now of the basis-dependent Shannon entropy. In Figure 6a, we plot the Peres lattice of the Shannon entropy for the eigenstates of the Dicke model in the Fock and efficient bases. The eigenstates range from the ground state until a maximal converged eigenstate with eigenenergies k ∈ [ GS = −2.125, T = 2.356]. The inset in Figure 6a contains the same data, but shows the exponential of the entropy. The values of the Shannon entropy computed in the Fock basis are larger than in the efficient basis. For the Fock basis, the entropy grows unboundedly with energy, while for the efficient basis, S Sh saturates at high energies. These features make evident the advantages of using the efficient basis, since it requires fewer basis states to build a given eigenstate than what is needed by the Fock basis. While in Figure 6, our system size was restricted to j = 30, when we employ the efficient basis, the same computational resources allow us to go up to j = 100, as used in Figures 2, 5 and 7 below.  (21)) scaled to the systemsize dependence of the density of states, ν( ) ∝ 2j 2 (see Appendix B), for eigenstates |E k from both parity sectors of the Dicke model written in the efficient basis (red dots) and the Fock basis (purple dots). The inset shows the exponential values of the Shannon entropy for both bases. Panel (b): Deviation of Shannon entropy with respect to its microcanonical value (see Equation (15)). The inset shows the extremal deviation of the same quantity in the chaotic energy regime (see Equation (16)). The system size in both panels (a) and (b) is j = 30.
The analysis of the fluctuations of the Shannon entropy in Figure 6b and its inset shows that, similarly to what we observed for the entanglement entropy in Figure 5, they approach zero in the chaotic region. The results are comparable for both the Fock and efficient bases. In what follows, we examine the fluctuations of the Shannon entropy for different system sizes obtained for eigenstates written in the efficient basis only.
In Figure 7a, we plot the Peres lattice of the Shannon entropy for the eigenstates of the Dicke model written in the efficient basis. These eigenstates are contained in the energy interval k ∈ [ GS = −2.125, T = 1.841] for the same values of the system size which were considered previously, j = 30, 60, 100. Similarly to what we observe for the entanglement entropy in Figure 5a, the regular patterns and large fluctuations are visible at low energies, where we do not expect the validity of ETH. The fluctuations decrease as we move to high energies and as the system size is increased. This behavior is quantified with the deviations from the microcanonical average shown in Figure 7b and its inset. The decay of ∆ mic [S En ] to values close to zero for the entanglement in Figure 5b is more abrupt than what we see for ∆ mic [S Sh ] in Figure 7b, but both cases signal the transition to chaos. We close this section mentioning that the regular patterns seen at low energies in Figures 2 and 5-7 reflect the existence of families of periodic orbits, which were studied in previous works [61,62,106]. We also note that for high energies, the behavior of the Shannon entropy computed in the efficient basis and of the von Neumann entanglement entropy are similar in the sense that both measures nearly saturate. In contrast, the Shannon entropy in the Fock basis grows rapidly with energy, because of the infinite bosonic Hilbert subspace of the Dicke model. The modified bosonic sector of the efficient basis is responsible for the smaller values of the Shannon entropy when compared to the values for the Fock basis.

Conclusions
In this work, we confirmed the validity of the ETH in the chaotic region of the Dicke model. This was done by analyzing the diagonal and off-diagonal elements of the number of photons and the number of excited atoms for different system sizes. We corroborated that the validity of the ETH stems from the presence of chaotic eigenstates, which we showed by analyzing their components and measures of entanglement and delocalization.
The Shannon entropy, used to quantify the level of delocalization of the eigenstates in a given basis, made evident the advantages of using the efficient basis over the Fock basis. For high energies, the first leads to a slower growth of the entropy than the Fock basis, allowing us to reach converged states for larger system sizes than accessible with the Fock basis. Institutional Review Board Statement: Not applicable.

Data Availability Statement:
All the data that support the results and plots showed within this work are available from the corresponding authors upon request.

Acknowledgments:
We acknowledge the support of the Computation Center-ICN, in particular to Enrique Palacios, Luciano Díaz, and Eduardo Murrieta.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: ETH eigenstate thermalization hypothesis RMT random matrix theory OTOC out-of-time-ordered correlator

Appendix A. Diagonalization Bases
Appendix A. 1

. Fock Basis
The Fock basis is the natural choice to diagonalize the Dicke Hamiltonian, since it is composed by the tensor product between bosonic Fock states |n and atomic pseudo-spin states |j, m z , |n; j, m z = |n ⊗ |j, m z . On the one hand, the eigenvalues n = 0, 1, 2, . . . of the number operatorn =â †â (with eigenvalue equationn|n = n|n ) provide an infinite bosonic Hilbert subspace. On the other hand, the eigenvalues m z = −j, −j + 1, . . . , j − 1, j of the pseudo-spin operatorĴ z (with eigenvalue equationĴ z |j, m z = m z |j, m z ) provide a finite atomic Hilbert subspace with dimension d FB A = 2j + 1. In order to diagonalize the Dicke Hamiltonian a truncation value n max of the bosonic Hilbert subspace has to be chosen, which allows to define a finite bosonic dimension d FB B = n max + 1. In this way, the global Hilbert-space dimension of the Dicke model is given by the product d FB

Appendix A.2. Efficient Basis
The truncation required by the Fock basis for the convergence of the high-energy eigenstates rapidly increases with the system size j. This can be circumvented by using the so-called efficient basis, which was originally obtained by studying the integrable limit ω 0 → 0 [107,108]. This basis can be written in terms of a displaced bosonic annihilation operatorÂ =â + GĴ x with G = 2γ/(ω √ N ) and the eigenstates |j, m x ofĴ x (Ĵ x |j, m x = m x |j, m x ). The efficient basis states are defined by where |α m x = D(α m x )|0 is a coherent state centered at α m x = α * m x = −Gm x , andD(α m x ) = exp(α m xâ † − α * m xâ ) is the displacement operator. By commutingÂ † with the displacement operator and using that GĴ x |j, m x = −α m x |j, m x , one may further write where |N m x =D(α m x )|N are displaced Fock states, also called generalized coherent states. The eigenvalues N = 0, 1, 2, . . . of the displaced number operatorN =Â †Â (with eigenvalue equationN|N m x = N|N m x ) label a modified bosonic Hilbert subspace. The eigenvalues m x = −j, −j + 1, . . . , j − 1, j of the pseudo-spin operatorĴ x label the same finite atomic Hilbert sector with dimension d EB A = 2j + 1. As for the Fock basis, the modified bosonic sector must be truncated to some N max for diagonalization. This yields a finite modified bosonic dimension d EB B = N max + 1 and a global Hilbert-space dimension In contrast to the Fock basis, where convergence of high-energy eigenstates is infeasible for large system sizes (usually j > 30), the efficient basis allows to get thousands of converged eigenstates in high-energy regimes, even beyond j = 100 [103,104].

Appendix A.3. Mapping from Efficient Basis to Fock Basis
To compute the entanglement entropies shown in Figures 6 and 7, we perform a partial trace over the usual atomic sector of the Fock basis, but we diagonalize in the efficient basis. Thus, we need to map the wave function of the eigenstates from the efficient basis to the Fock basis. This is done as follows.
A general pure quantum state |Ψ can be expanded in the efficient basis and the rotated Fock basis (|n; j, m x ), respectively where x = N defines the efficient basis and x = n the Fock basis. Moreover, C x,m x = x; j, m x |Ψ are the coefficients of the arbitrary state in each basis, which must satisfy the normalization condition ∑ x,m x |C x,m x | 2 =1. Note that C n,m x = n; j, m x |Ψ using that j, m x |j, m x = δ m x ,m x . The term n|D(α m x )|N for n > N is given by [109][110][111] n|D(α m x )|N = N! n! α n−N m x e −|α mx | 2 /2 L n−N N (|α m x | 2 ), where L n 1 n 0 (x) is an associated Laguerre polynomial given by the Rodrigues formula L n 1 n 0 (x) = x −n 1 e x n 0 ! d n 0 dx n 0 (e −x x n 0 +n 1 ).
We numerically found that, in order to ensure a correct convergence of the coefficients C n,m x given by Equation (A4), the truncation value must be chosen as n max ≈ 3N max . The associated Laguerre polynomials can be efficiently calculated to arbitrary precision with a package included in the Wolfram Mathematica software [112].

Appendix B. Classical Limit of the Dicke Model
The classical limit of the Dicke model can be obtained taking the expectation value of the quantum HamiltonianĤ D under the tensor product of Glauber and Bloch coherent states |x = |q, p ⊗ |Q, P , and dividing it by the system size j [47,54,58,75,82,113,114] where h F (x) and h A (x) represent the Hamiltonians of two classical harmonic oscillators, and h I (x) the coupling between them. The bosonic Glauber and the atomic Bloch coherent states, represented by the canonical variables (q, p) and (Q, P) respectively, are given explicitly by |q, p = e −(j/4)(q 2 +p 2 ) e where |0 is the photon vacuum and |j, −j is the state with all the atoms in the ground state. The classical Dicke Hamiltonian h D (x), obtained with the latter method, has an infinite four-dimensional phase space M in the canonical variables x = (q, p; Q, P), where the atomic variables are bounded (Q 2 + P 2 ≤ 4). A useful property of this phase space is that it can be partitioned into a family of classical energy shells with finite volume V ( ) < ∞, given by where = E/j is the classical energy of the shell scaled to the system size j, which defines an effective Planck constanth eff = 1/j [115]. The finite volume V ( ) of the classical energy shells M( ) is obtained with a semiclassical approximation to the quantum density of states ν( ), using the Gutzwiller trace formula [114,116,117]. The explicit expression is given by where the density of states is proportional to the system size ν( ) ∝ 2j 2 , and can be derived explicitly following Ref. [114].
such that, the complete Tavis-Cummings Hamiltonian is given bŷ where the HamiltoniansĤ F andĤ A are the same terms given in Equations (3) and (4). The particularity of the Tavis-Cummings Hamiltonian is that it commutes with the operator Λ (see Equation (7)), which is a conserved quantity that defines the number of excitations. The last feature allows to diagonalize the Tavis-Cummings Hamiltonian in finite subspaces of such operator [119].