Entanglement and U(D)-spin squeezing in symmetric multi-quDit systems and applications to quantum phase transitions in Lipkin-Meshkov-Glick D-level atom models

Collective spin operators for symmetric multi-quDit (namely, identical $D$-level atom) systems generate a U$(D)$ symmetry. We explore generalizations to arbitrary $D$ of SU(2)-spin coherent states and their adaptation to parity (multicomponent Schr\"odinger cats), together with multi-mode extensions of NOON states. We write level, one- and two-quDit reduced density matrices of symmetric $N$-quDit states, expressed in the last two cases in terms of collective U$(D)$-spin operator expectation values. Then we evaluate level and particle entanglement for symmetric multi-quDit states with linear and von Neumann entropies of the corresponding reduced density matrices. In particular, we analyze the numerical and variational ground state of Lipkin-Meshkov-Glick models of $3$-level identical atoms. We also propose an extension of the concept of SU(2) spin squeezing to SU$(D)$ and relate it to pairwise $D$-level atom entanglement. Squeezing parameters and entanglement entropies are good markers that characterize the different quantum phases, and their corresponding critical points, that take place in these interacting $D$-level atom models.


I. INTRODUCTION
The development of quantum technologies partially relies on the efficient preparation of nonclassical atomic states and the exploitation of many-body entanglement [1][2][3] and spin squeezing [4], specially to enhance the sensitivity of precision measurements like in quantum metrology. Such is the case of many-body entangled (and spin-squeezed) states of cold atoms generated for instance in atom-atom collisions in Bose-Einstein condensates (BECs) [2].
Indistinguishable particles are naturally correlated due to exchange symmetry and there has been a long-standing debate on whether identical particle entanglement is physical or merely a mathematical artifact (see e.g [5] and references therein). Recent work like [6] shows indeed entanglement between identical particles as a consistent quantum resource in some typical optical and cold atomic systems with immediate practical impact. It can also be extracted and used as a resource for standard quantum information tasks [7]. Moreover, multipartite entanglement of symmetric multi-qubit systems can add robustness and stability against the loss of a small number of particles [5].
Understanding the role of the indistinguishableness of identical bosons and quantum entanglement has been the subject of many recent work (see e.g. [8,9] and references therein). We know that, for N = 2 particles, any quantum state is either separable or entangled. However, for N > 2, one needs further classifications for multipartite entanglement [10]. Many different measurements have been proposed to detect and quantify quantum correlations [3]. We shall restrict ourselves to bipartite entanglement of pure states, where necessary conditions for separability in arbitrary dimensions exist.
In order to quantify entanglement between identical particles we shall follow Wang and Mølmer's [11] procedure, who wrote the reduced density matrix (RDM) of one-and two-qubits, extracted at random from a symmetric multi-qubit state ψ, in terms of expectation values S ψ of collective spin operators S. For pairwise entanglement, the concurrence C (an entanglement measure introduced by Wootters [12]) was calculated for spin coherent states (SCSs) [13,14], Dicke and Kitagawa-Ueda [15] spin squeezed states, together with mixed states of Heisenberg models. Kitagawa-Ueda [15] spin squeezed states are the spin version of traditional parity adapted CSs, sometimes called "Schrödinger cat states" since they are a quantum superposition of weakly-overlapping (macroscopically distinguishable) quasi-classical coherent wave packets. They where first introduced by Dodonov, Malkin and Man'ko [16] and later adapted to more general finite groups than the parity group Z 2 = {1, −1} [17]. In this article we shall introduce U(D) SCSs (denoted DSCSs for brevity) adapted to the parity symmetry Z 2 × D−1 . . . ×Z 2 , which are a D-dimensional generalization of U(2) Schrödinger cats, and we shall refer to them as DCATs for short. In general, parity adapted CSs are a special set of "nonclassical" states with interesting statistical properties (see [18][19][20] for several seminal papers). Parity adapted DSCSs arise as variational states reproducing the energy and structure of ground states in Lipkin-Meshkov-Glick (LMG) D-level atom models (see [21] and later in Sec. V).
In Ref. [22], the concurrence C was related to the spin squeezing parameter ξ 2 = 4(∆S n ⊥ ) 2 /N introduced by [15], which measures spin fluctuations in an orthogonal direction to the mean value n ∝ S with minimal variance. Spin squeezing means that ξ 2 < 1, that is, when the variance (∆S n ⊥ ) 2 is smaller than the standard quantum limit S/2 = N/4 (with S the spin) attained by (quasiclassical) SCSs. This study shows that spin squeezing is related to pairwise correlation for even and odd parity multi-qubit states. Squeezing is in general a redistribution of quantum fluctuations between two noncommuting observables A and B while preserving the minimum uncertainty product ∆ ψ A∆ ψ B ≥ 1 2 | [A, B] ψ |. Roughly speaking, it means to partly cancel out fluctuations in one direction at the expense of those enhanced in the "conjugated" direction. For the standard radiation field, it implies the variance relation (∆q) 2 < 1/4 for quadrature (position q and momentum p) operators. For general U(D) spin systems of identical D-level atoms or "quDits", the situation is more complicated and we shall extend the D = 2 definition of spin squeezing to general D.
Spin squeezing can be created in atom systems by making them to interact with each other for a relatively short time in Kerr-like medium with "twisting" nonlinear Hamiltonians like H = λS 2 x [15], generating entanglement between them [23,24]. This effective Hamiltonian can be realized in ion traps [25] and has produced four-particle entangled states [26]. There are also some proposals for two-component BECs [23]. Likewise, the ground state at zero temperature of Hamiltonian critical many-body systems possessing discrete (parity) symmetries also exhibits a cat-like structure. The parity symmetry is spontaneously broken in the thermodynamic limit N → ∞ and degenerated ground states arise. Parity adapted coherent states are then good variational states, reproducing the energy of the ground state of these quantum critical models in the thermodynamic limit N → ∞, namely in matter-field interactions (Dicke model) of two-level [27,28] and three-level [29,30] atoms, BEC [31], U(3) vibron models of molecules [32,33], bilayer quantum Hall systems [34] and (LMG) models for two-level atoms [35][36][37]. Quantum information (fidelity, entropy, fluctuation, entanglement, etc) measures have proved to be useful in the analysis of the highly correlated ground state structure of these many-body systems and the identification of critical points across the phase diagram. Special attention must be paid to the deep connection between entanglement, squeezing and quantum phase transitions (QPTs); see [1,4] and references therein.
In this article we want to explore squeezing and interparticle and interlevel quantum correlations in symmetric multi-quDit systems like the ones described by critical LMG models of identical D-level atoms (se e.g. [21,[38][39][40][41][42] for D = 3 level atom models). The literature mainly concentrates on two-level atoms, displaying a U(2) symmetry, which is justified when we make atoms to interact with an external monochromatic electromagnetic field. However, the possibility of polychromatic radiation requires the activation of more atom levels and increases the complexity and richness of the system (see e.g. [29]). In any case, this also applies to general interacting boson models [43] and multi-mode BECs with two or more boson species. Collective operators generate a U(D) spin symmetry for the case of D-level identical atoms or D boson species (quDits). Recently [21] we have calculated the phase diagram of a three-level LMG atom model. Here we want to explore the connection between entanglement and squeezing with QPTs for this symmetric multi-qutrit system. For this purpose, we extend to D levels the usual definition of Dicke, parity-adapted SCSs and NOON states, and propose linear and von Neumann entropies of certain reduced density matrices as a measure of interlevel and interparticle entanglement. We also introduce a generalization of SU(2) spin squeezing to SU(D). The organization of the paper is as follows. In Sec. II we introduce collective U(D) spin operators S ij , their boson realization, their matrix elements in Fock subspaces of N symmetric quDits, DSCSs and their adaptation to parity ("DCATs"), and a generalization of NOON states to D-level systems ("NODONs"). In Sec. III we give a brief overview on the concepts and measures of interlevel and interparticle entanglement, considering different bipartitions of the whole system, that we put in practice later in Sec. VI. As entanglement measures, we concentrate on linear (unpurity) and von Neumann entropies. We compute entanglement between levels and atoms for DSCSs, DCAT and NODON states. In Sec. IV we extend Kitagawa-Ueda's definition of SU(2) spin squeezing to SU(D) and we also connect it with two-quDit entanglement introduced in the previous Section. In Sec. V we introduce D-level Lipkin-Meshkov-Glick (LMG) atom models (we particularize to D = 3 for simplicity) and study their phase diagram and critical points. In Sec. VI we analyze the ground state structure of the three-level LMG atom model across the phase diagram with the quantum information measures of Sec. III and the SU(3)-spin squeezing parameters of Sec. IV, thus revealing the role of entanglement and squeezing as signatures of quantum phase transitions and detectors of critical points. We only compute linear entropy in Sec. VI, since it is easier to compute than von Neuman entropy and eventually provides similar qualitative information for our purposes; the interested reader can consult Refs. [44,45] for a more general study on the relation between both entropies. Finally, Sec. VII is devoted to conclusions.

II. STATE SPACE, SYMMETRIES AND COLLECTIVE OPERATOR MATRIX ELEMENTS
We consider a system of N identical atoms of D levels (N quDits in the quantum information jargon). Let us denote by E ij = |i j| the Hubbard operator describing a transition from the single-atom level |j to the level |i , with i, j = 1, . . . , D. These are a generalization of Pauli matrices for qubits (D = 2), namely E 12 = σ + , E 21 = σ − , E 11 − E 22 = σ 3 and E 11 + E 22 = σ 0 (the 2 × 2 identity matrix). The expectation values of E ij account for complex polarizations or coherences between levels for i = j and occupation probability of the level i for i = j. The E ij represent the D 2 step operators of U(D), whose (Cartan-Weyl) matrices l|E ij |k = δ il δ jk fulfill the commutation relations Let us denote by E µ ij , µ = 1, . . . , N the embedding of the single µ-th atom E ij operator into the N -atom Hilbert space; namely, The collective U(D)-spin (D-spin for short) operators are They are the generators of the underlying U(D) dynamical symmetry with the same commutation relations as those of E ij in (1). When focusing on two levels i > j, we might prefer to use (roman i denotes the imaginary unit throughout the article) with commutation relations [J (ij) (and cyclic permutations of x, y, z), which is an embedding of D(D − 1)/2 SU(2) subalgebras into U(D). Although the form (3) for D-spin operators could be more convenient to extrapolate all the D = 2 level machinery to arbitrary D, we shall still prefer the form (2) (at least in this paper), which allows for more compact formulas.
The D N -dimensional Hilbert space is the N -fold tensor product [C D ] ⊗N . The tensor product representation of U(D) is reducible and decomposes into a Clebsch-Gordan direct sum of mixed symmetry invariant subspaces. Here we shall restrict ourselves to the N +D−1 N -dimensional fully symmetric sector (see [21] for the role of other mixed symmetry sectors), which means that our N atoms are indistinguishable (bosons). Denoting by a † i (resp. a i ) the creation (resp. annihilation) operator of an atom in the i-th level, the collective D-spin operators (2) can be expressed (in this fully symmetric case) as bilinear products of creation and annihilation operators as (Schwinger representation) S ii is the operator number of atoms at level i, whereas S ij , i = j are raising and lowering operators. The fully symmetric representation space of U(D) is embedded into Fock space, with Bose-Einstein-Fock basis (| 0 denotes the Fock vacuum) when fixing n 1 + · · · + n D = N (the linear Casimir C 1 = S 11 + · · · + S DD ) to the total number N of atoms. There are other realizations of D-spin operators in terms of more than D bosonic modes (e.g. when each level j contains M degenerate orbitals), which describe mixed symmetries [46][47][48], but we shall not consider them here. Collective D-spin operators (4) matrix elements are given by The expansion of a general symmetric N -particle state ψ in the Fock basis will be written as where is a shorthand for the restricted sum. D-spin operator expectation values (EVs) can then be easily computed as where we have used (6) and where by n ij we mean to replace n i → n i + 1 and n j → n j − 1 in n. Among all symmetric multi-quDit states, we shall pay special attention to U(D) SCSs (DSCSs for short) which are labeled by complex points z = (z 1 , . . . , z D ) ∈ C D . To be more precise, there is an equivalence relation: |z ∼ |z if z = qz for any complex number q = 0, which means that |z is actually labeled by class representatives of complex lines in C D , that is, by points of the complex projective phase space A class/coset representative can be chosen asz = z/z i when z i = 0, which corresponds to a certain patch of the manifold CP D−1 . This is equivalent to chose i as a reference level and set z i = 1. For the moment, we shall allow redundancy in z to write general expressions, although we shall eventually take i = 1 as a reference (lower energy) level and set z 1 = 1 in Section V. DSCSs (multinomial) can be seen as BECs of D modes, generalizing the spin U(2) (binomial) coherent states of two modes introduced by [13] and [14] long ago. For z = e i (the standard/canonical basis vectors of C D ), the DSCS corresponds to a BEC of N atoms placed at level i. 1 If we order levels i = 1, . . . , D from lower to higher energies, the state |e 1 would be the ground state, whereas general |z could be seen as coherent excitations. Coherent states are sometimes called "quasi-classical" states and we shall see in Section V that |z turns out to be 1 Note the difference between Fock states |n 1 , . . . , n D and DSCSs |(z 1 , . . . , z D ) , which are placed inside parentheses to avoid confusion.
a good variational state that reproduces the energy and wave function of the ground state of multilevel LMG atom models in the thermodynamic (classical) limit N → ∞.
Expanding the multinomial (9), we identify the coefficients c n of the expansion (7) of the DSCS |z in the Fock basis as where we have written |z| = ( However, contrary to the standard CSs, they can be orthogonal when z · z = 0. EVs of D-spin operators S ij (coherences for i = j and mean level populations for i = j) in a DSCS are simply written as DSCS non-diagonal matrix elements of D-spin operators can also be compactly written as Similarly, EVs of quadratic powers of D-spin operators in a DSCS state can be concisely written as and their DSCS matrix elements as Note that, for large N , quantum fluctuations are negligible and we have z|S ij S kl |z z|S ij |z z|S kl |z . Otherwise stated, in the thermodynamical (classical) limit we have We shall use these ingredients when computing one-and two-quDit RDMs in the next Section. We shall see that DSCSs are separable and exhibit no atom entanglement (although they do exhibit level entanglement). The situation changes when we deal with parity adapted DSCSs, sometimes called "Schrödinger cat states" (commented at the introduction), since they are a quantum superposition of weakly-overlapping (macroscopically distinguishable) quasi-classical coherent wave packets, as we shall explicitly see below. These kind of cat states arise in several interesting physical situations. As we have already mentioned, they can be generated via amplitude dispersion by evolving CSs in Kerr media, with a strong nonlinear interaction, like the already commented spin-squeezed states of [15]. They exhibit statistical properties similar to squeezed states, with deviations from Poissonian (CS) distributions. Squeezing and multiparticle entanglement are important quantum resources that make Schrödinger cats useful for quantum enhanced metrology [2]. They are also good variational states [14], reproducing the energy of the ground state of quantum critical models in the thermodynamic limit N → ∞. To construct them, we require parity operators defined as They are conserved when the Hamiltonian scatters pairs of particles conserving the parity of the population n j in each level j = 1, . . . , D. It is easy to see that Π j (a † j ) nj | 0 = (−a † j ) nj | 0 , so that the effect of parity operations on number states (5) is Π j | n = (−1) nj | n . Likewise, using the multinomial expansion (9), it is easy to see that the effect of parity operators on symmetric DSCSs |z is then Note that Π −1 i = Π i and Π 1 . . . Π D = (−1) N , a constraint that says that the parity group for symmetric quDits is not In order to define a projector on definite parity (even or odd), we have to chose a reference level, namely i = 1. Doing so, the projector on even parity becomes where we denote the binary string Likewise, the projection operator on odd parity is Π odd = 1 − Π even . Choosing level i = 1 as a reference level is equivalent to choose a patch on the manifold CP D−1 where z 1 = 0; in this way, any coherent state |z is equivalent to the class representative |z/z 1 , due to equivalence relation |z ∼ |z if z = qz with q = 0. Let us simply denote by z = (1, z 2 , . . . , z D ) the class representative in this case. It will be useful, for later use as variational ground states, to define the (unnormalized) generalized Schrödinger even cat state where and we are using b as a shorthand for b∈{0,1} D−1 . It is just the projection of a DSCS on the even parity subspace. The state (20) is a generalization of the even cat state for D = 2 in the literature [16], given by for the class representative z = (z 1 , z 2 ) = (1, α). The shorthand |α = |(1, α) is used in the literature when a class representative (related to highest |e 1 or lowest |e 2 weight fiducial vectors) is implicitly chosen. The squared norm of |2CAT} is simply Note that the overlap 1, α|1, −α = ((1 − |α| 2 )/(1 + |α| 2 )) N N →∞ −→ 0, which means that |(1, α) and |(1, −α) are macroscopically distinguishable wave packets for any α (they are orthogonal for |α| = 1). Likewise, the unnormalized 3CAT is explicitly given by when setting z = (z 1 , z 2 , z 3 ) = (1, α, β) as a class representative. The squared norm is now These expressions can be generalized to arbitrary D as We shall use (23) and (24) in Sections V and VI, when discussing a LMG model of atoms with D = 3 levels. These 3CAT states have also been used in U(3) vibron models of molecules [32,33] and Dicke models of 3-level atoms interacting with a polychromatic radiation field [29,30].
The D-spin EVs on a DCAT state (20) can be now computed and the general expression is where we set (−1) bi = 1 = |z i | for i = 1 (reference level). Similarly, EVs of quadratic powers of D-spin operators in a DCAT state can be concisely written as To finish this Section, let us comment on a generalization to arbitrary D of another prominent example of quantum states that are useful for quantum-enhanced metrology and provide phase sensitivities beyond the standard quantum limit. We refer to Greenberger-Horne-Zeilinger (GHZ) or "NOON" (when considering bosonic particles) states. For D = 2 level systems, NOON states can be written in the Fock state notation (5) as (see e.g. [2]) Using the canonical basis vectors {e 1 , e 2 } of C 2 , we can write the NOON state as a linear superposition of U(2) SCSs with phases e iφ1,2 , which coincides with (28) (up to an irrelevant global phase e iφ1 ) for the relative phase φ = φ 2 − φ 1 . Multi-mode (or multi-level, in our context) generalizations of NOON states have been proposed in the literature (see e.g. [49,50]). In our scheme, this generalization of NOON states (29) to D level systems adopts the following form EVs of linear and quadratic powers of D-spin operators in NODON states can be easily calculated as The computations in this section will be necessary to discuss entanglement and squeezing properties of all these states in the Sections III and IV.

III. ENTANGLEMENT MEASURES IN MULTI-QUDIT SYSTEMS
In this Section we define several types of bipartition of the whole system, computing the corresponding RDMs and entanglement measures for different kinds of symmetric multi-quDit states ψ in terms of linear L and von Neumann S entropies. We start computing interlevel entanglement in Section III A and then (one-and two-quDit) interparticle entanglement in Section III B.

A. Entanglement among levels
For a general symmetric N -particle state ψ like (7), the RDM on the level i is i (ψ) = tr j =i n, n c n c n |n 1 , . . . , n D n 1 , . . . , n D | = n |c n | 2 |n i n i |.
Thus i (ψ) lies in a single boson Hilbert space of dimension N + 1. Its purity is then Here the superscript makes reference to "level", to distinguish it from "atom" purity P a in the next section. It can also make reference to entanglement between different boson species , like rotational-vibrational entanglement [32,33] in algebraic molecular models [51,52] such as the vibron model based on a bosonic U(3) spectrum-generating algebra [53,54]. For the case of the DSCS |z in (9), taking the coefficients c n in (10), and after a lengthy calculation, the RDM on level i turns out to be diagonal Note that the eigenvalues λ n can be expressed in terms of only two positive real coordinates (x i , y i ), except for the reference level z i = 1, for which x i = 1 and therefore there is only one independent variable y i = |z| 2 − 1. For example, for U(3) SCSs, choosing i = 1 as the reference level and using the parametrization z = (1, α, β) for the phase space CP 2 in (23), we have x 1 = 1, x 2 = |α| 2 , x 3 = |β| 2 and y 1 = |α| 2 + |β| 2 , y 2 = 1 + |β| 2 , y 3 = 1 + |α| 2 . The purity of i (z) is simply P i (x i , y i ) = N n=0 λ 2 n (x i , y i ). In Figure 1 we represent the Linear and von Neumann entanglement entropies for a general level i as a function of (x, y) [for the reference level i = 1, we have to restrict ourselves to the cross section x = 1]. We normalize linear and von Neumann entropies so that their interval range is [0, 1], the extremal values corresponding to pure and completely mixed RDMs, respectively. We shall see that both entropies, L and S, provide similar qualitative behavior for the bipartitions studied in this paper. Interlevel isentropic contours correspond to the straight lines y = mx (see Figure 1), and the maximum is attained for y = x.
The large N behavior of the interlevel linear entanglement entropy L i (z) around the maximum y = x is Figure 2 we represent contour plots of L 1,2 and S 1,2 for the RDM of a 3CAT of N = 20 qutrits on levels i = 1 and i = 2. Note that linear and von Neumann entropies display a similar structure. We omit L 3 and S 3 since they are just the reflection in a diagonal mirror line of L 2 and S 2 , respectively. L 1 attains its maximum at the isentropic circle |α| 2 + |β| 2 = 1, whereas L 2 attains its maximum at the isentropic hyperbola |α| 2 − |β| 2 = 1. The large N behavior of the interlevel linear entanglement entropy for a DCAT around the maximum Figure 2 also shows (in magenta color) the parametric curve (α(λ), β(λ)) obtained later in Section V and related to the stationary points (59) of the energy surface (57) in the quantum phase diagram of a LMG 3-level atom model, where λ is the atom-atom interaction coupling constant. For high interactions we have (α(λ), β(λ)) λ→∞ −→ (1, 1), which does not lie inside the maximum isentropic curve, although the difference between both entropies tends to zero in the limit N → ∞ (see later in Sections V and VI for more information).
For NODON states (30), the RDM on level i and its purity are

B. Entanglement among atoms
We compute the one-and two-particle RDMs for a single and a pair of particles extracted at random from a symmetric N -quDit state. The corresponding entanglement entropies are expressed in terms of EVs of collective D-spin operators S ij .

One-quDit reduced density matrices
Any density matrix of a single quDit can be written as a combination of Hubbard matrices E ij with commutation relations (1) as with r ij complex numbers (the generalized Bloch vector) fulfillingr ij = r ji and (the generalized Bloch sphere) We want to construct the one-quDit RDM for one quDit extracted at random from a symmetric N -quDit state ψ. The procedure consists of writing the one-quDit RDM entries in terms of expectation values (EVs) of collective D-spin operators (2). Remember the definition of E µ ij , µ = 1, . . . , N after (1) as the embedding of the single µ-th atom E ij operator into the N -atom Hilbert space. Atom indistinguishableness implies that E µ ij = 1 N S ij , for any µ = 1, . . . , N , and therefore the one-quDit RDM of any normalized symmetric N -quDit state ψ like (7) can be expressed as with D-spin EVs (8). Note that tr[ρ 1 (ψ)] = 1 since D i=1 S ii = N (total population of the D levels), which is related to the linear Casimir operator C 1 = D i=1 S ii of U(D). From the condition tr(ρ 1 (ψ) 2 ) ≤ 1 we obtain the general relation which could be seen as a measure of the fluctuations or departure from the linear C 1 and quadratic C 2 Casimir operators given by C 1 = N 1 and The quantum limit N 2 in (40) is attained for DSCSs. Indeed, for |ψ = |z in (9), the DSCS operator EVs were calculated in (12). Therefore, the purity of the corresponding one-quDit/atom RDM is simply [we denote interatom purity by P a to distinguish from the interlevel purity P discussed in the previous section] which means that there is not entanglement between atoms in a DSCS . This is because a DSCS is eventually obtained by rotating each atom individually. The situation changes when we deal with parity adapted DSCSs or "Schrödinger cat states" (20). Indeed, the one-quDit RDM ρ 1 (DCAT) does not correspond now to a pure state since, using the D-spin EVs on a DCAT state (26), the purity gives That is, unlike |z , the Schrödinger cat |DCAT is not separable in the tensor product Hilbert space [C D ] ⊗N . In Figure  3, we represent contour plots of linear and von Neumann   1), which means that highly coupled atoms are maximally entangled in a cat-like ground state (see later in Sections V and VI for more information).
To finish this Section, let us comment on one-quDit entanglement for NODON states (30). Taking into account the D-spin EV (31), the one-quDit RDM of a NODON is simply ρ 1 (NODON) = 1 D 1 D and its linear entropy L a 1 = 1, implying maximally mixed RDM.

Two-quDit reduced density matrices
Likewise, any density matrix of two quDits can be written as withr ijkl = r jilk complex parameters subject to tr[ρ 2 ] = 1 and 0 < tr[(ρ 2 ) 2 ] ≤ 1. Now we need to express the RDM on a pair of particles, extracted at random from a symmetric state of N D-level atoms, in terms of EVs of bilinear products of collective D-spin operators S. In particular, we have due to indistinguishableness. Therefore, the two-particle RDM of a symmetric state ψ of N > 2 quDits is written as Using the Casimir values (41), one can directly prove that tr[ρ 2 (ψ)] = 1 for any normalized symmetric state ψ. The case D = 2 was considered by Wang and Mølmer in [11]. The procedure is straightforwardly extended to ρ M for an arbitrary number M ≤ N/2 of quDits. The purity of ρ 2 (ψ) can be compactly written as In order to construct the two-particle RDM of a DSCS (9), we need the EVs of quadratic powers (14). With these ingredients, we can easily compute the two-particle RDM of a DSCS (9) which, for large N has the following asymptotic expression The purity of ρ 2 (z) is 1 since z is separable in the tensor product Hilbert space [C D ] ⊗N , as we have already commented. Moreover, one can see that ρ 2 (z) = ρ 1 (z) ⊗ ρ 1 (z). However, the Schrödinger cat (20) is non-separable and has an intrinsic pairwise entanglement. Taking into account the particular estructure of linear (26) and quadratic (27) D-spin operator EVs, The general formula (48) becomes In Figure 4, we represent contour plots of normalized linear and von Neumann entanglement entropies for the two-qutrit RDM ρ 2 (3CAT) of a U(3) Schrödinger cat (23) as a function of the phasespace CP 2 coordinates α, β [they just depend on the moduli]. As for the one-quDit case, they attain their maximum value at the phase-space point (α, β) = (1, 1); however, unlike the one-quDit case, pairwise entanglement entropies do not attain the maximum value of 1 at this point, but L a 2 = 5/6 and S a 2 0.623 for large N . As already commented, variational (spin coherent) approximations to the ground state of the LMG 3-level atom model [discussed later in Section V] recover this maximum entanglement point (α, β) = (1, 1) at high interactions λ → ∞ (as can be seen in the magenta curve).

IV. SU(D) SPIN SQUEEZING: A PROPOSAL
As we have already commented in the introduction, Wang and Sanders [22] showed a direct relation between the concurrence C, extracted from the two-qubit RDM (47) for D = 2, and the SU(2) spin J = (J x , J y , J z ) squeezing parameter introduced by [15], which measures spin fluctuations in an orthogonal direction to the mean value J with minimal variance. Actually, the definition (53) refers to even and odd symmetric multi-qubit states [remember the extension of this concept to multi-quDits after (17)] for which J = (0, 0, J z ) and therefore the orthogonal direction lies in the plane XY. This definition can be extended to even and odd symmetric multi-quDit states for which S ij ∝ δ ij [see e.g. (26) for the case of the even DCAT state]. Using the embedding (3) of D(D − 1)/2 SU(2) spin subalgebras into U(D), and mimicking (53), we can define D(D − 1)/2 spin squeezing parameters ξ ij , i > j for D-spin systems as: We have chosen the normalization factor 1 N (D−1) so that (54) reduces to (53) for D = 2 and so that the total D-spin squeezing parameter is one (no squeezing) for the DSCSs |z in (9). Actually, for DSCSs we have that ξ 2 ij = (|z i | 2 + |z j | 2 )/(|z| 2 (D − 1)) is written in terms of average level populations z|S ii |z = N |z i | 2 /|z| 2 of levels i and j, acording to (12). Therefore, the presence of D-spin squeezing means in general that ξ 2 D < 1. Using the EVs (31) for NODON states (30), the corresponding spin squeezing parameters are ξ ij = 2/[D(D − 1)], which gives ξ 2 D = 1, thus implying that NODON states do not exhibit spin squeezing.
Note that D-spin squeezing parameters ξ ij are constructed in terms of D-spin quadratic EVs, as the two-quDit RDM (47) and its purity (48) do. Therefore, the deep relation between pairwise entanglement and spin squeezing revealed by Wang-Sanders in [22] for symmetric multi-qubit systems is extensible to symmetric multi-quDits in the sense proposed here. In Figure 5 we show a contour plot of the total D-spin squeezing parameter ξ 2 D for the 3CAT. As in previous figures, the magenta curve represents the trajectory in phase space of the stationary points (59) of the energy surface (57) of the three-level LMG Hamiltonian (56) as a function of the control parameter λ. See later around Figure 9 in Section VI for further discussion.

V. LMG MODEL FOR THREE-LEVEL ATOMS AND ITS QUANTUM PHASE DIAGRAM
In this section we apply the previous mathematical machinery to the study and characterization of the phase diagram of quantum critical D-level Lipkin-Meshkov-Glick atom models. The standard case of D = 2 level atoms has already been studied in the literature (see e.g. [37]). We shall restrict ourselves to D = 3 level atoms for practical calculations, although the procedure can be easily extended to general D. In particular, we propose the following LMG-type Hamiltonian (17), since the interaction term scatters pairs of particles conserving the parity of the population n j in each level j = 1, . . . , D. Energy levels have good parity, the ground state being an even state. We divide the two-body interaction in (56) by the number of atom pairs N (N − 1) to make H an intensive quantity, since we are interested in the thermodynamic limit N → ∞. We shall see that parity symmetry is spontaneously broken in this limit.
As already pointed long ago by Gilmore and coworkers [14,55], coherent states constitute in general a powerful tool for rigorously studying the ground state and thermodynamic critical properties of some physical systems. The energy surface associated to a Hamiltonian density H is defined in general as the coherent state expectation value of the Hamiltonian density in the thermodynamic limit. In our case, the energy surface acquires the following form where we have used DSCS EVs of linear (12) and quadratic (14)  E 0 ( , λ) = min α,β∈C E (α,β) ( , λ) is attained at the stationary (real) phase-space values α ± 0 = ±α 0 and β ± 0 = ±β 0 with In Figures 2, 3, 4 and 5 we plotted (in magenta color) the stationary-point curve (α 0 (λ), β 0 (λ)) on top of level, oneand two-qutrit entanglement entropies, and squeezing parameter, noting that (α 0 (λ), β 0 (λ)) → (1, 1) for high λ → ∞ interactions. We will come to this later in Section VI. Inserting (59) into (57) gives the ground state energy density at the thermodynamic limit Here we clearly distinguish three different phases: I, II and III, and two second-order QPTs at λ are discontinuous. In the stationary (magenta) curve (α 0 (λ), β 0 (λ)), the phase I corresponds to the origin (α 0 , β 0 ) = (0, 0) (squared point), phase II corresponds to the horizontal part β 0 = 0 up to the star point, and phase III corresponds to β 0 = 0.
Note that the ground state is fourfold degenerated in the thermodynamic limit since the four U(3) SCSs |z ±± 0 = |1, ±α 0 , ±β 0 have the same energy density E 0 . These four U(3) SCSs are related by parity transformations Π j in (17) and, therefore, parity symmetry is spontaneously broken in the thermodynamic limit. In order to have good variational states for finite N , to compare with numerical calculations, we have two possibilities: 1) either we use the 3CAT (23) as an ansatz for the ground state, minimizing 3CAT|H |3CAT , or 2) we restore the parity symmetry of the U(3) SCS |1, α 0 , β 0 for finite N by projecting on the even parity sector. Although the first possibility offers a more accurate variational approximation to the ground state, it entails a more tedious numerical minimization than the one already obtained in (58) for N → ∞. Therefore, we shall use the second possibility which, despite being less accurate, it is straightforward and good enough for our purposes. That is, we shall use the 3CAT (23), evaluated at α = α 0 and β = β 0 and conveniently normalized (24), as a variational approximation |3CAT 0 to the numerical (exact) ground state |ψ 0 for finite N .

