Quantum phases and entanglement properties of an extended Dicke model
QQuantum phases and entanglement propertiesof an extended Dicke model
Michal Kloc a, ∗ , Pavel Str´ansk´y a , Pavel Cejnar a a Institute of Particle and Nuclear Physics, Faculty of Mathematics and Physics, Charles University,V Holeˇsoviˇck´ach 2, Prague, 18000, Czech Republic
Abstract
We study a simple model describing superradiance in a system of two-level atoms interact-ing with a single-mode bosonic field. The model permits a continuous crossover betweenintegrable and partially chaotic regimes and shows a complex thermodynamic and quantumphase structure. Several types of excited-state quantum phase transitions separate quantumphases that are characterized by specific energy dependences of various observables and bydifferent atom–field and atom–atom entanglement properties. We observe an approximaterevival of some states from the weak atom–field coupling limit in the strong coupling regime.
Keywords:
Single-mode superradiance model, Excited-state quantum phase transitions,Thermal and quantum phases, Entanglement properties of excited states
1. Introduction
Since its prediction in 1954 [1], the effect of superradiance has attracted a lot of theoreticaland experimental attention [2–4]. Its basic principle—the fact that a coherent interaction ofan unexcited gas with the vacuum of a common field can create a spontaneous macroscopicexcitation of both matter and field subsystems—appears in various incarnations in diversebranches of physics [5, 6].The Dicke model [1, 7, 8] of superradiance resorts to a maximal simplification of theproblem to capture the main features of the superradiant transition in the most transparentway. The model shows a thermal phase transition, analyzed and discussed in Refs. [9–12], aswell as a zero-temperature (ground-state) Quantum Phase Transition ( qpt ) [13–15], whichwas addressed experimentally and realized with the aid of a superfluid gas in a cavity [16–18].Recent theoretical analyses showed that the model exhibits also a novel type of criticality—the one observed in the spectrum of excited states [19–22]. These so-called Excited-StateQuantum Phase Transitions ( esqpt s) affect both the density of quantum levels as a functionof energy and their flow with varying control parameters, and are present in a wide varietyof quantum models with low numbers of degrees of freedom [23–28]. ∗ The corresponding author; email address: [email protected]ff.cuni.cz
Preprint submitted to Annals of Physics October 1, 2018 a r X i v : . [ qu a n t - ph ] M a r n this work, we analyze properties of a simple superradiance model interpolating betweenthe familiar Dicke [1] and Tavis-Cummings [7, 8] Hamiltonians. A smooth crossover betweenboth limiting cases is achieved by a continuous variation of a model parameter, which allowsone to observe the system’s metamorphosis on the way from the fully integrable (hence atleast partly understandable) to a partially chaotic (so entirely numerical) regime. One ofthe aims of our work is to survey the phase-transitional properties of the extended modeland to investigate the nature of its quantum phases —the domains of the excited spectrumin between the esqpt critical borderlines. A usual approach to characterize different phasesrelated to the ground state of a quantum system makes use of suitable “order parameters”,i.e., expectation values of some observables. We show that an unmistakable characterizationof phases involving excited states is not achieved through the expectation values alone butrather through their different smoothed energy dependences (trends).The second part of our analysis is devoted to the entanglement properties of excitedstates across the whole spectrum and their potential links to the esqpt s and quantumphases of the model. It is known that a continuous ground-state qpt in many models(including the present one) is characterized by a singular growth of entanglement within thesystem, which can be seen as a quantum counterpart of the diverging correlation length incontinuous thermal phase transitions [29–36]. A question therefore appears whether thereexist any entanglement-related signatures of esqpt s. The extended Dicke model is rathersuitable for a case study of this type since it allows one to analyze at once various types ofentanglement—that between the field and all atoms, and that between individual atoms.The plan of the paper is as follows: Basic quantum and classical features of the modelare described in Sec. 2. Thermal and quantum critical properties and a classification ofthermodynamic and quantum phases are presented in Sec. 3. The atom–field and atom–atom entanglement properties are investigated in Sec. 4. Conclusions come in Sec. 5.
2. Extended Dicke model
Consider single-mode electromagnetic field with photon energy ω (polarization neglected)interacting with an ensemble of N two-level atoms, all with the same level energies ± ω / λ and an additional interaction parameter δ (whoserole will be explained later), the Hamiltonian can be written as H = ω b † b + ω N (cid:88) k =1 12 σ kz + λ √ N (cid:34)(cid:0) b † + b (cid:1) N (cid:88) k =1 12 (cid:0) σ k + + σ k − (cid:1) − (1 − δ ) N (cid:88) k =1 12 (cid:0) b † σ k + + bσ k − (cid:1)(cid:35) = ω b † b + ω J z (cid:124) (cid:123)(cid:122) (cid:125) H free + λ √ N (cid:2) b † J − + bJ + + δ b † J + + δ bJ − (cid:3)(cid:124) (cid:123)(cid:122) (cid:125) H int , (1)where operators b † and b create and annihilate photons, while σ k • stands for the respectivePauli matrix with subscript { + , − , z } or { x, y, z } acting in the 2-state Hilbert space of the2 th atom. The part of H denoted as H free represents the free Hamiltonian of the fieldand atomic ensemble, while the part H int constitutes the atom–field interaction with aconveniently scaled strength λ/ √ N . Defining collective quasi-spin operators J • = (cid:80) k σ k • for the whole atomic ensemble (which is possible due to the long wavelength assumption),we rewrite the whole Hamiltonian in the simplified form given in the second line of Eq. (1).The interaction part of the Hamiltonian (1) contains a parameter δ . For δ = 1 weobtain the standard Dicke Hamiltonian [1], in which the interaction is written in the dipoleapproximation (the dipole operator for the k th atom is proportional to σ k + + σ k − = 2 σ kx and thecoupling strength in units of energy is given by λ = ω dN / /(cid:15) ωV / , with the electric dipolemoment matrix element d , vacuum permitivity (cid:15) and cavity volume V ). This Hamiltonianis sometimes simplified by omitting the terms b † σ k + and bσ k − that for very small λ yieldnegligible contributions to the transition amplitudes [7, 8]. The reduced model with δ = 0,in case of N > δ in Eq. (1) varies smoothly within the interval δ ∈ [0 , J = J x + J y + J z with eigenvalues j ( j +1), where j is integer for N even or half-integer for N odd[1]. The full atomic Hilbert space H A is the span of all 2 N possible configurations of atoms,but due to the conservation of J the dynamics can be investigated separately in any of thesingle- j subspaces H j,l A with dimension 2 j +1. The decomposition reads as follows [8] H A = N (cid:78) k =1 H k A (cid:124)(cid:123)(cid:122)(cid:125) C = N (cid:76) j =0 or (cid:18) R j (cid:76) l =1 H j,l A (cid:19) , (2)where l enumerates replicas (their number is R j = [ N !(2 j + 1)] / [( N + j + 1)!( N − j )!]) ofthe space with given j differing by the exchange symmetry of atomic components. In thefollowing, we will investigate thermodynamic properties of the model in the full atomicspace H A , as well as quantum properties in a single- j space H j,l A . The most natural choicein the latter case is the unique ( R j =1) subspace with maximal value j = N/
2, which is fullysymmetric under the exchange of atoms and therefore emphasizes the collective character ofthe superradiance phenomenon. A general- j subspace has a mixed exchange symmetry suchthat only a number N ∗ = 2 j ≤ N of atoms can be excited independently, while excitations ofthe remaining N − N ∗ atoms have to compensate each other (in Ref. [1] the quantum number j is called a “cooperation number of the atomic gas”). The reduced single- j model has onlytwo effective degrees of freedom f , one associated with the bosonic field, the other with theSU(2) algebra of collective quasi-spin operators, hence f = 2 [13, 14]. In contrast, in thefull (all- j ) model the SU(2) algebra of Pauli matrices for each atom brings an independentdegree of freedom, so the whole atom–field system has f = N +1.The Tavis-Cummings Hamiltonian with δ = 0 conserves the sum M (cid:48) = b † b + J z [8]. For3 E j λ (a) (b) (c) Figure 1: Spectra of quantum energies of the j = N/ λ for ω = ω = 1and N = 6. The panels correspond to (a) the Tavis-Cummings limit δ = 0, (b) intermediate regime δ = 0 . δ = 1. any fixed j , the conserved quantity can be written as M = M (cid:48) + j = b † b (cid:124)(cid:123)(cid:122)(cid:125) n + J z + j (cid:124) (cid:123)(cid:122) (cid:125) n ∗ , (3)where n is the number of field bosons and n ∗ the number of excited atoms (taking values n ∗ = m + j ∈ [0 , N ∗ ], where m is the J z quantum number). The solutions of the δ = 0 modelare therefore restricted to any fixed- M subspace H j,lM of the full Hilbert space H = H A ⊗ H F ⊃ H j,l A ⊗ H F = ∞ (cid:76) M =0 H j,lM . (4)Quantity (3) is not conserved in δ (cid:54) = 0 cases, but the full Hamiltonian (1) always conserves“parity” Π = ( − ) M as the sum of atomic and field excitation quanta is varied only in pairs.Figure 1 shows quantum spectra of Hamiltonian (1) with j = N/ λ for three different values of δ . The evolution of the ground state,namely its sudden drop to lower energies above a certain value of λ , indicates the superradi-ant transition at zero temperature. We note that the eigensolutions of the model for a given j are in general obtained by a numerical diagonalization of Hamiltonian (1) in the space H j,l A ⊗ H F from Eq. (4), using the basis | m (cid:105) j,l A ≡ | m (cid:105) A in H j,l A (with m = − j, − j +1 , ..., + j ) and | n (cid:105) F in H F (with n = 0 , , , ... ). Since the diagonalization procedure requires a truncationof H F , the eigensolutions must be checked for convergence.A general state vector for any fixed j is expressed as | ψ (cid:105) = + j (cid:88) m = − j ∞ (cid:88) n =0 α mn | m (cid:105) A | n (cid:105) F , (5)4
20 -10 0 10 200123456 -20 -10 0 10 20-20 -10 0 10 200123456 -20 -10 0 10 20-20 -10 0 10 20 -20 -10 0 10 20 φ x (a)(b) (c)(d) (e)(f) φ x - - - - (g)(h) λ =2.5 λ =1.1 Figure 2: Squared wave functions | ψ ( φ, x ) | of selected eigenvectors of the Hamiltonian (1) with N = 40, ω = ω = 1 for δ = 0 . λ × E values indicated by dots in Fig. 3d. Panels (a)–(f), respectively,correspond to the vertical row of points at the right of Fig. 3d from bottom to top [panel (a) shows theground state], panels (g),(h) to the left pair of points. Note that variable φ is 2 π periodic. where α mn are expansion coefficients satisfying the normalization condition (cid:80) m,n | α mn | = 1.Note that for even- and odd-parity states, respectively, the sums in Eq. (5) go either overeven or odd values of M = n + m + j . Any vector (5) can be visualized as a wave function ψ ( φ, x ), where the angle φ ∈ [0 , π ) comes from the quasi-spin representation of the atomicsubsystem and variable x ∈ ( −∞ , + ∞ ) from the oscillator representation of the bosonicfield. For integer j (even N ), the wave function can be obtained by substitutions (up tonormalization constants) | m (cid:105) A (cid:55)→ e imφ and | n (cid:105) F (cid:55)→ H n ( ω / x ) e − ωx / , with H n standing forthe Hermite polynomial. For half-integer j (odd N ), however, the wave function ψ ( φ, x )acquires a spinorial character. A possible representation on the interval φ ∈ [0 , π ) can beobtained by coupling an integer angular momentum j (cid:48) = j − / N − j (cid:48)(cid:48) = 1 / j as in the case of spinor spherical harmonics.This leads to a mapping of | m (cid:105) A to a 2-valued function ∝ C ± e i ( m ± / φ , which is 2 π -periodic,with C ± denoting the corresponding Clebsch-Gordan coefficients. Examples of squared wavefunctions of the j = N/ N in the medium- δ regime aredepicted in Fig. 2. They are taken at the parameter and energy values indicated by dots inFig. 3d ( δ = 0 . φ (cid:55)→ ( φ − π ) mod 2 π and x (cid:55)→ − x .In the following, we will need the classical limit of the model. It is constructed by a5 E w j -10123-2-3-1012-20 1 2 3 0 1 2 3 0 1 2 3 ( )a =0 d ( )b =0.1 d ( )c =0.2 d ( )d =0.3 d ( )e =0.5 d ( )f =1 d SNTC SNTCD SNTCDSNTCD SND SND l l c 0 = l c l l c l l c l l c l l c f reg Figure 3: (Color online) The λ × E maps of the single- j model with ω = ω = 1 for six values of parameter δ .Tones of blue indicate the classical regular fraction (9). The curves correspond to the ground state energy E (full line) and the critical esqpt energies E c (dotted), E c (dashed), and E c (dash-dotted) from Eqs. (29)–(31). Acronyms d (Dicke), tc (Tavis-Cummings), n (normal) and s (saturated) indicate different “quantumphases” of the model (Sec. 3.2). Dots in panel (d) mark places where wave functions in Fig. 2 are taken. simple operator to c-number mappings according to (see, e.g., Ref. [37])( J x , J y , J z ) (cid:55)→ j (sin θ cos φ, sin θ sin φ, − cos θ ) , (6)( b, b † ) (cid:55)→ √ (cid:0) x + ip, x − ip (cid:1) . (7)The spherical angles φ ∈ [0 , π ) and θ ∈ [0 , π ] determine the orientation of the classicalquasi-spin vector (cid:126)j = ( j x , j y , j z ). Note that angle θ is measured here from the south to thenorth pole, thus a fully de-excited state of atoms ( j z = − j ) corresponds to θ = 0 whilea maximally excited state of atoms ( j z = + j ) is associated with θ = π . It can be shownthat the pair of variables φ and j z = − j cos θ , respectively, represents canonically conjugatecoordinate and momentum associated with the collective states of the atomic subsystem.Similarly, quantities x ∈ ( −∞ , + ∞ ) and p ∈ ( −∞ , + ∞ ) in Eq. (7) are suitable coordinateand momentum of the field subsystem. Let us point out that while the atomic phase spacefor each N is finite (it is a ball with radius j ), the field phase space covers the whole plane,expressing the fact that the number of bosons is unlimited. Note that an alternative routeto the classical limit of the Dicke model was taken in Refs. [13, 14].The classical Hamiltonian resulting from application of Eqs. (6) and (7) in (1) reads as H cl = ω p + x ω j z + λ (cid:114) N ( j − j z ) (cid:20) (1+ δ ) x cos φ − (1 − δ ) p sin φ (cid:21) (8)6see also Ref. [22]). Since the classical limit coincides with j, N → ∞ , the above expres-sion needs an appropriate scaling. It is achieved by introducing size-independent quantities H cl / j and j z / j together with ( x, p ) / √ j . Their application in Eq. (8) leads to an ex-pression for scaled energy containing an effective coupling strength λ eff ≡ λ (cid:112) j/N , whichincreases with j and reaches the bare λ for the maximal j = N/ δ drives the system away from the integrable Tavis-Cummings limit,it is relevant to measure the degree of chaos present at various stages of the transition tothe Dicke limit. We note that chaos in the basic and extended Dicke model was studied alsoin Refs. [13, 14, 38–40]. Classical chaos in a general Hamiltonian (8) can be quantified by aso-called regular fraction f reg . It is defined as the volume of the phase space domain filledwith regular orbits relative to the total phase space volume available at given energy E : f reg ( E ) = (cid:82) dφ dj z dx dp δ ( E − H cl ) χ reg (cid:82) dφ dj z dx dp δ ( E − H cl ) . (9)Here, χ reg ( φ, j z , x, p ) stands for a characteristic function defining the regular part of thephase space, i.e., the domain where classical dynamics is regular in the sense of vanishingLyapunov exponents ( χ reg = 1 in the regular part and χ reg = 0 in the chaotic part). Theway how the regular fraction is obtained from numerical simulation of classical dynamics isdescribed in Ref. [41].The evolution of the classical regular fraction for Hamiltonian (8) with δ increasing from0 to 1 is shown in Fig. 3. Each panel displays a map of f reg for a given δ in the plane λ × E .The value of f reg is encoded into the tones of blue: the full color indicates perfectly regularareas ( f reg = 1) and white completely chaotic ones ( f reg = 0). Not surprisingly, we observethat the degree of chaos exhibits an overall increase with δ . However, even in the Dickelimit δ = 1 the model is not entirely chaotic. The most chaotic domain in all δ > λ above the superradiant transition and energy E exceeding a certainvalue above the ground-state E . Besides the main trends in the dependence of f reg we alsoobserve some surprising fine structures—for instance the “ribs” in panels (d)–(f). M for δ = 0Finding eigensolutions of Hamiltonian (1) is much simpler in the δ = 0 limit [8]. In thatcase, the additional conserved quantity (3) splits the Hilbert space H j,l A ⊗ H F from Eq. (4)into a sum of dynamically invariant subspaces H j,lM . These are spanned by states | m (cid:105) A | n (cid:105) F satisfying constraint n + m = M − j and therefore have dimension d ( j, M ) = min { j, M } + 1.As seen in Fig. 1a, the full δ = 0 spectrum is comprised of a number of mutually non-interacting (crossing each other) spectra with different values of M . Each of these fixed- M spectra is obtained by truncation-free diagonalization in a finite dimension. Using theunperturbed basis | i (cid:105) ≡ | m = − j + i (cid:105) A | n = M − i (cid:105) F enumerated by i = 0 ,..., min { j, M } , weexpress the δ = 0 Hamiltonian in a single M -subspace in the following tridiagonal form: (cid:104) i | H | i (cid:48) (cid:105) = [ ω ( M − i ) + ω ( i − j )] (cid:124) (cid:123)(cid:122) (cid:125) E i (0) δ ii (cid:48) + λ (cid:114) ( i +1)(2 j − i )( M − i ) N (cid:2) δ i ( i (cid:48) − + δ ( i − i (cid:48) (cid:3)(cid:124) (cid:123)(cid:122) (cid:125) H (cid:48) ii (cid:48) . (10)7 igure 4: The j = N/ N = 40) in the M = 2 j subspace for the δ = 0 model with ω = ω = 1(panel a), ω = 2 , ω = 1 (panel b) and ω = 0 . , ω = 1 (panel c). The wave-function entropies S of individualeigenstates in the unperturbed basis | m (cid:105) A | n (cid:105) F are encoded in the varying shade of the corresponding lines. Examples of single- M spectra obtained from this Hamiltonian are given in Fig. 4.The simplest spectrum shown in Fig. 4a corresponds to the tuned case with ω = ω . Theeigenvalues of the tuned Hamiltonian (10) read as E i = E (0)+ λE (cid:48) i , where E (0) = ω ( M − j )and E (cid:48) i are eigenvalues of H (cid:48) ii (cid:48) . Since the latter come in pairs with opposite signs (for odddimensions there is an additional unpaired eigenvalue E (cid:48) i = 0), the spectrum linearly expandswith increasing λ around E (0). The full spectrum in Fig. 1a results just from a pile up ofspectra similar to that from Fig. 4a with all values of M . Moreover, as the eigenvectors ofmatrix (10) with ω = ω coincide with those of H (cid:48) ii (cid:48) , they do not depend on λ . This is alsoverified in Fig. 4a, where the wave-function entropy S ( λ ) of individual eigenstates in theunperturbed basis is encoded in the shade of the respective line. The constancy of S ( λ ) foreach state confirms the absence of its structural evolution. The wave-function entropy for ageneral state (5) with respect to the basis | m (cid:105) A | n (cid:105) F is defined by the expression S ( ψ ) = − j +1) (cid:88) m,n | α mn | ln | α mn | (11)(cf. Ref. [42]), so it measures the degree of delocalization of the given state in the fixed basis.The minimal value S = 0 is assigned to any of the basis states (no delocalization), whilethe maximum S = ln d ( j, M ) / ln(2 j +1) is taken by a uniformly distributed superposition ofall basis states (maximal delocalization). In Sec. 4 we will show that Eq. (11) quantifies theatom–field entanglement contained in the given state for δ = 0.The spectrum of a detuned Hamiltonian with ω (cid:54) = ω is less trivial. Examples of leveldynamics in a single- M subspace and the corresponding wave-function entropies for ω > ω and ω < ω are shown in panels (b) and (c) of Fig. 4. We see that the eigenvalues are no morestraight lines and that the eigenvectors change with λ . The peculiar shapes of these spectrawill be explained in Sec. 3.2, here we only point out that the tuned situation is restoredin an approximate sense for very large λ , when the last term in the detuned Hamiltonian8 = ω ( n + J z ) − λH int / √ N − ( ω − ω ) J z becomes negligible. Therefore, investigating propertiesof the model for large coupling strengths, one can make an assumption that ω ≈ ω .Let us finally return to the classical analysis, now focused specifically on the δ = 0 system.The conservation of quantity M from Eq. (3) introduces particular correlations between thedegrees of freedom associated with atomic and field subsystems. As a consequence, thesystem with f = 2 degrees of freedom is for any fixed value of M reduced to a system withjust a single effective degree of freedom. To make this explicit, we employ a transformation xpφj z (cid:55)→ x (cid:48) = x cos φ − p sin φp (cid:48) = p cos φ + x sin φφ (cid:48) = φ + j z + ( p + x ) / M (cid:48) = j z + ( p + x ) / , (12)where M (cid:48) = M − j is conserved. It can be checked that this transformation is canonical,so ( x (cid:48) , p (cid:48) ) and ( φ (cid:48) , M (cid:48) ) are new pairs of conjugate coordinates and momenta. The ( x (cid:48) , p (cid:48) )variables for each M (cid:48) satisfy a constraint ( x (cid:48) + p (cid:48) ) / ∈ [max { , M (cid:48) − j } , M (cid:48) + j ], so theyform a disc (for M (cid:48) ≤ j ) or an annulus (for M (cid:48) > j ). The classical Tavis-Cummings ( δ = 0)Hamiltonian in terms of the new variables becomes H TCcl = ( ω − ω ) p (cid:48) + x (cid:48) ω M (cid:48) + λ x (cid:48) (cid:118)(cid:117)(cid:117)(cid:116) N (cid:34) j − (cid:18) M (cid:48) − p (cid:48) + x (cid:48) (cid:19) (cid:35) , (13)which for a constant M (cid:48) depends only on ( x (cid:48) , p (cid:48) ), so has effectively f = 1 degree of freedom.This not only guarantees the full integrability of the δ = 0 model, but also results in someemergent critical phenomena which will be discussed below.
3. Phases and phase transitions
The extended Dicke Hamiltonian (1) has a complex phase structure [21, 22, 43]. In thissection, we analyze thermal phase transition of the full model with atomic Hilbert space (2)and the number of degrees of freedom f = N +1. Following the standard approach describedin Ref. [9], we use the Glauber coherent states | α (cid:105) ∝ e αb † | (cid:105) , where α ≡ ( α (cid:48) + iα (cid:48)(cid:48) ) ∈ C is aneigenvalue of the b operator and | (cid:105) denotes the field vacuum. An average number of bosons (cid:104) n (cid:105) α = (cid:104) α | b † b | α (cid:105) = | α | in the state | α (cid:105) represents an indicator of the transition from normalto superradiant phase: (cid:104) n (cid:105) α = 0 in the normal phase and (cid:104) n (cid:105) α ∼ O ( N ) in the superradiantphase. On the way to the thermodynamic limit N → ∞ , it is convenient to perform thescaling transformation α (cid:55)→ ¯ α = α/ √ N , so that the scaled coherent-state variable ¯ α ≡ ¯ α (cid:48) + i ¯ α (cid:48)(cid:48) becomes a suitable complex order parameter of the superradiant phase transition.We start by introducing two special values of the atom–field coupling strength: λ c = √ ωω δ , λ = √ ωω − δ . (14)9hile the first one will be shown to represent a critical strength for the occurrence of thesuperradiant phase, the second value will turn out to be a strength at which a specific formof the superradiant phase occurs for δ ∈ (0 , T (measured in energy units), we minimize in variable ¯ α a scaledfree energy F ¯ α ( T ) /N = − T ln Z ¯ α ( T ) /N , where Z ¯ α ( T ) = (cid:104) ¯ α | Tr A e − H/T | ¯ α (cid:105) is the partitionfunction, calculated for a given ¯ α from a partial trace of the Hamiltonian exponential overthe full atomic Hilbert space (2). For Hamiltonian (1) this procedure yields Z ¯ α ( T ) = e − ω | ¯ α | T Tr N exp (cid:26) − T (cid:18) ω λ [(1+ δ ) ¯ α (cid:48) + i (1 − δ ) ¯ α (cid:48)(cid:48) ] λ [(1+ δ ) ¯ α (cid:48) − i (1 − δ ) ¯ α (cid:48)(cid:48) ] ω (cid:19)(cid:27) . (15)The two eigenvalues of the matrix in the exponential are ± (cid:113)(cid:0) ω (cid:1) + λ (cid:2) (1+ δ ) ¯ α (cid:48) +(1 − δ ) ¯ α (cid:48)(cid:48) (cid:3) ≡ ± e ¯ α , (16)which leads to a simple expression of the scaled free energy: N F ¯ α ( T ) = ω | ¯ α | − T ln (cid:0) e ¯ α T (cid:1) . (17)Stationary points of F ¯ α ( T ) /N can be obtained from a simple analysis. The first derivativein variable ¯ α (cid:48) or ¯ α (cid:48)(cid:48) , respectively, vanishes at the points satisfying the first or second line ofthe following array: ¯ α (cid:48) = 0 or 2 ω e ¯ α = λ (1 − δ ) tanh e ¯ α T , (18)¯ α (cid:48)(cid:48) = 0 or 2 ω e ¯ α = λ (1+ δ ) tanh e ¯ α T . (19)The point ¯ α = 0 is always a trivial solution of both lines, but additional solutions appearfor λ > λ c . This coupling strength represents a critical point where the normal phase ¯ α = 0becomes unstable at the lowest temperatures and a new, superradiant equilibrium is createdat some non-zero values of ¯ α . The critical temperature for the superradiant transition, i.e.,the upper temperature limit for the existence of the ¯ α (cid:54) = 0 solution in the region λ > λ c , is T c = ω − λ c λ . (20)For T > T c , the stable equilibrium of the system is again only the normal solution ¯ α = 0.For δ = 0 (the Tavis-Cummings limit), both right-side conditions in Eqs. (18) and (19)become identical and yield a solution | ¯ α | = const that grows from zero with increasingdifference λ − λ c . Therefore, the superradiant minimum of the free energy forms a circlearound ¯ α = 0. For δ >
0, however, a simultaneous solution of both right-side conditions in(18) and (19) is no more possible. The superradiant equilibrium is then represented by apair of points on the real axis, ( ¯ α (cid:48) , ¯ α (cid:48)(cid:48) ) = ( ± ¯ α (cid:48) , α (cid:48) > δ < λ > λ , a pair of saddle points appears for low temperatures on theimaginary axis, ( ¯ α (cid:48) , ¯ α (cid:48)(cid:48) ) = (0 , ± ¯ α (cid:48)(cid:48) ), where ¯ α (cid:48)(cid:48) ≥ ¯ α (cid:48) solves the right-side equation in (19).These unstable solutions exist for temperatures below T = ω − λ λ . (21)10 F T=10 ' -3-2-10120 3 0 3 -3-2-10120 -4-202 -8-6-4 T=0.2< T T=3< T c T=4.95 T c '''''''' -2 ' -303 '' -3 0 3 ' -303 '' -3 0 3 ' -303 '' FF ' '' -303 -- -3 0 3 ' -303 '' F Figure 5: Thermal phase diagram of the all- j model ( ω = ω = 1, δ = 0 .
3) and sample landscapes of the scaledfree energy at the λ = 2 . T c for the transition to the normal ( n ) phase, while the dotted curve indicates thetemperature T separating the tc and d superradiant phases (with or without saddles of free energy). Thefree energy landscapes in the complex- α plane are, for selected temperatures, visualized by contour-shadeplots (darker areas indicate lower values) and by the cuts along real and imaginary axes. A possibility of thermodynamic quasi-equilibrium states associated with the saddle pointsof free energy (17) was recently discussed in Ref. [43].A thermal phase diagram in the λ × T plane for δ = 0 . T = T c . The critical tem-perature T c from Eq. (20) determines the phase transition between the superradiant andnormal (acronym n ) phases of the model. The superradiant phase exists in two forms: theTavis-Cummings phase (acronym tc ) with the saddles in the free energy landscape andthe Dicke phase (acronym d ) without the saddles. While the tc phase is the only type ofsuperradiant phase in the Tavis-Cummings limit (where λ → λ c and T → T c ), the d phaseis exclusive in the Dicke limit (in which λ → ∞ ). For intermediate δ both phases coexist,being separated by temperature T from Eq. (21). Let us stress that the tc → d transition(in contrast to d → n ) is not a phase transition in the standard sense since it does not affectthe global minimum of the free energy. 11 .2. Ground-state and excited-state quantum phase transitions We return now to the analysis of the extended Dicke model with the restriction to asingle- j collective subspace of atomic states. The number of relevant degrees of freedom ofthe single- j model is f = 2, independently of the size parameter N , which implies that theinfinite-size limit, N → ∞ , coincides with the classical limit [19, 26]. Both ground-stateand excited-state quantum phase transitions can be predicted from the classical version ofthe model, namely from the behavior of stationary points of the classical Hamiltonian (8).In particular, the esqpt s result from singularities (non-analyticities) in the semiclassicaldensity of states (cid:37) cl ( E ) = ∂∂E π ) (cid:90) dφ dj z dx dp Θ( E − H cl ) , (22)where Θ( x ) is the step function (Θ = 0 for x < x ≥ ∇ H cl = 0 (with ∇ standing for the gradient in the phase space).It has been shown [28] that the esqpt s caused by non-degenerate stationary points—those which are locally quadratic as their Hessian matrix of second derivatives has onlynon-zero eigenvalues—can be classified by a pair of numbers ( f, r ), where a so-called indexof the stationary point r is a number of negative eigenvalues of the Hessian matrix. Inparticular, for f = 2, the first derivative of the level density in a vicinity of a stationary-point energy E c behaves as ∂(cid:37) cl ∂E ∝ (cid:26) ( − ) r/ Θ( E − E c ) for r = 0 , , − ) ( r +1) / ln | E − E c | for r = 1 , r even) or a logarithmic divergence (for r odd) at thecritical energy E c .In the present case, the evaluation of the level density according to Eq. (22) can besimplified with the aid of the volume-preserving substitution (cid:18) xp (cid:19) (cid:55)→ ξ = x + (1+ δ ) λω (cid:113) N ( j − j z ) cos φη = p − (1 − δ ) λω (cid:113) N ( j − j z ) sin φ , (24)which transforms the classical Hamiltonian (8) to the form H (cid:48) cl = ω η + ξ ω j z − λ ω j − j z N (cid:0) δ cos 2 φ + δ (cid:1)(cid:124) (cid:123)(cid:122) (cid:125) h cl ( φ,j z ) (25)with fully separated variables of the field and atomic subsystems. Although the mapping (24)does not represent a canonical transformation (so ξ, η is not a new coordinate-momentumpair—otherwise the model would be separable and thus fully integrable for any δ ), it simpli-fies the analysis of the level density. In particular, the integration in Eq. (22) over variables12 - - 10 -
00 - 2- 10 -
00 - 4- 3- 2- 10 -3-2-100 -3-2-100 θ cos φ sin θ φ =0 =0.77 c =1.43 =2 h cl j θ cos φ sin θ φ θ cos φ s i n θφ θ cos φ -3-2-100 θ cos φ sin θ φ - 00 θ cos φ h cl j -3-2-100 θ cos φ sin θ φ θ cos φ s i n θφ s i n θφ h cl jh cl j Figure 6: The function h cl from Eq. (25) in the phase space of the atomic subsystem defined by sphericalangles θ and φ for ω = ω = 1 and δ = 0 .
3. Various panels correspond to the indicated values of thecoupling strength λ . Each contour-shade polar plot (with φ ≡ angle and θ ≡ radius) is accompanied by thecorresponding horizontal and vertical cuts. ξ, η , on which the transformed Hamiltonian depends quadratically, can be performed ex-plicitly (cf. Ref. [26], where an analogous calculation is performed for a Hamiltonian with aquadratic kinetic term). The result is a simplified expression (cid:37) cl ( E ) = 12 πω (cid:90) dφ dj z Θ( E − h cl ) , (26)in which the integration, involving function h cl ( φ, j z ) defined in Eq. (25), goes only over the2-dimensional collective phase space of the atomic subsystem—a ball with radius j . Eq. (26)is proportional to an area of the ball region where h cl takes values less than (or equal to)the chosen energy E , hence it can be visualized as flooding of a landscape with profile h cl on a globe. The h cl function for selected values of λ and δ is depicted in Fig. 6.It is clear from the expression (25) that the stationary points of H cl in all four variablescorrespond to ξ, η = 0 and φ, j z determined as stationary points of the function h cl . Sincethe quadratic minimum in variables ξ, η has a null index, any stationary point of H cl has anindex r equal to that of the corresponding stationary point of h cl . For the determination of esqpt s it is therefore sufficient to find and classify stationary points of the function h cl on aball. Non-degenerate stationary points of h cl with index r = 0, 1, and 2 cause, respectively,an upward jump, logarithmic divergence and a downward jump of ∂(cid:37) cl /∂E .Taking into account that the effective coupling parameter λ eff in the scaled classicalHamiltonian is reduced with respect to actual λ by a factor (cid:112) j/N , see the text belowEq. (8), we obtain from Eq. (14) the following j -dependent values of the critical couplings: λ c ( j ) = (cid:115) N j √ ωω δ , λ ( j ) = (cid:115) N j √ ωω − δ . (27)13 - 4- 3- 2- 1 ε d ρ d ε
010 1.5 3 d ρ d ε -4 -3 -2 -1 0 1 -4 ε (a) d ρ d ε ε (a) (b) ε Figure 7: Energy derivative of the semiclassical level density, ∂(cid:37) cl /∂E , for the single- j model with ω = ω = 1and δ = 0 . λ and ε ≡ E/ ( ω j ). The shade diagram in the left panel was obtained fromEq. (26); darker areas represent larger values and vice versa. Panels (a) and (b) show cuts (full curves)at (a) λ = 1 and (b) λ = 2 . N = 40 and Gaussian smoothening of the spectrum ( σ = 0 .
04 and0 .
07 for λ = 1 and 2 .
5, respectively). Scales in panels (a),(b) are arbitrary and not the same.
These will play the roles of λ c and λ in individual subspaces of states with fixed values of j .Note that in the highest subspace with j = N/ r = 0 represent the global minimum of both h cl and H cl functions, demarcating the ground state of the N → ∞ system. The minimum appearsat j z = − j ( φ arbitrary) for λ < λ c ( j ), and at j z = − jλ c ( j ) /λ , φ = 0 and π (a pair ofdegenerate minima) for λ ≥ λ c ( j ). The ground-state energy is given by the formula E ( j ) = (cid:40) − ω j for λ ∈ [0 , λ c ( j )) , − ω j (cid:104) λ c ( j ) λ + λ λ c ( j ) (cid:105) for λ ∈ [ λ c ( j ) , ∞ ) , (28)which exhibits a second-order qpt from normal to superradiant ground-state phase at thecritical coupling λ c ( j ), where d E /dλ has a discontinuity.(ii) Stationary points with r = 1 represent saddles of h cl . They are located at j z = − j ( φ arbitrary) for λ c ( j ) ≤ λ < λ ( j ), and at j z = − jλ ( j ) /λ , φ = π/ π/ λ ≥ λ ( j ).These stationary points correspond to an esqpt (a logarithmic divergence of ∂(cid:37) cl /∂E ) atthe critical energy E c ( j ) = (cid:40) − ω j for λ ∈ [ λ c ( j ) , λ ( j )) , − ω j (cid:104) λ ( j ) λ + λ λ ( j ) (cid:105) for λ ∈ [ λ ( j ) , ∞ ) . (29)14 -101234 λ ε =0' c =0.5=1.5
10 5 5 1024
10 5 5 1024
10 5 5 1024 TC ω j H cl -10 -5 0 5 10 x' ρ Figure 8: Semiclassical level density in the M = 2 j subspace of the j = N/ δ = 0)with ω = 2 , ω = 1 as a function of λ and ε ≡ E/ω j (shade plot on the left, with dark areas indicatinglarger values and vice versa) and three dependences of the corresponding Hamiltonian function (13) on x (cid:48) for p (cid:48) = 0 (on the right). An esqpt due to an inflection point of the classical Hamiltonian above the qpt critical point (32) results in a logarithmic divergence of the semiclassical level density (the dark band in theshade plot), cf. Fig. 4b. For δ → λ ( j ) → ∞ and Eq. (29) is reduced to its first line.(iii) Stationary points with r = 2 are maxima of h cl at j z = − j ( φ arbitrary) for λ ≥ λ ( j )and at j z = + j ( φ arbitrary) for λ ≥
0. Related esqpt s (downward jumps of ∂(cid:37) cl /∂E ) appearat critical energies E c ( j ) = − ω j for λ ∈ [ λ ( j ) , ∞ ) , (30) E c ( j ) = + ω j for λ ∈ [0 , ∞ ) . (31)The second maximum of h cl at energy (31) is the global one, so for E ≥ E c the formula(26) yields a constant (saturated) value of the level density equal to (cid:37) cl = 2 j/ω .All esqpt critical borderlines E c , E c and E c from Eqs. (29)–(31) for various valuesof δ are demarcated in Fig. 3 above. Their existence is numerically verified in Fig. 7, whichdepicts the λ × E dependence of ∂(cid:37) cl /∂E for j = N/ f = 2 character of our model, the esqpt singularities occur in the first derivative of the level density with respect to energy.However, the conservation of quantity (3) in the δ = 0 limit and the corresponding reductionof the number of effective degrees of freedom to f = 1 (see Sec. 2.2) leads to a possibility togenerate an esqpt singularity in the level density itself .Such an effect can really be identified within a particular single M -subspace of the full15odel. The lowest-energy state of any of such subspaces for N → ∞ coincides with the globalminimum of classical Hamiltonian (13) within the available (for a given M ) domain of thephase space. It turns out that the most interesting M -subspace is the one with M = 2 j .This particular subset of states exhibits for ω > ω a second-order ground-state qpt at thecoupling strength λ equal to a critical value λ (cid:48) c ( j ) = (cid:115) N j ω − ω . (32)At this coupling, the main minimum of function (13) moves from ( x (cid:48) , p (cid:48) ) = (0 ,
0) away to x (cid:48) < ,
0) becomes an inflection point. The inflection point is present for λ ∈ ( λ (cid:48) c ( j ) , ∞ ) andgenerates an esqpt (a logarithmic divergence of the level density in the M = 2 j subspace)at energy E c from Eq. (31). Note that for ω < ω the esqpt appears for λ > | λ (cid:48) c ( j ) | at theupper energy of the unperturbed spectrum.This effect, addressed already in Ref. [19], is demonstrated in Fig. 8. It displays a shadeplot of the semiclassical density of levels with M = 2 j (obtained by the f = 1 phase–spaceintegration) and the underlying forms of the Hamiltonian (13). A finite-size sample of thepredicted singularity in the level density was seen in panel (b) of Fig. 4 above. A moredetailed analysis of Eq. (13) shows that the present type of criticality is absent in the M (cid:54) = 2 j subspaces. This can be intuitively understood from the evaluation of matrix elements H (cid:48) ii (cid:48) ofthe interaction Hamiltonian with δ = 0 in the unperturbed basis, see Eq. (10). The interactionmatrix elements quantify the mixing induced by H int in the unperturbed eigenbasis of H free .Neglecting the trivial state with M = 0 which does not mix at all, the matrix elements H (cid:48) ii (cid:48) are particularly small in the M = 1 subspace (dimension 2) and for the i = 2 j state of the M = 2 j subspace (dimension 2 j +1). Only the latter state can develop a singularity in the j, N → ∞ limit. It is the state | m = + j (cid:105) A | n = 0 (cid:105) F , which for ω > ω represents the loweststate of the M = 2 j subspace, while for ω < ω it is the highest state; cf. panels (b) and (c)of Fig. 4. Although the M (cid:54) = 2 j subspaces show no quantum critical effects, it is interestingto realize that a pile up of all M -subspaces produces the downward jump of ∂(cid:37) cl /∂E (with (cid:37) cl standing for the total semiclassical level density), as observed in the E = E c esqpt ofthe δ = 0 system. The critical borderlines E c , E c and E c in the λ × E plane separate spectral domainsthat we consider to constitute different “quantum phases” of the model (in analogy tothermodynamic phases). They are denoted by acronyms d (Dicke), tc (Tavis-Cummings), n (normal) and s (saturated), see Fig. 3. Both phases d and tc , which coexist for δ inbetween 0 and 1, contain quantum states of a superradiant nature because their energy islowered with respect to the minimal non-radiant value E = − ω j . In the limits of δ , themodel shows only one type of the superradiant phase: tc for the Tavis-Cummings limit δ = 0 and d for the Dicke limit δ = 1. In contrast, both phases n and s contain states thatresemble excitations in the fully non-radiant regime at λ = 0. The s phase above E c3 yieldsa constant, saturated value of level density (cid:37) cl .16 igure 9: Peres lattices of quantum expectation values of observables J z (panel a) and n (panel b) inindividual eigenstates of the j = N/ N = 40, ω = ω = 1, λ = 2 .
5. The upper panelsdepict lattices for the limiting values of δ = 0 and 1, whereas the lower panels correspond to δ = 0 .
3. Theinsets of the δ = 0 . While the definition of esqpt s is obvious from the behavior of quantities like level densityat the critical energies, the meaning of quantum phases in between the critical energies isnot a priori clear. It should be looked for in the structure of the energy eigenstates inthe corresponding energy domains. However, sample wave functions in Fig. 2 [where panels(a), (b), (d), (f) and (h) correspond to phases d , d , tc , n and n , respectively] indicatethat individual eigenstates do not show sufficient clues for the identification of phases. Wetherefore resort to a more efficient visualization tool allowing us a bulk inspection of theeigenstate properties, namely to the method of so-called Peres lattices [44]. It was provenuseful in various models (see Ref. [45] and references therein) including the Dicke model[22, 38]. A general Peres lattice shows expectation values (cid:104) P (cid:105) i = (cid:104) ψ i | P | ψ i (cid:105) of a selectedobservable P in individual energy eigenstates | ψ i (cid:105) (enumerated by integer i ) arranged intolattices with energy E i on one of the axes. Examples of Peres lattices for observables (a) J z = n ∗ − j and (b) n = b † b are shown in the respective panels of Fig. 9. Their comparisonwith the map of chaos in Fig. 3d shows—in agreement with the original conjecture [44]—anoverall correlation of orderly arranged lattice domains with more regular regions of classicaldynamics and vice versa.The δ = 0 . D TC N S
E/( j) j NH int Figure 10: Average slopes of the bunches of neighboring 20 levels for the j = N/ N = 40, ω = ω = 1 and δ = 0 . λ = 2 .
5. Quantum phases d , tc , n and s are distinguished by differentsmoothed energy dependences of the average slopes—see the piecewise fits indicated by full black lines. the model. There exist apparent similarities between parts of both δ = 0 . tc and d phases, respectively, and the corresponding Tavis-Cummings and Dickelattices displayed in the δ = 0 and 1 upper panels. However, a more specific distinction ofquantum phases results from an averaging of both lattices over the neighboring eigenstates,as shown in the insets of both lower panels of Fig. 9. Points in these smoothed lattices represent averages (cid:104) J z (cid:105) i and (cid:104) n (cid:105) i of the respective expectation values (cid:104) J z (cid:105) i and (cid:104) n (cid:105) i over 20neighboring eigenstates. Smoothed dependences of the averaged quantities on energy aregiven by lines, resulting from piecewise fits within the four quantum phases. We observe thatvarious quantum phases are recognized by different characters of the energy dependences,namely: (i) the phase d by slowly descending dependences of both (cid:104) J z (cid:105) i and (cid:104) n (cid:105) i averages,(ii) the phase tc by roughly constant dependences, (iii) the phase n by linearly increasingdependences, and (iv) the phase s by saturated (cid:104) J z (cid:105) i ≈ (cid:104) n (cid:105) i .Peres lattices for both quantities J z and n are connected with the lattice of the interactionHamiltonian H int through the energy conservation. Using the Hellmann-Feynman formula, dE i /dλ = (cid:104) dH/dλ (cid:105) i , we derive the following relation for the slope of individual energy levels dE i dλ = (cid:104) H int (cid:105) i √ N = E i − ω (cid:104) J z (cid:105) i − ω (cid:104) n (cid:105) i λ . (33)The slopes averaged over bunches of neighboring 20 levels, as in the insets of Fig. 9, arepresented in Fig. 10 for the same values of control parameters. Piecewise fits of the energydependences in individual quantum phases are again indicated by lines. We see that any ofthe d , tc , n and s quantum phases is characterized by a specific, roughly invariant energy18ependence of the averaged level slope dE i /dλ within the corresponding energy domain.Abrupt changes of the character of these dependences coincide with critical esqpt energies.This is in accord with a general relation between the density and flow signatures of esqpt s,namely with the fact that a typical esqpt generates the same type of non-analyticity in theenergy dependences of both quantities (cid:37) cl and dE i /dλ [26, 28]. Indeed, the discontinuitiesof ∂ ( dE i /dλ ) /∂E at the critical energies E c and E c , as observed in Fig. 10, are consistentwith the analogous behavior of the ∂(cid:37) cl /∂E (Fig. 7). On the other hand, the anticipatedpoint of a singular growth of the level slope (logarithmic divergence of its energy derivative)at E c is smoothed out in the finite spectrum for a moderate system’s size.
4. Atom–field and atom–atom entanglement
We turn to the study of quantum entanglement properties of individual eigenstate of themodel Hamiltonian and their links to the esqpt s. Our analysis includes two types of bipar-tite entanglement: (a) that between the bosonic field and the set of all atoms (atom–fieldentanglement), and (b) that between any pair of individual atoms (atom–atom entangle-ment). The entanglement of type (a) is an important ingredient of superradiance since theinteraction term of the Hamiltonian (1) carries a direct coupling between the atomic andfield subsystems. In contrast, the entanglement of type (b) appears only due to an indi-rect coupling of individual atoms via the bosonic field, so it may be expected to be just a“higher-order” effect. We start with a brief description of the measures used to quantifyboth types of entanglement.In the following, the atom–field and atom–atom entanglement will be evaluated in indi-vidual eigenstates | ψ i (cid:105) ∈ H j,l A ⊗H F of Hamiltonian (1). Some examples of the eigenstate wavefunctions were shown in Fig. 2. Since the wave-function arguments φ and x , respectively,correspond directly to the atomic and field coordinates, a compound state | ψ (cid:105) is factorizedwith respect to the atom–field partitioning of the system if it has a product wave function ψ ( φ, x ) = ψ (cid:48) A ( φ ) ψ (cid:48)(cid:48) F ( x ), where ψ (cid:48) A and ψ (cid:48)(cid:48) F are arbitrary atomic and field wave functions. Themethod to quantify a departure of a given pure state | ψ (cid:105) from exact factorization, i.e., anamount of atom–field entanglement involved in | ψ (cid:105) , makes use of the von Neumann entropy corresponding to the reduced density operators ρ A and ρ F of the atomic and field subsys-tem, respectively [46]. These operators are obtained by partial tracing of the total densityoperator ρ = | ψ (cid:105)(cid:104) ψ | over the irrelevant part of the total Hilbert space, that is ρ A = Tr F ρ (trace goes over H F ) for the atomic subsystem and ρ F = Tr A ρ (trace goes over H j,l A ) for thefield subsystem. Von Neumann entropy of the pure compound state ρ is by definition zero,but entropies of both reduced density operators ρ A and ρ F satisfy: S A = S F ≥
0. The case S A = S F = 0 implies a separable compound state, while the case S A = S F > ρ to a single subsystem leads to a loss of information on mutual entangle-ment between subsystems. A maximally entangled state yields S A = S F = ln(2 j +1), where2 j +1 is a dimension of H j,l A , the smaller of both subspaces. We therefore define a normalized19tom–field entanglement entropy in a compound state | ψ (cid:105) as S ( ψ ) = − Tr [ ρ A ln ρ A ]ln(2 j +1) = − Tr [ ρ F ln ρ F ]ln(2 j +1) . (34)It changes between S = 0 for separable states and S = 1 for maximally entangled states.Quantifying the atom–atom entanglement, i.e., quantum correlations between a ran-domly chosen pair { k, l } of atoms, is a more complex problem. The use of the above en-tropic approach is disabled by the fact that for any state | ψ (cid:105) of the whole atom–field system,an arbitrary pair of atoms generically occurs in a mixed quantum state. It is known thatmutual entanglement of a pair of objects in a mixed compound state cannot be recognizedby a non-zero entropy of the reduced density operators [47]. Indeed, a mixed compoundstate of the pair generates mixed reduced states of the objects even in absence of entangle-ment. In such cases, the evaluation of an entropy-based measure of entanglement (so-calledentanglement of formation) has to be performed with respect to all possible decompositionsof the compound density operator into statistical mixtures of pure states [48], which is fromthe computational viewpoint a difficult task [49].A way to bypass this obstacle for two-qubit systems was proposed in Refs. [50, 51] interms of a quantity called concurrence . The idea was applied and further elaborated [29–31, 52, 53] for multi-qubit systems in fully symmetric states , like our ensemble of N two-levelatoms with j = N/
2. The atom–atom entanglement in this case is characterized by a scaledconcurrence C ( ψ ) = ( N −
1) max (cid:110)(cid:112) λ − (cid:112) λ − (cid:112) λ − (cid:112) λ , (cid:111) , (35)where λ ≥ λ ≥ λ ≥ λ are eigenvalues (real and non-negative) of a non-Hermitianmatrix ρ kl A ˜ ρ kl A , with ρ kl A denoting a reduced density matrix of the pair of atoms { k, l } and˜ ρ kl A ≡ (cid:0) σ ky ⊗ σ ly (cid:1) ρ kl ∗ A (cid:0) σ ky ⊗ σ ly (cid:1) the corresponding “spin-flipped” (conjugate under the timereversal) state (with ⊗ denoting the tensor product of operators for the selected atoms and ρ kl ∗ A a complex conjugated density matrix in the σ kz , σ lz eigenbasis). The argument ψ in thedefinition (35) again reminds the compound state | ψ (cid:105) in which the concurrence is calculated.It has to be stressed that the full symmetry of | ψ (cid:105) under the exchange of atoms, which isguaranteed in the j = N/ ρ kl A , as well as˜ ρ kl A and C , are the same for any pair of atoms, hence do not in fact depend on k, l .Matrix elements of the reduced two-atom density matrix ρ kl A for any symmetric stateof an N -qubit system can be expressed through expectation values (cid:104) J z (cid:105) , (cid:104) J z (cid:105) and (cid:104) J + (cid:105) ofcollective quasi-spin operators in a given state—see formulas (4) and (11) in Ref. [52]. Thismakes it possible to perform a straightforward numerical computation of C in the j = N/ ξ though therelation C = max { − ξ , } , so it varies between C = 0 for separable states and C = 1for the states that exhibit the maximal entanglement allowed for a given N . An entropycorresponding to the entanglement of formation of a single atomic pair reads as [51] s ( ψ ) = − (cid:88) σ = ± a σ ln a σ , a ± = (cid:20) ± (cid:113) − C ( ψ ) ( N − (cid:21) . (36)20 C δ S λ (a) (b) Figure 11: (Color online) Entanglement properties of the ground state of the j = N/ λ and δ (for ω = ω = 1, N = 40). The red curve indicates the qpt critical coupling λ c .Left: atom–field entanglement measured by the entropy (34). Right: atom–atom entanglement measuredby the concurrence (35). The darker areas indicated larger entanglement. For C = 0 we have s = 0, while for C = 1 we obtain s ≈ ln 2 / ( N − +ln( N − / N − ,which decreases from s = ln 2 for N = 2 to zero in the classical limit N → ∞ . With N (cid:29) s ∼ ln N/N . Note that the N (cid:29) per atom decreases as S/N ∼ ln N/N , see Eq. (34), so it exceedsthe atom–atom entanglement entropy about N times. These scaling properties are verifiedtheoretically and by large- N numerical calculations [32].The atom–field and atom–atom entanglement properties have been studied for the groundstate of the Dicke model [32, 33, 36]. As in other models of similar nature, see e.g. Refs. [29,30, 34, 35], the second-order qpt was shown to induce a singularity in both entanglementmeasures (34) and (35). Figure 11 verifies this conclusion in our extended model withvariable parameter δ . We observe that both S and C measures exhibit an increase atabout the critical coupling λ c from Eq. (27). For δ >
0, the atom–atom entanglement dropsback to nearly zero values with λ > λ c . The atom–field entanglement saturates at a value S ≈ ln 2 / ln( N + 1) for λ (cid:29) λ c , which is due to an irreducible atom–field coupling in thelowest positive parity state in the strong coupling limit [32]. In the δ ≈ λ c . This is due to a specific mechanism, in which the ground state at each λ is formed via unavoided crossings of levels with different values of the conserved quantumnumber M , see Fig. 1a. We will analyze the integrable case in more detail below. δ = 0 case In the following, we present results of a numerical study of the atom–field entanglement inindividual eigenstates of Hamiltonian (1) with j = N/
2. We start by analyzing the integrable21 SE j λ M=8, dim = 9M=0, dim = 1
Figure 12: (Color online) The full energy spectrum of the δ = 0 model with ω = ω = 1 for j = N/ N = 40(left panel), and the atom–field entanglement entropies S in individual eigenstates corresponding to λ = 1 . M -subspaces, the second chain to thesecond states etc. The points corresponding to two selected M -subspaces ( M = 0 and 8) are highlighted. Tavis-Cummings limit δ = 0. This simple setting will allow us to obtain some insight into theentanglement properties from quasi-analytic solutions, which will serve as a useful referencefor the less trivial δ > δ = 0 the Hilbert space splits into thesubspaces H M with fixed values of quantum number M , see Eq. (4). The reduced densitymatrix of the atomic subsystem within each H M reads as ρ A = M (cid:88) n =0 min { M − j,j } (cid:88) m,m (cid:48) = − j α mn M ( m ) α ∗ m (cid:48) n M ( m (cid:48) ) | m (cid:105) A (cid:104) n | n M ( m ) (cid:105) F (cid:104) n M ( m (cid:48) ) | n (cid:105) F (cid:104) m (cid:48) | A = min { M − j,j } (cid:88) m = − j | α mn M ( m ) | | m (cid:105) A (cid:104) m | A , (37)where n M ( m ) = M − j − m is a number of bosons associated with quasi-spin projection m for a given value of M . This expression implies that for δ = 0 the entanglement entropy (34)of any eigenstate is equal to the wave-function entropy (11) corresponding to its expansionin the non-interacting basis.Let us focus first on the tuned case ω = ω . The wave-function entropies for the M = 2 j spectrum of a tuned δ = 0 Hamiltonian were shown in Fig. 4a. We notice that the entropyof individual eigenstates does not change with λ and that it has an apparent symmetrywith respect to the vertical reflection of the spectrum around the centroid energy E (0). The22ormer feature was explained by constancy of eigenstates of the simple Hamiltonian (10), thelatter follows from a recursive relation for the eigenstate components which yields the samedistributions | α mn | , hence the same entropies, for the pair of levels with opposite slopes.The full spectrum of the Tavis-Cummings model is obtained by combining the spectra forall values of M . This is for ω = ω seen in Fig. 12. It shows the energy spectrum as a functionof λ (in the left panel) and the atom–field entanglement entropy of individual eigenstatesfor one particular value of the coupling strength (the right panel). The qpt critical coupling λ c coincides with the point where the energy of the M = 0 ground state is crossed by thelower state from the M = 1 space. A further increase of λ above λ c leads to a sequenceof consecutive level crossings in which the lowest states from subspaces with increasing M become instantaneous ground states of the system. A similar mechanism applies also inthe spectrum of excited states, where we observe a sequence of separated caustic structuresformed by states with increasing ordinal numbers within each M -subspace (the n th causticstructure in the vertical direction represents an envelope of lines corresponding to the n thstates from individual M -subspaces).The entropic spectrum on the right of Fig. 12 arises in a similar way. It is composed ofseveral mutually shifted V-shaped chains of points, each of them containing a collection ofstates from different M -subspaces. For example, the lowest V-chain of points is formed by thelowest-energy states from all M -subspaces, the second lowest chain by second lowest-energystates and so on. To indicate a contribution of a single M -subspace to the entropic spectrum,we highlighted the points corresponding to two values of M . For M >
0, the contributionhas a seagull-like form (see the M = 8 example in the figure), whose reflection symmetryaround E (0) results from the above-discussed properties of Eq. (10). The instantaneousground state at the given λ = 1 .
5, which is well above the critical point λ c , has a high valueof M and therefore carries a considerable amount atom–field entanglement (as in generalexpected in superradiance). On the other hand, the trivial one-dimensional subspace with M = 0 formed by the state | m = − j (cid:105) A | n = 0 (cid:105) F , which was the ground state of the δ = 0Hamiltonian at λ < λ c and becomes an excited state with energy E = E c = − ω j for λ > λ c ,has apparently zero atom–field entanglement. We will see below that some remnants of thisstate remain present in entropic spectra also in the non-integrable δ > detuned case ω (cid:54) = ω . The atom–field entropic spectra of aHamiltonian with δ = 0 and ω = 2 ω are displayed in Fig. 13 for increasing values of thecoupling strength λ (the critical coupling is λ c = √ λ (cid:29) λ c , the detuned situation becomes very similar to the tunedone. This is obvious from a comparison of the λ = 3 . λ (cid:46) λ c , the detuning∆ ω plays an important role. A detuned system shows generally lower values of the atom–field entanglement than the tuned one, which can be attributed to a smaller perturbationefficiency of the interaction Hamiltonian due to larger energy gaps between the unperturbed(factorized) states, see Eq. (10).It is clear that the factorized state | m = − j (cid:105) A | n = 0 (cid:105) F with M = 0 at critical energy E c = − ω j is an invariant eigenstate of any δ = 0 Hamiltonian, irrespective of ω, ω and λ .So the S = 0 point at the place of the unperturbed ground state is present in all entropic23 igure 13: The atom–field entanglement entropy of individual eigenstates of the δ = 0 model with j = N/ N = 40 in the detuned case ( ω = 2 , ω = 1) for several values of the coupling strength λ . spectra of Fig. 13. However, for an ω (cid:54) = ω system, a local decrease (converging to S = 0with j, N → ∞ ) of the atom–field entanglement is present also at the highest critical energy E c = + ω j . This is due to quantum critical properties of the M = 2 j subspace discussedin Sec. 3.2. We consider here the ω > ω case (see Fig. 4b), but the same phenomenon takesplace also in the reversed detuning hierarchy ω < ω (Fig. 4c). For λ < λ (cid:48) c , where λ (cid:48) c is thecritical coupling (32), the lowest state of the M = 2 j subspace at energy E c keeps thefactorized form | m = + j (cid:105) A | n = 0 (cid:105) F (mind the sign difference in m from the other factorizedstate, which means that now a maximal possible, for a given j , number of atoms is in the upper level). Distinct lowering of the atom–field entanglement entropy at E = E c is observedin the λ = 0 . λ (cid:48) c = 0 . qpt -related effect becomessharper with increasing j and N , and in the j, N → ∞ limit it results in a similar S = 0cusp as that at energy E c . For ω < ω , the above factorized state, causing the same effect,is the highest state of the M = 2 j subspace at λ = 0.An analogous lowering of the atom–field entanglement entropy at the upper criticalenergy E c appears also for the coupling strengths above the qpt critical value λ (cid:48) c . The effectin this domain is caused by the f = 1 type of esqpt in the M = 2 j subspace. Indeed, whenexcited states within this subspace cross the critical energy, they get temporarily localized inthe coordinate region around the point x (cid:48) = 0, i.e., at the position of the factorized state. Thisis due to “dwelling” of the E = E c classical trajectory and a semiclassical wave functionat an inflection point of the classical Hamiltonian (13). An analogous phenomenon wasobserved also in other quasi f = 1 models, see e.g. Ref. [54]. The lowering of the eigenstateentropy in the E ≈ E c domain was visible already in panels (b) and (c) of Fig. 4, but it isnot apparent in Fig. 13. To see the effect in the pile up of all- M spectra, one would have togo to higher j, N values, when the localization becomes stronger and the entropy decreases24
01 0 1.5 3 λ E ω j λ S -101 λ (a) (b) (c) log( j ) l o g ( S ) -0.19-0.181 2 3 Figure 14: The atom–field entanglement entropy for three states from the M = 2 j subspace of the δ = 0detuned ( ω = 2 , ω = 1) model with j = N/
2: the ground state (panel a) and excited states which at λ = 0have 15 % (panel b) and 24 % (panel c) of the maximal excitation energy. The dashed curves are for N = 40,the full ones for N = 1000. The energy of the corresponding level (for N = 40) is shown at the bottom ofeach panel. The log-log plot in panel (b) inset shows the minimum value of S in the dip as a function of j for the level at 15 % of the spectrum (for this level, the dependence is roughly S min ∝ /j . ). to its lower limit S = 0.The esqpt -based localization effect is nevertheless well documented in Fig. 14, whichdepicts the evolution of the atom–field entanglement for three selected levels (ground stateand two excited states) from the M = 2 j = N subspace for N = 40 and 1000. We see thatthe passage of the selected excited states through the esqpt energy E c is correlated witha local decrease of the entanglement entropy [sharp dips in the dependences in panels (b)and (c)]. As demonstrated by the curves for two values of N , the drop of entropy becomessharper and deeper with increasing size of the system. The decrease of S at the lowest pointof the dip with j = N/ S min ∝ /N a .The exponent a is rather small and depends on the selected level (e.g., for the displayedlevel at 15 % of the spectrum we obtain a = 0 . a = 0 . λ saturates at avalue close to S ≈ ln(0 . d max ) / ln d max , where d max ≡ d ( j, j ) = 2 j +1 is the dimension ofthe M = 2 j subspace. This expression follows from the use of the random matrix theoryfor estimating an average wave-function entropy for d max orthonormal states in the limit ofstrong mixing [42]. Note also that the critical coupling λ (cid:48) c of the ground-state qpt coincideswith the sharpest increase of the ground-state entanglement entropy (Fig. 14a).25 igure 15: Evolution of the j = N/ δ ∈ [0 ,
1] at λ = 2 . ω = ω , N = 40. The esqpt energies are indicated by vertical lines. The spectrumcontains 5000 well converged eigenstates. - 2 - 1 0 1 20123456 x φφφ φ xxx (a) (b) Figure 16: (Color online) A detailed view of wave functions for (a) the eigenstate closest to the transitionbetween the tc and n phases at λ = 2 . δ = 0 . E = E c eigenstate for δ = 0. The other parameters are the same as in Fig. 15. δ > case Now we are ready to turn our attention to the non-integrable δ > δ ∈ [0 ,
1] for ω = ω and λ = 2 .
5, well above the critical coupling λ c for any δ . The vertical lines demarcate the esqpt energies, the zones between them corresponding to the d , tc , n and s quantum phases. Onecan see that the departure from integrability (an increase of δ ) gradually destroys regularpatterns in the entropic spectra, in analogy with Peres lattices. Orderly patterns in theentropic spectrum occur in the domains of regular dynamics, while disorderly scatteredpoints indicate chaotic domains. Note that the accumulation of points in a narrow bandclose to (slightly below) the maximal entropy S ≈ S ≈ ln 2 / ln(2 j +1), corresponding to states with a good parity quantum number [32].An interesting aspect of the entropic spectra in Fig. 15 is the evolution of the atom–field entanglement in transition between the superradiant and normal phases at the criticalenergy E = − ω j . For δ >
0, this energy corresponds to E c for λ ∈ ( λ c , λ ] and to E c for λ ∈ ( λ , ∞ ). We know that in the δ = 0 limit, there always exists the fully factorizedtrivial eigenstate with M = 0 located right at this energy and a group of weakly entangledeigenstates with low values of M located nearby (see Figs. 12 and 13). To what extent dothese structures survive an increase of parameter δ ? The answer can be read out fromFig. 15. We see that the eigenstates with lowered entropy and even a single state with S ≈ tc phase is present for the chosen valueof the coupling strength λ . This is so if λ > λ , hence δ < − (cid:112) N ωω / jλ , see Eq. (27).If λ ∈ [ λ c , λ ], that is for δ large enough to avoid the existence of the tc phase for a given λ , the decrease of entropy no longer takes place. These observations allow us to say thatthe occurrence of states with lowered atom–field entanglement entropy is a signature of the esqpt between the n and tc phases, but not of that between the n and d phases.The wave function corresponding to the E ≈ E c eigenstate with the lowest entanglemententropy for δ = 0 . =8, dim = 9M=0, dim = 1M=20, dim = 21 C E/( j) Figure 17: (Color online) The scaled concurrence for the δ = 0 model with ω = ω = 1, λ = 2 . j = N/ N = 40. The contribution of selected M -subspaces (the U-shaped chains of points) is highlighted. comparison with the fully factorized state with M = 0 shown in panel (b). Although the finestructures of the δ > δ = 0 case), a distant view shows a great deal of similarity to the M = 0 state.Following our findings for δ = 0, one might expect an analogous decrease of the atom–field entanglement entropy also at the critical energy E c (in transition between the s and n phases) in a detuned ( ω (cid:54) = ω ) system with δ >
0. A numerical evidence of this phenomenon ishowever hindered by the above-presented (cf. the inset in Fig. 14b) slowness of the decreaseof the entanglement entropy dip to zero. To see an effect in the δ > j, N values, which is computationally demanding.
We turn now to a brief numerical analysis of the atom–atom entanglement, which islimited to the symmetric j = N/ δ = 0 and λ > λ c is plotted in Fig. 17. The highlighted points denote states belonging toselected M -subspaces. Like in the entropic spectrum (Fig. 12), the states corresponding thesame M -subspace form a pattern which is always symmetric with respect to the reflectionaround the central energy E (0). Trivially, the M = 0 state | m = − j (cid:105) A | n = 0 (cid:105) F has zeroconcurrence as the atomic state with minimal quasi-spin projection is fully factorized (allatoms in the lower level). So this state shows no entanglement in either atom–field oratom–atom sense. In contrast, some states with low values of M , which all have only weak(increasing with M ) atom–field entanglement, show relatively large (decreasing with M )atom–atom entanglement. 28 igure 18: An evolution of the atom–atom entanglement with increasing value of parameter δ for λ = 2 . ω = ω = 1, j = N/ N = 40. The concurrence spectra contain 5000 well converged states. Consider as an example a doublet of states with M = 1. The scaled concurrence for thesestates is C = ( N − /N , which changes between C = 1 / N = 2 and the maximal value C = 1 for N → ∞ (see the upper pair of points in Fig. 17). As follows from Eq. (10), the M = 1eigenstates are expressed by the superpositions ∝ | m = − j (cid:105) A | n = 1 (cid:105) F ± | m = − j +1 (cid:105) A | n = 0 (cid:105) F ,whose first term is fully factorized in the atomic subspace (as in the M = 0 case) while thesecond term is maximally entangled (a symmetric state with one atom in the upper level29 d =0.3 d =0.5 l c l l c l E w j l Figure 19: Spectra of even-parity states of the tuned j = N/ N = 12 for two values of δ . Thenumber of principal components of an actual eigenstate in the unperturbed basis is encoded in the shades ofgray of the respective line (black means perfect localization, ν i = 1, white complete delocalization, ν i → ∞ ).Revival of some λ < λ c eigenstates at λ ∼ λ is seen in both spectra. and the rest of atoms in the lower level). It turns out that for increasing M >
1, only thestates with the largest positive and negative slopes within the given M -subspace yield anon-vanishing atom–atom entanglement. These states form the “antennae” of the U-shapeddependences of C in each M -subspace. The scaled concurrence of these states decreaseswith M , forming together the left and right chains of C > δ = 0 eigenstates has C = 0.Figure 18 shows an evolution of the concurrence spectrum with increasing δ for a fixed λ > λ c . Apparently, an overall trend of the atom–atom entanglement in individual eigenstatesis a decrease with increasing δ . We stress that the concurrence of an absolute majority ofstates for any δ is zero (the spectrum contains 5000 states) and both the number of stateswith C > C in these states further decrease with δ . The only states inwhich the atom–atom entanglement remains significant even for a relatively large values of δ appear in a narrow energy interval around the esqpt critical energy E c demarcating thetransition between the tc and n phases. These states partly coincide with those in whichwe previously observed a reduced atom–field entanglement, see Fig. 15. For δ large enoughto avoid the existence of the tc phase, the atom–atom entanglement of all states vanishes.The states around the critical energy E c with anomalous atom–field and atom–atomentanglement properties form only a small subset of all states, their appearance is never-theless surprising. To make this clear, let us look at these states from a slightly differentperspective. So far we have fixed a sufficiently large value of λ and changed δ , but now let30s do it the other way round. If δ is fixed at a finite value and λ increases from 0 above λ c and λ , the states with anomalous entanglement first (not later than at λ c ) disappear fromthe spectrum, but above the critical coupling (somewhere in the λ ∼ λ region) some ofthem reappear again! In fact, a systematic analysis of the spectrum discloses approximaterevivals of a number of weak-coupling ( λ (cid:46) λ c ) eigenstates (not only those with anomalousentanglement) in the strong-coupling ( λ (cid:38) λ ) regime. This is illustrated in Fig. 19 where weagain show spectra of the tuned model with j = N/ δ . The fullness (shadeof gray) of each level encodes the number of principal components ν i = 1 (cid:80) i (cid:48) |(cid:104) ψ i (cid:48) (0 + ) | ψ i ( λ ) (cid:105)| (38)of the corresponding Hamiltonian eigenstate | ψ i ( λ ) (cid:105) in the basis of Hamiltonian eigenstates | ψ i (cid:48) (0 + ) (cid:105) taken at λ = 0+ dλ ≡ + infinitesimally displaced from the degeneracy point λ = 0(this basis is therefore different from the | m (cid:105) A | n (cid:105) F basis used above in the evaluation ofthe wave-function entropy). The number of principal components varies from ν i = 1 for aperfectly localized state (identical with one of the λ = 0 + eigenstates) to ν i → dim H = ∞ for totally delocalized states. The revival of several λ = 0 + eigenstates is seen in Fig. 19 asreappearance of certain dark lines in the spectrum for large λ in the tc and n phases. Themechanism underlying this phenomenon is unknown.
5. Conclusions
We have analyzed properties of a generalized Dicke model of single-mode superradianceallowing for a continuous (governed by parameter δ ) crossover between integrable Tavis-Cummings and partly chaotic Dicke limits. We have considered three versions of the model,differing by constraints on the Hilbert space H = H A ⊗ H F . In particular: (i) the all- j modelwith the entire atomic space H A and f = N +1 classical degrees of freedom, (ii) a single- j model with a single quasispin subspace H j,l A ⊂ H A and f = 2, and (iii) a single- j, M modeltaken for δ = 0 in a subspace H M ⊂ H j,l A ⊗ H F ; in this case f = 1.Our first aim was to determine all types of thermal and quantum phase transitions andto characterize various phases of the atom–field system for intermediate values of δ . We haveshown that the thermal phase diagram of the all- j model exhibits a coexistence of the tc and d types of superradiant phase (with and without saddles in the free energy landscape),which contract to a single type tc or d in the limits δ = 0 or 1, respectively. The quantumphase diagram of a single- j model contains the superradiant qpt and three types of esqpt s,characterized by singularities (jumps and a logarithmic divergence) in the first derivativeof the semiclassical level density as a function of energy. The quantum phase diagram ofa critical single- j, M model with M = 2 j shows for ω > ω another second-order qpt andfor ω (cid:54) = ω also another esqpt . The latter is characterized by a singularity (logarithmicdivergence) directly in the level density in the M = 2 j subspace.We have associated four spectral domains in between the esqpt critical borderlines ofthe single- j model with the tc , d , n and s quantum phases of the model. These phasesare characterized by different shapes of a smoothed energy dependence of some expectation31alues, which show abrupt changes at the esqpt energies. A natural choice of the phase-defining quantity is the interaction Hamiltonian, whose expectation value (cid:104) H int (cid:105) i determinesthe slope dE i /dλ of the given level. According to the density–flow relation of the esqpt signatures [28], a smoothed level slope should exhibit the same type of non-analyticity asthe semiclassical level density.Our second aim was to analyze atom–field and atom–atom entanglement properties ofthe model in a wide excited domain. This was first done in the integrable single- j, M model, which allowed for a qualitative explanation of results, and then numerically in thesingle- j model with j = N/
2. We have found that an absolute majority of eigenstates inthe spectrum for any choice of model parameters has large atom–field entanglement butvanishing atom–atom entanglement. Exceptional in this sense for δ = 0 are the eigenstateswith low values of M , which yield a weak atom–field entanglement and simultaneously anincreased atom–atom entanglement. We showed that remnants of these states exist relativelyfar in the non-integrable regime with δ >
0. This concerns in particular the state with M = 0located at the esqpt energy E c , which has zero atom–field entanglement. Another statewith reduced atom–field entanglement entropy, originating in the critical M = 2 j subspace,appears at the esqpt energy E c for ω (cid:54) = ω .Although the above-discussed states with anomalous entanglement properties have beenlocated near the esqpt critical borderlines, we do not regard this connection as system-atic. The esqpt singularity in the spectrum has strong consequences on the structure ofeigenvectors—hence possibly also on their entanglement properties—in the f = 1 case, asindeed observed in the M = 2 j subspace of the present model. However, for f > esqpt sincetransitions between quantum phases affect (as we have seen) the trends of eigenstate varia-tions rather than the eigenstates themselves. Nevertheless, it should be stressed that robust,systematic studies of entanglement in excited states are still rather scarce. We hope thatresults of our analysis will initiate similar studies in other relevant systems. Acknowledgments
We acknowledge discussions with T. Brandes, N. Lambert, J. Hirsch and M.A. Bastarrachea-Magnani, and support of the Czech Science Foundation, project no. P203-13-07117S.
References [1] R.H. Dicke, Phys. Rev. , 99 (1954).[2] M. Gross and S. Haroche, Phys. Rep. , 301 (1982) .[3] M.G. Benedict (ed.), Super-radiance: Multiatomic Coherent Emission , Taylor & Francis, New York,1996.[4] T. Brandes, Phys. Rep. , 315 (2005).[5] J.D. Bekenstein and M. Schiffer, Phys. Rev. D , 064014 (1998).[6] N. Auerbach and V. Zelevinsky, Rev. Prog. Phys. , 89 (1963).[8] M. Tavis and F.W. Cummings, Phys. Rev. , 379 (1968); , 692 (1969).[9] Y.K. Wang and F.T. Hioe, Phys. Rev. A , 831 (1973).
10] K. Hepp and E.H. Lieb, Phys. Rev. A , 2517 (1973).[11] K. Hepp and E.H. Lieb, Ann. Phys. (N.Y.) , 360 (1973).[12] K. Rzazewski, K. Wodkiewicz and W. Zakowicz, Phys. Rev. Lett. , 432 (1975).[13] C. Emary and T. Brandes, Phys. Rev. Lett. , 044101 (2003).[14] C. Emary and T. Brandes, Phys. Rev. E , 066203 (2003).[15] J. Vidal and S. Dusuel, Europhys. Lett. , 817 (2006).[16] F. Dimer, B. Estienne, A.S. Parkins and H.J. Carmichael, Phys. Rev. A , 013804 (2007).[17] K. Baumann, C. Guerlin, F. Brennecke and T. Esslinger, Nature , 1301 (2010).[18] K. Baumann, R. Mottl, F. Brennecke and T. Esslinger, Phys. Rev. Lett. , 140402 (2011).[19] P. P´erez-Fern´andez, P. Cejnar, J.M. Arias, J. Dukelsky, J. E. Garc´ıa-Ramos and A. Rela˜no, Phys. Rev.A , 033802 (2011).[20] P. P´erez-Fern´andez, A. Rela˜no, J.M. Arias, P. Cejnar, J. Dukelsky and J. E. Garc´ıa-Ramos, Phys. Rev.E , 046208 (2011).[21] T. Brandes, Phys. Rev. E , 032133 (2013).[22] M.A. Bastarrachea-Magnani, S. Lerma-Hern´andez and J.G. Hirsch, Phys. Rev. A , 032101, 032102(2014).[23] P. Cejnar, M. Macek, S. Heinze, J. Jolie and J. Dobeˇs, J. Phys. A , L515 (2006).[24] M.A. Caprio, P. Cejnar and F. Iachello, Ann. Phys. (N.Y.) ,1106 (2008).[25] P. Cejnar and P. Str´ansk´y, Phys. Rev. E , 031130 (2008).[26] P. Str´ansk´y, M. Macek and P. Cejnar, Ann. Phys. , 73 (2014).[27] P. Str´ansk´y, M. Macek, A. Leviatan and P. Cejnar, Ann. Phys. , 57 (2015).[28] P. Str´ansk´y and P. Cejnar, Phys. Lett. A , 2637 (2016).[29] A. Osterloh, L. Amico, G. Falci and R. Fazio, Nature , 608 (2002).[30] J. Vidal, G. Palacois and R. Mosseri, Phys. Rev. A , 022107 (2004).[31] J. Vidal, R. Mosseri and J. Dukelsky, Phys. Rev. A , 054101 (2004).[32] N. Lambert, C. Emary and T. Brandes, Phys. Rev. Lett. , 073602 (2004).[33] N. Lambert, C. Emary, and T. Brandes, Phys. Rev. A , 053804 (2005).[34] T. Barthel, S. Dusuel and J. Vidal , Phys. Rev. Lett. , 220402 (2006).[35] J. Vidal, S. Dusuel and T. Barthel, J. Stat. Mech. , P01015 (2007).[36] L. Bakemeier, A. Alvermann and H. Fehske, Phys. Rev. A , 043821 (2012).[37] L. Bakemeier, A. Alvermann and H. Fehske, Phys. Rev. A , 043835 (2013).[38] M.A. Bastarrachea-Magnani, B. L´opez-del-Carpio, S. Lerma-Hern´andez and J.G. Hirsch, Phys. Scr. A , 068015 (2015).[39] M.A. Bastarrachea-Magnani, B. L´opez-del-Carpio, J. Ch´avez-Carlos, S. Lerma-Hern´andez and J.G.Hirsch, Phys. Rev. E , 022215 (2016).[40] J. Ch´avez-Carlos, M.A. Bastarrachea-Magnani, S. Lerma-Hern´andez and J.G. Hirsch, Phys. Rev. E ,022209 (2016).[41] P. Cejnar and P. Str´ansk´y, AIP Conf. Proc. , 23 (2014).[42] P. Cejnar and J. Jolie, Phys. Rev. E , 387 (1998).[43] M.A. Bastarrachea-Magnani, S. Lerma-Hern´andez, J.G. Hirsch, J. Stat. Mech. , 093105 (2016).[44] A. Peres, Phys. Rev. Lett. , 1711 (1984).[45] P. Str´ansk´y, P. Hruˇska and P. Cejnar, Phys. Rev. E , 066201 (2009).[46] V. Vedral, Rev. Mod. Phys. , 197 (2002).[47] A. Peres, Phys. Rev. Lett. , 1413 (1996).[48] C.H. Bennett, D. DiVincenzo, J. Smolin and W.K. Wooters, Phys. Rev. A , 3824 (1996).[49] Y. Huang, New J. Phys. , 033027 (2014).[50] S. Hill and W.K. Wootters, Phys. Rev. Lett. , 5022 (1997).[51] W.K. Wootters, Phys. Rev. Lett. , 2245 (1998).[52] X. Wang and K. Mølmer, Eur. Phys. J. D
385 (2002),[53] X. Wang and B.C. Sanders, Phys. Rev. A , 012101 (2003).[54] S. Heinze, P. Cejnar, J. Jolie and M. Macek, Phys. Rev. C014306 (2006).