{'fraction_non_alphanumeric': 0.04242843396680298, 'fraction_numerical': 0.011907625691604522, 'mean_word_length': 4.336489088575096, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Exact spectrum of the Lipkin-Meshkov-Glick model in the thermodynamic limit and finite-size corrections

The spectrum of the Lipkin-Meshkov-Glick model is exactly derived in the thermodynamic limit by means of a spin coherent states formalism. In the first step, a classical analysis allows one to distinguish between four distinct regions in the parameter space according to the nature of the singularities arising in the classical energy surface; these correspond to spectral critical points. The eigenfunctions are then analyzed more precisely in terms of the associated roots of the Majorana polynomial, leading to exact expressions for the density of states in the thermodynamic limit. Finite-size effects are also analyzed, leading in particular to logarithmic corrections near the singularities occuring in the spectrum. Finally, we also compute expectation values of the spin operators in a semi-classical analysis in order to illustrate some subtle effects occuring in one region of the parameter space.


I. INTRODUCTION
The Lipkin-Meshkov-Glick (LMG) model was proposed in 1965 to describe shape phase transitions in nuclei 1,2,3 . This model is often used to describe the magnetic properties of molecules such as Mn 12 acetate 4 . However, it also captures the physics of interacting bosons in a double-well-like structure 5,6 and is thus relevant to (two-mode) Bose-Einstein condensates 7 as well as Josephson junctions. It has also been recently used in optical cavity quantum electrodynamics in its dissipative version 8,9 for studying the decoherence of a single spin coupled to a spin bath 10,11 or quench dynamics 12 . Note also that, in recent years, the entanglement properties of its ground state 13,14,15,16,17,18,19,20,21,22 as well as the finite-size behavior 23,24,25,26 have focused much attention on this model.
An exact solution of this model has been derived 27,28,29 , but it requires the solution of Bethe-like equations, which is more costly in terms of computational effort than exact diagonalization. Although the low-energy physics of the model has been widely studied through different approaches (variational 1,30,31 , bosonization 23,32,33 , and coherent states 33,34 ), its high-energy properties have only been very recently investigated numerically 35,36,37 and several interesting features have been revealed. More precisely, for special values of the energy, the spectrum has been shown to display singularities which are reminiscent of the critical point responsible for the well-known quantum phase transition at zero temperature.
In a recent Letter 38 , we proposed a theoretical framework which allows for an exact computation, in the thermodynamic limit, of the LMG model spectrum for the whole range of parameters and leads to a precise description and understanding of these so-called exceptional points. The present paper is an extension of that work, in which we will detail some of its main results and ex-tend the analysis along several directions (relation to the semiclassical treatment, first-order finite-size corrections and expectation values of observables).
The paper is organized as follows. In Sec. II, we introduce the LMG model and the spin-coherent-state formalism 39 , which is the key ingredient of our approach. We then derive the classical energy surface 37 , whose extrema lead to a qualitative phase diagram; these extrema are related to the density-of-states singularities. In a second step, most importantly, we analyze this phase diagram quantitatively. In Sec. III, we introduce the Majorana polynomial and map the timeindependent Schrödinger equation onto a first-order nonlinear (Riccatti-like) differential equation. In Sec. IV, we give solutions of this equation in the thermodynamic limit. This leads to simple expressions of the density of states in the whole phase diagram. In Sec. V, we go beyond this limit and compute the leading finite-size corrections to the density of states. Finally, in Sec. VI, we compute the expectation values (throughout the spectrum) of some spin observables, paying particular attention to one region for which spectral subtleties prevent us from using the Hellmann-Feynman theorem. Some technical details are given in the Appendix.

II. COHERENT-STATE REPRESENTATION AND CLASSICAL ENERGY SURFACE
A. Lipkin-Meshkov-Glick model The LMG model describes a set of N spin-1 2 particles mutually interacting through an (anistropic) XYlike Hamiltonian and coupled to an external transverse magnetic field h. The Hamiltonian of this system can be expressed in terms of the total spin operators S α = N i=1 σ i α /2 where the σ α 's are the Pauli matrices: In the following, for simplicity, we only consider the maximum spin sector s = N/2 with N even. Given the symmetry of the spectrum of H, we focus on the parameter range h 0; |γ y | γ x . Note also that H, S 2 = 0 and H, e iπ(Sz−s) = 0 (spin-flip symmetry). In the standard eigenbasis {|s, m } of S 2 and S z , this latter symmetry implies that odd-and even-m states decouple. In the thermodynamic limit, both subspaces are isospectral so that we further limit the following analysis to the (s + 1)dimensional sector with m even. It is known that H exhibits a quantum phase transition for h = γ x or h = γ y .

B. Coherent-state representation of the spin operators
To determine the spectrum of the Hamiltonian H, it is convenient to use a spin-coherent-state representation 39 . Let us denote by {|s, m } the standard eigenbasis of S 2 , S z with eigenvalues s(s + 1) and m, respectively. The unnormalized spin coherent state |α is then defined as The scalar product of two such states is whereᾱ is the complex conjugate of α. These coherent states obey the following closure relation: dᾱdα π (2s + 1) where dᾱdα = dRe(α) dIm(α). In this representation, a quantum state Ψ(α) = α|Ψ is a polynomial in α, and the action of the spin operators on Ψ translates into differential operators: where S ± = S x ± iS y . We shall discuss below the representation of Ψ(α) in terms of its zeros (the Majorana representation).

C. Classical energy surface
In the thermodynamic limit, a variational description of the ground state 1,30,31 , built with respect to the |α states, leads to the dominant behavior of the model and, in particular, the location of the quantum phase transition. The latter can obtained from an analysis of the minima of the variational energy H 0 : Note that, in this limit, a classical spin description is valid, such that the correspondence between a state |α and a classical vector is simply obtained via a stereographic map from the complex plane onto the S 2 sphere [with α = e iθ tan(φ/2)], leading to the parametrization Here we shall first be interested in the geometrical properties of the whole classical energy surface H 0 (ᾱ, α). Its extrema, obtained by imposing ∂ᾱH 0 = ∂ α H 0 = 0, are given in Table I together with the corresponding energy. When one further imposes that α andᾱ be complex conjugate, the configuration space (spanned by the Hamiltonian parameters) is split into distinct regions characterized by the number of extrema and saddle points in H 0 (ᾱ, α).
This phase diagram coincides with that derived from the analysis of density of states singularities, as done in the next section. We shall describe below how far the classical analysis can help in understanding the spectral results. Note that a related analysis of the classical energy surface, including comparisons to numerically derived spectra, has already been proposed by Castaños et al. 37 in terms of the (θ, φ) angles instead of the present (ᾱ, α).

D. Classical description of the phase diagram
The zero-temperature phase diagram of the LMG model is usually discussed in terms of its groundstate properties. In this case, only two phases are distinguished 1,24,31 . For h > γ x (symmetric phase), the ground state is unique and lim s→∞ S z /s = 1, whereas for h < γ x (broken phase), the ground state is twofold degenerate and lim s→∞ S z /s = h/γ x . Note that the degeneracy in the broken phase arises only in the thermodynamic limit, where the gap between the ground and first excited states vanishes exponentially with s. The quantum phase transition at h = γ x is of second order and characterized by mean-field critical exponents 31 as well as nontrivial finite-size scaling behavior 23,24,25 .
We have shown 38 that, when considering the full spectrum, four different zones arise instead of two, corresponding to a splitting of the broken phase region into three distinct parts characterized by different singularities in the density of states (see Fig. 1).
Note that such singularities have already been pointed out in the numerical study of the special case γ x = −γ y 35,36 and were called "exceptional points." We emphasize in the present study that these exceptional points are associated with saddle points of the energy surface. Of course, the absolute minimum (maximum) gives the lower (upper) bound of the spectrum. Note that these bounds may be degenerate.
In the thermodynamic limit, to a given energy in the spectrum corresponds a level set on H 0 (ᾱ, α). At that energy, the Husimi function local maxima (defined in the next section) are known to concentrate along this level set, which forms the classical orbit. Singularities of the surface (maxima, minima, or saddle points) translate into singularities of the level sets (a main ingredient in Morse surface theory). This, in turn, affects the density of states computation, as illustrated in the next section, and explains why the singularities in the H 0 (ᾱ, α) surface and in the density of states are in close correspondance.
As an illustration, we display in Fig. 2 the classical energy surface for (γ x = 5, γ y = 3, h = −1), which is precisely the point of zone IV whose density of states is shown in Fig. 1. As can be seen, the density of states contains two different types of singular points, being the locus of either a divergence or discontinuity. The analysis of the classicalenergy surface allows one to qualitatively understand all these features. Indeed, it contains two absolute minima (denoted m) which provide the lower bound of the spectrum (twofold-degenerate ground-state energy); two saddle points (denoted s) corresponding to the singular behavior of density of states; one local maximum (denoted M ) which is associated with the discontinuity, and one absolute maximum, not shown here, giving the upper bound of the spectrum.
The same geometrical analysis can be performed throughout the configuration space. A typical classical surface in zone I displays one minimum and one maximum, which, respectively, signal the lower and upper edges of the spectrum. A zone-II surface has two absolute minima (corresponding to the broken-phase degenerate ground states), a saddle point (corresponding to the density-of-states singularity), and one maximum (the upper spectrum edge). Finally, a generic zone-III surface has (again) two absolute minima, two saddle points (corresponding to the two singularities in the spectrum, arising at different energies), and two absolute maxima (corresponding to a degenerate upper state). Note that, when displayed on the sphere, one recovers the standard result for surfaces singularities, which states that the number of maxima plus the number of minima minus the number of saddle points equals the genus of the sphere, i. e., 2.
Thus, the analysis of the classical energy surface allows us to qualitatively describe the phase diagram shown in Fig. 1. However, it does not give any quantitative information concerning the density of states. The aim of what follows is to develop a reliable method to exactly compute the full spectrum of the LMG model.

A. Majorana polynomial and Majorana sphere
The first step consists in analysing the eigenstates in the spin-coherent-state formalism. Any |Ψ can be represented by its Majorana polynomial 40 defined as where d 2s is the degree of this polynomial in α (d = 2s for a generic state). The roots α k of Ψ(α) fully characterize the normalized quantum state |Ψ up to a global phase.
It may be more convenient to represent such a state |Ψ on the so-called Majorana sphere, which can be seen as a generalization of the celebrated Bloch sphere used for spin-1 2 states. To do so, one first complements the d roots Ψ(α) with 2s − d roots at infinity in the complex plane. Next, the resulting set of 2s complex numbers α k is mapped onto 2s points on the unit sphere by an inverse stereographic map. For instance, the basis states |s, m are represented by s − m points on the north pole and s + m points on the south pole. Less trivial examples can be found in Fig. 3 for eigenstates of H in the zone III.
Let us also introduce G(α), the logarithmic derivative of Ψ(α) The 1/2s factor is here to ensure that G is well behaved at the (infinite-s) thermodynamic limit. Let us also define the Husimi function associated with a general state Ψ(α), We shall further need to locate the Husimi function extrema, which are easily found to satisfy As explained above, for the Hamiltonian eigenstates, these maxima converge at the thermodynamic limit, toward the semiclassical orbits, which are the level sets of the classical energy surface H 0 .

B. From Schrödinger to Riccati
Let us now write the time-independent Schrödinger equation H|Ψ = E|Ψ in the coherent-state representation. Using relations (5), (6), and (7), one transforms the Schrödinger equation into the linear differential equation where ε = E/s and The next step consists in converting the linear secondorder differential equation (18) for Ψ into a nonlinear first-order differential equation for its logarithmic derivative G(α), which satisfies the following Riccati-like equation C. Density of states and poles of G The density of states is then obtained from the analysis of the poles of the function G. To illustrate the poles location, several typical states are displayed in Fig. 3 on the Majorana sphere. Each dot represents one pole of G, i. e., one Majorana zero α k , which is mapped from the complex plane to the sphere by an inverse stereographic projection.
The cornerstone of this study is that, for the LMG model, the α k 's spread over two curves C 0 and C 1 in the complex plane. In addition, the nth excited state of H has 2n poles on C 1 and 2(s − n) on C 0 (thus defining both curves). This remarkable property stems mainly from existing maps (which may differ between parameter space regions) between the LMG model and the problem of a particle in an effective one-dimensional potential (see the Appendix). In the latter case, the oscillation theorem indexes the excited states by the number of wave-function nodes on the real axis. This leads here to (at least one set of) zeros lying on simple lines in the complex plane, where the pole density varies monotoneously with energy.
Let us consider the normalized integrated density of states, N (ε) ∈ [0, 1]. We shall enumerate by n the eigenstates of increasing energy, starting from n = 0 for the ground state to n = s for the highest-energy state. The special location of the G poles leads to a simple relation between N (ε) and p, the number of poles lying in C 1 ,  (39) and (44), respectively] in the thermodynamic limit.
which reads where C 1 is a contour that surrounds C 1 and oriented such that N ≥ 0. For the sake of simplicity, we shall further consider the density of poles in C 1 , called I ∈ [0, 1], which simply reads In general, Eqs. (25) and (22) cannot be exactly solved for arbitrary s. The main goal of this paper is to solve these in the thermodynamic limit (s → ∞) and to capture the leading finite-size corrections in a 1/s expansion.

