Exact solvability, non-integrability, and genuine multipartite entanglement dynamics of the Dicke model

In this paper, the finite size Dicke model of arbitrary number of qubits is solved analytically in an unified way within extended coherent states. For the $N=2k$ or $2k-1$ Dicke models ($k$ is an integer), the $G$-function, which is only an energy dependent $k \times k$ determinant, is derived in a transparent manner. The regular spectrum is completely and uniquely given by stable zeros of the $G$-function. The closed-form exceptional eigenvalues are also derived. The level distribution controlled by the pole structure of the $G$-functions suggests non-integrability for $N>1$ model at any finite coupling in the sense of recent criterion in literature. A preliminary application to the exact dynamics of genuine multipartite entanglement in the finite $N$ Dicke model is presented using the obtained exact solutions.


Introduction
The Dicke model [1] describes the interaction of a number (N > 1) of two-level atoms (qubits) with a single bosonic mode and has been a paradigmatic example of collective quantum behavior. It is closely related to many fields in quantum optics, quantum information science, and condensed matter physics. The superradiant behavior by an ensemble of quantum dots [2], Bose-Einstein condensates [3,4], and coupled arrays of optical cavities has been used to simulate the strongly correlated systems [5].
Quantum integrability is an often-mentioned concept in the description of a quantum system, however it is not well defined to date [6], in sharp contrast with the classical integrability, mostly due to the markedly different ways to count the degrees of freedom (d.o.f) in the quantum and classical mechanics. This situation becomes more serious for the Dicke model, where the classical limit is lacking because of the quantum nature of the finite number of levels, one cannot define integrability from any classical limit. In the semi-classical limit, where the boson mode is treated as a classical field while the discrete level is kept as a quantum entity, the signatures of chaos emerge in the Dicke model [7], even in the Rabi model (single qubit) [8]. One appealing way to address the quantum integrability is looking at the energy level statistics. The corresponding Berry-Tabor criterion [9] states that a quantum system is integrable if the level statistics is Possionian. However it relies on the semiclassical arguments and only concerns quantum systems with continuous d.o.f. The level statistics of the finite but large N Dicke model is shown to be Possionian below and Wigner-Dyson above the critical coupling [10,11]. While the level statistics of Rabi model is neither Poissonian nor Wigner [12], which is characteristics of genuinely integrable models. Recently, Braak proposed that if the eigenstates of a quantum system can be uniquely labeled by f = f 1 + f 2 quantum numbers, where f 1 and f 2 are the numbers of the discrete and continuous d.o.f, then it is integrable [13]. By this novel criterion, the N = 3 Dicke model is non-integrable [14], while the Rabi model is integrable [13]. Implications of the (non)-integrability in the Dicke model with arbitrary number of qubits should be generally interesting.
Despite the simplicity of the full Dicke model, its solution is however highly nontrivial. Many approaches have been developed and extended to the this model [11,15,16,17,18,19,20,21,22,23,24]. Among them, the exact eigensolutions in the finite N Dicke model are very crucial to get some deep insights, but they are mainly limited to numerical diagonalization in truncated Bosonic Fock space [11,17,19]. Recently, Chen et al. [24] have presented numerical exact solutions to this model using extended coherent states (ECS) [25], where the truncation of the Hilbert space can be alleviated systematically. It has been extensively shown in Refs. [26,27,28] that as compared with the photonic number (Fock) basis, ECS is exhibited to be valid for a large region of the Hamiltonian parameter space by analyzing the converged energy eigenvalues and eigenfunctions.
More recently, an analytical exact solution, which is mathematically well defined, to the Rabi model [29] has been discovered by Braak [13] using the representation of bosonic creation and annihilation operators in the Bargmann space of analytical functions [30]. A so-called G-function with a single energy variable was derived yielding exact eigensolutions. Alternatively, using ECS, this G-function was recovered in a simpler, yet physically more transparent manner by Chen et al. [31]. The continuous extensions to the analytical treatments for the Dick-type model with less than 3 qubits have also achieved progress [14,32,33,34,35,36,37,38,39,40,41]. Among them, complicated 6×6 determinant G-functions were derived in the Bargmann representation for the N = 3 Dicke model [14]. By using the ECS, G-function resembling the most simple one without a determinant in the Rabi model was given for the N = 2 Dicke model [39] and 2 × 2 determinant G-function was derived in the Rabi model with two arbitrary qubits [41]. A generally concise analytic solutions to the Dicke model of arbitrary number of qubits, which is a counterpart to the previous numerical exact solution [24], should be highly called for.
On the other hand, quantum correlations among qubits, viewed as quantum information resources, have also attracted extensive interest. A dozen qubits are sufficient in the applications to quantum information technology. However, only bipartite quantum correlations have been studied in the Dicke model [17,18,19,20,21,22,42], and multi-partite quantum correlations have not been explored, to the best of our knowledge. The multi-partite quantum correlations can be used to implement operations in the measurement-based quantum computation [43] and high-precision metrology [44], thus is most important from the experimental viewpoint. However, it was difficult to characterize multi-partite quantum correlations in many previous efforts [45]. Recently, the theory of multipartite entanglement has been achieved much progress, especially the computable entanglement monotone [46]. On the experimental side, multi-qubits have been constructed in several solid devices recently, such as the well-known Greenberger-Horne-Zeilinger (GHZ) states [47] of more than two qubits using the complete circuit and full set of gates [48,49,50]. The path to scalability is prerequisite to realize the quantum computation and quantum information processes.
The goal of this paper is twofold. First, we demonstrate a successful derivation of a concise G-function for the arbitrary N Dicke model for by means of ECS [31,39,41]. For N = 2k − 1 or 2k, the G-function is a k × k determinant, which are very feasible to calculate all eigensolutions for the Dicke model with a dozen qubits. The non-integrability of the finite N Dicke model is then discussed in terms of the level distribution dictated by the G-functions. Second, using the exact eigensolutions, we study the dynamics of the genuine multipartite entanglement (GME) from the maximum entangled states, such as Bell states for N = 2 and GHZ states for N = 3 and 4.
The paper is organized as follows. In Sec. 2, we introduce the model Hamiltonian and briefly review the previous numerical exact techniques based on the ECS. Section 3 describes the detailed analytic scheme for the solutions to the finite N Dicke model. The G-functions are derived explicitly. In Sec. 4, we illustrate that the exact spectra can be easily obtained by the stable zeros of the G-function and discuss the implications of the present derivation with the non-integrability of the finite N Dick model. As an example of the applications of the present method, the exact dynamics of the GME of the Dicke model is studied using the obtained exact eigensolutions in Sec. 5. A summary is given in the final section.

