The Critical Point Entanglement and Chaos in the Dicke Model

Ground state properties and level statistics of the Dicke model for a finite number of atoms are investigated based on a progressive diagonalization scheme (PDS). Particle number statistics, the entanglement measure and the Shannon information entropy at the resonance point in cases with a finite number of atoms as functions of the coupling parameter are calculated. It is shown that the entanglement measure defined in terms of the normalized von Neumann entropy of the reduced density matrix of the atoms reaches its maximum value at the critical point of the quantum phase transition where the system is most chaotic. Noticeable change in the Shannon information entropy near or at the critical point of the quantum phase transition is also observed. In addition, the quantum phase transition may be observed not only in the ground state mean photon number and the ground state atomic inversion as shown previously, but also in fluctuations of these two quantities in the ground state, especially in the atomic inversion fluctuation.


Introduction
It is well known that the Dicke model [1] exhibits a "superradiant" quantum phase transition (QPT) in the thermodynamic limit as shown by Hepp and Lieb [2].To date, many aspects of the model have been studied, such as the dynamical properties of quantum entanglement and decoherence [3][4][5][6][7][8], quantum phase transitions (QPT) [2,[9][10][11] and chaotic properties [3,5,9,10,[12][13][14][15].Furuya et al. [3] did the initial work on the entanglement process in a varied version of the model and concluded that for short times, there is a faster increase in decoherence for chaotic initial conditions as compared to regular ones, which have an oscillatory increase.Hou and Hu studied this problem further in [5], concluding that the entanglement makes a distinct change at the QPT and that for regular initial conditions, the entanglement is smaller than for chaotic ones in the case of weak coupling and fluctuates with a small amplitude for strong coupling and chaotic initial conditions.They paid special attention to the relation between chaos and the QPT, e.g., in [10], they show that the system undergoes a transition from quasi-integrable to quantum chaotic and that the transition is linked to precursors of the QPT.The relation between entanglement and the QPT has also been investigated by others.For example, the entanglement and the QPT in the rotating-wave approximation (RWA) is reported in [6], which suggests that the entanglement only occurs in the strong coupling limit.The model has also been extensively studied without the RWA, with its quantum-chaotic properties discussed by several other authors [12][13][14][15], without reference to the QPT.
In the thermodynamic limit, the ground state in the model is exactly solvable in the whole coupling range, as shown in [9,10] based on the Holstein-Primakoff realization of the su(2) Lie algebra , which only provides a method to get exact results for the ground state, but not for the whole spectrum.However, the level statistics only becomes meaningful when a sufficient number of excited states in the model is known.On the other hand, for finite N cases, the system is non-integrable in general [16], with the finite-size effect shown to be crucial in understanding the properties of entanglement [5,7,8,[17][18][19].As has been shown in [20] recently, the only exception is the model with the j = 1/2 case, which is not only exactly solvable, but also integrable.Moreover, since the photon number is not conserved, a procedure that considers only a truncated number of Dicke states is used together with standard diagonalization methods [7,8].Such straightforward diagonalization yields useless results, especially when the atom-field coupling becomes strong, unless a huge subspace spanned by the Dicke states is considered in the diagonalization.Usually, convergence is assumed to be achieved if the ground-state energy is determined within small relative errors.Within this approximation, one must typically diagonalize a relatively large sparse energy matrix in the strong coupling regime, which becomes prohibitive for large finite system sizes, which carries the further complication that the convergence of the ground-state energy is very slow.As a remedy, instead of the Dicke basis, the Dicke Hamiltonian can be diagonalized in a shifted boson basis where convergence is reached with a smaller number of shifted boson states across the whole coupling regime [21].This approach was also used to study the quantum criticality, the finite size effects, fidelity and the order parameter in the Dicke model [22].
Alternatively, a progressive diagonalization scheme (PDS) for the Dicke model was offered [23] based on a similar proposed methodology for solving the nuclear pairing problem [24] and the Rabi Hamiltonian [25], where the latter is just a special case of the Dicke model with N = 1.This approach shows that the lower part of eigenstates and the corresponding eigenvalues can be obtained progressively within a finite number of steps.In each step, only a single variable finite-order polynomial equation needs to be solved due to the fact that only a finite number of the shifted boson states correlates among each other, which is mainly due to the fact that the overlaps |α k nµ | of the shifted boson states |nµ) with the Dicke states |k; jµ diminish rapidly with the increasing |n − k|.An analysis of outcomes suggests the effectiveness of this diagonalization method for the model proposed in [21], in which the shifted boson basis was also adopted in the diagonalization.The convergence is tied to the number of steps that are needed, which, in turn, depends on the coupling parameter λ/ω and the number of the atoms N = 2j.Very recently, a detailed numerical analysis of the solutions of the Dicke model obtained by using the approach shown in [21] has been made [26][27][28], and a comparative quantum and semi-classical analysis of the model has also been shown [29,30], which all justify the effectiveness of the approach proposed in [21].
An application of the PDS to the Dicke model can be used for a determination of ground state properties, such as photon number, the atomic inversion, fluctuations of the latter two quantities, an entanglement measure, the Shannon information entropy, etc., and, in addition, for the study of level statistics with a larger, but finite number of atoms.The purpose of this paper is to show that, for a finite number of atoms, the entanglement measure defined in terms of the normalized von Neumann entropy of the reduced density matrix of the atoms reaches its maximum value at the critical point of the QPT, at which the system is most chaotic.Noticeable change in the Shannon information entropy near or at the critical point of the QPT is also observed.In addition, the QPT may be observed not only in the ground state mean photon number and the ground state atomic inversion, as shown in [10], but also in fluctuations of these two quantities in the ground state, especially in the atomic inversion fluctuation.

The PDS for the Dicke Model
In this section, we briefly review the PDS for the model that we will use to calculate level statistics and other quantities of the ground state.The single-mode Dicke Hamiltonian without the RWA, which describes N two-level atoms interacting with a quasi-monochromatic field, may be written as: where J z , J ± = J x ± iJ y are the usual angular momentum operators for a pseudo-spin of j = N/2, and a, a † are the bosonic operators of the field.The atomic level-splitting is given by ω 0 ; ω is the field frequency; and λ is the atom-field coupling.The parity with the operator P = e iπ(a † a+N/2+Jz) is conserved, since it commutes with the Hamiltonian ĤD .
In order to diagonalize Hamiltonian (1), we make a π/2 rotation around the y-axis to obtain pseudo-spin operators: with R = e −i π 2 Jy , which was also done in [21].After the rotation, Hamiltonian (1) becomes: Let |Ψ D be an eigenstate of (1) satisfying: After the rotation (2), we have: Therefore, once the eigenvalue Problem (3) with Ĥ|Ψ = E Ψ |Ψ is solved, all eigenvalues of (1) are thus obtained with the corresponding eigenstate |Ψ D = R −1 |Ψ .Hence, we only focus on the diagonalization procedure for Hamiltonian (3).Hamiltonian (3) can also be written as: which can easily be diagonalized when ω 0 = 0 with eigen-energies: where The corresponding eigenstates in the shifted boson basis can be expressed as: where |0; jµ for µ = −j, − j + 1, • • • , j − 1, j are the boson vacuum state satisfying a|0; jµ = 0 and the eigenstates of the pseudo-spin operators J 2 and J z , simultaneously.Like the Dicke states |n; jµ ≡ |n |jµ , where |n is the eigenstate of the boson number operator satisfying a † a|n = n|n , the basis vectors {|nµ)} also form a complete set of orthonormal basis vectors for the Hamiltonian (6).Thus, (6) can be written as: According to the analytical step-by-step diagonalization procedure proposed in [25], the procedure starts with the first term and the k = 0 sector in the last term of (9) with: For given j, the matrix elements (J x ) µµ can be expressed as: Let: be the eigenstate of Ĥ(0) , where ξ is used to label the ξ-th eigenvalue of Ĥ(0) , and: is the normalization constant.The eigenequation Ĥ(0) results in the following relation for the expansion coefficients c ξ nµ : where the overlap α k nµ is given by: It should be noted that there are some special eigenstates of (10) with nonzero expansion coefficient c ξ nµ for some fixed n, either when α 0 nµ = 0 or when the overlap: for some µ values according to (10).These eigenstates are not correlated among the shifted boson states with different n and only linear combinations of |nµ) for fixed n, so that they can be easily determined: (a) When j is an integer, since α 0 nµ=0 = 0 for n ≥ 1 according to (15), Eigenequation ( 14) provides the special solutions for this case with: for n = 1, 2, • • • , of which the corresponding eigen-energies are given by: with the parity P = (−1) n .
(b) Besides the case studied in (a), for given n, one can construct the states: satisfying (16).One can verify that there are only possible solutions (19) satisfying ( 16) when j = 1.In this case, since α 0 nµ=−1 = (−1) n α 0 nµ=1 and α 0 nµ=0 = α 0 0µ=0 δ n0 , ( 14) results in the solutions with: For given n, the eigen-energies of (20) are independent of ω 0 with: with parity P = −1.It should be emphasized that there are special solutions shown in (a) only when j is an integer and those shown in (b) only when j = 1.
In the following, we only focus on solutions other than those shown in Cases (a) and (b).Generally, these solutions are linear combinations of the shifted boson states |nµ) with different n, which are more complicated than those special cases shown in (a) and (b).When α 0 nµ = 0, let: By substituting ( 22) into ( 14), Equation ( 14) results in the linear relations among {γ µ (E ξ )} with: for µ = −j, −j + 1, • • • , j − 1, j, where: The linear relations ( 23) can be used to get expressions of γ µ (E and the equation for determining the eigenvalue ξ )} being arbitrary as long as j is finite, in which we set γ j (E 23) is solved, Hamiltonian (9) may be expressed as: where the sum in the first term should run over those provided by ( 18), ( 21) and ( 22) with (23), respectively.
Then, we do the next diagonalization with: with: where: |Ψ (1) N In this step, the special solutions in (b) will not be possible when j = 1.However, there are still solutions in (a) that exist when j is an integer with: for n = 2, 3, • • • , of which the corresponding eigen-energies are still given by: Besides these special solutions, similar to the initial step, other solutions can be obtained by setting: in which: The parameters γ ξ 1 ) in ( 31) should satisfy: Similar to the initial step, (33) provides solutions to the expansion coefficients {γ ξ 1 )} up to an overall factor and the corresponding eigenvalue E (1) Similarly, in the k-th step, we take a diagonalized part of the Hamiltonian H (k−1) and the k-th projection in the second term of ( 9) with: to do the next diagonalization for k = 1, 2 • • • , where: with The special solutions in (a) still exist when j is an integer with: for n = k + 1, k + 2, • • • , of which the corresponding eigen-energies are still given by: Besides these special solutions, other solutions can be obtained by setting: in which and the parameters γ ξ k ) should satisfy: Then, after the k-th step of the diagonalization, Hamiltonian (9) can be written as: In the k-th step, one can verify that the eigenstates (35) for any j case are mutually orthogonal.Moreover, the eigenstates of the special solutions (37) are also orthogonal with (35), which is obvious, because (35) only involves |n, 0) for n = 0, 1, • • • , k, while these states are excluded in the special solutions (37).
As shown in [23,25], the advantage of the PDS is that only a few steps are needed for lower excited states, especially when λ/ω ≤ 1, which is mainly due to the fact that the overlap |α k nµ | diminishes rapidly with the increasing of |n − k|.As a result, not only do the lower excited energies remain unchanged after a number of steps, with which the overlaps of the Dicke states in the second term of (42) with the lower excited states determined in that step are almost negligible, but also, only a finite number of states {|Ψ (k−1) ξ k−1 } correlate among each other in the k-th step, which, in turn, effectively truncates the infinite sum in (35) and ( 41) into a finite sum.Thus, other uncorrelated states and the corresponding eigenvalues in the k-th step diagonalization keep unchanged as those in the (k − 1)-th step, which has already been demonstrated in [25] for j = 1/2 case.Actually, the above feature in the PDS applies to other j cases, as well.For fixed j or λ/ω, |n − k|, with which the overlap |α k nµ=±j | is non-negligible, increases with the increasing of λ/ω or j.Namely, the more the atoms involved with the larger the coupling parameter, the more the shifted boson states |n, µ) will be correlated, which also answers why the number of steps required in order to get an accurate solution increases with the increasing of λ/ω and j.
Once all solutions in the (k − 1)-th step are known, they can be used to construct the next k-th step solutions.For example, using the results shown in (49), (50) and (46) of the initial step for j = 1/2 case, we can use ( 27)- (33) to do the k = 1 step diagonalization.In this case, according to (35), we have: where: with which (51) can be expressed as: ξ 1 ) where: Since ξ 1 ), Equation ( 53) results in the following equation for determining the eigenvalues E (1) where ν 1 = 1 or −1, from which we get: Similar linear relations of {γ ξ k )} can easily be obtained according to (41) for other j cases in the k-th step.The expressions of {γ ξ k )} become cumbersome only when j is really large.In order to demonstrate the convergence of the PDS, Figure 1a shows the first 20 eigenenergies of the model with j = 1/2 for ω = ω 0 and λ/ω 0 = 0.5 starting from k = 0 in the first step, which is only efficient for the ground state.In particular, the ground state energy remains unchanged up to the fourth decimal place after k = 3 steps of the diagonalization.More steps are needed for higher excited states.For example, the 10th excited energy remains unchanged up to the fourth decimal place after k = 10 steps.Similarly, the 20th excited energy remains unchanged up to the fourth decimal place after k = 15 steps.With the same precision, the number of steps needed increases slightly with increasing strength of the coupling parameter λ/ω 0 .For example, in the model with j = 1/2 for ω = ω 0 and λ/ω 0 = 1.0, of which the first 20 eigenenergies are shown in Figure 1b, k = 7 steps are needed for the ground state energy, while k = 20 steps are needed for the 20th excited energy.Moreover, with the same precision, the number of steps needed also slightly increases with the increasing of j.As shown in Figure 1c for the model with j = 20 for ω = ω 0 and λ/ω 0 = 0.5, k = 10 steps are needed for the ground state energy, while k = 20 steps are needed for the 20th excited energy.Therefore, the PDS starting by including the k = 0 sector in the last term of (9) and then by adding other terms successively with the increasing order of k arranged is only efficient for lower part of the spectrum with finite j, especially when j and λ/ω 0 are small.The fact that the number of steps needed is always finite is just the consequence that the overlap |α k nµ | diminishes rapidly with the increasing of |n − k|.Actually, instead of k = 0, we can use the k = t sector in the second term of (9) as an initial value for an excited state.Then, by including a few terms in the second part of ( 9) with k = t, t ± 1, • • • , t ± q with a finite desired q, the number of steps, Q ∼ 2q, needed for the excited state is similar to that for the ground state.In each step, only a single variable finite-order polynomial equation needs to be solved, because only a finite number of the shifted boson states correlate among each other, which is mainly due to the fact that the overlaps |α k nµ | of the shifted boson states |nµ) with the Dicke states |k; jµ diminish rapidly with the increasing of |n − k|.Hence, the infinite sum involved, for example that on the right-hand side of (46) or those involved in (55), becomes a finite sum.The CPU time needed for solving the equation depends on the order of the polynomial and the solver used.In the j = 20 case, as an example, the CPU time needed in each step is always about v seconds with v ∼ 2.3-15.0when λ/ω 0 ≤ 2 by using Wolfram Mathematica, which is typical not only for ground state starting from the first step with k = 0, but also for an excited state starting from the first step with k = t, where t is chosen close to the quantum number of the excited state determined by the first term of (9).Thus, besides the time needed for reorganizing the coefficients of the next step polynomial, the total CPU time needed for getting the desired result for each excited state is about vQ seconds.As can be observed from Figure 1c, one needs 10 steps in order to get the ground state energy converged, which is also typical for other excited energies as described above; the CPU time needed for each excited energy is about 10v seconds.Therefore, the PDS is quite efficient if only a few excited states are concerned, especially for the ground state.However, the PDS is still time consuming if more excited states are needed.