Leading-order expansion for G
Let us assume that G and ε can be expanded in the form At leading order (1/s) 0 , Eq. (22) becomes a second-order polynomial equation for G 0 whose solutions are The four roots of Q, ±r ± , are branch points of G 0 . The integrated density of states in the thermodynamic limit, N 0 (ε 0 ), now reads A natural choice for the G 0 branch cuts is given by the curves C 0 and C 1 , on which the G poles accumulate as s increases. It indeed corresponds to the direction, in the complex plane, for which the quantity computed in Eq. (32) is real at each (infinitesimal) step of the integration. This latter condition was in fact implemented to draw the curves C 0 and C 1 in the different figures.
In the next section, we analyze in detail the four abovementioned different regions in the phase diagram, in terms of N 0 (ε 0 ), its derivative, and the density of states ρ 0 (ǫ 0 ) = ∂ ε0 N 0 (ε 0 ). These quantities are, in most cases, computed as indicated in Eq. (32). It may happen, as noted below, that the C 1 curve has a complex shape, while C 0 is simple. Since the integral over all branch cuts, corresponding to C 0 and to C 1 , sums to unity, we can safely consider the integral over C 0 , instead of the nontrivial one over C 1 , and write N 0 (ε 0 ) as one minus this integral. We also face the case of state degeneracies, with corresponding symmetric or nonsymmetric classical orbits. Each such orbit is considered separately by imposing the analyticity of G 0 in the region containing this orbit, bounded eventually by a closed branch cut on the sphere. The related Ψ(α) is zero along this line and can be considered as vanishing outside the considered region. This corresponds quite well to the (numerically derived) eigenstate in the nonsymmetric case. However, in the symmetric case this description fails to reproduce the exact eigenstates since the latter is generically a linear combination of states located close to the classical orbits.