VI. ENTANGLEMENT AND SQUEEZING AS SIGNATURES OF QPTS
The objective in this Section is to use level and particle entanglement and squeezing measures as signatures of QPTs in these LMG models, playing the role of order parameters that characterize the different phases or markers of the corresponding critical points. We restrict ourselves to linear entropy which, as already shown, gives qualitative information similar to von Neumann entropy for this study, with the advantage that it requires less computational resources. As already commented, Refs. [44,45] contain more general information about the relation between both entropies. Linear entropies of one-and two-qutrit RDMs turn also to provide similar qualitative information, although pairwise entanglement shows a more direct relation to spin squeezing.
We have numerically diagonalized the Hamiltonian (56) for N = 50 3-level atoms, and several values of λ (in units), and we have calculated level, one-and two-qutrit entanglement linear entropies for the ground state |ψ 0 = n c n | n , plugging the coefficients c n into (8,33,48). We have also calculated level and atom entanglement linear entropies for the variational approximation |3CAT 0 to the ground state |ψ 0 discussed in the previous Section for N = 50. In Figure  7 we compare numerical with variational ground state entanglement measures between levels i = 1, 2, 3. According to Figure 7, we see that, in phase I, 0 ≤ λ ≤ 1/2, variational results indicate that there is no entanglement between levels, whereas numerical results show a small (but non-zero) entanglement for N = 50. In phase II, 1/2 ≤ λ ≤ 3/2, levels i = 1 and i = 2 get entangled, but level i = 3 remains almost disconnected. In phase III, λ ≥ 3/2, level i = 3 gets entangled too. Interlevel entanglement grows with λ attaining the maximum value of 0.84 at the limiting point (α 0 (∞), β 0 (∞)) = (1, 1) for N = 50. This behavior of the interlevel entropy for the 3CAT variational state can be also appreciated by looking at the stationary (magenta) curve in Figure 2 with relation to the isentropic curves.
Concerning atom entanglement, Figure 8 shows a better agreement between variational and numerical results, showing a rise of entanglement when the coupling strength λ grows across the three phases, attaining values close to the large N maximum values L a 1 = 1 and L a 2 = 5/6 at the limiting point (α 0 (∞), β 0 (∞)) → (1, 1). The entanglement growth is more abrupt between phases I and II than between phases II and III. We see that both, level and atom entanglement measures capture differences between the three phases, even for finite N , and therefore they can be considered as precursors of the corresponding QPT. The main features of the inter-atom entanglement entropy for the 3CAT variational state are also captured by the trajectory of the stationary (magenta) curve in Figures 3 and 4 through the isentropic curves.
In Figure 9 we represent the D = 3 spin total squeezing parameter ξ 2 D (55) of the variational and numerical ground states for N = 50 atoms, as a function of the control parameter λ (in units). The results reveal a clear growth of squeezing at the critical points, the change being more abrupt at these points for the variational (parity adapted mean field) than for the numerical ground state. Note that the variational ground state only shows squeezing (ξ 2 D < 1) at the critical points, whereas the numerical ground state exhibits squeezing for any λ = 0. Looking at the stationary (magenta) curve of Figure 5 we appreciate that it practically lies in regions of no squeezing (in red color) except near the critical points, where squeezing suddenly increases (yellow color regions).

