{'fraction_non_alphanumeric': 0.04135811094033813, 'fraction_numerical': 0.012668250197941409, 'mean_word_length': 4.145458902468248, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spontaneous symmetry breaking and Husimi Q-functions in extended Dicke model

We study parity breaking transition in the ground state of a quantum array of N two-level systems (represented by pseudospins 1/2) coupled to photon mode. The spin-photon coupling occurs via a shift of the photon oscillator equilibrium position by an amount linear in one projection of the total (pseudo)spin. We demonstrate analytically, that the system is unstable to emergence of a finite shift of photonic equilibrium that causes upturn from concave to convex dependence of the ground state energy on the total spin projection involved in light-matter coupling. Simultaneously, parity symmetry of the ground state wave function gets broken in strong coupling limit even for finite N. We show that previously suggested subradiant domain of the phase diagram of the system acquires a new feature: photonic shift is predicted inside the subradiant area along the strong coupling interval on the coupling strength axis. Husimi Q-functions of the photon and spin subsystems for different values of light-matter coupling are presented.


I. INTRODUCTION
The Dicke model in cavity QED describes interaction of a single photonic mode with an ensemble of the two level systems [1,2]. In thermodynamic limit there is a quantum phase transition described by parity symmetry breaking of the ground state wave function, though it is preserved by the Dicke Hamiltonian [3,4]. The extended Dicke model with a quadratic term (Ŝ i are pseudospin operators): (the Plank constant is put to one) arises for various light-matter interaction models [5][6][7][8], with photon momentump and coordinateq operators being defined below. When there is ε = 0 in the coefficient in front of the quadratic term the full square with respect top is completed and the spin-projection S y acts as a shift of the equilibrium position of the photonic oscillator.
The phase transition properties for ε = 0 were studied in [6,7]. It was shown, that for ε < 0 arXiv:2004.13772v1 [cond-mat.stat-mech] 28 Apr 2020 there is a transition to superradiant state, while a subradiant state realizes for ε > 0 and in narrow region of ε < 0 domain. In this work we focus on superradiant transitions on the ε = 0 line of the ε − g phase diagram.
The case of ε = 0 was investigated previously [8], where a first order quantum phase transition (QPT) was found in thermodynamic limit N 1 using rotating Holstein-Primakoff transformation. The parity symmetry of the ground state wave function gets broken as a result of formation of the superradiant photonic condensate. However, the existence of a phase transition in the Dicke model with quadratic term is subjected to debate [9][10][11][12]. Besides, it is obvious, that for finite total spin the parity symmetry of the Hamiltonian is always preserved [13,14]. This contradiction is studied in the present paper. Beyond a critical coupling g c the system with the Hamiltonian (1) gets spontaneous symmetry break along the ε = 0 g-axis.
Though the Hamiltonian symmetry is formally conserved for any finite values of parameters, the parity symmetric phase becomes unstable to a formation of a shift of the equilibrium position of momentum p of the cavity mode. We prove this by a straightforward analytical derivation, that demonstrates how an emergence of photonic shift upturns from concave to convex the dependence of the ground state energy on the total spin-projection S y involved in light-matter coupling in Hamiltonian (1).
The plan of the paper is as follows. First, we briefly dwell on the superradiant and subradiant quantum phase transitions discussed in the earlier works [3,6,7], then focus on the case of ε = 0.
The Born-Oppenheimer potentials will be used to illustrate the difference between different types of transitions. An analysis of spontaneous symmetry breaking is then made using idea [15,16] of non-commuting ultrastrong coupling limit and limit of non-zero symmetry breaking term (finite superradiant density). For this purpose a perturbation theory is used to obtain ground state energy as a function of the spin-projection S y involved in light-matter coupling and of a finite photonic equilibrium shift α, that breaks the symmetry of the ground state in the strong coupling limit. The analytic results demonstrate that emergence of photonic shift upturns the ground state energy dependence on the total spin projection S y from concave to convex. We also demonstrate numerically the occurrence of the phase transition in the strong coupling limit in the case ε = 0 using trial wave-function technique. In Appendix A the Husimi Q-functions of the photon and spin subsystems are presented that directly depict changes in behaviour of the subsystems under an increase of spin-photon coupling strength g at different values of parameter ε.

II. MODEL
We consider an extended Dicke Hamiltonian of an ensemble of N two-level systems interacting with bosonic mode: A single two-level system is described by spin 1/2 with projections operatorsσ x,y,z j . The ∝ g 2 ε term arises from dipole-dipole interaction [6,7]. Depending on ε one obtains the usual Dicke model (ε = −1) for strong ferroelectric inter-dipole interaction (i.e. inter-spin interaction in the spin-1/2 representation of the two-level systems), or extended Dicke model for repulsive interdipole interaction (ε > 0). The Hamiltonian of non-interacting dipoles with ε = 0 is obtained for Josephson junctions interacting merely with photonic mode [5,8] via gauge invariant shift of Cooper pairs phase differences across the junctions.
Introducing collective spin operatorsŜ x,y,z = jσ x,y,z j we rewrite (2) in the form: For qualitative investigation it is convenient to switch to momentum and coordinate operators constituting harmonic oscillator of the photon mode: The Hamiltonian then is expressed as: One can see that a full square can be completed: When ε = 0 and E J = 0 theŜ y can be treated as a gauge field and the Hamiltonian becomes gauge invariant because gŜ y only shifts equilibrium position of the oscillator, but does not effect its energy.

III. SUPERRADIANT AND SUBRADIANT PHASES
In previous works [3,6] it was shown, that depending on the sign of ε in the Hamiltonian (7), different types of quantum phase transitions exist in the system. The transition occurs at some critical coupling g c and manifests itself in a change of the system macroscopic observables, such as average photons number and spin projection.
To illustrate the different phases we calculate Born-Oppenheimer potential of the Hamiltonian in strong coupling limit (i.e. E J /g 2 → 0). It is a set of parabolas: with gaps opened at their crossings due to −E JŜz term, σ y is an eigenvalue ofŜ y . The multiwell potential possesses 2S +1 minima at p = −gσ y [6,17]. Depending on the sign of ε, different potential wells will be the lowest ones: at p = ±gS for negative ε , and at p = ±g/2 or p = 0 for positive ε. In the latter case the minimum is at p = 0 for integer spin and at p = ±g/2 for half-integer.

A. Superradiant phase
The superradiant quantum phase transition occurs when ε < 0, in this case the critical coupling is [6] g c ∼ E J −εN , average photonic oscillator momentum and spin projection on yaxis after the transition are p ∼ g Ŝ y ∼ gN . The Born-Oppenheimer potential has lowest wells at its ends, so the system finds itself in the corresponding state with non-zero averages. This is illustrated in fig. 1.

B. Subradiant phase
In subradiant phase (i.e. for ε > 0) the projection on y-axis is minimal, being ±1/2 for odd N and 0 for even one. Correspondingly the average photon-oscillator momentum is p ∼ g/2 or 0. In Born-Oppenheimer potentials one can see, that the lowest wells are the central ones.
If N is odd there is no central well with σ y = 0 and the averages are minimal, but non-zero. In opposite case, there is the σ y = 0 well where the averages are zero. The potentials are depicted in fig. 2. The subradiant phase transition occurs as a cross-over [6] in the limit of infinite coupling strength g, which is not the thermodynamic limit, because the ground state wave function does not depend on system's size.

IV. PARITY CONSERVATION AND SPONTANEOUS SYMMETRY BREAKING
The Hamiltonian (3) commutes with the parity operator so there is a corresponding conservation law. Neighbouring energy levels in the spectrum have different parity. The operatorN is the excitations number operator. Each term in the Hamiltonian changes the number of excitations by 2 because it is either quadratic inp,q andŜ y or is their productpŜ y , so the parity stays unchanged.
The Hamiltonian commutes withΠ for any values of parameters and formally the symmetry will be always preserved. However, there is still a possibility for a phase transition in the strong coupling limit. First of all, as coupling increases, the energy gap between states with different parity tends to zero, meaning they will mix and the conservation law could be broken. A clear sign of spontaneous symmetry breaking is non-commuting limits of infinite coupling and infinitesimal symmetry breaking terms. We consider this in detail in the next section.
The Hamiltonian with ε = 0 arises in the description of an array of Josephson junctions coupled to microwave cavity via gauge-invariant shift of the junctions phases [5,8]. The Hamiltonian expressed via photon oscillator coordinate and momentum is: One can clearly see, thatŜ y acts as a gauge field. For large couplings g, when theŜ z term is negligible, the Hamiltonian is invariant with respect to the change of projection σ y . This means that in Born-Oppenheimer potential all wells share the same level, see fig. 3. However, any small E J will destroy the degeneracy making preferable the spin state with maximal total spin S and with the minimal projection on y-axis [6]. At first glance, this will lead to subradiant phase transition as if ε > 0, but it is not the case. The transition has its own unique properties which we discuss below (see also [8]).

A. Spontaneous symmetry breaking
The system's wave function ψ depends on coupling constant g in the Hamiltonian (3) the system will undergo a spontaneous symmetry breaking in g → ∞ limit [15].
To study this phenomenon of not commuting limits analytically in the strong coupling limit of g we add small symmetry breaking term αp to the Hamiltonian and treat −E JŜz term as a The energy spectrum (counted from zero-point energy ω/2) and wave functions of unperturbed HamiltonianĤ 0,α are: ψ α n,σy = n gσy+α , σ y , where index n gσy+α indicates photon number n state of the photon mode in the Hamiltonian (13) with shifted equilibrium position −(gσ y + α) of the harmonic oscillator, and index σ y defines total spin projection. Now, to second order in E J , perturbation theory gives the following correction to the energy of the ground state E α,0 (compare [7]): where: γ 2 = g 2 /ω and δ = gα/ω. In the r.h.s limit in (12), γ 1, and δ < γ 2 , one finds from (16) and (14): First, minimize E 0 with respect to σ y : Now, substitute the resulting expression for α m back into (17) and find the minimal energy as function of σ y : On the other hand, if the l.h.s limit in (12) is implemented first, i.e. α → 0, one finds [7]: Remarkably, the sign of the curvature of the ground state energy E 0 as function of pseudo-spin projection σ y is different in (19) and (20). Hence, indeed, the limits in (12) do not commute.
Namely, the limit of conserved parity α = 0 produces concave dependence of E 0 on σ y , (20), as is the case in subradiant configuration of the Born-Oppenheimer effective potential in Fig. 2.
On the other hand, in the limit γ 1, with a finite parity breaking term αp in the Hamiltonian (13), the ground state energy E 0 is convex function of pseudo-spin projection σ y , (19), with α m = E 2 J S/g 3 → 0 in the limit g → ∞, as is the case in the superradiant configuration of the Born-Oppenheimer effective potential in Fig.1. Thus, a transition into parity broken "dipole ordered" state with σ y = ±S = ±N/2 takes place in this case.The latter state is stabilized by the condition g ≥ g c ,which ensures the energy of the state is the lowest possible one, where [8]: A few observations are in order here. First, we emphasise, that both dependencies E 0,rhs and E 0,lhs in (19) and (20) The shift of photonic equilibrium p in (21) is in accord with the result obtained in [8] for the amplitude of photonic superradiant condensate in the limit N 1. On the other hand, inequalities: g 2 ω and E 2 J S/g 4 << 1 leave space for occurrence of the symmetry breaking transition already in the case N ∼ 1, though, it would be not a thermodynamic phase transition in a strict sense, compare [13,14].

B. Variational approach at ε = 0
Case of ε = 0 is hard to investigate numerically, because even small numerical error will bring the system into ε > 0 or ε < 0 case. Therefore, here we use variational method to study spin configuration in the ground state. The test function is Coefficient b can be expressed via a by using normalization condition: and we obtain test function ψ(a, b(a)) = ψ(a), that depends on a single parameter a. The result now does not depend on spin being integer or half-integer. In Fig. 4 one can see that the minimum of E(a) changes its position from a 0 > 0 to a 0 = 0 when coupling constant g increases. It means that at g > g c the system behaves as it does at ε < 0 case having maximal absolute value of spin projection on the y -axis. We calculate average of momentum operator ψ(a 0 )|p|ψ(a 0 ) = −g ψ(a 0 )|S y |ψ(a 0 ) in the ground state as function of coupling constant g.
using the variational wave function ψ(a 0 ) found above. Then, there is a jump from p = 0 to p = 0 at the transition point when S y becomes non-zero. This jump signifies emergence of superradiant condensate in the system (compare [8]). We conclude that at ε = 0 the system behaves like in ε < 0 case, i.e. there is a transition to maximal absolute value projection of spin on y-axis. This is in agreement with analysis presented in [8] and with the results presented above. We have summarised our findings in the diagram in Fig. 6.

C. Ground state wave function
We now focus on the ground state wave function. The energy (14) is minimal for n = 0 and σ y = S if g, α > 0. Then, the ground state wave function is: From Born-Oppenheimer potential point of view this may be thought of as if the symmetry breaking term has added a slope to the initially "horizontal" bottom line in Fig. 3, thus, driving the system to the well at the edge, corresponding to σ y = S. Since the lowest potential well is now positioned at p = gS, in the limit of g → ∞ the system occupies state far from the central potential well (p = 0) in Fig.3. Simultaneously the parity operator (9) does not commute with Hamiltonian (13).
The ground state energy in this case does not depend on σ y so the ground state is highly degenerate, but symmetric. It means, that the state of the system is a superposition of states symmetric with respect to p = 0 wells.
Let us start from the state ψ αβ 0,S . If we first take the limit of zero symmetry breaking, i.e. α → 0, the wave function will restore the superposition with σ y = −S state via tunnelling to symmetric potential well and acquiring the form (25). It will now also remain symmetric in the limit of g → ∞.
Alternatively, if we take the limits in the different order, i.e the limit g → ∞ will be taken first, then the system will not be able to restore the superposition state (25) after α → 0 limit, because the symmetric counterpart of the system's state will be localized infinitely far at p ∼ −gS. One can see, that depending on the order of the limits the system reaches the different states: either symmetric or not. This means that in the infinite coupling limit (or in the zero photonic mode limit ω → 0 [13]) the symmetric state is unstable and a spontaneous symmetry breaking will occur (compare [14]).

VI. CONCLUSION
The phase transition at ε = 0, see Fig. 6, that was previously considered using Goldstein-Primakoff large N approximation in [8], has properties of both superradiant and subradiant transitions.
After the transition the system finds itself in the state with non-zero macroscopic observables, i.e p ∼ g Ŝ y ∼ gN like in the superradiant phase. But in contrast to superradiant transition and similar to subradiant one, the ultrastrong critical coupling condition The transition is caused not by a direct change of the Hamiltonian symmetry, but due to spontaneous symmetry breaking in the limit of ultrastrong coupling g, that leads to exponentially small energy gap separating states of different parity. In paper [6] the phase diagram in It is instructive to look into the qualitatively different behaviour of the model in the cases ε < 0 and ε > 0 using the method of Husimi Q-functions [19]. The parity breaking causes the change of the states of spin and photon subsystems. In order to study them we numerically calculate corresponding Husimi Q-functions. Given an eigenstate |Ψ of the whole system one can obtain reduced density matrix of a constituting subsystem by a partial trace operation: where Tr S,Ph denotes partial trace over spin or photon subsystem respectively. Husimi Qfunction is then the reduced density matrix averaged over a coherent state, i.e. if the coherent state is parametrised by a set of variables α 1 , . . . , α n Q(α 1 , . . . , α n ) = α 1 , . . . , α n | ρ |α 1 , . . . , α n .
Its value gives the probability of the subsystem being in the state |α 1 , . . . , α n .

Spin subsystem
Spin coherent state is parametrised by spherical angles φ and θ and can be obtained by an action of rotation operator [20]R (θ, φ) = e iθ(Ŝz sin φ−Ŝy cos φ) (A4) on the eigenstate ofŜ z , i.e. |θ, φ =R(θ, φ) |σ z = S . Husimi Q-function is given then by We will analyse spin configuration that arises from Husimi Q-function, which is calculated using density matrix of the ground state. The latter is distinguished by the minimal eigenvalue of the Hamiltonian matrix. At zero coupling the energy is minimal when spin-projection on z-axis has maximal value. In the opposite case of strong coupling, the spin orientation differs depending on ε. According to (8), projection on y-axis may increase or decrease the system's energy depending on the sign of ε. The energy is minimal at maximal σ y for ε < 0, so that spin aligns along y-axis. For ε > 0 projection modulus on y-axis tends to be minimal. Husimi Q-functions posses different symmetries in the two limiting cases. In the limit of zero coupling there is only one preferable direction along z-axis while the system is invariant with respect to rotation in x-y plane. So, the Q-function is oriented in the positive direction of z-axis and is rotationally symmetric around it. At large coupling the system's state is characterised by σ y and the energy does not depend on x and z projections. Due to energy being quadratic in σ y there is no preferred y-direction and the functions have a corresponding symmetry. As a result, g < g c S = 3 2 , ε < 0, g > g c S = 3 2 , ε > 0, g > g c S = 2, ε < 0, g > g c S = 2, ε > 0, g > g c Figure 7: Husimi Q-function for spin subsystem in the ground state the Q-function is symmetric in x-z plane and orients itself along y-axis in both directions or conversely minimizes its value along y-axis depending on the sign of ε as discussed previously.
One can see in Fig. 7, that the symmetry of the Q-function of the system in the ground state is very different in two limiting cases of weak and strong coupling g. Hence, a symmetry breaking transition must occur during the gradual increase of the coupling constant g.

Photon subsystem
One can write Q-function for photons as well: Q P h (p, q) = p, q| ρ P h |p, q , where ρ P h is the photon density matrix, |p, q is a coherent photon state with p = p and q = q. This function demonstrates probability density of photon having coordinates p and q in the phase space.
At coupling below critical the distribution is peaked at p = 0, q = 0 meaning there is no photon condensate. As coupling grows above critical the situation differs depending on spin configuration -the position of the distribution peak is proportional to projection of spin on y-axis. So, for integer spin and ε > 0 when in the ground state σ y = 0 there is still no photon condensate and a single peak remains. If the spin is not integer and σ y = 1/2 peaks appear at p ∼ g/2. For ε < 0 the spin projection in ground state is maximal and p ∼ gS. Due to a b Figure 8: Husimi Q-function for photon subsystem in q ω/2, p/ √ 2ω normalized coordinates. a) State without photon condensate at g < g c . It remains also at g > g c in the ground state of integer-spin Hamiltonian for ε > 0. b) Two-peeked distribution state at g > g c with photon condensate, in case of ε < 0, or ε > 0 and half-integer spin.
symmetry of spin with respect to changing sign of σ y there are two peeks centred along p-axis and placed on it symmetrically. The two latter cases look qualitatively similar and at first glance differ only quantitatively i.e. by position of Q-function maximum and value of critical coupling. But one should remember that these two transitions differ by spin configuration and also the ε < 0 case is not sensitive to spin being integer or half-integer. One-peeked and two-peeked Husimi Q-functions are shown in Fig. 8.