Analytical expressions of the densities of states
A precise study of the branch cuts C 0 and C 1 allows one to distinguish between five different forms of the density of states (labeled (a), (b), (c), (d), and (e) below) that can be expressed in terms of a complete elliptic integral of the first kind, an incomplete elliptic integral of the third kind, and a complete elliptic integral of the third kind, Π(n|m) = Π(n, π/2|m).
Depending on the Hamiltonian parameters, we have already distinguished between four different zones, following the classical surface singularities. We will now show how these zones are characterized in terms of the density-of-states behavior. Indeed, each time a classical surface singularity (maximum, minimum, or saddle point) is crossed, the level sets (classical orbits or Husimi function local maxima) experience topological changes, as well as the integration contours, leading to a new expression for the integrated density of states. We now detail these different expressions, by describing each zone.
• Zone I: |γ y | < γ x < h. Within this range of parameters (which coincides to the "symmetric phase" discussed in Sec. II D) the spectrum lies in the interval −h ε 0 h and the density of states is a smooth decreasing function of the energy as can be seen in Fig. 1. The distribution of Majorana polynomial roots in this zone is similar to that displayed in Fig. 3(b). In the complex plane, C 0 and C 1 lie in the imaginary and real axes respectively. The integrated density of states is given by • Zone II: |γ y | < h < γ x .
In this region, one must distinguish between two cases − II (a): − ε 0 −h. C 0 coincides with the whole imaginary axis while C 1 is made of two disconnected segments in the real axis as depicted in Fig. 3(a). Here, the integrated density of states reads − II (b): −h ε 0 h. C 0 and C 1 are the same as in zone I and the analytic expression of the density of states is given by Eq. (36).
These two branches (a) and (b) of the density of states diverge at ε 0 = −h. Indeed, the integrated density of states can be simplified into the form and one can check that ρ II 0 (−h) = ∂ ε0 N II 0 (ε 0 )| −h diverges. One can further extract the leading behavior of the density of states near this point to obtain • Zone III: h < −γ y < γ x .
In this region, one must distinguish between three cases − III (a): − . C 0 is made of two disconnected segments on the imaginary axis while C 1 coincides with the whole real axis as depicted on the Majorana sphere in Fig. 3(c). The integrated density of states simply reads where N (a) 0 is given in Eq. (38). In this zone III, the density of states has two singularities at ε 0 = ±h. The integrated density of states for these energies is given by (44) As done in zone II, one can compute the leading behavior of the density of states near these points and one gets For γ x = −γ y , the spectrum is symmetric with respect to ε 0 = 0 and the above expression gives the exact location, in the thermodynamic limit, of the so-called exceptional point observed in Refs. 35,36 where a more complex diverging behavior was conjectured.
• Zone IV: h < γ y < γ x . In this zone the density of states presents three different regions, of type (d), (e), and (b). The curve C 1 is more complex here, while C 0 always lies on a straight line in the complex plane. This is why we choose to integrate . C 0 coincides with the whole imaginary axis while C 1 has two disconnected branches lying symmetrically on the unit circle with respect to the imaginary axes. We are here facing a case where the classical orbits are related by symmetry [see Fig. 4(d)]. One finds, for this region, where This region shows two disconnected classical trajectories not related by symmetry (see Fig. 4), corresponding to two qualitatively different kinds of states which alternate in the spectrum. C 0 comprises two disconnected components lying in the imaginary axis, while C 1 is still complex and, moreover, is different for the two kinds of states. One finds For the critical energy, at the boundary between IV(d) and IV(e), the integrated density of states simplifies to In addition, the density-of-states singular behavior is not symmetrical and reads h. C 0 is simply connected and lies on the imaginary axes. Like in the previous case, C 1 is nontrivial (see Fig. 4). Nevertheless, the expression found for N 0 in this region coincides with that given by Eq. (36).
We now discuss the particular features found in the spectral region IV(e). At ε 0 = −h, the density of states is discontinuous (see Fig. 1), a fact which can be understood already from the topological analysis of the classical surface H 0 . Indeed, the transition from zone (e) to zone (b) corresponds to leaving a local maximum of H 0 (see Fig. 2); therefore, a family of classical orbits no longer contributes to the density of states.
In addition, as opposed to all other regions, the energy difference between two consecutive levels, ∆ (i) = E (i+1) − E (i) , computed for increasing s, does not converge towards the analytical result and, actually, does not converge at all. In region IV(e), ∆ (i) spreads over two branches (+) and (−), depending on the parity of the i, which oscillate without converging as s increases, as can be seen in Fig. 5. In this case, the gap we compute, in the thermodynamic limit, is actually the average gap, namely ∆ 0 (ε 0 ) = 1 2 ∆ (+) (ε 0 ) + ∆ (−) (ε 0 ) . This is clearly to be understood in relation to the existence of two kinds of states alternating in the spectrum. Indeed, when analyzed separately within each set of states (e + or e − ), the computed energy gaps (between levels j and j + 2 in the energy spectrum) converge as s → ∞. In addition, both such gaps converge to twice the value of ∆ 0 (ε 0 ) (otherwise the two kind of states would not alternate as observed numerically). The oscillatory behavior noted in Fig. 5 signals an energy drift (with s) of one set of energy levels with respect to the other. In the central region, one sees a real lack of convergence toward the red line when increasing s, which is the average gap as computed in the thermodynamic limit.