Model
The Hamiltonian of the finite N Dicke model can be written as ( = 1) where ∆ is the qubit splitting, d † creates one photon in the common single-mode cavity with frequency ω, λ describes the atom-cavity coupling strength, (we usually is the Pauli operator of the i-th qubit, J is the usual angular momentum operator, and J + and J − are the angular raising and lowing operators and obey the SU(2) Lie algebra [J + , J − ] = 2J z , [J z , J ± ] = ±J ± . This Hamiltonian posses a Z 2 -symmetry like the Rabi model and exhibits a second-order quantum phase transition at λ c = 1/2 in the thermodynamic limit [11]. In this work, the Hilbert space is spanned by the Dicke state N 2 , m , which is the eigenstate of J 2 and J z with the eigenvalues N 2 ( N 2 + 1) and m. Throughout this paper, the unit is taken of ω = 1.
Previously, we have proposed a numerical exact approach for Dicke model of arbitrary N by ECS [24]. The wave function is described in the basis j = N 2 , m |ϕ m d where the photonic state is given by where |0 is the photonic vacuum state, g m = 2gm (m = j, ..., −j) ,and N c is the truncated bosonic number in the Fock space of the new operator A † m = d † + g m . The coefficient u m,n can be determined through the exact Lanczos diagonalization. In that work, we have solved the Dicke model for more than one thousand qubits numerically, and obtained some critical exponents with high accuracy. It was numerically shown later that increasing the number of atoms imposes a strong limit to the states in the Fock basis, but not in the ECS scheme [26,27]. Interestingly, the similar or same ansatz for the wavefunction like (2) have been used in many recent works for the Dick-type models.
In the next section, we will propose an analytic scheme to the exact solution for the finite N Dicke model.

Analytical scheme to exact solutions
We will demonstrate that the N = 2k − 1 and N = 2k Dicke models can be treated in an unified way. In the matrix form of the Dicke state basis, the Hamiltonian takes a tridiagonal type. The diagonal element is and δ is the Kronecker delta. We introduce k pairs of displaced operators as where m = 0 for N even case is excluded here, because it is just the original photonic operators without displacement and will be independently treated later.
Specially, the diagonal element H m,m for m ′ = m can be transformed into the form with free particle number operators.

Then we propose the wavefunction in Fock space of
where c m ′ ,n are coefficients and |n A +m are the number states in A † +m . The latter is just called as ECS [31] with the following property where the vacuum state |0 A +m is the eigenstate of the one-photon annihilation operator d with eigenvalue −2mg.
The Schrödinger equation leads to the recurrence relation for the expanded coefficients as and We can see that c For the same non-degenerated states, |A +m should be proportional to |A −m . Comparing the Dicke states |j, m and |j, j − m gives Left multiplying 0| and with the use of where +(−) in the l.h.s denotes the corresponding states of positive (negative) parity. Therefore we have k sets of linear equations for k pairs of displaced operators. For the most simple Rabi model, we only have one linear equation with either parity, and the relevant coefficients c ,n can be determined recursively from the initial value c , because Eqs. (7) and (8) in this case can be reduced to a linear threeterm recurrence relation. For N > 1, such linear three-term recurrence relations are not available, and the coefficients cannot be obtained from a single initial one, but instead, from N initial parameters for each k. Then one may naively think that all coefficients in the k sets of linear equations (10) should be determined by k × N initial parameters. Fortunately, we will show below that k independent initial parameters would determine all coefficients recursively.
The wavefunction can be also expressed in original Fock space as where a −m ′ ,n = ± (−1) n a m ′ ,n with +(−) positive (negative) parity. Note that this auxiliary expansion itself exists intrinsically for m = 0 in the N even case, c.f. Eq. (3). By the Schrödinger equation, all coefficients a m ′ ,n can be determined by k initial ones a m ′ =0,n=0 recursively and Note that a 0,n is only present for the N even Dicke model. Similarly, the same non-degenerated states in Eqs. (5) and (11) yields Exact solvability, non-integrability, and genuine multipartite entanglement dynamics of the Dicke model7 Projecting onto A +m 0| gives the initial coefficient c where the use has been made of To this end, through Eqs. (7), (12), and (15) , all coefficients c For the mth linear equation, set a m ′ ,n=0 = 1, and all other initial variables to be zero, the matrix element G ± (m ′ , m) is given by the summation in Eq. (10). The G-function for the Dicke model is thus defined as the following k × k determinant For nonzero k unknown variables a m =0,n=0 , G ± (E) = 0 is required. The zeros thus give all eigenvalues of the Dicke model, which in turn give the eigenstates using Eq. (5) or Eq. (11). It is interesting to note from the coefficients in Eqs. ( 7) and (15) that the present G-function is a well defined transcendental function. Thus analytical exact solutions have been formally found for the Dicke model. Exceptional solutions: Similar to the Rabi model [13,31], for some special model parameters, there are exceptional solutions that do not correspond to the zeros of the G-functions. They can be also obtained from the pole singularities in Eq. (8) for m > 0. It is just the mth kind exceptional solution. So we totally have k kinds of exceptional solutions, which are just the eigenvalues in the atomic degenerate limit ex into the k sets of linear equations (10), we can note that the equation corresponding to A † +m is not available due to the singularity, which can be replaced by i.e. the numerator in Eq. So there are still k sets of linear equations available for the mth kind exceptional eigenvalue. The zeros of the corresponding k ×k determinant will yield the condition for the occurrence of the exceptional eigenvalues. For a given ∆, we can obtain the coupling constants g n,m ± , which can also be located by the crossing points of curves for Eq. ( 17) and the corresponding energy levels in the spectral graph. Generally, g n,m + = g n,m − , owing to the different parity dependent conditions. It follows that exceptional eigensolutions Exact solvability, non-integrability, and genuine multipartite entanglement dynamics of the Dicke model8 are non-degenerate, as noted in the N = 3 Dicke model [14]. It should be stressed here that the level crossing in two parity subspace does not occur at an exceptional eigenvalue, and therefore the doubly degenerate state at the level crossing has the regular spectrum with fixed parity for the N > 1 Dicke model, in sharp contrast with the Rabi model.
These exceptional solutions result in k kinds of pole structure in the G-function for both the N = 2k − 1 and the 2k Dicke models. For N = 2k, there exists additional pole at E = n even (odd) for the positive (negative) parity, as can be seen from Eq. ( 13). It is just the eigenvalue for j = 0 subspace, and isolated from other Dicke states with nonzero j. It is not an exceptional solution but yields divergence as well, as demonstrated in the N = 2 Dicke model [39]. It will be shown later that all these poles are very helpful to analyze the distribution of the zeros, namely eigenvalues, in the G-function.

Numerical zeros of G-functions and non-integrability
Analogous to the Rabi model [13,31], zeros of the G -functions in the Dicke model cannot be obtained analytically in the closed form either. Searching for the zeros numerically should be required practically. The curves of G ± (E) defined in Eq. (16) is plotted in Fig. 1(a) for ∆ = 0.7 and g = 0.25 for the N = 3 Dicke model. The stable zeros give all eigenvalues in the plotted energy regime, which has been confirmed with the numerical exact solutions. Typically, the convergence is assumed to be achieved if zeros (i.e. E) are determined within relative errors |(E Nc − E Nc+1 ) /E Nc+1 | < 10 −8 , where N c is the truncated number in the series expansions for the displaced operators. We also calculate R Nc = ln (E Nc /E Nc+1 ) for all zeros in Fig. 1(b). The stable zeros, where R Nc is almost the same as the relative errors, should be on the R Nc = 0 line within error 10 −8 , much smaller than the symbol size.
One can also notice that very few unstable zeros emerge in practical calculations because of inevitable finite truncations. They do not belong to the true eigenvalues and can be excluded very easily. They are very sensitive to the truncated number N c and cannot converge with N c , because the corresponding coefficients c (m) m,n oscillate with increasing magnitudes as n increases. In sharp contrast, for the stable zeros, the coefficients c (m) m,n converge to zero rapidly with increasing n. The positions of unstable zeros must change if increasing N c by 1. So they are easily figured out in Fig. 1 (b) for apparent deviation from R Nc = 0 line. It should be pointed out here that unstable zeros are absolutely not the true zeros of G-function if the summations are really performed infinitely. As demonstrated in Fig. 1 (b) that their positions shift to higher energy regime (E Nc+1 > E Nc ) with increasing N c . Theoretically, the unstable zeros disappear in the finite energy spectra if N c → ∞.
The baselines shown in Fig. 1 (a) are close to m = 3/2 and m = 1/2 exceptional eigenvalues due to the divergence in the G-functions. For the present model parameters, no exceptional solution is missed in the stable zeros of the G-function. In principle, the condition for the occurrence of the exceptional eigenvalues is hardly satisfied for given rational model parameters. Note that for the N = 3 Dicke model, the present G-function within ECS is only a 2 × 2 determinant, much simpler than G-functions with 6 × 6 determinant derived by Braak using the Bargmann representation [14]. The G-curves are different for the same model parameters, but the stable zeros should be the same. What is more, we can study the large N Dicke model straightforwardly. The only thing we need do is changing the value of N in these formulae. We like to stress here that increasing the number of qubits N would almost not bring much more additional effort to search for the zeros of the present G-function numerically. By 6 × 6 determinant, we can study the N = 12 Dicke model, which is sufficient for the recent multi-qubits constructed on the superconducting circuits for scalability. The corresponding G-function is exhibited in Fig. 2 for ∆ = 0.7 and g = 0.15, equivalently λ ≈ 0.52, slightly higher than the critical coupling. The unstable zeros can be also easily excluded by the same procedure as stated above. To reduce the magnitude of the value of G-function in the large N Dicke model, in practice, one can re-scale the elements in Eq. (16) in each row by a same real number. The shape of the G function will be modified, but the positions of the stable zeros remain unchanged.
It is interesting to link coefficients in wavefunction (2) in previous numerically exact techniques [24] and the present wavefunctions (5), (9), and (11) in the same Dicke state The third one only exists for N even. So the convergency with the truncated number N c of the ECS in both approaches should display the similar behavior in the practical calculations, although the previous coefficients are obtained from numerical diagonalization and the present ones by the zeros of the G-functions.   To show the convergence, we present the relative difference of the energy η Nc = |(E(N c ) − E(N c − 1)) /E(N c )| as a function of N c in Fig. 3 for N = 12 model. It is observed that the energy can converge rapidly with N c to any desired accuracy. Interestingly, η Nc curves almost saturate for the high excited states, demonstrating the common advantage of the ECS technique. Generally N c = 20 is sufficient to achieve the eigenvalues with very high accuracy for a large region of the model parameter space, similar to the numerical diagonalization in the ECS basis [24,26,27,28]. The magnitude of the value of G-function will increase rapidly with N in the low energy region, so rescaling and/or high precision computations are needed in some cases. Actually, we do not have to plot the G-functions to obtain the exact solutions and only locate all zeros numerically for both N c and N c + 1 instead, followed by checking the stability of all zeros.
Non-integrability: Both Dicke and Rabi models possess the Z 2 symmetry and have therefore two conserved quantities, energy and parity. The dimension of the finitedimensional factor of the Hilbert space is N + 1 in the Dicke case, thus greater than two for N > 1, while in the Rabi model it is two. The parity symmetry provides two labels and this suffices to label each state uniquely for N = 1 but not for N > 1. It follows that the Dicke model is non-integrable [14]. It is, however, analytically solvable, as shown in the previous section. This demonstrates again (another case is the driven Rabi model [13]) that integrability and solvability are different concepts, at least concerning systems with constituents in the deep quantum limit. Braak's criterion is consistent with both the manifestation of chaos in the semiclassical counterpart of the Dicke model [7] and the "picket-fence" character of spectrum in the Rabi model [12].
By the structure of the derived G-functions, we will discuss the non-integrability of the finite Dicke model in terms of the energy level distribution in some detail below.
First, we review two limits of the Dicke model. In the limit of g = 0, the atoms are decoupled from the field, the eigenvalues are where n ′ (= 0, 1, 2, ...∞) is the photonic number in the original Fock space d † , m ′ (= −j, ..., +j) is the eigenvalues of J x . In this case, the system is certainly integrable. The symmetry is not weaker than that in the rotating-wave approximation [51]. In the strong coupling limit g → ∞, the off-diagonal elements in (1) can be neglected, the system can be described by the direct product of the Dicke state N 2 , m and number states in the displaced operators A † m , the eigenvalues are then given by E m,n = n − 4m 2 g 2 , where n (= 0, 1, 2, ...∞) is the photonic number in A † m , m (= −j, ..., +j) is the eigenvalue of J z . They are doubly degenerate for m = 0. The whole spectrum must show a regular behavior. Interestingly they have the same expressions as the poles of the G-functions at finite coupling. Based on the closed-form solution (21), the finite N Dicke model in this limit is considered to be integrable by Emary and Brandes [10,11]. The degenerate atomic limit ∆ = 0 is actually equivalent to the strong coupling limit.  Recently, Batchelor and Zhou argued that both Dicke and Rabi models are integrable in this limit [52] in the sense of Yang-Baxter concept [53].
Generally, there are N + 1 eigenvalues in the unit energy interval in the weak coupling limit and N +2 2 N +1 2 eigenvalues for N even (odd) in the strong coupling limit in the E > 0 region. The nontrivial eigenvalues at finite coupling should interpolate from two limits given by Eqs. (20) and (21) which can be located by the zeros of G-function. Without loss of generality and also for brevity, we demonstrate the G + (x)-functions, where x = E + N 2 g 2 , for N = 3 and 5 in Fig. 4 for two typical coupling strengths below and above λ c = 1/2 to explain the non-trivial level statistics. The baselines due to divergence are at the poles x (m) n = n + (N 2 − 4m 2 ) g 2 . The G-function is not analytic at the baselines, so the positions of the zeros are dictated by the pole structure of Gfunction. It is clearly exhibited in Fig. 4 that all segments of the G-curves within the same parity start from one baseline and end at the adjacent baseline, indicating that the zeros must be restricted in the subintervals given by the nearest neighbor baselines. In this way, the level distribution is dominantly controlled by the pole structure of G-functions.
For small coupling, the baselines x (m) n within the same n are very close, so zeros in each parity subspace in the interval [n, n + 1] prefer to stay in the wider subintervals on the right side rather than the narrow subintervals on the left side. As N increases, one can conjecture that the zeros are squeezed to few wider subintervals, leading to relatively more very small avoided crossings exhibited in the spectral graph. Interestingly, more small avoiding crossings for N = 5 than N = 3 are obviously seen in Fig. 4 at λ = 0.2. The general trend can be confirmed numerically. Note that the present G-function is single valued for the energy variable, so only the level repulsion is allowed and the level crossing can be excluded in each parity subspace. Although we could not rule out the rare events that the G curve is exactly tangent to G = 0 line if the extremal condition is also met, leading to the level "collisions" in the spectral graph, they are however only related to fragile degeneracies and by no means the true level crossing. There is no visible difference between a true crossing and a small avoided crossing in the level statistics, so the Poissonian statistics observed at the weak coupling [11] is a consequences of rather many very small avoided crossings than many true level crossings. The level crossing is forbidden by non-integrability. The Berry-Tabor criterion may only be suited to the quantum models having the classical limits, which comprise many important systems including the infinite N Dicke model [11], but not the finite Dicke model.
As the coupling increases, the largest subinterval between x (N/2) n and x (N/2−1) n can exceed the main interval [n, n + 1]. In this way, the crossover coupling constant can be estimated as Very interestingly, for large N, λ (N ) c is close to 1/2, the critical point of the quantum phase transitions. For finite but large N, many wide subintervals emerge above λ (N ) c , and then provide the substantial ways to assign the zeros, allowing an irregular distribution of the energy level. This trend can be also seen in Fig. 4 at λ = 0.8. Alleviating the jam of the zeros in few subintervals would reduce the avoiding crossing and then lead to the Wigner-like statistics at strong coupling, consistent with the observation in the level statistics [11]. Specially, for the Rabi model, there is no subinterval and the only interval [n, n + 1] remains unchanged for arbitrary coupling, so the energy distribution is more regular than that in the Dicke model [12,13].
To this end, we believe that the Dicke model is exactly solvable but not integrable for any finite λ and N > 1. Besides these fundamental issues, we may also find many practical applications of the present analytical exact technique in the physical problems. In the next section, we will apply it to the exact dynamics of GME of this model, which has not been explored before.

Dynamics of genuine multipartite entanglement
First, we briefly review the scheme to the calculation of the GME. A state ρ for multipartite systems is called biseparable if it can be written as a mixture of states separable with respect to different bipartitions. It is well-known that separable states are always positive partial transpose (PPT) [54], indicating that a set of separable states with respect to some bipartition is contained in a larger set of states with PPT for the same partition. Since any biseparable state is a PPT mixture, a state which is not PPT mixture implies genuine multipartite entanglement.
Recently, a new approach to characterize genuine multipartite entanglement has been proposed [46] by considering the PPT mixture. PPT mixtures can be characterized by the method of semidefinite programming [55]. Given a state ρ for multipatite systems, one can search the minimization of the trace of matrix W ρ, where W is a decomposable entanglement witness for any bipartition. It has been proved that a state is PPT mixture if minimum of T r(W ρ) is positive. If the minimum is negative, then ρ is not a PPT mixture hence is genuinely multipartite entangled. The optimization problem can be solved by using YALMIP [56], SEDUMI [57] or SDPT3 [58], a readyto-use freely available implementation. Moreover, the absolute value of such minimum (denoted byE(ρ)) was proved to be an entanglement monotone for genuine multipartite entanglement [46].
In the finite N Dicke model, we start from the maximum entangled states in the N-qubit basis which are just the Bell states for N = 2 and GHZ states for N > 2. In the basis of Dicke states, it can be rewritten as The initial field state is the vacuum state |0 d . Since [H, J 2 ] = 0, the state will only evolve in the subspace spanned by the basis N 2 , m |l d . Using the n-th eigenvector of the Dicke model given in Sec. II we can expand initial state (22) as The evolution of wave function thus is given by: To calculate the dynamics of the GME, each Dicke state | N 2 , m should be expressed in N-qubit basis |j q (j = 1, ..., 2 N ). The atomic reduced density matrix ρ can be calculated by tracing out the photonic degrees of freedom ρ(t) =Tr photon |ψ(t) ψ(t)| with dimension 2 N .
The evolution of the GME initiated from the qubit maximum entangled state and field vacuum state can be studied within the above definition and the exact eigensolutions. Fig. 5 presents the GME dynamics for different coupling constant g at resonance ∆ = 1 for N = 2, 3, and 4 respectively. Periodic or quasi periodic behavior for the GME dynamics is observed at the very weak coupling, which is basically the same as the results in the RWA where the photonic number is conserved. With the increasing coupling, the regular behavior is gradually destroyed, mostly due to the activation of more photons. It is worth further study whether the "chaotic" behavior at the strong coupling is relevant with the quantum chaos in this model [11]. Very interestingly, we find that at the same coupling strength (g) between the single-mode cavity and the individual qubit, the GME for more qubits are more stronger, which will be an advantage if GME for more qubits act as the resource in the quantum information processing.

Summary
In this work, we have derived a concise G-function by a k × k determinant for the arbitrary N = 2k − 1 and 2k Dicke model. This G-function is a well defined transcendental function, so formally analytical exact solutions to the Dicke model is found in a strict mathematical sense. Without built-in truncations, the present method is essentially different from the previous numerical exact ones, therefore of more academic values. The mathematics behind the G-function is very interesting and may be worth further exploration in the future. This work is to extend the methodology of G-function in the Rabi model to the identical multi-qubit cases, thereby allowing not only in-depth studies in some fundamental issues but also practically feasible treatment to energy spectra and eigenstates.
The stable zeros of this G-function will lead to exact eigensolutions of the model, which are useful for many applications in the Dicke model. The GME dynamics has been calculated for Dicke model with a few number of qubits as a example in this paper. It is shown that GME dynamical behavior is strongly N dependent. The GME becomes stronger with the increasing qubit number for the same coupling.
In the Dicke model, we find that the exact regular spectrum can be described in terms of infinite polynomials, and can be simply located numerically. In this sense, one may argue that these exact solutions are not fully analytic. While the isolated exact exceptional ones can be expressed algebraically, actually a parabolic function of g. In any case, the exact solvability has been demonstrated without doubt. The eigenstates cannot be uniquely specified by the quantum numbers representing the continuous and discrete d.o.f, suggesting non-integrability of the N > 1 Dicke model at any finite λ in terms of Braak's criterion. The level distribution is mainly controlled by the pole structure of the G-functions, which is determined by the exceptional spectrum. We argue that the Poissonian statistics at the weak coupling observed in literature is the consequence of rather many very small avoided crossings than many true level crossings. The absence of the more level crossings than allowed by parity symmetry is a strong indication of non-integrability. The present work may add a new example to the nonintegrable but exactly solvable models.