Critical Phenomena in Light–Matter Systems with Collective Matter Interactions

Abstract

We study the quantum phase diagram and the onset of quantum critical phenomena in a generalized Dicke model that includes collective qubit–qubit interactions. By employing semiclassical techniques, we analyze the corresponding classical energy surfaces, fixed points, and the smooth Density of States as a function of the Hamiltonian parameters to determine quantum phase transitions in either the ground (QPT) or excited states (ESQPT). We unveil a rich phase diagram, the presence of new phases, and new transitions that result from varying the strength of the qubits interactions in independent canonical directions. We also find a correspondence between the phases emerging due to qubit interactions and those in their absence but with varying the strength of the non-resonant terms in the light–matter coupling. We expect our work to pave the way and stimulate the exploration of quantum criticality in systems combining matter–matter and light–matter interactions.

1. Introduction

Quantum phase transitions (QPT) are generally defined as the sudden change in the properties of the ground state of a quantum system as a function of a control parameter. They possess an essential role in modern physics, especially in studying many-body quantum systems, quantum information, and quantum control [1,2]. The active interest in quantum critical phenomena during the last two decades stems from their impact on the spectral features and dynamics of complex quantum systems, leading, e.g., to the development of new concepts such as that of Excited-State Quantum Phase Transition (ESQPT), which is meant to explain the consequences of the propagation of critical behavior from the ground state to the rest of the spectrum of a quantum system [3,4,5,6]; and that of Dynamical Quantum Phase Transition (DQPT), seeking to fathom the onset of criticality exhibited in non-equilibrium phenomena [7,8,9]. Typically, understanding Quantum Phase Transitions depends on the specific system of study. Thus, the field remains an open challenge with exciting avenues, striving to reach a general framework to describe the interplay between many-body properties, strong interactions, and critical phenomena.

A paradigmatic example of a QPT is the Dicke Hamiltonian’s superradiant phase transition [10,11]. The Dicke model describes a collection of atoms within the two-level approximation interacting with a single-mode radiation field inside a cavity [12]. The superradiant QPT is characterized by a non-zero expectation value of the photon number when the light–matter strength reaches a critical value in the thermodynamic limit. As it describes the collective degrees of freedom of a set of two-level systems (qubits), the Dicke Hamiltonian offers a general description of the spin–boson interaction. Additionally, it constitutes a paradigmatic example for the study of the ultra-strong coupling (USC) regime [13,14,15] Consequently, the model has found a great reception in the description of several setups, mainly in the context of quantum information [16,17,18,19]. In recent years, it has been experimentally realized in a broad range of tunable systems, from Bose–Einstein condensates in optical lattices [7,20,21,22,23], superconducting qubits [24,25,26] to cavity-assisted Raman transitions [27,28]. Not only formal derivations of Dicke-like Hamiltonians have been found in the framework of ultracold atoms in optical lattices [29,30,31], and for superconducting qubits [32,33,34,35], but superradiant effects have also been proposed in nuclear physics [36], solid-state physics [37], bidimensional materials [38,39], and quantum dots [40], among others. Moreover, its algebraic simplicity allows one to employ it as a test bed for studying critical features in the spectrum, several topics relevant to quantum information such as quantum chaos and quantum correlations [41,42,43], and the quantum-classical correspondence [44,45,46,47]. The latter is possible because the Hamiltonian can be mapped to only two relevant degrees of freedom (those of the boson and the collective spin), so there is a well-defined classical limit. This feature has raised questions about the nature of the superradiant QPT, being deemed as a mean-field QPT due to its classicality and smallness of quantum fluctuations, which are the ones driving phase transitions at low temperatures [48].

An additional feature in the systems where the Dicke model finds application is the possibility to build up collective qubit–qubit interactions. Previous works have addressed this problem by considering, e.g., dipolar interactions in atomic systems [49], shifts due to the Stark effect in optomechanical setups [50,51], Josephson dynamics in a two-component BEC [52,53,54], or the onset of chaos [55,56]. Two general results stand due to the presence of matter interactions: The prediction of a first-order phase transition [53,55,57,58,59], the shift of the critical coupling of the standard superradiant phase transition [32,49]—including its possible suppression—and a richer phase diagram [55,60]. Despite previous studies of the ground-state properties in this system, the interplay of these collective interactions on the spectral properties of the Hamiltonian, as well as the understanding of the different energy domains marked by the presence of quantum phase transitions, has not been investigated exhaustively.

In this work, we are interested in studying the critical behavior of a generalized Dicke Hamiltonian that includes collective qubit–qubit interactions. It will constitute an example of the intriguing combination between matter–matter interactions and (ultra) strong light–matter ones and the rich phase diagrams they can produce. Unlike previous works, here, we add a general combination of non-linear interactions in the form of quadratic terms in the collective pseudo-spin operators to the standard Dicke Hamiltonian in all the x-, y-, and z-directions. Then, we perform a standard semiclassical analysis to obtain the behavior of energy surfaces, the ground-state energy, and the Density of States (DoS) as a function of the Hamiltonian parameters. This exploration allows us, from a unified perspective, to obtain indicators of critical quantum behavior, i.e., both QPT and ESPQT, as has been done previously in other works [55,61,62]. We offer a general overview that unifies some results of previous works and found new behavior unlocked by the unique interplay between the different directions of the interactions.

The article is organized as follows. In Section 2, we present the generalized Dicke Hamiltonian, including qubit–qubit interactions. Next, in Section 3, we discuss the corresponding classical Hamiltonian obtained via coherent states, the Hamilton equations of motion, and the fixed points, commonly associated with critical behavior in the ground state of the related quantum system. Furthermore, we present an overview of the classical energy surfaces and classify the different phases according to the ground-state properties. In Section 5, we calculate the semi-classical density of states (DoS) to identify the ESQPT and the spectral domains. Finally, in Section 6, we present our conclusions. We include several appendices with details on the calculations.

2. Generalized Dicke Hamiltonian

We study a generalized Dicke Hamiltonian that includes collective qubit–qubit interactions

$$\begin{array}{c}{\widehat{H}}_D={\widehat{H}}_0+{\widehat{H}}_{\mathrm{I}}+{\widehat{H}}_{\mathrm{qq}},\end{array}$$

(1)

where

$$\begin{array}{c}{\widehat{H}}_0=\omega {\widehat{a}}^{†}\widehat{a}+{\omega }_0{\widehat{J}}_z,\\ {\widehat{H}}_{\mathrm{I}}=\frac{\gamma }{\sqrt{N}}\left[(\widehat{a}{\widehat{J}}_{+}+{\widehat{a}}^{†}{\widehat{J}}_{-})+\xi ({\widehat{a}}^{†}{\widehat{J}}_{+}+\widehat{a}{\widehat{J}}_{-})\right],\\ {\widehat{H}}_{\mathrm{qq}}=\frac{1}{N}\left({\eta }_x{\widehat{J}}_x^2+{\eta }_y{\widehat{J}}_y^2+{\eta }_z{\widehat{J}}_z^2\right).\end{array}$$

The first term denotes the non-interacting Hamiltonian ${\widehat{H}}_0$, the second one is the ${\widehat{H}}_{\mathrm{I}}$ usual spin–boson interaction, and the last one contains the qubit–qubit interactions ${\widehat{H}}_{\mathrm{qq}}$. Here, ${\widehat{a}}^{†}$ ($\widehat{a}$) is the creation (annihilation) boson operator, and ${\widehat{J}}_{z,x,y}$ are the pseudo-spin operators representing the collective degrees of freedom of the set of N qubits, which follow the rules of the su(2)-algebra. Generally, a set of N qubits is spanned into a $2^N$ dimensional Hilbert space; however, as we are interested in describing the collective degrees of freedom, it suffices to work in the totally symmetric subspace corresponding to $j=N/2$, where $j(j+1)$ is the eigenvalue of the pseudo-spin length operator ${\widehat{\mathbf{J}}}^2$. Thus, the dimension of the Hilbert space is reduced to only $N+1$, where the collective ground state lies. The Hamiltonian parameter set is given by $\omega$, ${\omega }_0$, and $\gamma$, which are the boson frequency, the qubit energy splitting, and the spin–boson interaction. Additionally, we have ${\eta }_i$ with $i=x,y,z$, the collective qubit–qubit couplings in each direction. Depending on the setup, one can grant a specific meaning to the interactions in the z and its perpendicular directions. An intuitive approach comes from interacting Bose–Einstein condensates in a two-site trap and Josephson effects. There, ${\widehat{J}}_z$ is related to a relative population of particles in the condensates and ${\widehat{J}}_x$ (${\widehat{J}}_y$) to ladder operators and relative phases between them. Thus, ${\eta }_z$ and ${\eta }_x$ (${\eta }_y$) represent the strength of collective on-site and between neighboring sites interactions (hopping effects), respectively [52]. Otherwise, interactions from the x and y directions arise from dipolar coupling in atomic setups [49,63] or interactions between superconducting qubits [32,34,35].

The Hamiltonian in Equation (1) possesses several well-known limits. In the absence of qubit–qubit interactions (${\eta }_i=0$ for $i=x,y,z$), one recovers the standard light–matter interaction. The parameter $\xi$ takes the system from the integrable Tavis–Cummings model ($\xi =0$) [64], which describes a system in the strong coupling regime under the Rotating-Wave Approximation (RWA), to the standard, non-integrable Dicke model ($\xi =1$) typically describing the USC [12]. In both limits, the superradiant QPT takes place when the light–matter coupling crosses the critical value ${\gamma }_{\xi +}=\sqrt{\omega {\omega }_0}/(1+\xi )$. For values below the coupling ($\gamma <{\gamma }_{\xi +}$), the system is in a normal phase, characterized by a zero-average of photon population in the thermodynamical limit ($\overline{n}=\langle {\widehat{a}}^{†}\widehat{a}\rangle /N=0$), while for ${\gamma }_{\xi +}>\gamma$ one finds a finite photon number $\overline{n}\ne 0$, thus called superradiant phase. Additionally, the Hamiltonian exhibits two ESQPTs [61,65,66,67], which are identified as non-analyticities in the derivative of the smooth DoS as a function of energy in the thermodynamic limit. One at a critical energy $E_{-}^{\left(c1\right)}/{\omega }_{0}j={ϵ}_{-}^{\left(c1\right)}=-1$ that only appears in the superradiant phase (characterized by a logarithmic divergence in the derivative of the DoS as energy increases), and a second one at ${ϵ}_{+}^{\left(c2\right)}=+1$ that appears for every coupling as a jump singularity (a step function in the DoS derivative) [68] and is related to the saturation of the collective qubit Hilbert space.

A finite value of $\xi \in (0,1)$ leads to the generalized or extended Dicke model instead [6,45,61,69]. There, a new superradiant phase appears whose critical point occurs at ${\gamma }_{\xi -}=\sqrt{\omega {\omega }_0}/(1-\xi )$ [62,69,70]. While in the TC model ${\gamma }_{0+}={\gamma }_{0-}$ and the new phase is equal to the standard one; in the Dicke model, ${\gamma }_{1-}\to \infty$, so it becomes not observable. The ESPQTs predicted in the Dicke model persist in the generalized one; the only difference is that the ESQPT changes its type from a logarithmic singularity in the derivative of the smooth DoS to a step function with a downward jump from lower to higher energies in the interval ${\gamma }_{\xi +}<\gamma <{\gamma }_{\xi -}$ [4,6,70]. On the other hand, in the absence of light–matter interaction, the boson is decoupled from the collective spin. Then, one arrives at a version of Lipkin–Meshkov–Glick Hamiltonian (LMG) [71,72,73], a well-known model with one degree of freedom, originally stemming from nuclear physics, that nowadays is connected to Josephson junctions and cold atoms in optical lattices [74]. Critical phenomena in the LMG model have been extensively studied, and the system exhibits both a first-order and a second-order QPT, as well as ESQPTs [75,76,77,78,79,80,81,82,83], to cite some works. Naturally, we expect that the generalized Dicke Hamiltonian inherits critical features from the LMG model by including the qubit–qubit interactions.

3. Classical Corresponding Hamiltonian

A classical Hamiltonian can be obtained by taking the expectation value of Equation (1) over a tensor product of Glauber $|z\rangle$ and Bloch $|w\rangle$ coherent states as trial states [45,47,61,67,84], where $|0\rangle$ and $|j,-j\rangle$ are the boson and pseudo-spin vacuum states, respectively [85],

$$\begin{array}{c}|z\rangle \otimes |w\rangle =\frac{e^{-{\left|z\right|}^2/2}}{{(1+|w|}^2{)}^j}{e}^{z{\widehat{a}}^{†}}{e}^{w{\widehat{J}}_{+}}|0\rangle \otimes |j,-j\rangle .\end{array}$$

(2)

By dividing over j we obtain

$$\begin{array}{c}{H}_{\mathrm{cl}}^{\left(\xi \right)}(z,w)=j^{-1}\langle z,w|{\widehat{H}}_D|z,w\rangle ={\omega \left|z\right|}^2-\left(\frac{1-{\left|w\right|}^2}{1+{\left|w\right|}^2}\right)\left[{\omega }_0-\frac{{\eta }_z}{2}\left(\frac{1-{\left|w\right|}^2}{1+{\left|w\right|}^2}\right)\right]+\\ \frac{1}{2{\left(1+{\left|w\right|}^2\right)}^2}\left[\left({\eta }_x-{\eta }_y\right)\left(w^2+{\overline{w}}^2\right)+\left({\eta }_x+{\eta }_y\right)2w\overline{w}\right]+\frac{\gamma \left(z+\overline{z}\right)\left(w+\overline{w}\right)}{\sqrt{2}{(1+|w|}^2)}.\end{array}$$