V. FINITE-SIZE CORRECTIONS
In the previous section, we have analyzed the thermodynamic limit of the LMG model spectrum by considering the leading terms in the expansion (26) [order (1/s) 0 ]. We now express the next-order corrections, which have already been shown, at least for the ground state, to display nontrivial scaling properties 23,24,25,31 . For the sake of simplicity, we limit the present analysis to the case γ x = 1, and γ y = 0.

A.
First-order expansion for G Identifying terms of order 1/s in Eq. (22), one obtains the following form for the first-order term of G: G 1 is thus an analytic function of α with poles at ±r − and ±r + whileG 1 has the same branch cuts as G 0 . I(ε) reads, recalling Eq. (25) and developing up to first order, where I 0 (ε) is given in Eq. (32) and where one can rewrite the 1 4 coming from the integration over the poles. For γ x = 1, γ y = 0, one has only zones I and II to consider, which focuses the analysis on only two energy regions. In zones I and II(b) one obtains whereas in region II(a) one finds Now, for all s, we expect that I(ε) = I 0 (ε 0 ), which implies, at order 1/s, I 1 (ε − ε 1 /s) = I 1 (ε 0 ) = 0. This condition allows one to compute the first-order correction to the energy, ε 1 , which is displayed in Fig. 6 (lower left) and compares nicely with the numerical values, already for small values of s (here s = 50).