VII. CONCLUSIONS AND OUTLOOK
We have extended the concept of pairwise entanglement and spin squeezing for symmetric multi-qubits (namely, identical two-level atoms) to general symmetric multi-quDits (namely, identical D-level atoms). For it, we have firstly computed expectation values of U(D) spin operators S ij in general symmetric multi-quDit states like: U(D)spin coherent states, their adaptation to parity (Schrödinger DCAT states), and an extension of NOON states to D levels (NODON states). The reduced density matrices to one-and two-quDits extracted at random from a symmetric multi-quDit state exhibit atom entanglement for DCAT states, but not for U(D)-spin coherent states. We have used entanglement to characterize quantum phase transitions of LMG D-level atom models (we have restricted to D = 3 for simplicity), where DCAT states (as an adaptation to parity of mean-field spin coherent states) turn out to be FIG. 8. One-qutrit L a 1 and two-qutrit L a 2 entanglement linear entropies of the ground state of the three-level atom LMG model Hamiltonian (56) as a function of the control parameter λ (in units). Critical points, at which a QPT takes place, are marked with vertical grid lines, whereas horizontal grid lines label the asymptotic values L a 1 → 1 and L a 2 → 5/6 of the entropies. We compare numerical with variational results for N = 50 atoms. a reasonable good variational approximation to the exact (numerical) ground state. We have also proposed an extension of standard SU(2)-spin squeezing to SU(D)-spin operators, which recovers D = 2 as a particular case. We have evaluated SU(3)-spin squeezing of the ground state of the LMG 3-level atom model, as a function of the control parameter λ, and we have seen that squeezing grows in the neighborhood of critical points λ c , therefore serving as a marker of the corresponding quantum phase transition. A deeper study and discussion of squeezing in these models requires a phase space approach in terms of a coherent (Bargmann) representation of states, such as the Husimi and Wigner functions, and it will be the subject of future work. all thank Octavio Castaños for his valuable collaboration in the early stages of this work.