(3)

Instead of employing the complex numbers z and w, it is more convenient to use canonical classical variables $(q,p)$ and $(j_z,\varphi )$ for the boson and spin spaces, respectively. Here, $z=\sqrt{j/2}\left(q+ip\right)$ and $w=\sqrt{(1+j_z)/(1-j_z)}{e}^{-i\varphi }$. Additionally, due to the fixed value of the pseudospin lenght $j=N/2$, we have $j_x=\sqrt{1-j_z^2}cos\varphi$ and $j_y=\sqrt{1-j_z^2}sin\varphi$. In this manner, we obtain a classical generalized Dicke Hamiltonian that reads

$$\begin{array}{c}{H}_{cl}^{\left(\xi \right)}=\frac{\omega }{2}(q^2+p^2)+j_z\left({\omega }_0+\frac{{\eta }_z{j}_z}{2}\right)+\frac{1}{2}\left(1-j_z^2\right)\left({\eta }_x{cos}^2\varphi +{\eta }_y{sin}^2\varphi \right)+\\ +\gamma \sqrt{1-j_z^2}\left[(1+\xi )qcos\varphi -(1-\xi )psin\varphi \right].\end{array}$$

(4)

To characterize the energy surfaces and identify the critical behavior, we will need the equations of movement. The Hamilton equations are

$$\begin{array}{c}\dot{q}=\frac{\partial H_{cl}^{\left(\xi \right)}}{\partial p}=\omega p-\gamma \sqrt{1-j_z^2}(1-\xi )sin\varphi .\end{array}$$

(5)

$$\begin{array}{c}\dot{p}=-\frac{\partial H_{cl}^{\left(\xi \right)}}{\partial q}=-\omega q-\gamma \sqrt{1-j_z^2}(1+\xi )cos\varphi ,\end{array}$$

(6)

$$\begin{array}{c}\dot{\varphi }=\frac{\partial H_{cl}^{\left(\xi \right)}}{\partial j_z}={\omega }_0+{\eta }_z{j}_z-j_z({\eta }_x{cos}^2\varphi +{\eta }_y{sin}^2\varphi )-\frac{\gamma j_z}{\sqrt{1-j_z^2}}\left[(1+\xi )qcos\varphi -(1-\xi )psin\varphi \right],\end{array}$$

(7)

$$\begin{array}{c}\dot{j_z}=-\frac{\partial H_{cl}^{\left(\xi \right)}}{\partial \varphi }=\left(1-j_z^2\right)({\eta }_x-{\eta }_y)cos\varphi sin\varphi +\gamma \sqrt{1-j_z^2}\left[(1+\xi )qsin\varphi +(1-\xi )pcos\varphi \right].\end{array}$$

(8)

In Appendix A we present the Hamilton equations for $\xi =0$ and $\xi =1$.

4. Energy Surfaces and Their Extrema

In this section, we obtain the fixed, stationary, or equilibrium points $(q_s,p_s,j_{zs},{\varphi }_s)$ of the energy surface $H_{cl}(q_s,p_s,j_{zs},{\varphi }_s)$ from Hamilton equations. They ease characterizing the system’s different quantum phases and transitions as a function of the Hamiltonian parameters, employing that the minimum of the energy surface can be identified with the ground-state energy [86]. In this case, the thermodynamic limit coincides with the classical limit because we have an effective Planck’s constant given by ${\hslash }_{eff}=\hslash /N$. Therefoe, as $N\to \infty$, ${\hslash }_{eff}\to 0$ [87]. Thus, to find the fixed points, we make Hamilton Equations (5)–(8) equal to zero. From the first two, we obtain a pair of equations defining the quadratures

$$\begin{array}{c}{p}_s=\frac{\gamma }{\omega }\sqrt{1-j_{zs}^2}(1-\xi )sin{\varphi }_s,and{q}_s=-\frac{\gamma }{\omega }\sqrt{1-j_{zs}^2}(1+\xi )cos{\varphi }_s,\end{array}$$

(9)

Next, we can insert them into Equations (7) and (8) to obtain a second pair of equations that set the atomic (collective spin) variables

$$\begin{array}{c}{\omega }_0+j_{zs}\left\{{\eta }_z-\left({\eta }_x{cos}^2{\varphi }_s+{\eta }_y{sin}^2{\varphi }_s\right)+\frac{{\gamma }^2}{\omega }\left[{(1+\xi )}^2{cos}^2{\varphi }_s+{(1-\xi )}^2{sin}^2{\varphi }_s\right]\right\}=0,\end{array}$$

(10)

$$\begin{array}{c}(1-j_{zs}^2)\left\{({\eta }_x-{\eta }_y)-\frac{{\gamma }^2}{\omega }\left[{(1+\xi )}^2-{(1-\xi )}^2\right]\right\}cos{\varphi }_{s}sin{\varphi }_s=0.\end{array}$$

(11)

We observe that Equations (10) and (11) are enough to determine the general conditions to find the fixed points. To better visualize the energy surfaces we study throughout this work, we use Equations (9) and a new set of atomic variables, as described in Appendix B. Then, the energy surface is restricted to the atomic space, simplifying the identification of fixed points.

4.1. Deformation of the Normal Phase

Two fixed points exist for every value of the Hamiltonian parameters. They come from Equation (11), when one makes $j_{zs}=\pm 1$. From Equation (9), at these values, one automatically gets that $p_s=q_s=0$, where ${\varphi }_s$ is left indeterminate given that they coincide with the poles of the unitary Bloch sphere. The coordinates of these stationary points are

$$\begin{array}{c}(p_s,q_s,j_{zs},{\varphi }_s)=(0,0,\pm 1,\mathrm{indeterminate})\end{array}$$

(12)

and their energy is given by

$$\begin{array}{c}{ϵ}_{\pm }=\pm 1+\frac{{\eta }_z}{2{\omega }_0}\end{array}$$

(13)

In the normal phase, the stationary point at $j_{zs}=-1$ is a stable, absolute minimum [61]. It corresponds to the lowest energy value of the system, marking the quantum ground-state energy. As the expectation value of the photon number at the ground state $\langle \mathrm{g}.\mathrm{s}.\left|{\widehat{a}}^{†}\widehat{a}\right|\mathrm{g}.\mathrm{s}.\rangle /N$ is in general proportional to ${\left|z\right|}^2=q^2+p^2$, it characterizes the quantum features of each phase. At the fixed point $j_z=-1$, we have that $q_s=p_s=0$, so we speak of a normal phase, distinguished by the absence of a strong-correlated light–matter quantum state that could lead, e.g., to a collective emission of photons. On the other hand, the point at $j_{zs}=+1$, which always belongs to a higher energy domain, is typically an unstable fixed point [88]. As this point signals the maximum energy of the pseudospin (given that $|j_{zs}|\le 1$), the entire phase space associated with the Bloch sphere becomes available for the pseudospin dynamics. Any additional energy will only increase the boson energy. Thus, it marks the onset of an ESQPT [61,66].

The existence of a single global minimum in the standard normal phase of both the TC and the Dicke models is followed by energy surfaces that are invariant under $\varphi$ rotations. This is connected to the conservation of the total number of excitations operator $\widehat{\Lambda }={\widehat{a}}^{†}\widehat{a}+{\widehat{J}}_z+j\widehat{\mathbb{I}}$, as in the normal phase, the system is virtually decoupled. In this case, the shape of the potential corresponds to a single well, as shown in Figure 1(c3). However, even though the nature of the fixed points does not change, when ${\eta }_x-{\eta }_y\ne 0$, the energy surface becomes deformed thanks to the influence of interactions in x and y directions, and the rotational symmetry is broken at higher energies. We call this situation the deformed normal phase. The energy surfaces corresponding to this situation are shown, for example, in Figure 1(c4),(c7), where either the interactions in the y or x directions are present. Later, once we identify the parameter domains for the other phases, we will explain how the surface stretches depending on the qubit interaction directions.

Finally, we notice that the energy of the points at $j_z=\pm 1$ is invariant concerning the qubit interactions in the x and y directions. Still, it is uniformly shifted by the z-interactions, as was noted before in previous works [53,55,89]. The normal phase is thus identified by the presence of only these two fixed points. Additional stationary points will emerge, and the point at $j_z=-1$ will change its type of extrema according to the onset of the other phases. Next, we solve Equations (9)–(11) to find those points and the phases for three different situations: (1) The Tavis–Cummings limit ($\xi =0$), (2) the Dicke limit ($\xi =1$), and (3) for an arbitrary value of $\xi$.

4.2. Tavis–Cummings Limit

By setting $\xi =0$, we cancel the non-resonant terms in the Hamiltonian (${\widehat{a}}^{†}{\widehat{J}}_{+}$ and $\widehat{a}{\widehat{J}}_{-}$). Hence, we recover a Tavis–Cummings model modified by the qubit–qubit interactions. It corresponds to the situation where the Rotating-Wave Approximation (RWA) holds [90]. The Hamiltonian becomes

$$\begin{array}{c}{H}_{cl}^{\left(0\right)}=\frac{\omega }{2}(q^2+p^2)+j_z\left({\omega }_0+\frac{{\eta }_z{j}_z}{2}\right)+\frac{1}{2}\left(1-j_z^2\right)\left({\eta }_x{cos}^2\varphi +{\eta }_y{sin}^2\varphi \right)\\ +\gamma \sqrt{1-j_z^2}\left(qcos\varphi -psin\varphi \right),\end{array}$$

(14)

and Equations (10) and (11) are in this case

$$\begin{array}{c}{\omega }_0+j_{zs}\left({\eta }_z-\left({\eta }_x{cos}^2{\varphi }_s+{\eta }_y{sin}^2{\varphi }_s\right)+\frac{{\gamma }^2}{\omega }\right)=0,\end{array}$$

(15)

$$\begin{array}{c}(1-j_{zs}^2)({\eta }_x-{\eta }_y)cos{\varphi }_{s}sin{\varphi }_s=0.\end{array}$$

(16)

We find five different solutions for the stationary points (including $j_{zs}=\pm 1$). The other three conditions are given by the cases: $cos{\varphi }_s=0$ ($sin{\varphi }_s=\pm 1$), $sin{\varphi }_s=0$ ($cos{\varphi }_s=\pm 1$), and ${\eta }_x={\eta }_y$.

4.2.1. Superradiant-Symmetric Phase

We start studying the situation when the interactions in the x and y directions are the same ${\eta }_x={\eta }_y=\eta$. Equation (15) becomes

$$\begin{array}{c}{\omega }_0+j_{zs}\left({\eta }_z-\eta \right)+\frac{{\gamma }^2}{\omega }=0.\end{array}$$

(17)

Thus, we can obtain $j_{zs}$ immediately as

$$\begin{array}{c}{j}_{zs}=-\frac{{\omega }_0}{{\eta }_z-\eta +\frac{{\gamma }^2}{\omega }}=-\frac{1}{f_0},f_0=\frac{\Delta {\eta }_{zs}}{{\omega }_0}+f_{0+},\end{array}$$

(18)

where $f_{0+}={\gamma }^2/{\gamma }_{0+}^2$, ${\gamma }_{0+}=\sqrt{\omega {\omega }_0}$ is the critical coupling of the superradiant phase in the standard TC model, and $\Delta {\eta }_{zs}={\eta }_z-\eta$. Substituting the value of $j_{zs}$ in the definitions Equation (9), we obtain the stationary points

$$\begin{array}{c}(p_s,q_s,j_{zs},{\varphi }_s)=\left(\frac{\gamma }{\omega }\sqrt{1-\frac{1}{f_0^2}}sin{\varphi }_s,-\frac{\gamma }{\omega }\sqrt{1-\frac{1}{f_0^2}}cos{\varphi }_s,-\frac{1}{f_0},\mathrm{indeterminated}\right).\end{array}$$

(19)

In other words, there is a continuum of fixed points, associated with the conservation of the total number of excitations which makes the standard TC Hamiltonian integrable. This leads to the standard result from the TC model where the energy surface takes the form of the Mexican hat potential [61], as shown in Figure 1(c3). These points are valid when $f_0\ge 1$, so the value of $j_{zs}$ remains real. Thus, there is a critical coupling given by

$$\begin{array}{c}{\gamma }_0^c={\gamma }_{0+}\sqrt{1-\frac{\Delta {\eta }_{zs}}{{\omega }_0}},\end{array}$$

(20)

where we obtain the standard critical coupling of the TC model modified by a factor that depends on the qubit–qubit interactions (for $\Delta {\eta }_{zs}\ge 0$ the critical coupling ${\gamma }_0^c$ becomes zero). This set of points has an energy

$$\begin{array}{c}{ϵ}_{s0\varphi }=-\frac{1}{2}\left(f_0+\frac{1}{f_0}\right)+\frac{{\eta }_z}{2{\omega }_0}.\end{array}$$

(21)