B. Energy gaps
The gap between two successive levels has already been discussed above in the zone-IV case. At the thermodynamic limit, it generically reads With the analysis done in the previous section, we can now compute finite size corrections to the gap. To first order, we obtain The above derived values of ε 1 allow us to get a closed form for ∆ 1 , which nicely compares to the numerical values, as can be seen in Fig. 6 (lower right) for s = 50.
The ∆ 1 correction is singular at the exceptional points, which are, as discussed in Sec. IV, located at ε 0 = −h. Note that Leyvraz and Heiss numerically found a logarithmic singularity at the exceptional points 25 . A related feature was already observed for the gap between the ground state and the first-excited state 23,24 . In the latter case, a scaling hypothesis led to a derivation of the first-order correction, showing a N −1/3 behavior. Unfortunately, the scaling hypothesis cannot be used here at the exceptional points. We have determined the behavior of the gap in their vicinity; setting η = |h + ε 0 |, one gets Note that the leading term is simply the inverse of ρ 0 , which is given in Eq. (42) and vanishes when η goes to zero.

VI. OBSERVABLE EXPECTATION VALUES
In this section, we discuss the expectation values of spin observables for generic eigenstates of the LMG model. The simplest way to perform such a calculation is to use the Hellmann-Feynman theorem, which relates these expectation values to the partial derivative of the eigenenergies with respect to Hamiltonian parameters. For instance As an illustration, we compare in Figs. 7 and 8 three cases, computed numerically (at finite s) and via the Hellmann-Feynman theorem in the thermodynamic limit, i. e. replacing E by s ε 0 . As expected, one can see an almost perfect agreement, except for zone IV(e) discussed below.
Let us still make use of the semiclassical analysis discussed in previous sections. The expectation value Ψ|Ô|Ψ for an observable O reads 41 , at leading order, where T is the period of the classical orbit with energy ε 0 and α(t) the solution of the classical dynamics equation 42 . Let us focus on the S z case. In zone I, it is maximal for the ground state. Indeed, in that region, H 0 is minimum for α = 0, where the classical orbit degenerates to a single point at which the ground-state amplitude |Ψ(α)| 2 is concentrated. As a result, although this true ground state differs from the simple fully polarized state 21 , S z reaches its maximum value s.
This also occurs in regions II and III, for energies corresponding to the exceptional points. Here, classical orbits display a characteristic "figure-8" shape, with the values of α therefore differing from zero. The saturation effect results in that case from the fact that the period of the orbit diverges, with a vanishingly small classical velocity near α = 0, forcing the expression in Eq. (69) to saturate. In all cases except zone IV(e), this latter computation leads to the same result as that simply obtained from the Hellmann-Feynman theorem.
In zone IV(e), the numerically computed expectation values alternate along two distinct curves, differing from the Hellmann-Feynman result. This corresponds to the already discussed existence, for the same cenergy ε 0 , of two kinds of classical trajectories nonrelated by symmetry (see Fig. 4). For each numerically derived eigenstate, the associated |Ψ(α)| 2 concentrates alternatively near one of the two classical orbits. Integrating separately along each orbit precisely gives the two branches that are observed numerically (Fig. 8), while the Hellmann-Feynman computation leads to an averaged value.