The QPT and Entanglement
It is well known that the Dicke model undergoes the QPT at zero temperature [2,[9][10][11], for which many ground state quantities, such as the ground state mean photon number, the atomic inversion and their fluctuations can be taken as order parameters to signify the QPT.Specifically, there is a noticeable change in these quantities at the critical point of the QPT, which is specified by the coupling parameter λ/ω.Therefore, the coupling parameter λ/ω serves as the control parameter of the QPT in the model.
By using the PDS, these order parameters at the resonance point with ω = ω 0 for some finite j cases were calculated.In Figure 2, we plot the ground state mean photon number n = Ψ g |a † a|Ψ g and the ground state atomic inversion J z = Ψ g |J z |Ψ g as functions of λ/ω 0 for the j = 10, 15 and 20 cases, respectively, which clearly illustrates the nature of the QPT.When λ/ω 0 < (λ/ω 0 ) c , the system is only microscopically excited, whereas both the field and the atomic ensemble acquire macroscopic excitations when λ/ω 0 ≥ (λ/ω 0 ) c .These results are all consistent with those shown in [10].In addition, the QPT may also be observed in fluctuations of the photon number and atomic inversion in the ground state defined by: and: These two quantities as functions of λ/ω 0 for j = 10, 15 and 20, respectively, are shown in Figure 3.In Figures 2 and 3, the vertical dashed line indicates the critical point value, (λ/ω 0 ) c = 0.5, determined in the large-N limit of the model shown in [9,10].Abrupt change in these fluctuations at or near the critical point is noticeable, especially in ∆J z .Though noticeable change in the atomic inversion fluctuation was also observed in the large-N limit of the Dicke model calculated from the exact steady-state solution through the equivalent Fokker-Planck equation of the system [31,32], there is no sharp peak emerging in the atomic inversion fluctuation at the critical point from the calculation, while, as clearly shown in the right panel of Figure 3, a sharp peak develops in the atomic inversion fluctuation near the critical point in the finite-N cases studied in this paper.Moreover, the rescaled concurrence, which can be used as an entanglement measure and also related to the atomic inversion fluctuation, in the dispersive limit of the model with ω {ω 0 , λ} was also studied [33].It should be noted that the Dicke Hamiltonian (1) in the dispersive limit is reduced to the Lipkin-Meshkov-Glick (LMG) model.As noticed in [34], though the Dicke model shares many features with the LMG model, the presence of a bosonic mode coupled to a set of two-level systems in the Dicke model may give rise to interesting challenges.Actually, the QPT in the Dicke model is mainly driven by the coupling between the bosonic field and the atoms, while the QPT in the LMG model is driven by interactions among atoms.In this sense, the QPT in the two models is different.The connection between entanglement and particle number statistics in the ground state of quantum many-body systems has been noticed.For instance, the variance of the particle number in a subgroup is an entanglement measure of the corresponding mode-bipartition for any pure state [35].Moreover, coherence between subspaces of fixed particle number in a subgroup immediately implies mode entanglement [36,37].It is clearly shown in this paper that, besides the ground state mean photon number and the ground state atomic inversion, as shown previously [10], the QPT in the Dicke model can also be observed in fluctuations of these two quantities in the ground state, especially in the atomic inversion fluctuation.In a recent paper [38], Brandes has found that there is an abrupt change in observables, such as the mean atomic inversion and the boson (photon) number and its fluctuations at arbitrary energies, which demonstrate that the QPT occurs not only at the ground state, but also in excited states.However, it should be noted that the QPT occurring in excited states is not the same as that in the ground state.The former results in a singularity in the spectrum, while the latter results in a singularity in the properties of the ground state.
To show the ground state entanglement quantitatively, we adopt the simple entanglement measure defined by: where ρ a is the reduced density matrix obtained by taking a partial trace over the field.We use the logarithm to the base N + 1 instead of base two to ensure that the maximal measure is normalized to one [39,40].In addition, Koashi entangled webs [41] were used to measure entanglement in the atomic part of the ground state of finite-size models related to the extended Dicke models [42,43], which provide another good measure of average qubit-qubit entanglement to show the highest entanglement occurring at the critical point.Moreover, the Shannon information entropy of the ground state of quantum many-body systems is also a good measure of correlations among local states [44], which, in the Dicke model, is defined by: where {w nµ } is the expansion coefficients of the ground state |Ψ g in terms of the Dicke states |n; jµ .Obviously, I H = 0 when λ = 0.In this case, all Dicke states are uncorrelated, and the ground state of the system is in the normal phase.When λ > 0, the system moves from the uncorrelated normal phase toward the correlated "superradiant" phase with I H > 0. In [34], the mutual information entropy, which is related to the reduced von Neumann entropy (59), but different from the Shannon information entropy (60), was calculated for the LMG model.Since the LMG model is equivalent to the Dicke model in the dispersive limit, the analysis provided in [34] also applies to the Dicke model in the ω {ω 0 , λ} limit.The ground state entanglement and the Shannon information entropy of the model at the resonance point with ω = ω 0 for j = 10, 15 and 20 as functions of λ/ω 0 are shown in Figure 4. Since the system studied is finite, rigorously speaking, only a cross-over occurs at or near the critical point.Nevertheless, as shown in Figure 4, with the increasing of the number of the atoms, a peak emerges in the entanglement measure near λ/ω 0 ∼ 0.5, near which noticeable changes in the photon number, the atomic inversion and their fluctuations are observed.This is consistent with what is called critical point entanglement [45].Moreover, a noticeable change in the ground state Shannon information entropy near or at the critical point can also be observed when j is large.Anyway, the above analysis justifies the connection of particle number statistics with entanglement observed in [35][36][37].

Level Statistical Properties
We use two typical statistical measures to analyze the level fluctuation properties (that is, departures from uniformity in the spectrum).One is the nearest-neighbor level spacing distribution P (S) [46], which was also calculated in [10] for several finite j cases, and the other is the spectral rigidity ∆ 3 (L) [47], which was not discussed in [10].These measures are determined through the unfolded levels { Ẽξ } that are obtained from the mapping N (E ξ ) → Ẽξ , where {E ξ } are obtained from the PDS after a number of steps until its value is unchanged in the next step, where ξ labels the ξ-th excited level, and the mapping is determined by fitting a smooth polynomial function to the number staircase function, N (E ξ ), that counts the number of levels below E ξ .In our analysis, we use a six-degree polynomial in E ξ that has been shown to be sufficient [46,48].By construction, the unfolded spectrum, { Ẽξ }, is dimensionless with an average level spacing of unity.
The distribution of the nearest-neighbor level spacing P (S) is defined as the probability of two nearest-neighbor energy levels to have a spacing S, which is calculated from the unfolded spectrum as: The distribution P (S) is shown by normalized histograms in our analysis.
As is commonly accepted, Poisson statistics with: characterizes a regular system (uncorrelated level spacings), while the Wigner distribution, which is almost identical to the GOE prediction [46,47], characterizes a chaotic system (nearby levels are likely to repel each other).The spectral rigidity, that is the departure from uniformity (even spacings of a rigid spectrum) over a given span of levels, is measured by ∆ 3 (α,L) as defined by Dyson and Metha [47], which is the average of the least-squares deviations between the number staircase function N ( Ẽ) for the unfolded spectrum and its best linear fit (A Ẽ + B) over the energy interval [α, α + L].For a given L, smaller values of ∆ 3 imply stronger long-term correlations between the levels.The average of ∆ 3 (α, L) over n α intervals (α, α + L), which overlap by L/2 successively, yields a smoother ∆ 3 (L) measure, It is shown that: for a regular spectrum, while, in the large-L limit, for a chaotic system with GOE statistics.
In our calculations, 8000 excited levels with positive parity were considered, which are obtained from the PDS according to Section 2. Typically, the highest level energy is 390 ω 0 for λ/ω 0 = 0, 402.302 ω 0 for λ/ω 0 = 0.5, 862.495 ω 0 for λ/ω 0 = 5.0, 11,940.60ω 0 for λ/ω 0 = 100 and 58,789.7 ω 0 for λ/ω 0 = 500, respectively.We have checked that there are no obvious changes in P (S) and ∆ 3 (L) when more excited levels are taken in the calculations.For small j values, a non-generic distribution consisting of several isolated peaks in P (S) emerges, which is similar to the Rabi model case with j = 1/2.These cases are unusual, because the number of the atoms is too small to form a universal ensemble.Therefore, we only study the level statistics of the model with j = 20 as an example.Since the parity is conserved, only energy levels with positive parity are included in the statistics.Using the PDS, we calculate the positive parity level energies, which are then used to obtain both the nearest-neighbor level spacing distribution P (S) and the spectral rigidity ∆ 3 (L), of which the results are shown in Figures 5 and 6, respectively.Instead of the positive parity level energies, one can also use all negative parity level energies to calculate these two quantities, with results being quite similar to those obtained from the positive parity level energies.As shown in the upper left two panels in Figures 5 and 6, when λ/ω 0 is at or near zero, the model is obviously integrable or quasi-integrable, especially when λ/ω 0 = 0, for which the spectrum is quite similar to that of a harmonic oscillator.The P (S) reaches the unique maximal value at S = 0 when λ/ω 0 = 0, while ∆ 3 (L) in this case is far beyond the Poisson type, indicating much weaker correlation among the levels.As pointed out in [49], harmonic oscillator systems with equidistant excitation energies are unusual and non-generic, though they are integrable, which explains why the P (S) and ∆ 3 (L) for the case with λ/ω 0 ∼ 0 behave unusually, as shown in the upper left two panels in Figure 5.With increasing λ/ω 0 , both the nearest-neighbor level spacing distribution P (S) and the spectral rigidity ∆ 3 (L) gradually develop towards GOE-type statistics, which coincide with the corresponding GOE values at the critical point with λ/ω 0 = 0.5, as shown in the first lower left panel in Figures 5 and 6.By further increasing λ/ω 0 , the P (S) and ∆ 3 (L) tend to the Poisson type, which almost coincide with the corresponding Poisson values when λ/ω 0 is large enough.In fact, the Dicke model with j ≥ 1 for nonzero energy splitting between the two levels of the single atom is formally not integrable when λ/ω 0 → ∞, though the model is exactly solvable and integrable when j = 1/2, as recently pointed out by Braak [16,20].Braak's exact solution to the Dicke model with j = 1/2 was also recovered [50,51] by using the Bogoliubov operators and the method proposed in [22], respectively, while the finite size Dicke model was also studied within the shifted boson basis in [52], which confirms Braak's observations on the exact solvability of the Dicke model, though the model is not integrable in general.However, the Poissonian limit is reached as λ/ω 0 → ∞ at the resonance point, which is also demonstrated in [10].Generally, the GOE statistics is mainly due to the fact that no level crossing occurs among excited levels of the model, which is suggestive of a departure from regularity and emergence of quantum chaos [53,54].Actually, obvious level repulsion among excited levels at or near λ/ω 0 = 0.5, e.g., when j = 20, can be observed, though the spectrum is not shown here.In [38], Brandes observed the singular behavior of the derivative of the density of states at the critical point.Furthermore, it is shown in [29,30] that the density of states could be approximated semiclassically, and this approximation could be used as an exact unfolding to avoid fitting a polynomial function by hand in this model.Since the model can be studied semiclassically [30], a connection between the QPT in excited states of the model [29,38] and the onset of chaos was identified.It is also found in [30] that the onset of chaos is related to the breaking of the quadratic approximation of the Hamiltonian that is obtained by considering small oscillations around the global energy minimum.It was confirmed in the quantum and classical versions that chaos is present, both in the normal and superradiant phase.Anyway, our results clearly show that the system at the resonance point is most chaotic at the critical point with λ/ω 0 = 0.5.

Conclusions
In summary, ground state properties and the level statistics of the Dicke model for a finite number of atoms are investigated based on the progressive PDS.Particle number statistics, the entanglement measure, the Shannon information entropy and the level statistics at the resonance point for cases with a finite number of atoms as functions of the coupling parameter are calculated.The ground state mean photon number, the ground state atomic inversion and the distribution of the nearest-neighbor level spacing P (S) obtained from the results of the PDS agree with those shown in previous studies, e.g., those shown in [10].In addition, the QPT may also be observed in fluctuations of these two quantities in the ground state, especially in the atomic inversion fluctuation.A noticeable change in the Shannon information entropy near or at the critical point of the QPT is also observed.It is clearly shown that the entanglement measure defined in terms of the normalized von Neumann entropy of the reduced density matrix of the atoms reaches its maximum value at the critical point of the QPT, when the system is most chaotic.When a quantum many-body system undergoes a QPT, there is always a noticeable change in many ground state quantities, such as particle number statistics and the entanglement measure at the critical point [35][36][37].However, it is not necessarily true that the system is most entangled at the critical point.For example, the entanglement measure increases with increasing value of the control parameter in the Jaynes-Cummings model [6], which is also the Dicke model with the RWA, the spin chain models [55], the Bose-Hubbard model with on-site repulsion [56], etc., though saturation in the measure will reach beyond the critical point in the strong coupling regime.On the other hand, in general, chaotic systems tend to produce larger entanglement than for regular systems, but there are also exceptions for classically regular systems, as shown in [57,58].Therefore, the critical point entanglement with the concurrent onset of chaos seems unique in the Dicke model without the RWA.

Figure 2 .
Figure 2. The ground state mean photon number n (left panel) and the atomic inversion J z (right panel) as functions of λ/ω 0 at the resonance point with ω = ω 0 for 10 (red solid dots), 15 (open circles) and 20 (solid line), where the vertical dashed line indicates the critical point position determined in the thermodynamic limit [10].

Figure 3 .
Figure 3.The same as Figure 2, but for the ground state photon number fluctuation ∆n (left panel) and the atomic inversion fluctuation ∆J z (right panel) as functions of λ/ω 0 at the resonance point with ω = ω 0 for 10 (red solid dots), 15 (open circles) and 20 (solid line).

Figure 4 .
Figure 4.The same as Figure 2, but for the ground state entanglement (left panel) and the Shannon information entropy (right panel) as functions of λ/ω 0 at the resonance point with ω = ω 0 for 10 (red solid dots), 15 (open circles) and 20 (solid line).

Figure 5 .
Figure 5. Level spacing distribution P (S) of the model with j = 20 at the resonance point with ω = ω 0 for various coupling strengths λ/ω 0 .In all panels, the (blue) dashed line describes the Poisson statistics, while the (red) dash-dotted line describes the GOEstatistics.

Figure 6 .
Figure 6.The spectral rigidity ∆ 3 (L) of the model with j = 20 at the resonance point with ω = ω 0 for various coupling strengths λ/ω 0 .In all of panels, the (blue) dashed line describes the Poisson statistics, while the (red) dash-dotted line describes the GOE statistics.