This value is obtained from $ϵ=H_{cl}(q_s,p_s,j_{zs},{\varphi }_s)/{\omega }_0$. As it can be straightforwardly seen, for this parameter domain, ${ϵ}_{s0\varphi }<{ϵ}_{-}$. Thus, if one calculates the Hessian matrix (see Appendix C), one can identify these points as a set of minima. Instead, for $f_0\ge 1$ the point at $j_{zs}=-1$ becomes a local maximum, while the one at $j_{zs}=+1$ remains the absolute maximum. We notice that the ground-state energy is continuous, so ${ϵ}_{s\varphi }={ϵ}_{-}$ at $f_0=1$. According to the form of ground-state energy, we have $q_s\ne 0$ and $p_s\ne 0$. As a result, the domain where $f_0\ge 1$ is recognized as a superradiant phase. One of the major effects of the interactions with respect to the standard TC Hamiltonian is the shift in the critical coupling: both the interactions in the z and x (y) directions change it. Moreover, we notice that the critical coupling ${\gamma }_0^c$ becomes zero when ${\eta }_z-\eta ={\omega }_0$. This means that, for interacting values where ${\eta }_z-\eta \ge {\omega }_0$, there is only a superradiant phase, but not a normal phase for every value of the coupling! Thus, the onset of the superradiant phase can be suppressed or stimulated by choosing the right value of the relative interactions in the z and perpendicular directions.

Interestingly, we observe that the rotational symmetry is not broken in this case because the qubit–qubit interactions are balanced in x and y directions (${\eta }_x={\eta }_y$). Hence, we call the quantum phase existing for $f_0\ge 1$ and ${\eta }_x\ne {\eta }_y$ the superradiant-symmetric phase. Next, we consider the imbalanced case (${\eta }_x\ne {\eta }_y$), which leads to two different, but symmetric to each other, superradiant phases.

4.2.2. Superradiant-x Phase

Now, we consider ${\eta }_x\ne {\eta }_y$ ($\Delta {\eta }_{zx}\ne \Delta {\eta }_{zy}$) and $cos{\varphi }_s=\pm 1$ ($sin{\varphi }_s=0$). Here, we obtain $p_s=0$, and from Equation (10) we obtain

$$\begin{array}{c}{j}_{zs}=-{\left(\frac{{\eta }_z-{\eta }_x}{{\omega }_0}+\frac{{\gamma }^2}{\omega {\omega }_0}\right)}^{-1}=-\frac{1}{f_{0x}},f_{0x}=\frac{\Delta {\eta }_{zx}}{{\omega }_0}+f_{0+},\end{array}$$

(22)

which looks exactly as in the previous case with a critical coupling given by

$$\begin{array}{c}{\gamma }_{0x}^c={\gamma }_{0+}\sqrt{1-\frac{\Delta {\eta }_{z,x}}{{\omega }_0}}.\end{array}$$

(23)

where $\Delta {\eta }_{zx}={\eta }_z-{\eta }_x$. Similarly, for $\Delta {\eta }_{zx}\ge {\omega }_0$, the critical coupling ${\gamma }_{0x}^c$ becomes zero. Unlike the previous case, however, there is not an infinite set of stationary points, but only two degenerated ones. This is the most common case for all the superradiant phases we will see in the following for arbitrary $\xi$. Substituting in Equation (9), we derive $q_s$, so the fixed points are given by

$$\begin{array}{c}(p_s,q_s,j_{zs},{\varphi }_s)=\left(0,\mp \frac{\gamma }{\omega }\sqrt{1-\frac{1}{f_{0x}^2}},-\frac{1}{f_{0x}},\pi or0\right),\end{array}$$

(24)

with energy

$$\begin{array}{c}{ϵ}_{s0x}=-\frac{1}{2}\left(f_{0x}+\frac{1}{f_{0x}}\right)+\frac{{\eta }_z}{2{\omega }_0},\end{array}$$

(25)

which lies always below ${ϵ}_{-}$, as in the superradiant-symmetric case. Again, the expectation value of the photon number operator becomes different from zero in the thermodynamic limit, so the phase where these points exist corresponds to a superradiant one. Its emergence is determined by ${\eta }_x$, independently of ${\eta }_y$, however. For this reason, we call it superradiant-x phase. Here, the point at $j_{zs}=-1$ becomes saddle point, as can be observed in Figure 1(c4)–(c6) when increasing $\gamma$.

4.2.3. Superradiant-y Phase

This is identical to the x case, but in the y direction. If we consider ${\eta }_x\ne {\eta }_y$ ($\Delta {\eta }_{zx}\ne \Delta {\eta }_{zy}$) and $cos{\varphi }_s=0$ ($sin{\varphi }_s=\pm 1$), now $q_s=0$ and

$$\begin{array}{c}{j}_{zs}=-{\left(\frac{{\eta }_z-{\eta }_y}{{\omega }_0}+\frac{{\gamma }^2}{\omega {\omega }_0}\right)}^{-1}=-\frac{1}{f_{0y}},f_{0y}=\frac{\Delta {\eta }_{zy}}{{\omega }_0}+f_{0+}\end{array}$$

(26)

where the critical is coupling given by

$$\begin{array}{c}{\gamma }_{0y}^c={\gamma }_{0+}\sqrt{1-\frac{{\eta }_z-{\eta }_y}{{\omega }_0}}.\end{array}$$

(27)

Substituting in Equation (9), we can obtain $p_s$. The two degenerated fixed points are

$$\begin{array}{c}(p_s,q_s,j_{zs},{\varphi }_s)=\left(\pm \frac{\gamma }{\omega }\sqrt{1-\frac{1}{f_{0y}^2}},0,-\frac{1}{f_{0y}},\pm \frac{\pi }{2}\right).\end{array}$$

(28)

Finally, their energy is

$$\begin{array}{c}{ϵ}_{s0y}=-\frac{1}{2}\left(f_{0y}+\frac{1}{f_{0y}}\right)+\frac{{\eta }_z}{2{\omega }_0}.\end{array}$$

(29)

Again, ${ϵ}_{s0y}\le {ϵ}_{-}$. The difference with respect to the previous superradiant phases is that now, the fixed points are rotated in phase space by $\pi /2$, as can be seen from Figure 1(c7)–(c9). As a result, the quadrature $q_s$ is the one that has become zero.

4.2.4. Quantum Phases in the Tavis–Cummings Limit

We have already seen that in the normal phase, there are only two extrema in the energy surface located at $j_{zs}=\pm 1$, whereas, in the superradiant phases described above, we have found four (or a continuous set, in the symmetric case). The energy of the ground-state across the different parameter domains can be expressed in a closed form as

$$\begin{array}{c}{ϵ}_0^{\mathrm{g}.\mathrm{s}.}=-\frac{1}{2}\left(F_0+F_0^{-1}\right)+\frac{{\eta }_z}{2{\omega }_0}\end{array}$$

(30)

with