VII. CONCLUSION
We have studied in detail the full spectrum of the Lipkin-Meshkov-Glick model by means of a coherent-states formalism. In a first step, we simply determined the main characteristics of the (zero temperature) phase diagram by analyzing extrema and saddle points of the classical energy surface. This leads us to distinguish between four zones in the phase diagram corresponding to various patterns of the density of states whereas the usual ground-state criterion leads to only two distinct phases.
In a second step, we analyzed more deeply the nature of the eigenstates in terms of their associated Majorana polynomial roots. This enabled us to exactly compute the integrated density of states in the thermodynamic limit as well as the first finite-size corrections. This remarkable result mainly stems from the fact that the roots of the Majorana polynomial lies on well-defined curves, where their density varies monotoneously with the energy. We also clarified the nature of the so-called "exceptional" points in the spectrum.
Finally, we addressed the question of computing generic observable expectation values, in particular when, owing to subtle spectral reasons, the Hellmann-Feynman theorem cannot be used.
In principle, the same type of analysis could be performed for any spin Hamiltonian expressed in terms of single-spin operators (so-called "collective models"). Preliminary investigations of such models with cubic or quartic interactions are currently under study. Another perspective, also presently under investigation, concerns the dynamical properties for evolutions under both fixed and variable Hamiltonian parameters. steps, given first for the case γ y < 0. 2. The unit circle being parametrized by an angle θ, we write Φ(θ) = e −isθ Ψ(e iθ ) for θ ∈ [0, 2π[. 3. Finally, we define a new function φ(x), which satisfies a one-dimensional Schödinger equation and such that part of its spectrum is put in one-to-one correspondance with the original spin spectrum. This is achieved by setting Φ(θ) = e f [x(θ)] φ[x(θ)] where f (x) and x(θ) are chosen to suppress the first-order derivative in the initial Equation (18) for Ψ(α) and to set the "mass" term equal to s. The resulting Schrödinger-like equation for φ(x), describing a particle in a one-dimensional periodic potential, reads Following this procedure, one obtains the effective potential V (x) = 1 2γ y − 2γ x sn(B|γ x /γ y ) 2 h(2s + 1) (γ x − γ y ) sn(B|γ x /γ y ) − h 2 s + (s + 1)γ x γ y cn(B|γ x /γ y ) 2 , (A2) with B = −γ y x + K γ x γ y .
Note that V is periodic with period L = 4 √ −γy K γx γy . The mapping onto a one-dimensional potential and the celebrated node oscillation theorem allows one to sort the eigenstates of increasing energy according to their number of nodes. Clearly, a φ(x) node leads to a Ψ(α) node for the corresponding LMG eigenstate. The first (2s + 1) eigenstates of this HamiltonianH correspond to the eigenstates of the LMG Hamiltonian with the same energy. Note that, since we focus in this paper on the (s+1)-dimensional "even-m" sector, this leads eventually to a node number inceasing by steps of 2 for each new eigenstate.
Typical potentials are shown in Fig. 9, with parameters associated with regions I, II and III of the LMG phase diagram. The LMG spectrum corresponds to the energies lying between the lower (blue) and the upper (red) lines. The qualitative differences between the three regions appear clearly here. Indeed, in region I the particle moves in a single-well potential whereas it is in a double-well potential in region II. In region III, a higher "allowed" energy region appears, with the extended (unbounded) states above the potential barrier. Crossing the latter corresponds to the upper density-of-states singularity discussed in the text. Note, however, that the extended or bounded nature of the eigenstates for this equivalent one-dimensional system does not have a direct translation into the nature of the corresponding eigenstates in the LMG problem. Similar transformations can be achieved for positive γ y but in this case, one must considerH = −e i π 2 Sy He −i π 2 Sy . Note the occurence of the minus sign which maps the high-energy states of the LMG model onto the low-energy states of the particle-problem (and reciprocally). Following steps (2)