$$\begin{array}{c}{F}_0=\left\{\begin{array}{cc}{f}_0& \mathrm{f}\mathrm{o}\mathrm{r}{\eta }_x={\eta }_z,\mathrm{a}\mathrm{n}\mathrm{d}\gamma \ge {\gamma }_0^c,\\ f_{0x}& \mathrm{f}\mathrm{o}\mathrm{r}{\eta }_x\ne {\eta }_y,\mathrm{a}\mathrm{n}\mathrm{d}\gamma \ge {\gamma }_{0x}^c,\\ f_{0y}& \mathrm{f}\mathrm{o}\mathrm{r}{\eta }_x\ne {\eta }_y,\mathrm{a}\mathrm{n}\mathrm{d}\gamma \ge {\gamma }_{0y}^c,\\ 1& otherwise\end{array}\right.\end{array}$$

(31)

Let us suppose that ${\eta }_y=0$, so there are no interactions in the y direction. Then, as a function of $\gamma$, the system undergoes a superradiant QPT at $\gamma ={\gamma }_{0x}^c$. As a result, the number of fixed points in the energy surface will change from two to four, and the minimum of the energy surface will be modified, reflecting the change in the ground-state energy associated with the QPT. The same will happen if we consider qubit–qubit interactions only in the y and z direction, but not in the x directions.

The situation is different when considering the combination of interactions in the x and y directions. In this case, there is a possibility that the fixed points arising from each direction simultaneously appear. Then, we will speak of a superposition of the phases. However, it is important to emphasize that even though the two phases appear superimposed, only one set of degenerate fixed points will correspond to the minimum of the energy surface, so the passage from a superradiant phase alone to a superposition of phases is not followed by a QPT (although the smooth DoS will abruptly change announcing the onset of new ESQPTs, hence the necessity to distinguish between a superradiant phase alone and a superimposed one). Here, we observe that, if ${\eta }_y\ge {\eta }_x$ ($\Delta {\eta }_{zx}\ge \Delta {\eta }_{zy}$), then we have that ${\gamma }_{0x}^c\le {\gamma }_{0y}^c$. Without loss of generality, we will take this as the standard case for most of our expressions (otherwise, the superradiant-x and y phases will exchange places). This condition separates the parameter domains in three zones: The normal phase $\gamma \in [0,{\gamma }_{0x}^c]$, the superradiant-x phase $\gamma \in [{\gamma }_{0x}^c,{\gamma }_{yx}^c]$, and a superposition of the superradiant-x and the superradiant-y phases $\gamma \in [{\gamma }_{0y}^c,\infty )$. Hence, depending on the values $\Delta {\eta }_{zx}$ and $\Delta {\eta }_{zx}$, the energy surface could have up to six stationary points, where only one (normal), two (superradiant), or an infinite set (superradiant-symmetric) can correspond to the ground-state. We also notice that for $\Delta {\eta }_{zx}\ge 0$ ($\Delta {\eta }_{zy}\ge 0$), the system enters into the superposition of phases for every value of $\gamma$, going from four to six stationary points in the energy surface, and with a ground state determined by the relationship between ${\gamma }_{0x}^c$ and ${\gamma }_{0y}^c$.

Thus, while in the normal phase there are only two relevant energies such that ${ϵ}_{-}\le {ϵ}_{+}$, and in the superradiant-symmetric phase—three ${ϵ}_{s0\varphi }\le {ϵ}_{-}\le {ϵ}_{+}$, in the superradiant-x and y, we will have four: ${ϵ}_{s0x}<{ϵ}_{s0y}\le {ϵ}_{-}\le {ϵ}_{+}$. This will become important later in Section 5 when discussing ESQPTs. Notice that the superradiant-symmetric phase exists for $\gamma \in [{\gamma }_0^c,\infty )$ because ${\eta }_x={\eta }_y$. Moreover, now we can explain the directions of the energy surface deformation in the normal phase. It turns out that if ${\gamma }_{0x}^c<{\gamma }_{0y}^c$, the deformation occurs in the x direction, as it is shown in Figure 1(c4) even if ${\eta }_x=0$. The opposite is true: if ${\gamma }_{0y}^c<{\gamma }_{0x}^c$, the deformation occurs in the y direction, as shown in Figure 1(c7).

The ground-state energy is continuous at the critical values of the light–matter coupling $\gamma$. Nevertheless, the derivatives are discontinuous. The order of the discontinuity allows us to classify the type of quantum phase transitions the system exhibits according to Ehrenfest’s classification of phase transitions. To do so, we calculate the gradient of the ground-state energy as a function of the interactions

$$\begin{array}{c}\nabla {ϵ}_0^{\mathrm{g}.\mathrm{s}.}=\left(\frac{\partial {ϵ}_0^{\mathrm{g}.\mathrm{s}.}}{\partial \gamma },\frac{\partial {ϵ}_0^{\mathrm{g}.\mathrm{s}.}}{\partial {\eta }_x},\frac{\partial {ϵ}_0^{\mathrm{g}.\mathrm{s}.}}{\partial {\eta }_y},\frac{\partial {ϵ}_0^{\mathrm{g}.\mathrm{s}.}}{\partial {\eta }_z}\right)=\\ \frac{1}{2{\omega }_0}\left\{\begin{array}{cc}\frac{1-f_{0x}^2}{f_{0x}^2}\left(\frac{2{\omega }_0}{\gamma }{f}_{0x},-1,0,1\right)+(0,0,0,1)& \mathrm{f}\mathrm{o}\mathrm{r}\gamma \ge {\gamma }_{0x}^c,\\ \frac{1-f_{0y}^2}{f_{0y}^2}\left(\frac{2{\omega }_0}{{\gamma }_{0+}}{f}_{0y},0,-1,1\right)+(0,0,0,1)& \mathrm{f}\mathrm{o}\mathrm{r}\gamma \ge {\gamma }_{0y}^c,\\ (0,0,0,1)& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.\end{array}$$

(32)

We need to evaluate the derivatives at three specific combinations of the parameters: $F_0^c=1$, $\Delta {\eta }_{zx}={\omega }_0$ ($\Delta {\eta }_{zy}={\omega }_0$) and $\Delta {\eta }_{zx}=\Delta {\eta }_{zy}$. At the critical light–matter coupling, we have $F_0^c=1$, so the ground-state energy from the normal to the superradiant phases as a function of the parameter $\gamma$ is continuous in the zeroth and first order, and only discontinuous at the second one. Thus, a second-order phase transition occurs from normal to superradiant, as expected from the standard TC model. As a function of ${\eta }_x$ and ${\eta }_y$, there is a first-order quantum phase transition at $\Delta {\eta }_x=\Delta {\eta }_y$ because it is the border between the superradiant-x and superradiant-y phases, i.e., the ground-state energy goes from being described by Equations (25)–(29), passing through Equation (21) (see Figure 1b). Finally, there are other two first-order QPTs where the system goes directly from the normal to the superradiant phase when $\Delta {\eta }_{zx}={\omega }_0$ or $\Delta {\eta }_{zx}={\omega }_0$ and $\gamma =0$. This first-order phase transition was identified in Refs. [53,91]. Here, we have found a generalization discovering that the relevant parameter is not just ${\eta }_z$, but $\Delta {\eta }_{zx}$ (or $\Delta {\eta }_{zx}$) instead. This is not surprising because given that the pseudospin length is conserved, a single direction can be expressed in terms of the others, i.e., ${\widehat{J}}_z^2=j(j+1)\widehat{\mathbb{I}}-{\widehat{J}}_x^2-{\widehat{J}}_y^2$, so the qubit-interacting Hamiltonian can be written as

$$\begin{array}{c}{\widehat{H}}_{qq}=-\frac{1}{N}\left(\Delta {\eta }_{zx}{\widehat{J}}_x^2+\Delta {\eta }_{zy}{\widehat{J}}_y^2\right)+\frac{1}{N}{\eta }_{z}j(j+1)\widehat{\mathbb{I}}\end{array}$$

(33)

Hence, the role of x or y interactions is to shift the critical coupling of the superradiant QPT and to shift the value of the interactions where the first-order phase transition emerges. In Figure 1a, we show the quantum phases in the $\Delta {\eta }_{zx}$ versus $\gamma$ space. There, we have included an artificial shift ${\omega }_0/2$ between ${\eta }_x$ and ${\eta }_y$ to exhibit the onset of the superradiant-x followed by a superposition of phases as the light–matter coupling increases, otherwise, the phases overlap in the diagram (because of the symmetry). In Figure 2b, we show the quantum phases in the $\Delta {\eta }_{zy}$ versus $\Delta {\eta }_{zx}$ space. This diagram depends on the value of $\gamma$; here, we choose an illustrative value to always be in the superradiant phases.

Further information about the system can be obtained by studying the energy surfaces as they are shown in Figure 1c–f. We employ representative values of the interaction strengths in each direction to explore and better show the evolution of the surface as the interacting parameters change. In the first row of Figure 1(c1)–(c3), we present a configuration ${\eta }_x={\eta }_y$, where we expect symmetry in the system. Increasing the light–matter coupling makes the surface go from a spherical well with a minimum at the center to the Mexican hat potential characteristic of the superradiant-symmetric phase, as shown in Figure 1(c1). Next, the symmetry of the surface in both the normal and superradiant phases is broken once we increase the interactions in a given direction other than z. In the middle row of Figure 1c, we make ${\eta }_y=0.9{\omega }_0$ and ${\eta }_x={\eta }_z=0$. As mentioned before, here, we identify a stretching of the energy surface in the x direction (because ${ϵ}_{s0x}<{ϵ}_{s0y}$). In Figure 1(c4), we select a coupling $\gamma <{\gamma }_{0x}^c$ that locates us within the normal phase. Eventually, increasing $\gamma$ leads us to the superradiant-x phase with two minima and a saddle point, as shown in Figure 1(c5). Then, in Figure 1(c6), we enter a regime of superposition between the two superradiant phases. Now, there are two minima points, two saddle points, and a maximum in the center. The situation is the same in the third row. From Figure 1(c7)–(c9), the stretching of the energy surface is in the y-direction, and we will recover the phenomenology of the last case, but rotated in $\pi /2$, characteristic of the superradiant-y phase. Here, for larger $\gamma$, the dominant phase is the superradiant-y as ${ϵ}_{s0y}<{ϵ}_{s0x}$. The reasoning is the same in the other examples we show in Figure 1d–f. In Figure 1d, we fix one of the directions to zero and tune the other two equal to 1. In the first two rows, we show the evolution of the surface as we increase $\gamma$, while making ${\eta }_x={\eta }_z=0.9$ (${\eta }_y={\eta }_z=0.9$). The situation becomes similar to that of Figure 1 because the relevant parameters to determine the phases are $\Delta {\eta }_{zx}$ and $\Delta {\eta }_{zy}$. In Figure 1e we choose a different combination of values for the interactions in each direction. Finally, in Figure 1f, we select negative values of the interactions. Given a set of interacting parameters ${\eta }_i$, we can identify how the energy surface will be deformed, and increasing $\gamma$ leads us always to the superposition of phases x-y. In each figure, we will have different couplings $\gamma$ to highlight the system’s phase.

Next, we will explore the other limit of the Hamiltonian, when the non-resonant terms are completely included, i.e., the Dicke limit.

4.3. Dicke Limit

We now fully include the counter-rotating terms, i.e., we set $\xi =1$. This corresponds to the Dicke limit, meaning we obtain a Dicke Hamiltonian plus the qubit–qubit interactions. It is given by

$$\begin{array}{c}{H}_{cl}^{\left(1\right)}=\frac{\omega }{2}(q^2+p^2)+j_z\left({\omega }_0+\frac{{\eta }_z{j}_z}{2}\right)+\frac{1}{2}\left(1-j_z^2\right)\left({\eta }_x{cos}^2\varphi +{\eta }_y{sin}^2\varphi \right)+\\ +2\gamma \sqrt{1-j_z^2}qcos\varphi .\end{array}$$

(34)

Equations (10) and (11) become

$$\begin{array}{c}{\omega }_0+j_{zs}\left({\eta }_z-\left({\eta }_x{cos}^2{\varphi }_s+{\eta }_y{sin}^2{\varphi }_s\right)+\frac{4{\gamma }^2}{\omega }{cos}^2{\varphi }_s\right)=0,\end{array}$$

(35)

$$\begin{array}{c}(1-j_{zs}^2)\left[({\eta }_x-{\eta }_y)-\frac{4{\gamma }^2}{\omega }\right]cos{\varphi }_{s}sin{\varphi }_s=0.\end{array}$$

(36)

Just as in the case of the TC limit, we obtain five possibilities when searching stationary points of the energy surface. The first four are identical: $j_{zs}=\pm 1$, $cos{\varphi }_s=0$ ($sin{\varphi }_s=\pm 1$), $sin{\varphi }_s=0$ ($cos{\varphi }_s=\pm 1$). However, the last one contains a different condition for the parameters, given by $4{\gamma }^2/\omega ={\eta }_x-{\eta }_y$. We immediately notice that the symmetry that was found in the TC limit ($\xi =0$) is now always broken because the relation ${\eta }_x={\eta }_y$ does not lead to the existence of stationary points anymore. This is an expected result, given that the standard Dicke model is non-integrable [47]. Moreover, from Hamilton equations, we observe that $p_s$ is always zero in this limit.

We have a normal phase in the Dicke limit too, where the fixed point at $j_{zs}=-1$ is an absolute minimum corresponding to the ground state, and the point at $j_{zs}=+1$ is an absolute maximum. Both points have an energy given by ${ϵ}_{\pm }=\pm 1+{\eta }_z/2{\omega }_0$ and the normal phase will be deformed if the interactions are privileged in either the x or y, as it happened in the $\xi =0$ limit. This case is shown in Figure 2(c1),(c4),(c7), where we encounter the same deformations as in the TC limit. Next, we will explore the other phases appearing for $\xi =1$.

4.3.1. Superradiant-x Phase

When we have $cos{\varphi }_s=\pm 1$ ($sin{\varphi }_s=0$), we encounter the same superradiant phase as in the $\xi =0$ limit, but with a modified critical coupling, as it would be expected [61]. Here, the value of the collective atomic variable becomes

$$\begin{array}{c}{j}_{zs}=-\frac{1}{f_{1x}},f_{1x}=\frac{\Delta {\eta }_{zx}}{{\omega }_0}+f_{1+}\end{array}$$

(37)

where $f_{1x}={\gamma }^2/{\gamma }_{1+}^2$ and ${\gamma }_{1+}=\sqrt{\omega {\omega }_0}/2$ is the critical coupling of the standard Dicke model. Using Equation (9) one gets the two degenerate minima typical of the superradiant phase

$$\begin{array}{c}\left(p_s,q_s,j{z}_s,{\varphi }_s\right)=\left(0,\pm \frac{2\gamma }{\omega }\sqrt{1-\frac{1}{f_{1x}^2}},-\frac{1}{f_{1x}},0\mathrm{or}\pi \right)\end{array}$$

(38)

The energy associated to these points is

$$\begin{array}{c}{ϵ}_{s1x}=-\frac{1}{2}\left(f_{1x}+\frac{1}{f_{1x}}\right)+\frac{{\eta }_z}{2{\omega }_0}.\end{array}$$

(39)

An almost identical result as that of the TC-like superradiant-x phase. The ground-state energy of the system is always below ${ϵ}_{-}$; thus, the energy of these points corresponds to the ground-state energy. This phase also is valid for values of the light–matter coupling $\gamma \ge {\gamma }_{1x}^c$, where

$$\begin{array}{c}{\gamma }_{1x}^c={\gamma }_{1+}\sqrt{1-\frac{\Delta {\eta }_{zx}}{{\omega }_0}}.\end{array}$$

(40)

Suppose the interactions ${\eta }_{x,z}$ vanish. In that case, we recover the standard result for the Dicke model, again, we observe that the role of the interactions is to shift the critical coupling. This effect has been observed before in several previous works for interactions in the x [32,34] and z [58,59,89,91] directions, and a combination of x and y directions [49]. Not only in the case of the Dicke limit but for arbitrary $\xi$, the critical coupling to attain superradiance can become zero with a suitable choice of the interactions in the z and x (y) directions, given that the relevant parameter is the difference between the $\Delta {\eta }_{zx}$ ($\Delta {\eta }_{zy}$) relative interactions, as we have seen before.

A specific feature of the $\xi =1$ limit is that breaking the rotational symmetry leads to the exclusion of the superradiant-y phase. Instead, we have the deformed phase, as we will immediately explain.

4.3.2. Deformed Phase

Instead of the superradiant-y phase, there is a new quantum phase emerging from the interactions. From the condition where $cos{\varphi }_s\pm 1$ ($sin{\varphi }_s=0$), we obtain

$$\begin{array}{c}{j}_{zs}=-\frac{{\omega }_0}{\Delta {\eta }_{zy}}=-\frac{1}{f_{1y}},\end{array}$$

(41)

This leads to two new fixed points

$$\begin{array}{c}\left(p_s,q_s,j{z}_s,{\varphi }_s\right)=\left(0,0,-\frac{1}{f_{1y}},\pm \frac{\pi }{2}\right),\end{array}$$

(42)

whose energy is given by

$$\begin{array}{c}{ϵ}_{s1y}=-\frac{1}{2f_{1y}}+\frac{{\eta }_y}{2{\omega }_0},\end{array}$$

(43)

However, if we add and subtract ${\eta }_z/2{\omega }_0$, it reads

$$\begin{array}{c}{ϵ}_{s1y}=-\frac{1}{2}\left(f_{1y}+\frac{1}{f_{1y}}\right)+\frac{{\eta }_z}{2{\omega }_0},\end{array}$$

(44)

The main difference with previous cases (the superradiant phases) is that $f_{1y}$ is independent of $\gamma$, so this phase exists for every value of the light–matter coupling given that $\Delta {\eta }_{zy}\ge {\omega }_0$. This phase was identified before in Ref. [55] in the absence of ${\eta }_y$. It is characterized by the two degenerate fixed points whose orientation in the atomic angle $\varphi$ is rotated by $\pi /2$ with respect to the superradiant-x phase’s fixed points (a result inherited from the superradiant-y phase). There, the photon number’s expectation value becomes zero because $q_s=p_s=0$. Furthermore, we have that because $|f_{1y}|\le 1$, ${ϵ}_{s1y}\le {ϵ}_{-}<{ϵ}_{+}$. Then, they mark the ground-state energy, and the point at $j_{zs}=-1$ becomes a saddle point, as it happens in the usual superradiant phases. However, because it is neither a normal, nor a superradiant phase, we deem it as a deformed phase, although one could name it a subrradiant phase. The energy surfaces in this phase are shown in Figure 2(e7),(f7).

Finally, we study the last condition for stationary points, given by the parameter relation

$$\begin{array}{c}\frac{{\eta }_x-{\eta }_y}{{\omega }_0}=\frac{\Delta {\eta }_{zy}-\Delta {\eta }_{zx}}{{\omega }_0}=\frac{{\gamma }^2}{{\gamma }_{1+}^2}=f_{1+}\end{array}$$

(45)

Since we can write Equation (35) as

$$\begin{array}{c}{\omega }_0+j_{zs}\left[{\eta }_z-{\eta }_y-{cos}^2\varphi \left({\eta }_y-{\eta }_x+\frac{4{\gamma }^2}{\omega }\right)\right]=0,\end{array}$$

(46)

it is clear that, when applying the condition in Equation (45), the factor multiplying the cosine function vanishes, so we obtain a stationary point. Here, we find the following value of $j_{zs}$ with energy

$$\begin{array}{c}\left(p_s,q_s,j{z}_s,{\varphi }_s\right)=\left(0,0,-\frac{1}{f_{1y}},\pm \frac{\pi }{2}\right),\\ {ϵ}_{s1y}=-\frac{1}{2}\left(f_{1y}+\frac{1}{f_{1y}}\right)+\frac{{\eta }_z}{2{\omega }_0},\end{array}$$

(47)

i.e., the stationary points from the deformed phase. Nevertheless, we can write Equation (45) as

$$\begin{array}{c}\frac{\Delta {\eta }_{zy}}{{\omega }_0}=\frac{\Delta {\eta }_{zx}}{{\omega }_0}+f_{1+}.\end{array}$$

(48)

Therefore, ${ϵ}_{s1x}={ϵ}_{s1y}$ and the stationary points coincide with those of the superradiant-x phase. This means that Equation (45) marks the frontier between the superradiant-x, the deformed phases, and the normal phase (for $f_{1y}=1$), a similar result to ${\eta }_x={\eta }_y$ for $\xi =0$.

4.3.3. Quantum Phases in the Dicke Limit

The existence of the deformed phase changes the quantum phase diagram in the Dicke limit. In the TC, there is a chance of a superposition of the two superradiant phases. We recall that, in this case, only one set of stationary points becomes the minimum, either those from the superradiant-x, or those from the superradiant-y phases. Instead, there is a stricter separation between phases because the deformed phase appears for $\Delta {\eta }_{zy}\ge {\omega }_0$. Then, we have only two stationary points at the normal phase and only four for both the superradiant-x and deformed phases. We can write the ground-state energy in closed form too:

$$\begin{array}{c}{ϵ}_1^{\mathrm{g}.\mathrm{s}.}=-\frac{1}{2}\left(F_1+F_1^{-1}\right)+\frac{{\eta }_z}{2{\omega }_0}\end{array}$$

(49)

with

$$\begin{array}{c}{F}_1=\left\{\begin{array}{cc}{f}_{1x}& \mathrm{f}\mathrm{o}\mathrm{r}\gamma \ge {\gamma }_{1x}^c\mathrm{a}\mathrm{n}\mathrm{d}\Delta {\eta }_{zx}\le {\omega }_0,\\ f_{1y}& \mathrm{f}\mathrm{o}\mathrm{r}\Delta {\eta }_{zy}\ge {\omega }_0,\\ 1& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.\end{array}$$

(50)

We also recall that $f_{1y}$ is independent from $\gamma$. Again, if we suppose ${\eta }_y={\eta }_z=0$, the ground state evolves as a function of $\gamma$ from the normal $\gamma \in [0,{\gamma }_{1x}^c]$ to the superradiant-x phase $\gamma \in [{\gamma }_{1x}^c,\infty )$. Then, the situation remains the same for $\Delta {\eta }_{zy}\ne 0$, but the deformed phase emerges above ${\omega }_0$. Similarly to the TC limit, we will have the following energy intervals ${ϵ}_{s1x}<{ϵ}_{s1y}<{ϵ}_{-}<{ϵ}_{+}$ in the superradiant-x phase and ${ϵ}_{s1y}<{ϵ}_{-}<{ϵ}_{+}$ in the deformed phase.

Next, we obtain the gradient of the ground-state energy as a function of the interactions, just like we did in the TC regime:

$$\begin{array}{c}\nabla {ϵ}_1^{\mathrm{g}.\mathrm{s}.}=\frac{1}{2{\omega }_0}\left\{\begin{array}{cc}\frac{1-f_{1x}^2}{f_{1x}^2}\left(\frac{2{\omega }_0}{\gamma }{f}_{1x},-1,0,1\right)+(0,0,0,1)& \mathrm{f}\mathrm{o}\mathrm{r}\gamma \ge {\gamma }_{1x}^c,\mathrm{a}\mathrm{n}\mathrm{d}\Delta {\eta }_{zx}\le {\omega }_0,\\ \frac{1-f_{1y}^2}{f_{1y}^2}\left(0,0,-1,1\right)+(0,0,0,1)& \mathrm{f}\mathrm{o}\mathrm{r}\Delta {\eta }_{zy}\ge {\omega }_0,\\ (0,0,0,1)& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.,\end{array}$$

(51)

There are three sets of parameter values to look for the presence of QPT. First, $F_1^c=1$, where we obtain a generalization of well-known second-order Dicke QPT from the normal to the superradiant-x phase. Next, at $\Delta {\eta }_{zy}^c={\omega }_0$, we obtain a first-order QPT from the normal and superradiant-x phases to the deformed one as a function of $\Delta {\eta }_{zx}$ or $\Delta {\eta }_{zy}$, recovering the results in Ref. [55], when ${\eta }_x={\eta }_y=0$. Finally, at $\Delta {\eta }_{zy}-\Delta {\eta }_{zx}={\omega }_0{f}_{1+}$, we have the corresponding behavior we found for the TC limit at ${\eta }_x={\eta }_y$: There is a first-order QPT signaling the border between the x and y sides. In Figure 2a,b, we show the different phases in the system as a function of the Hamiltonian parameters. The presence of first-order phase transitions in the Dicke model due to qubit–qubit interactions were predicted before considering interactions in the y direction [57], and later in the z [58,60,91], and are inherited from the LMG model as well.

The energy surfaces for the same parameter configurations we employed in the TC limit are shown in Figure 2c–f. In some cases, we slightly change the value of the interactions to highlight the phase in which the system is located. We observe the onset of the fixed points and what kind of extreme point they correspond to (minimum, maximum, saddle point) in each phase as a function of the interacting parameters. For ${\eta }_x={\eta }_y=0$ (Figure 2(c1)–(c3)), the system goes from the normal to the superradiant-x case, but the normal phase is not deformed. As usual, changing the balance between ${\eta }_x$ and ${\eta }_y$ breaks the overall symmetry of the normal phase but, in this case, it leaves unaffected the superradiant phase. In the Dicke limit, we have only a restriction for $\gamma$ given by ${\gamma }_{1x}^c$, as the deformed phase is independent of the light–matter coupling. As a direct consequence, the deformation tends to stretch in the horizontal direction, but we have cases as those in Figure 2(c7),(d7), where the energy surface in the normal phase is deformed in the y direction. We can observe that for most sets of parameters, the situation is similar to that of Figure 2(c4)–(c6). In (c4) $\gamma <{\gamma }_{1x}^c$, the system is in the normal phase. Increasing the parameter $\gamma$ makes the system enter the superradiant-x phase where the minima are in the horizontal direction, and a saddle point appears. We can identify the appearance of the deformed phase in (e7) and (f7), where the minima emerge in the vertical direction, and the point at the center becomes a saddle point. Finally, we stress that there is no situation where more than four stationary points appear, contrasting with the TC limit.

4.4. Arbitrary Coupling

We are ready to treat the general case, i.e., when $\xi \in (0,1)$. Here, we expect to find effects similar to those in the TC and Dicke limits. In the absence of qubit interactions, the main difference with those limits is the presence of two different superradiant domains that can coexist for some values of the light–matter coupling. Hence, three phases separated by the critical values of the light–matter coupling appear ${\gamma }_{\xi \pm }=\sqrt{\omega {\omega }_0}/(1\pm \xi )$, such that ${\gamma }_{\xi +}<{\gamma }_{\xi -}$ [6,62,69,70]. Given that ${\gamma }_{\xi +}<{\gamma }_{\xi -}$, in the interval ${\gamma }_{\xi +}<\gamma <{\gamma }_{\xi -}$, one finds the standard effect of the Dicke model that we will call here superradiance-$(+)$. Instead, for ${\gamma }_{\xi -}<\gamma$, there are two additional stationary points whose nature is that of saddle points that are attributed to a superradiance-$(-)$ effect. However, in this superposition of superradiant phases, the fixed points from the superradiance-$(+)$ are the minima, so there is no QPT [62]. We immediately notice the similarity with what we have observed in the TC model in the presence of qubit interaction. A result that will take importance later when we explore the energy domains. If we take $\xi \to 0$, then ${\gamma }_{0+}={\gamma }_{0-}$, and for the TC limit, the two superradiant phases become one. We can anticipate that the presence of the interactions ${\eta }_x$ and ${\eta }_y$ creates a similar result and produces two domains separated by the two critical couplings we have already described ${\gamma }_{0x}^c$ and ${\gamma }_{0y}^c$. On the other hand, when $\xi \to 1$ the critical coupling ${\gamma }_{1-}\to \infty$, the new superradiant-$(-)$ phase is pushed to larger values of the coupling until it vanishes. Hence, in the Dicke model, there is only one superradiant phase. In our case, this effect has been reflected on the onset of the deformed phase. As we will discuss below, by including the qubit–qubit interactions, we obtain a similar result for arbitrary $\xi$, where the superradiant-$(+)$ [($(-)$)] phase is modified by interactions in x (y) directions.

Once more, we obtain five conditions for fixed points from Equations (10) and (11): $j_{zs}=\pm 1$, $cos{\varphi }_s=\pm 1$ ($sin{\varphi }_s=0$), $sin{\varphi }_s=\pm 1$ ($cos{\varphi }_s=0$), and the special parameter relationship that now takes the form:

$$\frac{{\eta }_x-{\eta }_y}{{\omega }_0}=\frac{\Delta {\eta }_{zy}-\Delta {\eta }_{zx}}{{\omega }_0}=\frac{{\gamma }^2}{{\gamma }_{\xi +}^2}-\frac{{\gamma }^2}{{\gamma }_{\xi -}^2}=f_{\xi +}-f_{\xi -},$$

(52)

$j_{zs}=\pm 1$ leads to the two stationary points that mark the absolute minimum in the normal phase and the absolute maximum. Thus, we have again a (deformed) normal phase, as in the two previous cases, Next, we will recover the most general superradiant phases.

4.4.1. Superradiant-x and y Phases

If we evaluate the condition $cos{\varphi }_s=\pm 1$ ($sin{\varphi }_s=0$), we obtain the two degenerate stationary points corresponding to the superradiant ground state given by

$$\begin{array}{c}\left(p_s,q_s,j{z}_s,{\varphi }_s\right)=\left(0,\mp \frac{2\gamma }{\omega }\sqrt{1-\frac{1}{f_{\xi x}^2}},-\frac{1}{f_{\xi x}},0\mathrm{o}\mathrm{r}\pi \right),\end{array}$$

(53)

where now

$$\begin{array}{c}{f}_{\xi x}=\frac{\Delta {\eta }_{zx}}{{\omega }_0}+f_{\xi +}.\end{array}$$

(54)

At these points, the energy surface becomes

$$\begin{array}{c}{ϵ}_{s\xi x}=-\frac{1}{2}\left(f_{\xi x}+\frac{1}{f_{\xi x}}\right)+\frac{{\eta }_z}{2{\omega }_0}\end{array}$$

(55)

and they exist for $\gamma \ge {\gamma }_{\xi x}^c$ with

$$\begin{array}{c}{\gamma }_{\xi x}^c={\gamma }_{\xi +}\sqrt{1-\frac{\Delta {\eta }_{zx}}{{\omega }_0}}.\end{array}$$

(56)

Symmetrically, if we opt for the case $sin{\varphi }_s=\pm 1$ ($cos{\varphi }_s=0$), the stationary points are rotated by $\pi /2$, as expected,

$$\begin{array}{c}\left(p_s,q_s,j{z}_s,{\varphi }_s\right)=\left(\pm \frac{2\gamma }{\omega }\sqrt{1-\frac{1}{f_{\xi y}^2}},0,-\frac{1}{f_{\xi y}},\pm \frac{\pi }{2}\right),\end{array}$$

(57)

where

$$\begin{array}{c}{f}_{\xi y}=\frac{\Delta {\eta }_{zy}}{{\omega }_0}+f_{\xi -}.\end{array}$$

(58)

In the same way, these points appear only for $\gamma \ge {\gamma }_{\xi y}$ with

$$\begin{array}{c}{\gamma }_{\xi y}^c={\gamma }_{\xi -}\sqrt{1-\frac{\Delta {\eta }_{zy}}{{\omega }_0}}.\end{array}$$

(59)

Their energy is

$$\begin{array}{c}{ϵ}_{s\xi y}=-\frac{1}{2}\left(f_{\xi y}+\frac{1}{f_{\xi y}}\right)+\frac{{\eta }_z}{2{\omega }_0}.\end{array}$$

(60)

As anticipated, these phases correspond to a generalization of the superradiant$-x$ and y we have found in the TC and Dicke limits. Similar to the case of the Dicke model, the condition in Equation (52) lies at the border between the superradiant phases, as it is exhibited in Figure 3b.

4.4.2. Quantum Phases for Arbitrary Coupling

For arbitrary $\xi$, several of the results we have found before for $\xi =0,1$ stand. The major difference was the existence of the deformed phase in the Dicke. In fact, the phase diagram for $\xi \in (0,1)$ is very similar to the one for the TC limit, except for the absence of symmetry (for intermediate $\xi$ the Hamiltonian is non-integrable) and for the dependence of the border between the superradiant x and y phases given by Equation (52) that now is explicitly in terms of $\xi$ and $\gamma$.

What is truly new about the intermediate $\xi$ case is the correspondence between the interactions in the x direction and the superradiant-$(+)$ phase and those in the y direction and the superradiant-$(-)$ phase. This effect explains our previous findings. For $\xi =0$, both phases are completely analogous except that one depends on ${\eta }_x$ and the other on ${\eta }_y$. Instead, for $\xi =1$ the superradiant-y phase vanishes completely because ${\gamma }_{\xi -}$ goes to infinity. ($f_{\xi y}\to \Delta {\eta }_{zy}/{\omega }_0$). Therefore, in the Dicke limit, the superradiant-y transforms into the deformed phase. The correspondence is explicit once we recognize the dependence of the light–matter critical couplings on the interactions in the x and y direction.

Surprisingly, qubit interactions can shift the order in which the two superradiant phases x and y occur, a situation we have already encountered in the TC limit, but excluded in the standard case, i.e., for the superradiant-$(+)$ and $(-)$ phases. For ${\eta }_x={\eta }_y=0$, it holds that ${\gamma }_{\xi x}^c\le {\gamma }_{\xi y}^c$ (because ${\gamma }_{\xi +}^c\le {\gamma }_{\xi -}^c$). However, if we tune ${\eta }_x$ and ${\eta }_y$ independently, we can obtain that ${\gamma }_{\xi y}^c<{\gamma }_{\xi x}^c$. The condition to invert the order of the two critical couplings occurs for the set of parameters where

$$\begin{array}{c}\Delta {\eta }_{zy}=\frac{{\gamma }_{\xi +}^2}{{\gamma }_{\xi -}^2}\Delta {\eta }_{zx}-{\omega }_0\left[1-\frac{{\gamma }_{\xi +}^2}{{\gamma }_{\xi -}^2}\right]\end{array}$$

(61)

For $\xi =0$, this condition becomes ${\eta }_x={\eta }_y$, because ${\gamma }_{0+}^c={\gamma }_{0-}^c$. It is the same special parameter relation we find earlier, leading to the superradiant-symmetric phase. For $\xi =1$, we have that ${\left({\gamma }_{1-}^c\right)}^{-1}=0$. Then, the condition is now $\Delta {\eta }_{zy}={\omega }_0$, the limiting condition where the deformed phase emerges. This effect is shown in Figure 3a for $\xi =0.1$, where we have introduced an artificial shift $\Delta {\eta }_{zy}=\Delta {\eta }_{zx}+{\omega }_0/2$ (the equivalent of what we have done in Figure 1) to exhibit that the two curves of ${\gamma }_{\xi x}^c$ and ${\gamma }_{\xi y}^c$ cross as a function of $\gamma$ thanks to the interactions.

Once more, it is possible to express the ground-state energy in a general and simple form:

$$\begin{array}{c}{ϵ}_1^{\mathrm{g}.\mathrm{s}.}=-\frac{1}{2}\left(F_{\xi }+F_{\xi }^{-1}\right)+\frac{{\eta }_z}{2{\omega }_0}\end{array}$$

(62)

with

$$\begin{array}{c}{F}_{\xi }=\left\{\begin{array}{cc}{f}_{\xi x}& \mathrm{f}\mathrm{o}\mathrm{r}\gamma \ge {\gamma }_{\xi x}^c,\\ f_{\xi y}& \mathrm{f}\mathrm{o}\mathrm{r}\gamma \ge {\gamma }_{\xi y}^c,\\ 1& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.\end{array}$$

(63)

The gradient as a function of the interactions reads

$$\begin{array}{c}\nabla {ϵ}_{\xi }^{\mathrm{g}.\mathrm{s}.}=\frac{1}{2{\omega }_0}\left\{\begin{array}{cc}\frac{1-f_{\xi x}^2}{f_{\xi x}^2}\left(\frac{2{\omega }_0\gamma }{{\gamma }_{\xi +}^2},-1,0,1\right)+(0,0,0,1)& \mathrm{f}\mathrm{o}\mathrm{r}\gamma \ge {\gamma }_{\xi x}^c,\\ \frac{1-f_{\xi y}^2}{f_{\xi y}^2}\left(\frac{2{\omega }_0\gamma }{{\gamma }_{\xi -}^2},0,-1,1\right)+(0,0,0,1)& \mathrm{f}\mathrm{o}\mathrm{r}\gamma \ge {\gamma }_{1y}^c,\\ (0,0,0,1)& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.\end{array}$$

(64)

As a generalization of the TC and Dicke cases, we have a second-order QPT at $F_{\xi }^c=1$, a first-order QPTs between the superradiant-x and y phases at $\Delta {\eta }_{yz}=\Delta {\eta }_{zx}+{\omega }_0\left(f_{\xi +}-f_{\xi -}\right)$, and first-order QPTs from the normal to the superradiant-x (y) phases at the values of the parameters where the other superradiant-y (x) phase becomes prohibited. This is shown in Figure 3a,b. In Figure 3b, the region available for the superradiant-x and y phases change according to the value of $\gamma$. We are selecting the case for $\gamma /{\gamma }_{\xi +}=0.5$ and $\xi =0.1$ as an example. For other values of the light–matter coupling and the parameter $\xi$, we will have larger or smaller domains of validity for the superradiant phases. Similar to the two previous cases, in Figure 3c–f, we also illustrate the different behaviors and nature of the fixed points by showing the energy surfaces using the same interacting parameters as in Figure 1 and Figure 2, but for $\xi =0.5$.

5. Semiclassical Density of States

The study of quantum critical behavior is not limited to the ground state, but extends to excited energies. Next, we explore and classify the distinct energy domains emerging in each phase we have discussed in the previous section by considering the onset of Excited-State Quantum Phase Transitions (ESQPT) as a function of the light–matter and matter–matter interactions. To achieve this aim, we follow the standard methodology that has been developed for the study of ESQPTs, i.e., we analyze the energy dependence and singularities of a semi-classical approximation to the Density of States (DoS) ${\nu }_{\xi }\left(ϵ\right)$, obtained by calculating the available phase space volume at given energy using Weyl’s law [92]

$$\begin{array}{c}{\nu }_{\xi }\left(ϵ\right)=\frac{1}{{\left(2\pi \right)}^2}\int dqdpd{j}_{z}d\varphi \delta \left[ϵ{\omega }_0-H_{cl}^{\xi }(q,p,j_z,\varphi )\right].\end{array}$$

(65)

Although signatures of ESQPTs can be found in the smoothed level flow, the energy densities of some observables, and the oscillatory part of the DoS [4,6], the easiest way to identify them is via the smoothed DoS.

To integrate Equation (65), we need to eliminate first the bosonic degrees of freedom, following closely the methodology in Refs. [61,62], first, we clear the variable q from $H_{cl}^{\xi }(q,p,j_z,\varphi )=ϵ$ in terms of p, $j_z$ and $\varphi$. As one obtains a quadratic equation in q and p, it always yields two roots for every value of the parameters. The two solutions are

$$\begin{array}{c}{q}_{\xi \pm }=-\frac{\gamma }{\omega }\sqrt{1-j_z^2}(1+\xi )cos\varphi \pm \sqrt{-p^2+a_{\xi }p+b_{\xi }},\end{array}$$

(66)

where

$$\begin{array}{c}{a}_{\xi }=\frac{2\gamma }{\omega }\sqrt{1-j_z^2}\left(1-\xi \right)sin\varphi ,\end{array}$$

(67)

$$\begin{array}{c}{b}_{\xi }=-\frac{2}{\omega }{j}_z({\omega }_0+\frac{{\eta }_z{j}_z}{2})-\frac{1}{\omega }\left(1-j_z^2\right)\left({\eta }_x{cos}^2\varphi +{\eta }_y{sin}^2\varphi \right)+\frac{2ϵ{\omega }_0}{\omega }+\\ \frac{{\gamma }^2}{{\omega }^2}(1-j_z^2){(1+\xi )}^2{cos}^2\varphi .\end{array}$$

(68)

Then, we employ the properties of the Dirac delta function to obtain

$$\begin{array}{c}{\nu }_{\xi }\left(ϵ\right)=\frac{1}{{\left(2\pi \right)}^2}\int dqdpd{j}_{z}d\varphi \left[\delta (q-q_{\xi +}){\left|\frac{\partial H_{cl}^{\left(\xi \right)}}{\partial q}\right|}_{q_{\xi +}}^{-1}+\delta (q-q_{\xi -}){\left|\frac{\partial H_{cl}^{\left(\xi \right)}}{\partial q}\right|}_{q_{\xi -}}^{-1}\right].\end{array}$$

(69)

Evaluating the derivatives leads to

$$\begin{array}{c}{\left|\frac{\partial H_{cl}^{\left(\xi \right)}}{\partial q}\right|}_{q_{\xi \pm }}=\left|\omega \left[-\frac{\gamma }{\omega }\sqrt{1-j_z^2}(1+\xi )cos\varphi \pm \sqrt{-p^2+a_{\xi }p+b_{\xi }}\right]+\gamma \sqrt{1-j_z^2}(1+\xi )cos\varphi \right|=\\ \omega \sqrt{-p^2+a_{\xi }p+b_{\xi }}.\end{array}$$

(70)

Thus, the q integration yields

$$\begin{array}{c}{\nu }_{\xi }\left(ϵ\right)=\frac{1}{{\left(2\pi \right)}^2}\frac{2}{\omega }\int d{j}_{z}d\varphi \frac{dp}{\sqrt{-p^2+a_{\xi }p+b_{\xi }}}.\end{array}$$

(71)

The limits in the variables $j_z,\varphi$ and p are determined by the condition $-p^2+a_{\xi }p+b_{\xi }\ge 0$. The p integration is easily performed by writing

$$\begin{array}{c}-p^2+a_{\xi }p+b_{\xi }=(p_{\xi +}-p)(p-p_{\xi -}),\end{array}$$

(72)

where

$$\begin{array}{c}{p}_{\xi \pm }=\frac{1}{2}\left(-a_{\xi }\pm \sqrt{a_{\xi }^2+4b_{\xi }}\right)\end{array}$$

(73)

are the roots $(p_{\xi -}\le p_{\xi +})$ of the quadratic polynomial $-{\omega }^2{p}^2+a_{\xi }p+b_{\xi }=0$. Hence,

$$\begin{array}{c}{\nu }_{\xi }\left(ϵ\right)=\frac{2}{\omega {\left(2\pi \right)}^2}\int d{j}_z\int d\varphi {\int }_{p_{\xi -}}^{p_{\xi +}}dp\frac{1}{\sqrt{(p_{\xi +}-p)(p-p_{\xi -})}}=\frac{1}{2\pi \omega }\int d{j}_z\int d\varphi \end{array}$$

(74)

This result is valid, provided that the roots $p_{\xi \pm }$ are real, which, in turn, occurs only if the discriminant

$$\begin{array}{c}{a}_{\xi }^2+4b_{\xi }\ge 0\end{array}$$

(75)

is greater than or equal to zero. By substituting the values $a_{\xi }$ and $b_{\xi }$, this condition explicitly reads

$$\begin{array}{c}\frac{1}{2}\left(1-j_z^2\right)\left[\left(f_{\xi +}-\frac{{\eta }_x}{{\omega }_0}\right){cos}^2\varphi +\left(f_{\xi -}-\frac{{\eta }_y}{{\omega }_0}\right){sin}^2\varphi \right]\ge \frac{{\eta }_z}{2{\omega }_0}{j}_z^2+j_z-ϵ,\end{array}$$

(76)

or

$$\begin{array}{c}{cos}^2\varphi \ge g_{\xi }(j_z,ϵ),\end{array}$$

(77)

where

$$\begin{array}{c}{g}_{\xi }(j_z,ϵ)=\left\{\frac{2}{1-j_z^2}\left[\frac{{\eta }_z}{2{\omega }_0}{j}_z^2+j_z-ϵ\right]-\left(f_{\xi -}-\frac{{\eta }_y}{{\omega }_0}\right)\right\}{\left[\left(f_{\xi +}-f_{\xi -}\right)-\left(\frac{{\eta }_x}{{\omega }_0}-\frac{{\eta }_y}{{\omega }_0}\right)\right]}^{-1}.\end{array}$$

(78)

We observe these expressions only depend on the phase space volume over the region of the Bloch sphere covered at a given energy. They allow us to determine the limiting values for $(j_z,\varphi )$ in the Bloch sphere, given that $0\le cos{\varphi }_0\le 1$. If $f_{\xi }(j_z,ϵ)<0$, then the condition can be satisfied for all the values of $\varphi \in [0,2\pi )$, covering the Bloch sphere. Instead, if $f_{\xi }(j_z,ϵ)>1$, the condition cannot be fulfilled. It would be valid only within an interval of $\varphi$ given by the limiting angle

$$\begin{array}{c}{\varphi }_{\xi }=arccos\sqrt{g_{\xi }(j_z,ϵ)}=\\ arccos\left\{{\left[\frac{2}{1-j_z^2}\left[\frac{{\eta }_z}{2{\omega }_0}{j}_z^2+j_z-ϵ\right]-\left(f_{\xi -}-\frac{{\eta }_y}{{\omega }_0}\right)\right]}^{1/2}{\left[\left(f_{\xi +}-f_{\xi -}\right)-\left(\frac{{\eta }_x}{{\omega }_0}-\frac{{\eta }_y}{{\omega }_0}\right)\right]}^{-1/2}\right\}.\end{array}$$

(79)

such that $[-{\varphi }_{\xi },{\varphi }_{\xi }]$ or $[\pi -{\varphi }_{\xi },\pi +{\varphi }_{\xi }]$. In general, we can obtain limiting values for $j_z$ and $\epsilon$ where the condition is satisfied by taking into account the aforementioned limits of $cos\varphi$.

First, we consider the limits given by $cos{\varphi }_{\pm }=\pm 1$. It leads to a quadratic equation for $j_z$ which reads

$$\begin{array}{c}{j}_z^2\left\{\frac{{\eta }_z}{2{\omega }_0}+\frac{1}{2}\left(f_{\xi +}-\frac{{\eta }_x}{{\omega }_0}\right)\right\}+j_z-\left\{ϵ+\frac{{\eta }_z}{2{\omega }_0}+\frac{1}{2}\left(f_{\xi +}-\frac{{\eta }_z-{\eta }_x}{{\omega }_0}\right)\right\}=0,\end{array}$$

(80)

where we have inserted a zero by adding and substracting ${\eta }_z/2{\omega }_0$. We observe that the effect of the interactions in z direction is to shift the energy, while in the x direction, it is to shift the critical coupling. The resulting roots are

$$\begin{array}{c}{j}_{z\xi }^{(\pm )}\left(\epsilon \right)=-\frac{1}{f_{\xi x}}\left[1\mp \sqrt{2f_{\xi x}(ϵ-{ϵ}_{s\xi x})}\right]\end{array}$$

(81)

Second, we obtain the limits given by $cos{\varphi }_{1,2}=0$. Likewise, we obtain a quadratic equation which reads

$$\begin{array}{c}{j}_z^2\left[\frac{{\eta }_z}{2{\omega }_0}+\frac{1}{2}\left(f_{\xi -}-\frac{{\eta }_y}{{\omega }_0}\right)\right]+j_z-\left[ϵ+\frac{{\eta }_z}{2{\omega }_0}+\frac{1}{2}\left(f_{\xi -}-\frac{{\eta }_z-{\eta }_y}{{\omega }_0}\right)\right]=0.\end{array}$$

(82)

where we have also added and subtracted ${\eta }_z/2{\omega }_0$. We notice it is identical to Equation (80), but changing ${\eta }_x\to {\eta }_y$ and $f_{\xi -}\to f_{\xi +}$. Consequently, the solutions are given by

$$\begin{array}{c}{j}_{z\xi }^{(1,2)}\left(\epsilon \right)=-\frac{1}{f_{\xi y}}\left[1\mp \sqrt{2f_{\xi y}(ϵ-{ϵ}_{s\xi y})}\right].\end{array}$$

(83)

In the following, we will study the particular cases of the TC ($\xi =0$) and Dicke ($\xi =1$) to understand the effects of the interactions on the emergence of energy domains and critical energies. Finally, we will comment on the arbitrary $\xi$ case and the general typology of ESQPTs.

5.1. Energy Domains in the TC Limit

In the TC limit, the key functions determining the integral in Equation (74) are given by

$$\begin{array}{c}{\varphi }_0(j_z,\epsilon )=arccos\left\{{\left[\frac{2}{1-j_z^2}\left(\frac{{\eta }_z}{2{\omega }_0}{j}_z^2+j_z-\epsilon \right)-\left(f_{0+}-\frac{{\eta }_y}{{\omega }_0}\right)\right]}^{1/2}{\left(\frac{{\eta }_y}{{\omega }_0}-\frac{{\eta }_x}{{\omega }_0}\right)}^{-1/2}\right\},\end{array}$$

(84)

$$\begin{array}{c}{j}_{z0}^{(\pm )}\left(\epsilon \right)=-\frac{1}{f_{0x}}\left[1\mp \sqrt{2f_{0x}(ϵ-{ϵ}_{s0x})}\right],\end{array}$$

(85)

$$\begin{array}{c}{j}_{z0}^{(1,2)}\left(\epsilon \right)=-\frac{1}{f_{0y}}\left[1\mp \sqrt{2f_{0y}(ϵ-{ϵ}_{s0y})}\right].\end{array}$$

(86)

We observe the expressions for the critical coupling and the ground-state energy of the superradiant-x and y phases of the TC limit are recovered when one considers either $j_{z0}^{(\pm )}$ or $j_{z0}^{(1,2)}$, respectively. We note that $j_{z0}^{\left(1\right)}\le j_{z0}^{\left(2\right)}$. We must compare these values with those coming from the stationary points $j_z=\pm 1$, i.e., the ones that are present in all the phases and whose energy is given by ${\epsilon }_{\pm }=\pm 1+{\eta }_z/2{\omega }_0$. Subsequently, we recognize four different energy phases and three critical energies using the comparison between the values of $j_z$, the conditions in Equation (84), and what we learned in Section 3. These energy domains correspond to various behaviors of the function $g_{\xi }(j_z,ϵ)$. In turn, they determine the intervals of the variable $j_z$. Without loss of generality, let us assume that we select ${\eta }_x$ and ${\eta }_y$ such that ${ϵ}_{s0x}<{ϵ}_{s0y}$. Then, the energy phases are:

(a)

The upper interval where ${ϵ}_{+}<ϵ$. Here, the function $g_0(j_z,ϵ)$ is always less than one. The whole pseudospin sphere is available: $j_z\in [-1,1]$ and ${\varphi }_0\in [0,2\pi )$. Consequently, the available phase space volume (per j) saturates to its limiting value ${\nu }_0\left(ϵ\right)=2/\omega$.

(b)

The interval where ${ϵ}_{-}<ϵ<{ϵ}_{+}$. Here, $j_z$ takes values only in the interval $\left[-1,j_{z0}^{(+)}\right]$ with $\left|j_{z0}^{(+)}\right|\le 1$. This interval is always present for all values of parameters and corresponds to available phase space from the absolute minimum point at $j_z=-1$ and the absolute maximum at $j_z=+1$.

(c)

The interval that is only present in the superradiant-y phase, ${ϵ}_{0sy}<ϵ\le {ϵ}_{-}$.

(d)

The interval arising in presence of both the superradiant-x and y phases, ${ϵ}_{0sx}\le ϵ\le {ϵ}_{0sy}$. Here, the south pole of the pseudospin sphere ($j_z=-1$) is inaccessible and the variable $j_z$ is restricted to the interval $j_{z0}^{(-)}\le 1\le j_{z0}^{(+)}$. Considering that ${ϵ}_{0sx}<{ϵ}_{0sy}$ is the ground-state energy in the superradiant-x phase.

Clearly, we have three critical energies given by ${ϵ}_0^{\left(c1\right)}={ϵ}_{s0y}$, ${ϵ}_0^{\left(c2\right)}={ϵ}_{-}$, and ${ϵ}_0^{\left(c1\right)}={ϵ}_{+}$. All of them correspond to stationary points of the energy surface and to what we have already found in Section 3. The semiclassical approximation to the DoS in the TC model becomes

$$\begin{array}{c}\frac{\omega }{2}{\nu }_0\left(ϵ\right)=\left\{\begin{array}{cc}\frac{1}{\pi }{\int }_{j_{z0}^{(-)}}^{j_0^{(+)}}{\varphi }_0(j_z,ϵ)d{j}_z,& ϵ\in [{ϵ}_{0sx},{ϵ}_{0sy}]\mathrm{a}\mathrm{n}\mathrm{d}\gamma \in [{\gamma }_{0x}^c,{\gamma }_{0y}^c],\\ \frac{1}{\pi }\left[{\int }_{j_{z0}^{(-)}}^{j_{z0}^{\left(1\right)}}{\varphi }_0(j_z,ϵ)d{j}_z+{\int }_{j_{z0}^{\left(2\right)}}^{j_{z0}^{(+)}}{\varphi }_0(j_z,ϵ)d{j}_z\right]& ϵ\in [{ϵ}_{s0y},{ϵ}_{-}]\mathrm{a}\mathrm{n}\mathrm{d}\gamma \in [{\gamma }_{0y}^c,\infty ),\\ +\frac{1}{2}\left(j_{z0}^{\left(2\right)}-j_{z0}^{\left(1\right)}\right),& \\ \frac{1}{\pi }{\int }_{j_{z0}^{\left(1\right)}}^{j_{z0}^{+}}{\varphi }_0(j_z,ϵ)d{j}_z+\frac{1}{2}\left(j_{z0}^{\left(1\right)}+1\right),& ϵ\in [{ϵ}_{-},{ϵ}_{+}],\mathrm{a}\mathrm{n}\mathrm{d}\gamma \in [0,\infty ),\\ 1,& {ϵ}_{+}\le ϵ,\mathrm{a}\mathrm{n}\mathrm{d}\gamma \in [0,\infty ).\end{array}\right.\end{array}$$

(87)

The onset of these energy domains depends on the three intervals of $\gamma$ that we discussed in Section 3. The boundary between each energy domain signals the existence of an ESQPT, as the DoS has a critical change characterized by a singularity in its derivative, even though the DoS is continuous in the energy variable. The type of ESQPT is encoded in the first derivative of the DoS, $d{\nu }_0\left(\epsilon \right)/dϵ$, which in turn is in terms of the derivative

$$\begin{array}{c}\frac{\partial {\varphi }_0}{\partial ϵ}=\frac{1}{1-j_z^2}{\left\{\left[1-g_0(j_z,\epsilon )\right]g_0\left(j_z,\epsilon \right)\left(\frac{{\eta }_y-{\eta }_x}{{\omega }_0}\right)\right\}}^{-1/2}.\end{array}$$

(88)

It can be shown that those ESQPTs at ${ϵ}_{0sy}$ corresponds to a logarithmic-type discontinuity, and the ones at ${ϵ}_{-}$ and ${ϵ}_{+}$ are of the jump type. This is not the typical behavior of the standard TC model: we recover it only when ${\eta }_x={\eta }_y$. There, the symmetry leads to two jump-type singularities at ${ϵ}_{0sy}={ϵ}_{0sx}$ and ${ϵ}_{+}$ [61,66]. This is because ${ϵ}_{0sy}$ becomes the ground state, so only the ESQPTs corresponding to the fixed points at $j_{zs}=\pm 1$ remain. We will offer a unified explanation of this behavior later, when we discuss the arbitrary $\xi$ case. The volume of the available phase space for the TC model, encoded in the form of the semiclassical DoS, is shown for three different sets of interacting parameters, as a function of the energy, in the top row of Figure 4a–c, where we have chosen interacting parameters to highlight the different domains.

5.2. Energy Domains in the Dicke Limit

Following the same reasoning as in the previous section, we derive the expressions for $\xi =1$:

$$\begin{array}{c}{\varphi }_1(j_z,\epsilon )=arccos\left\{{\left[\frac{2}{1-j_z^2}\left[\frac{{\eta }_z}{2{\omega }_0}{j}_z^2+j_z-ϵ\right]+\frac{{\eta }_y}{{\omega }_0}\right]}^{1/2}{\left(f_{1+}-\left(\frac{{\eta }_x}{{\omega }_0}-\frac{{\eta }_y}{{\omega }_0}\right)\right]}^{-1/2}\right\},\end{array}$$

(89)

$$\begin{array}{c}{j}_{z1}^{(\pm )}\left(\epsilon \right)=-\frac{1}{f_{1x}}\left[1\mp \sqrt{2f_{1x}(ϵ-{ϵ}_{s1x})}\right],\end{array}$$

(90)

$$\begin{array}{c}{j}_{z1}^{(1,2)}\left(\epsilon \right)=-\frac{1}{f_{1y}}\left[1\mp \sqrt{2f_{1y}(ϵ-{ϵ}_{s1y})}\right].\end{array}$$

(91)

where $j_{z1}^{\left(1\right)}\le j_{z1}^{\left(2\right)}$. We recall that $f_{1y}$ is independent of the light–matter coupling. For the Dicke model, the range of the $j_z$ variable is given by the same expressions as in the TC model. Thus, we obtain the following intervals:

(a)

The interval ${ϵ}_{+}<ϵ$, where, as in the TC model the whole pseudo-spin sphere is available $j_z\in [-1,1]$, ${\varphi }_1\in [0,2\pi )$, and ${\nu }_1\left(ϵ\right)=2/\omega$.

(b)

The interval ${ϵ}_{-}<ϵ<{ϵ}_{+}$. Here, the $j_z$ variable takes values only in the interval $[-1,j_{z1}^{(+)}]$ and ${\varphi }_1$ is restricted. When $j_z\in [-1,ϵ]$, ${\varphi }_1$ takes values in the whole interval $[0,2\pi )$, but if $ϵ<j_z\le j_{z1}^{(+)}$, $0<{\varphi }_1<\pi$.

(c)

The interval ${ϵ}_{1sy}<ϵ\le {ϵ}_{-}$. It only appears in the deformed phase.

(d)

The lower interval ${ϵ}_{1sx}\le ϵ\le {ϵ}_{1sy}$. Here, the south pole of the Bloch sphere $(j_{zs}=-1)$ is inaccessible and the $j_z$ variable becomes restricted to the interval $j_z\in [j_{z1}^{(-)},j_{z1}^{(+)}]$.

The expression for the semiclassical DoS in the Dicke limit becomes

$$\begin{array}{c}\frac{\omega }{2}{\nu }_1\left(ϵ\right)=\left\{\begin{array}{cc}\frac{1}{\pi }{\int }_{j_{z1}^{(-)}}^{j_{z1}^{(+)}}{\varphi }_1(j_z,ϵ)d{j}_z,& ϵ\in [{ϵ}_{s1x},{ϵ}_{s1y}]\mathrm{a}\mathrm{n}\mathrm{d}\gamma \in [{\gamma }_{1x}^c,\infty ),\\ \frac{1}{\pi }\left[{\int }_{j_{z1}^{(-)}}^{j_{z1}^{\left(1\right)}}{\varphi }_1(j_z,ϵ)d{j}_z+{\int }_{j_{z1}^{\left(2\right)}}^{j_{z1}^{(+)}}{\varphi }_1(j_z,ϵ)d{j}_z\right]& ϵ\in [{ϵ}_{s1y},{ϵ}_{-}],\gamma \in [{\gamma }_{1x}^c,\infty ],\\ +\frac{1}{2}\left(j_{z1}^{\left(2\right)}-j_{z1}^{\left(1\right)}\right),& \mathrm{a}\mathrm{n}\mathrm{d}\Delta {\eta }_{zy}\ge {\omega }_0,\\ \frac{1}{\pi }{\int }_{j_{z1}^{\left(1\right)}}^{j_{z1}^{(+)}}{\varphi }_1(j_z,ϵ)d{j}_z+\frac{1}{2}\left(j_{z1}^{\left(1\right)}+1\right),& ϵ\in [{ϵ}_{-},{ϵ}_{+}]\mathrm{a}\mathrm{n}\mathrm{d}\gamma \in [0,\infty ],\\ 1,& {ϵ}_{+}\le ϵ\mathrm{a}\mathrm{n}\mathrm{d}\gamma \in [0,\infty ].\end{array}\right.\end{array}$$

(92)

If we cancel the interactions in z and y directions (${\eta }_x={\eta }_y$), we recover Equation (19) in [55], where ${ϵ}_{s1y}=-{\omega }_0/2{\eta }_z$ and

$$\begin{array}{c}{\varphi }_1(j_z,\epsilon )=arccos\left\{{\left[\frac{2}{1-j_z^2}\left(\frac{{\eta }_z}{2{\omega }_0}{j}_z^2+j_z-ϵ\right)\right]}^{1/2}{f}_{1+}^{-1/2}\right\}.\end{array}$$

(93)

The volume of the available phase space for Dicke model for three different couplings, as a function of the energy, is shown in the middle row of Figure 4d–f for representative values of the parameters. The singular behavior of the DoS is encoded in the derivative

$$\begin{array}{c}\frac{\partial {\varphi }_1}{\partial ϵ}=\frac{1}{1-j_z^2}{\left\{\left[1-g_1(j_z,\epsilon )\right]g_1(j_z,\epsilon )\left(f_{1+}-\frac{{\eta }_x-{\eta }_y}{{\omega }_0}\right)\right\}}^{-1/2}.\end{array}$$

(94)

Even though the fixed points belonging to the deformed phase do not appear as extrema in the superradiant-x phase, they still impact the energy domains, given that they are related to the interactions in the z directions via $j_{z1}^{(1,2)}$. Then, one can still find four different energy domains and three ESPQT. As it is shown in Figure 4 (middle row), the DoS curves for the Dicke and TC limits are very similar as a function of the energy $ϵ$ for small couplings. The behavior becomes analogous to the TC case. Although this result is similar to that of an extended Dicke model, it differs from the standard Dicke model, where the singularity at ${ϵ}_{-}$ is of the logarithmic type [61,66].

5.3. Energy Domains for Arbitrary Couplings and Typology of ESQPTs

Finally, we offer some considerations about the most general case. When an arbitrary value of $\xi \in (0,1)$ is chosen, we obtain a general expression by combining Equations (79), (81), and (83). It reads

$$\begin{array}{c}\frac{\omega }{2}{\nu }_{\xi }\left(ϵ\right)=\left\{\begin{array}{cc}\frac{1}{\pi }{\int }_{j_{z\xi }^{(-)}}^{j_{z\xi }^{(+)}}{\varphi }_{\xi }(j_z,ϵ)d{j}_z,& {ϵ}_{s\xi x}\le ϵ\le {ϵ}_{s\xi y},\mathrm{a}\mathrm{n}\mathrm{d}\gamma \in [{\gamma }_{\xi x}^c,{\gamma }_{\xi y}^c],\\ \frac{1}{\pi }\left[{\int }_{j_{z\xi }^{(-)}}^{j_z^{\left(1\right)}}{\varphi }_0(j_z,ϵ)d{j}_z+{\int }_{j_{z\xi }^{\left(2\right)}}^{j_{z\xi }^{(+)}}{\varphi }_0(j_z,ϵ)d{j}_z\right]& {ϵ}_{s\xi y}<ϵ\le {ϵ}_{-},\mathrm{a}\mathrm{n}\mathrm{d}\gamma \in [{\gamma }_{\xi y}^c,\infty ],\\ +\frac{1}{2}\left(j_{z\xi }^{\left(2\right)}-j_{z\xi }^{\left(1\right)}\right),& \\ \frac{1}{\pi }{\int }_{j_{z\xi }^{\left(1\right)}}^{j_{z\xi }^{(+)}}{\varphi }_0(j_z,ϵ)d{j}_z+\frac{1}{2}\left(j_{z\xi }^{\left(1\right)}+1\right),& {ϵ}_{-}<ϵ\le {ϵ}_{+},\mathrm{a}\mathrm{n}\mathrm{d}\gamma \in [0,\infty ),\\ 1,& {ϵ}_{+}<ϵ\mathrm{a}\mathrm{n}\mathrm{d}\gamma \in [0,\infty ).\end{array}\right.\end{array}$$

(95)

Given that ${ϵ}_{s\xi x}<{ϵ}_{s\xi y}$. This expression reunites the effects of the qubit–qubit interactions and the arbitrary light–matter coupling $\xi$. As discussed before, the phenomenology of this result is similar to what is found in absence of interactions but for arbitrary $\xi$, as shown in Refs. [6,62,69,70]. The main difference is the modification of the critical coupling of the superradiant-$(+)$ phase by ${\eta }_x$, the critical coupling of the superradiant-$(-)$ phase by ${\eta }_y$ and the direct shift of the energy and $j_z$ intervals of validity by ${\eta }_z$. Additionally, in this case, we have as a general expression

$$\begin{array}{c}\frac{\partial {\varphi }_{\xi }}{\partial ϵ}=\frac{1}{1-j_z^2}{\left\{\left[1-g_{\xi }(j_z,\epsilon )\right]g_{\xi }\left(j_z,\epsilon \right)\left[\left(f_{\xi +}-f_{\xi -}\right)-\frac{{\eta }_x-{\eta }_y}{{\omega }_0}\right]\right\}}^{-1/2}.\end{array}$$

(96)

ESQPTs can be classified using two numbers: The index of the transition r, which denotes the number of negative eigenvalues of the Hessian matrix of the classical Hamiltonian at the stationary points (where the phase space volume changes), and the number of relevant degrees of freedom f of the system, determining in which derivative of the smooth DoS the discontinuity associated to the ESQPT appears [68]. This classification is tied to the properties of the Hessian matrix describing the local dependence of the Hamiltonian around the fixed points of the energy surface, so it is valid only if the Hessian does not have zero or singular eigenvalues. r determines the type of singularity: for even f systems, $r=1$ a logarithmic-type singularity, while $r=2$ indicates a jump-type one [6,68,70]. As the Dicke model has only two degrees of freedom (the collective spin and the boson), it has $f=2$. The integrability of the standard TC model reduces to $f=1$ instead [70]. For arbitrary $\xi$, it has been shown that there are three ESPQT marked by three critical energies given by ${ϵ}_{\xi }^{c1}={ϵ}_{s\xi +}$, ${ϵ}_{\xi }^{c2}={ϵ}_{-}$, and ${ϵ}_{\xi }^{c3}={ϵ}_{+}$. Their indices correspond $r=1$, $r=2$, and $r=2$, respectively, corresponding to saddle points ($r=1$) or maxima $(r=2)$, as discussed in Section 3.

A major result of our exploration is that the interactions in ${\eta }_x$ and ${\eta }_y$ play a similar role to an arbitrary value of $\xi$ in both the ground-state and excited-state properties. ${\eta }_x$ modifies the magnitude of ${\gamma }_{\xi +}^c$, and ${\eta }_y$ does the corresponding for ${\gamma }_{\xi -}^c$. As a result, the ESQPT coming from the transition between the superradiant-x and superradiant-y domains (and vice versa) are analogous to the transition between the superradiant-$(+)$ and $(-)$ phase at ${\eta }_i=0$. This ESQPT has an index $r=1$ [69]. Likewise, for an arbitrary value of matter–matter interactions and $\xi$ both ESQPTs, the one at the stationary point $j_{zs}=-1$ and the one due to the saturation of the Bloch sphere ($j_{zs}=+1$) have $r=2$ [61,69,70]. This explains our findings for the TC and Dicke limits. In the first case, the extended TC model, including interactions, is not integrable anymore. In the Dicke case, we could still have a finite ${\eta }_x$ modifying the DoS even though the critical coupling ${\gamma }_{1-}^c$ goes to infinity. Therefore, as it is revealed in Figure 4g–i, where we use $\xi =0.5$ and various values of the qubit interactions as representative examples, we will find a similar typology to the cases we have mentioned in the absence of qubit–qubit interactions. Finally, in Figure 4j, we illustrate the possible energy domains and the location of the critical energies as a function of $\gamma$ for $\xi =0.2$, ${\eta }_x={\eta }_y=1$ and ${\eta }_z=2$.

6. Discussion and Conclusions

In this work, we have investigated the quantum phases emerging in an extended Dicke model that involves qubit–qubit interactions. We have also included the possibility of varying the strength of the non-resonant terms so that the system can go from the Tavis–Cummings to the Dicke regimes. To this end, we used standard semiclassical techniques, whose central element is considering the expectation value of the quantum generalized Hamiltonian over a tensor product of Bloch and Glauber coherent states. By studying the shape of the energy surfaces, their stationary points, and the behavior of the semiclassical approximation to the Density of States, one can identify and characterize the QPT, ESQPTs, quantum phases, and energy domains resulting from the combination of light–matter (spin–boson) and matter–matter (collective spin) interactions.

We have found general expressions for the ground-state energy and analyzed the QPTs as a function of the Hamiltonian parameters in three cases: for the Tavis–Cummings limit ($\xi =0$), the Dicke limit ($\xi =1$), and for an arbitrary interaction strength in between these two. We have considered a general combination of collective qubit interactions represented by operators ${\widehat{J}}_i^2$ with strengths ${\eta }_i$ and $i=x,y,z$. This is the most general case for two-body interactions between the collective degrees of freedom of the qubits. Each direction has a particular role in modifying the critical phenomena of the standard light–matter system for both the ground and excited states. To start, we examine the results for interactions in the z direction. As mentioned before, three main results have been discovered before due to a finite ${\eta }_z$: shifting of the ground-state by ${\eta }_z/2$, the onset of first-order phase transitions, and the modification of the critical value of the light–matter interaction where the superradiant QPT appears (see, e.g., Ref. [55]). Indeed, we have confirmed these results and generalized them as we have found that the same phenomena occur in the presence of interactions in the x and y directions. Furthermore, we have noted that the relevant parameters of the system are the differences $\Delta {\eta }_{z,x}$ and $\Delta {\eta }_{z,y}$, which is a natural result due to the conservation of the pseudospin length. Tuning these quantities allows for the stimulation and suppression of superradiance via manipulating the light–matter interaction.

However, this is not the only effect of the x and y interactions. In terms of the $\Delta {\eta }_{zi}$ ($i=x,y$) parameters, they produce two new quantum phases, the superradiant-x and superradiant-y phases. If we assume the interactions in the x and y directions are balanced, we recover the distinctive rotational symmetry of the standard TC model. Regardless of the value of $\xi$, the normal phase would be symmetric, and, in the case of the TC limit, the superradiant phase will correspond to the well-known Mexican hat potential. Thus, we call it the superradiant-symmetric phase. In the imbalanced case, we observe new effects. The integrability of the TC Hamiltonian breaks down, and the two superradiant phases appear. Moreover, the energy surface of the normal phase is deformed for every $\xi$ stretching in the x or y direction depending on the relationship between $\Delta {\eta }_{zx}$ and $\Delta {\eta }_{zy}$. Additionally, new effects appear, such as first-order QPT between the x, y, and normal phases and the existence of parameter domains where the fixed points of both the x and y phases coexist. We refer to this situation as a superposition of phases, where one of them can be dominant [62]. The passage between a single superradiant phase to one in a superposition does not imply a QPT because the ground state remains the same. Still, it will affect the energy domains and ESQPT present for that specific parameter set. On the other hand, the Dicke limit becomes a situation where the superradiant-y phase vanishes and leaves a deformed or subradiant phase first identified in Ref. [55]. It only occurs for $\Delta {\eta }_{zy}\ge 1$ independently of $\gamma$. The onset of this phase produces the development of a new first-order QPT between it, the superradiant-x, and normal phases. Additionally, it suppresses any superposition between phases.

Notoriously, one can understand these results from a unified point of view by looking at the arbitrary $\xi$ case in general. In the absence of interactions, an intermediate value of the light–matter interaction leads to the existence of two phases, the superradiant-$(+)$ and $(-)$. Their position in the quantum phases landscape is fixed, depending on the relationship between the critical couplings ${\gamma }_{\xi \pm }^c$. As a result, for a light–matter interaction larger than ${\gamma }_{\xi -}^c$, the two phases are superimposed [6,62,69,70]. It turns out the superradiant-x (y) phase is a generalization of the superradiant-$(+)$ [$(-)$] phase. Therefore, the phenomenology of critical phenomena for both the ground and excited states is similar. This has been confirmed by analyzing the semiclassical Density of States in the three regimes of the light–matter coupling and for the various cases of qubit interaction strengths. We have obtained general expressions for the DoS and the limiting values of the atomic variables $(j_z,\varphi )$ in the Bloch sphere that allow to identify energy domains and critical energies tied to ESQPTs. Finally, we have unveiled a unique feature due to the qubit interactions. Unlike the superradiant-$(\pm )$, the landscape of the superradiant x and y phases can be modified at will by independently tuning the qubit interaction strengths. This specific feature is left to be studied in the near future.

Our study provides a broad perspective of critical phenomena in collective models combining strong light–matter and matter–matter interactions. Future directions, such as the exploration of the existence and robustness of Goldstone and Higgs modes in quantum optical setups [89], may benefit from the general description of the quantum phases our results provide. Moreover, as experimental progress promises to make individually controlled interactions in each direction feasible soon, we expect our work to be a reference for exploring critical quantum phenomena in quantum information, atomic physics, quantum optics, and condensed matter systems involving collective qubits interactions.