Phase and Amplitude Modes in the Anisotropic Dicke Model with Matter Interactions

Abstract

Phase and amplitude modes, also called polariton modes, are emergent phenomena that manifest across diverse physical systems, from condensed matter and particle physics to quantum optics. We study their behavior in an anisotropic Dicke model that includes collective matter interactions. We study the low-lying spectrum in the thermodynamic limit via the Holstein–Primakoff transformation and contrast the results with the semi-classical energy surface obtained via coherent states. We also explore the geometric phase for both boson and spin contours in the parameter space as a function of the phases in the system. We unveil novel phenomena due to the unique critical features provided by the interplay between the anisotropy and matter interactions. We expect our results to serve the observation of phase and amplitude modes in current quantum information platforms.

1. Introduction

Amplitude and phase modes are collective excitations of an order parameter usually studied in connection to a continuous symmetry breaking. For the $U\left(1\right)$ symmetry, they are known as massive Anderson–Higgs (AH) [1,2] and massless Nambu–Goldstone (NG) modes [3,4], respectively, and are present in several systems, from condensed matter [5,6,7,8,9], quantum antiferromagnets [10,11,12] to particle physics [13,14,15]. Phase and amplitude modes help to describe critical phenomena such as quantum phase transitions (QPTs) [3], i.e., a parametric sudden change in the ground state’s properties of a quantum system [16,17]. The paradigmatic case is the quartic potential $V\left(\alpha \right)=r{\alpha }^2+g{\alpha }^4$, expressed in terms of the scalar complex parameter $\alpha$ that possesses the $U\left(1\right)$ continuous symmetry [16]. The parameter space is separated into two phases: the normal or non-ordered phase ($r,g>0$), where the potential has the form of a single well with a minimum at $\alpha =0$, and the ordered phase ($r<0$), where the potential has a Mexican hat form. The ground state is a degeneration of states with $\left|\alpha \right|>0$. In the ordered phase, fluctuations of the parameter quantify the two possible excitations: (1) the phase mode (NG), which occurs at a constant amplitude and has vanishing excitation energy, thus being massless; and (2) the amplitude mode (AH), which has a constant phase and possesses finite excitation energy, thus being massive. In general, both modes are hard to observe experimentally. The NG is particularly sensitive to symmetry deviations, as it gains finite excitation energy [18,19]. Likewise, the AH mode is challenging to observe because of its massive character. Ultracold atoms setups are particularly suitable for their observation and control due to their high tunability [19,20,21,22,23,24].

A standard example of a mean-field quantum phase transition (QPT) is the superradiant-to-normal transition of the Dicke Hamiltonian [25,26,27]. This model represents a collection of atoms, each simplified to a two-energy level transition, interacting with a single-mode radiation field within a cavity [28]. When the light–matter coupling strength exceeds a critical threshold, the QPT is indicated by a non-zero photon number expectation value in the thermodynamic limit. The Dicke Hamiltonian provides a comprehensive framework for describing spin–boson interactions by capturing the collective behavior of a set of two level systems or qubits in both equilibrium and non-equilibrium setups [29,30,31,32]. Its algebraic simplicity has allowed the investigation of fundamental topics such as quantum chaos [33,34,35,36], the quantum–classical correspondence [37,38,39,40], Excited-State Quantum Phase Transition (ESQPT) [41,42,43,44], non-equilibrium phenomena [30,45,46,47], and the physics of the ultrastrong light–matter coupling regime [48,49,50]. Thanks to this, the model has become a general framework to describe collective qubit effects in a vast array of systems, including ultracold atoms [51,52,53], superconducting qubits [54,55,56,57], nuclear physics [58], solid-state systems [59,60,61], and quantum dots [62]. Also, superradiance has been realized in superconducting qubits [63,64,65], cavity-assisted Raman transitions [66,67], bose-Einstein condensates in optical lattices [68,69,70,71,72], and unitary Fermi gases [73,74]. It is precisely in these setups where the amplitude and phase modes have been studied in connection to supersolid states exhibiting self-organization [22,75,76].

Under the rotating-wave approximation, the Dicke model becomes the Tavis–Cummings (TC) model [77], which is integrable in a classical sense, conserving the total number of excitations. This is reflected in the semi-classical limit, where a Mexican hat potential emerges in the superradiant phase [78]. In this symmetry-broken phase, many states cluster together to form a quasi-continuous band of energies, where the phase can take any value in $[0,2\pi )$. Hence, there is no cost of energy to displace the state or change the order parameter along the potential’s minima. Because of this, the TC model exhibits both a gapless phase mode and a gapped amplitude mode. Instead, the Dicke model has only a ${\mathbb{Z}}_2$ symmetry. As a result, it has a doubly degenerate spectrum in the superradiant phase corresponding to a double-well semi-classical potential [78]. The pair of (symmetry-breaking) ground states correspond to those at the bottom of a potential, reflecting two possible values for the phase (0 and $\pi$). Hence, the phase mode becomes gapped (Higgs-like) and corresponds to what has been called a roton-type mode [79,80,81] in analogy to the roton excitations in superfluid helium. The phase and amplitude modes in the Dicke model are called polariton modes [82], as well. A polariton is a hybrid quantum state emerging from strong light–matter interactions [83]. Because the phase and amplitude modes in the Dicke model represent the collective light–matter quantum state at low energy, they correspond to effective upper and lower polaritons [84,85]. The differences between the Dicke and TC models can be identified by studying the gap between the ground and first excited states. It has been studied for finite size [19,86], and in the thermodynamical limit [35,87]. Also, the phase and amplitude modes have been observed for large light–matter coupling [81,88]. A way to do this is by exploiting their singular behavior at the superradiant QPT [18]. The symmetry-breaking effect in the Dicke model has been measured experimentally [27,70], and employed for quantum sensing [89].

The so-called anisotropic, generalized, or unbalanced Dicke model [18,38,78,90,91,92] is a useful tool to explore the parametric change in the gap and understand the polariton modes as one passes from a continuous to a discrete symmetry, where the phase model can gain mass. There, one assumes a relative strength between the rotating (excitation conserving) and counter-rotating terms, allowing to tune the Hamiltonian between the TC and Dicke limits. Both the equilibrium [93,94,95,96] and non-equilibrium properties [72,97,98,99,100] of the anisotropic Dicke model have been studied, exhibiting a rich phase space and the presence of novel critical phenomena. Phase modes have been found in other extensions of the Dicke model, such as generalizations of the spin–boson coupling [101], Jahn–Teller–Dicke models [102], two-mode Dicke models [88,103,104,105], Jaynes–Cummings lattices [106], the three-level Dicke model [107,108], and also including matter–matter interactions [19].

In light of experimental progress in the combination of strong matter–matter interactions and ultrastrong light–matter systems, such as the recent realization of the Dicke model in a solid-state system with tunable spin–magnon interactions [109,110,111], in this work, we investigate the effects of collective matter interactions over the phase and amplitude or polariton modes in the anisotropic Dicke model. Collective matter interactions mediated by light introduce unique behavior in the Dicke model, such as first-order QPT [112,113,114,115,116,117,118], parametric shifts over the superradiant QPT [54,119,120], new quantum phases [117,121,122,123], and novel perspectives on quantum chaos [117,124], such as the amplification of the regularity-to-chaos transition [125]. Also, recently, the matter and light linear response has been investigated for this model [126,127] in the context of the interplay between cavity QED and magnetic systems.

To achieve our end, we employ the standard approach of applying the Holstein–Primakoff (HP) transformation to obtain a low-energy approximation of the excitation modes around the minima [79]. Given the limitations of the HP transformation [128], we interpret the results under the light of the semi-classical corresponding energy surfaces that can be obtained via coherent states [122]. Additionally, we calculate the geometric or Berry phase around the critical points and relate its behavior to the phase and amplitude modes [129,130]. The geometric phase possesses non-analytic behavior along a QPT, so it has been employed as a signature of criticality [131,132,133]. It has been observed experimentally in Heisenberg chains [134] and qubit systems [135,136]. It has been studied for the standard Dicke model [137,138,139] and extensions such as including two-photon processes [140], impurity coupling [141], and collective matter interactions in the z direction [142].

The article is organized as follows. In Section 2, we present the anisotropic Dicke Hamiltonian with matter interactions and details over its classical limit. Next, in Section 3, we obtain exact expressions for the low-lying polariton branches using the Holstein–Primakoff approximation. In Section 4, we discuss the gap behavior as a function of the different phases in the system. We calculate the geometrical phase for arbitrary circulations generated by the photon number or the relative population operator in Section 5. Finally, in Section 6, we present our conclusions and perspectives.

2. Anisotropic Dicke Hamiltonian with Matter Interactions

The anisotropic Dicke Hamiltonian including collective matter interactions is given by

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

(1)

The first two terms denote the non-interacting parts of the Hamiltonian, where $\omega$ is the boson frequency, ${\omega }_0$ is the qubit energy splitting, ${\widehat{a}}^{†}$ ($\widehat{a}$) is the creation (annihilation) boson operator, and ${\widehat{J}}_z$ is the relative population of excited qubits, given that ${\widehat{J}}_{z,x,y}$ are pseudospin operators representing the collective degrees of freedom of the set of N qubits, which follow the rules of the su(2)-algebra. The third term is the spin–boson interaction, with $\gamma$ as the coupling strength and $\xi$ as the anisotropy parameter, that allows tuning the relative strength of the counter-rotating term, so $\xi =0$ and $\xi =1$ are the TC and Dicke limits, respectively. Being a system made of the tensor product of two Hilbert spaces, the Hamiltonian follows both the algebraic properties of the su(2) and the Heisenberg–Weyl algebras. Like the standard Dicke model [78], it commutes with the parity operator $\widehat{\Pi }=exp\left(i\pi {\widehat{n}}_e\right)$, where $n_e={\widehat{a}}^{†}\widehat{a}+{\widehat{J}}_z+j\widehat{\mathbb{I}}$ is the total number of the excitations operator. Because this is a discrete parity, the Hamiltonian is non-integrable in the classical sense, except for the TC limit, when the Hamiltonian commutes directly with ${\widehat{n}}_e$. Moreover, the Hamiltonian also commutes with the pseudospin length ${\widehat{\mathbf{J}}}^2=J_x^2+J_y^2+J_z^2$, so the totally symmetric subspace, where $j=N/2$, contains the ground-state of the collective system. The last term in Equation (1) accounts for the matter interactions, where ${\eta }_i$ with $i=x,y,z$ are the couplings in each direction, whose meaning depends on the particular setup realizing the Dicke model, including Josephson dynamics [115,143,144], atomic dipolar couplings [119,145], optomechanical setups [146,147], or interactions between superconducting qubits [54,56,57,148,149]. A standard approach is regarding ${\eta }_z$ as the strength of collective on-site interaction and ${\eta }_x$ (${\eta }_y$) as a strength of the collective inter-qubit interactions. It is worth emphasizing that the relevant interacting parameters are $\Delta {\eta }_{z\mu }={\eta }_z-{\eta }_{\mu }$, with $\mu =x,y$, given that the pseudospin length is conserved, and one can express one direction in terms of the others ${\widehat{J}}_z^2=j(j+1)\widehat{\mathbb{I}}-{\widehat{J}}_x^2-{\widehat{J}}_y^2$ [122].

The Hamiltonian in Equation (1) commutes with the squared pseudospin length operator $[{\widehat{H}}_{\xi },\widehat{{\mathbf{J}}^2}]$, so the collective ground state lies in the totally symmetric subspace that corresponds to $j=N/2$, and the Hilbert space of the system is effectively reduced to $N+1$ states. As a result, a classical corresponding Hamiltonian can be obtained by mapping the boson and collective pseudospin degrees of freedom to the classical limit via coherent states. To do so, one takes the value of the quantum Hamiltonian over the tensor product of Glauber $|z\rangle$ and Bloch $|w\rangle$ coherent states as trial states [38,40,78,150,151], where $|0\rangle$ and $|j,-j\rangle$ are the boson and pseudospin vacuum states, respectively [152]. In canonical classical variables for the boson $(q,p)$, $z=\sqrt{j/2}\left(q+ip\right)$ and the pseudospin $(j_z,\varphi )$, $w=\sqrt{(1+j_z)/(1-j_z)}{e}^{-i\varphi }$, the classical corresponding Hamiltonian reads [122]

$$\begin{array}{c}{H}_{\xi }^{\left(cl\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}$$

(2)

A standard method to determine quantum phases in this class of systems with collective degrees of freedom is to analyze the extreme points of the corresponding classical energy surfaces under catastrophe theory [153,154], here, the one in Equation (2). Because the ground state is well-described by coherent states as trial states, dramatic changes in the energy surface’s minima point out the existence of ground-state QPT. This method is extended to other extreme points on the energy surface to identify excited-state QPT or ESQPT [78,93,117]. In the following, we will focus only on the absolute minima of the energy surfaces, as we are interested in polariton modes. We employ the classical energy surfaces to interpret the behavior of these modes across the parameter space. We picture the surfaces in a set of variables $(u,v)=arccos(-j_z)(cos\varphi ,sin\varphi )$ that, by using Hamilton equations to eliminate the bosonic quadratures q and p, allows one to visualize the fixed points in the collective pseudospin space alone [78,122]. In this picture, the energy surface reads as

$$\begin{array}{c}E(\xi ,u,v)={\omega }_0{sin}^2\sqrt{u^2+v^2}\frac{1}{2\left(u^2+v^2\right)}\left[u^2\left(\frac{{\eta }_x}{{\omega }_0}-f_{\xi +}\right)+v^2\left(\frac{{\eta }_y}{{\omega }_0}-f_{\xi -}\right)\right]-\\ {\omega }_{0}cos\sqrt{u^2+v^2}\left(1-\frac{{\eta }_z}{2{\omega }_0}cos\sqrt{u^2+v^2}\right).\end{array}$$

(3)

where $f_{\xi \pm }={\gamma }^2/{\gamma }_{\xi \pm }^c$ and ${\gamma }_{\xi \pm }^c=\sqrt{\omega {\omega }_0}/(1\pm \xi )$ are the critical couplings for the superradiant-$(\pm )$ phases as detailed below.

The anisotropy introduces new quantum phases to the standard Dicke model that matter–matter interactions then modify. For zero matter interactions (${\eta }_i=0$), the Dicke ($\xi =1$) and TC ($\xi =0$) models exhibit only two phases, the normal (non-ordered) and the superradiant (ordered), separated by the critical coupling ${\gamma }_{\xi +}^c$. The anisotropy modifies the energy landscape such that two different superimposed phases appear, the superradiant-$(+)$ and superradiant-$(-)$ phases, the latter with a critical coupling ${\gamma }_{\xi -}^c$ [44,93,94]. These phases are also called the electric and magnetic superradiance, with a different coupling combination for the counter and counter-rotating terms [18,92]. Each phase has its extreme points, but only one predominates as the ground state, depending on which of the two terms, the rotating or the counter-rotating, is stronger. For our case of study $\xi \in [0,1]$, the superradiant-x phase defines the ground state.

In the presence of interactions (${\eta }_i\ne 0$) but in the absence of light–matter coupling ($\gamma =0$), one obtains a Lipkin–Meshkov–Glick (LMG) Hamiltonian [155,156,157]. Hence, the anisotropic Dicke Hamiltonian inherits the algebraic properties and critical phenomena from the LMG [158,159,160,161,162,163,164,165,166]. There are three major modifications to the ground-state quantum phase diagram of the anisotropic Dicke model due to matter interactions [122]. First, interactions in the x and y direction produce a deformation, creating a deformed normal phase [92]. Second, the superradiant-$(+)$ (superradiant-$(-)$) phase transforms into the superradiant-x (superradiant-y), and the critical coupling is modified to ${\gamma }_{\xi x}^c={(1+\Delta {\eta }_{zx}/{\omega }_0)}^{1/2}{\gamma }_{\xi +}^c$ (${\gamma }_{\xi y}^c={(1+\Delta {\eta }_{zy}/{\omega }_0)}^{1/2}{\gamma }_{\xi -}^c$) [122]. Finally, while ${\eta }_z$ produces an overall energy shift, in the case of the Dicke limit ($\xi =1$), a new phase emerges, the deformed phase, which suppresses superradiance for $\Delta {\eta }_{zy}\ge 1$ and is the source of a first-order QPT for arbitrary coupling [117]. In the next section, we discuss the behavior of the polariton modes in each of these quantum phases.

3. Phase and Amplitude Modes in the Anisotropic Dicke Model

The Holstein–Primakoff transformation (HPT) [167] provides a quadratic approximation to the lower energy modes of Equation (1) in the thermodynamic limit ($N\to \infty$). After diagonalizing the quadratic, approximated Hamiltonian, the resulting branches correspond to the polariton modes or mean-field low-lying excitations of the photon–matter quantum superposition [18,80]. The low-lying spectrum can also be obtained by linearization of the classical equations of motion, i.e., calculating the small oscillations around the minima [168,169]. The HPT reads

$$\begin{array}{c}{\widehat{J}}_z={\widehat{b}}^{†}\widehat{b}-j,{\widehat{J}}_{+}=\sqrt{2j}{\widehat{b}}^{†}\sqrt{1-\frac{{\widehat{b}}^{†}\widehat{b}}{2j}},{\widehat{J}}_{-}=\sqrt{2j}\sqrt{1-\frac{{\widehat{b}}^{†}\widehat{b}}{2j}}\widehat{b}.\end{array}$$

(4)

where $\widehat{b}$ (${\widehat{b}}^{†}$) is an annhilation (creation) boson operator such that $[\widehat{b},{\widehat{b}}^{†}]=\mathbb{I}$. First, by substituting Equation (4) in Equation (1), one obtains

$$\begin{array}{c}\widehat{H}=\omega {\widehat{a}}^{†}\widehat{a}+{\omega }_0\left({\widehat{b}}^{†}\widehat{b}-j\right)\\ +\gamma \left[\left(\widehat{a}{\widehat{b}}^{†}\sqrt{1-\frac{{\widehat{b}}^{†}\widehat{b}}{2j}}+{\widehat{a}}^{†}\sqrt{1-\frac{{\widehat{b}}^{†}\widehat{b}}{2j}}\widehat{b}\right)+\xi \left({\widehat{a}}^{†}{\widehat{b}}^{†}\sqrt{1-\frac{{\widehat{b}}^{†}\widehat{b}}{2j}}+\widehat{a}\sqrt{1-\frac{{\widehat{b}}^{†}\widehat{b}}{2j}}\widehat{b}\right)\right]\\ \frac{1}{4}\left[{\eta }_x{\left({\widehat{b}}^{†}\sqrt{1-\frac{{\widehat{b}}^{†}\widehat{b}}{2j}}+\sqrt{1-\frac{{\widehat{b}}^{†}\widehat{b}}{2j}}\widehat{b}\right)}^2-{\eta }_y{\left({\widehat{b}}^{†}\sqrt{1-\frac{{\widehat{b}}^{†}\widehat{b}}{2j}}-\sqrt{1-\frac{{\widehat{b}}^{†}\widehat{b}}{2j}}\widehat{b}\right)}^2\right]+\frac{{\eta }_z}{2j}{\left({\widehat{b}}^{†}\widehat{b}-j\right)}^2\end{array}$$

(5)

Next, to study the excitation modes, we employ a mean field approximation by displacing both $\widehat{a}$ and $\widehat{b}$ bosonic operators

$$\begin{array}{c}{\widehat{a}}^{†}\to {\widehat{c}}^{†}+\alpha \sqrt{2j},\mathrm{and}{\widehat{b}}^{†}\to {\widehat{d}}^{†}-\beta \sqrt{2j}.\end{array}$$

(6)

where $\alpha ,\beta \in \mathbb{C}$ are scaled displacements. The square roots in the HPT become

$$\begin{array}{c}\sqrt{1-\frac{{\widehat{b}}^{†}\widehat{b}}{2j}}=\sqrt{1-\frac{{\widehat{d}}^{†}\widehat{d}-\beta \sqrt{2j}({\widehat{d}}^{†}+\widehat{d})+2j{\beta }^2}{2j}}=\sqrt{k}\sqrt{{\xi }_0},\\ k=1-{\beta }^2,\sqrt{{\xi }_0}=\sqrt{1-\frac{{\widehat{d}}^{†}\widehat{d}-\beta \sqrt{2j}({\widehat{d}}^{†}+\widehat{d})+2j{\beta }^2}{2jk}}\end{array}$$

(7)

Next, we make use of the thermodynamic limit $N\to \infty$. Here, the customary assumption is that $\langle {\widehat{b}}^{†}\widehat{b}\rangle \ll 2j$, so $\sqrt{1-{\widehat{b}}^{†}\widehat{b}/2j}\simeq 1-{\widehat{b}}^{†}\widehat{b}$. In particular,

$$\begin{array}{c}\sqrt{{\xi }_0}\simeq 1-\frac{{\widehat{d}}^{†}\widehat{d}-\beta \sqrt{2j}\left({\widehat{d}}^{†}\widehat{d}\right)}{4jk}-\frac{{\beta }^2{({\widehat{d}}^{†}+\widehat{d})}^2}{16j{k}^2}\end{array}$$

(8)

By keeping only terms up to the quadratic order, the Hamiltonian becomes

$$\begin{array}{c}\widehat{H}=\omega {\widehat{c}}^{†}\widehat{c}+\left[{\omega }_0+\gamma {\displaystyle \frac{\alpha \beta }{\sqrt{k}}}(1+\xi )+{\eta }_z(2{\beta }^2-1)-{\beta }^2{\eta }_x\right]{\widehat{d}}^{†}\widehat{d}\\ +\left[\omega \alpha -\gamma \sqrt{k}\beta (1+\xi )\right]\sqrt{2j}({\widehat{c}}^{†}+\widehat{c})\\ +\left[-{\omega }_0\beta +\gamma \alpha \left({\displaystyle \frac{k-{\beta }^2}{\sqrt{k}}}\right)(1+\xi )-{\eta }_{x}k\beta \left(1-{\displaystyle \frac{{\beta }^2}{k}}\right)-{\eta }_z(2{\beta }^2-1)\beta \right]\sqrt{2j}({\widehat{d}}^{†}+\widehat{d})\\ +\left[\gamma {\displaystyle \frac{\alpha \beta }{4\sqrt{k}k}}(2-{\beta }^2)(1+\xi )+{\eta }_x{\displaystyle \frac{k}{4}}\left(1-{\displaystyle \frac{4{\beta }^2}{k}}\right)+{\eta }_z{\beta }^2\right]{({\widehat{d}}^{†}+\widehat{d})}^2-{\eta }_y{\displaystyle \frac{k}{4}}{({\widehat{d}}^{†}-\widehat{d})}^2\\ -\gamma \left[{\displaystyle \frac{{\beta }^2}{2\sqrt{k}}}(1+\xi )\right]({\widehat{c}}^{†}+\widehat{c})({\widehat{d}}^{†}+\widehat{d})+\gamma \sqrt{{\xi }_0}\left[(\widehat{c}{d}^{†}+{\widehat{c}}^{†}\widehat{d})+\xi ({\widehat{c}}^{†}{\widehat{d}}^{†}+\widehat{c}\widehat{d})\right]\\ +\left[\omega {\alpha }^2+{\omega }_0{\beta }^2-{\displaystyle \frac{{\omega }_0}{2}}-\gamma \sqrt{k}2\alpha \beta (1+\xi )+{\eta }_{x}k{\beta }^2+{\eta }_z{\left({\beta }^2-{\displaystyle \frac{1}{2}}\right)}^2\right]2j\\ -\gamma {\displaystyle \frac{\alpha \beta }{2\sqrt{k}}}(1+\xi )+{\eta }_x{\displaystyle \frac{{\beta }^2}{2}},\end{array}$$

(9)

3.1. Deformed Normal Phases

The deformed normal phases are characterized by a zero expectation value of photon and matter excitations. Thus, by setting $\alpha ,\beta \to 0$, we obtain

$$\begin{array}{c}\widehat{H}=\omega {\widehat{c}}^{†}\widehat{c}+{\omega }_0\left({\widehat{d}}^{†}\widehat{d}-j\right)+\gamma \left[\left(\widehat{c}{\widehat{d}}^{†}+{\widehat{c}}^{†}\widehat{d}\right)+\xi \left({\widehat{c}}^{†}{\widehat{d}}^{†}+\widehat{c}\widehat{d}\right)\right]+\\ \frac{1}{4}\left[{\eta }_x{\left({\widehat{d}}^{†}+\widehat{d}\right)}^2-{\eta }_y{\left({\widehat{d}}^{†}-\widehat{d}\right)}^2\right]+\frac{{\eta }_z}{2j}{\left({\widehat{d}}^{†}\widehat{d}-j\right)}^2.\end{array}$$

(10)

Next, one writes the Hamiltonian in terms of the boson quadratures [79]

$$\begin{array}{c}\widehat{c}=\sqrt{\frac{\omega }{2}}\left(\widehat{x}+\frac{i}{\omega }{\widehat{p}}_x\right),\widehat{d}=\sqrt{\frac{{\omega }_0}{2}}\left(\widehat{y}+\frac{i}{{\omega }_0}{\widehat{p}}_y\right),\end{array}$$

(11)

By substituting Equation (11), the Hamiltonian is now represented as follows:

$$\begin{array}{c}\widehat{H}=\frac{1}{2}\left[{\omega }^2{\widehat{x}}^2+{\omega }_{zx}^2{\widehat{y}}^2+{\widehat{p}}_x^2+{\omega }_{zy}^2{\widehat{p}}_y^2+\frac{\gamma }{{\gamma }_{\xi +}}2\omega {\omega }_0\widehat{x}\widehat{y}+\frac{\gamma }{{\gamma }_{\xi -}}2{\widehat{p}}_x{\widehat{p}}_y-{ϵ}_0\right].\end{array}$$

(12)

where

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

(13)

We consider the following dimensionless frequencies containing the matter interaction strengths

$$\begin{array}{c}{\tilde{\omega }}_{zx}=\sqrt{1-\frac{\Delta {\eta }_{zx}}{{\omega }_0}},\mathrm{and}{\tilde{\omega }}_{zy}=\sqrt{1-\frac{\Delta {\eta }_{zy}}{{\omega }_0}}.\end{array}$$

(14)

Next, two rotations are performed to eliminate the cross terms $\widehat{x}\widehat{y}$ and ${\widehat{p}}_x{\widehat{p}}_y$. We will apply a rotation with angle ${\theta }_1$ for the variables $\widehat{x}$ and $\widehat{y}$, and another with angle ${\theta }_2$ for ${\widehat{p}}_x$ and ${\widehat{p}}_y$:

$$\begin{array}{c}\left(\begin{array}{c}\widehat{x}\\ \widehat{y}\end{array}\right)\left(\begin{array}{cc}cos{\theta }_1& sin{\theta }_1\\ -sin{\theta }_1& cos{\theta }_1\end{array}\right)\left(\begin{array}{c}{\widehat{q}}_1\\ {\widehat{q}}_2\end{array}\right),\left(\begin{array}{c}{\widehat{p}}_x\\ {\widehat{p}}_y\end{array}\right)\left(\begin{array}{cc}cos{\theta }_2& sin{\theta }_2\\ -sin{\theta }_2& cos{\theta }_2\end{array}\right)\left(\begin{array}{c}{\widehat{p}}_1\\ {\widehat{p}}_2\end{array}\right).\end{array}$$

(15)

The following conditions arise:

$$\begin{array}{c}tan2{\theta }_1=2f_{\xi x}^{1/2}\frac{\omega {\omega }_0{\tilde{\omega }}_{zx}}{{\omega }_0^2{\tilde{\omega }}_{zx}^2-{\omega }^2},\mathrm{and}tan2{\theta }_2=2f_{\xi x}^{1/2}\left(\frac{1-\xi }{1+\xi }\right)\frac{{\tilde{\omega }}_{zx}}{{\tilde{\omega }}_{zy}^4-1},\end{array}$$

(16)

where $f_{\xi x}={(\gamma /{\gamma }_{\xi x})}^2$, and the critical coupling indicating the onset of superradiance is

$$\begin{array}{c}{\gamma }_{\xi x}^c={\gamma }_{\xi +}^c{\tilde{\omega }}_{zx}.\end{array}$$

(17)

Substituting these conditions yields a two-mode decoupled Hamiltonian, determining the energies ${ϵ}_{1\pm }^N$ that represent the low-lying excitation modes:

$$\begin{array}{c}\widehat{H}=\frac{1}{2}\left[{\left({ϵ}_{1-}^{\mathrm{N}}\right)}^2{\widehat{q}}_1^2+{\left({ϵ}_{1+}^{\mathrm{N}}\right)}^2{\widehat{q}}_2^2+{\left({ϵ}_{2-}^{\mathrm{N}}\right)}^2{\widehat{p}}_1^2+{\left({ϵ}_{2+}^{\mathrm{N}}\right)}^2{\widehat{p}}_2^2-{ϵ}_0\right],\end{array}$$

(18)

with energies given by

$$\begin{array}{c}{ϵ}_{1\pm }^N=\sqrt{\frac{1}{2}\left[({\omega }^2+{\omega }_0^2{\tilde{\omega }}_{zx}^2)\pm \sqrt{{({\omega }^2-{\omega }_0^2{\tilde{\omega }}_{zx}^2)}^2+4{\omega }^2{\omega }_0^2{\tilde{\omega }}_{zx}^2{f}_{\xi x}}\right]},\end{array}$$

(19)

$$\begin{array}{c}{ϵ}_{2\pm }^N=\sqrt{\frac{1}{2}\left[(1+{\tilde{\omega }}_{zy}^2)\pm \sqrt{{(1-{\tilde{\omega }}_{zy}^2)}^2+4{\tilde{\omega }}_{zy}^2{g}_{\xi y}}\right]}\end{array}$$

(20)

where $g_{\xi y}={\left(\gamma /{\gamma }_{\xi y}\right)}^2$, with

$$\begin{array}{c}{\gamma }_{\xi y}={\gamma }_{\xi -}^c{\tilde{\omega }}_{zy}\end{array}$$

(21)

the critical coupling of the superradiant-y behavior. It is convenient to rewrite Equation (18) in a second quantized form in order to obtain a decoupled Hamiltonian with new bosonic excitations ${\widehat{a}}_i$ (${\widehat{a}}_i^{†}$) with $i=1,2$ following $[{\widehat{a}}_i,{\widehat{a}}_i^{†}]=\mathbb{I}$. They read

$$\begin{array}{c}{\widehat{q}}_1=\frac{({\widehat{a}}_1^{†}+{\widehat{a}}_1)}{\sqrt{2{\omega }_{{\mathrm{N}}_{-}}}},{\widehat{p}}_1=i\sqrt{\frac{{\omega }_{{\mathrm{N}}_{-}}}{2}}({\widehat{a}}_1^{†}-{\widehat{a}}_1),\end{array}$$

(22)

$$\begin{array}{c}{\widehat{q}}_2=\frac{({\widehat{a}}_2^{†}+{\widehat{a}}_2)}{\sqrt{2{\omega }_{{\mathrm{N}}_{+}}}},{\widehat{p}}_2=i\sqrt{\frac{{\omega }_{{\mathrm{N}}_{+}}}{2}}({\widehat{a}}_2^{†}-{\widehat{a}}_2),\end{array}$$

(23)

where ${\omega }_{{\mathrm{N}}_{-}}={ϵ}_{1-}^{\mathrm{N}}/{ϵ}_{2-}^{\mathrm{N}}$, and ${\omega }_{{\mathrm{N}}_{+}}={ϵ}_{1+}^{\mathrm{N}}/{ϵ}_{2+}^{\mathrm{N}}$. Thus, we obtain the quantum oscillation modes and the low-energy spectrum. It is essential to highlight that these expressions depend on the $\Delta {\eta }_{xz}$ and $\Delta {\eta }_{yz}$ parameters, allowing us to modulate the behavior of the modes both in the normal phase and in the superradiant phase as we will see later on. Finally, the Hamiltonian reads

$$\begin{array}{c}\widehat{H}={ϵ}_0^{\mathrm{N}}+{ϵ}_{-}^{\mathrm{N}}{\widehat{a}}_1^{†}{\widehat{a}}_1+{ϵ}_{+}^{\mathrm{N}}{\widehat{a}}_2^{†}{\widehat{a}}_2,\end{array}$$

(24)

with

$$\begin{array}{c}{ϵ}_0^{\mathrm{N}}=\frac{1}{2}\left({ϵ}_{-}^{\mathrm{N}}+{ϵ}_{+}^{\mathrm{N}}-{ϵ}_0\right)\end{array}$$

(25)

The phase ${ϵ}_{-}^{\mathrm{N}}$ and amplitude ${ϵ}_{+}^{\mathrm{N}}$ modes in the deformed normal phases become

$$\begin{array}{c}{ϵ}_{-}^{\mathrm{N}}={ϵ}_{2-}^{\mathrm{N}}{ϵ}_{1-}^{\mathrm{N}},{ϵ}_{+}^{\mathrm{N}}={ϵ}_{2+}^{\mathrm{N}}{ϵ}_{1+}^{\mathrm{N}}.\end{array}$$

(26)

3.2. Superradiant Phases

Now, we consider the case where $\alpha ,\beta \ne 0$. To obtain a quadratic Hamiltonian, we need to eliminate the linear terms in $\widehat{c},{\widehat{c}}^{†}$, and $\widehat{d},{\widehat{d}}^{†}$ from Equation (3). Hence, we obtain the following conditions determining the values of $\alpha$ and $\beta$,

$$\begin{array}{c}\omega \alpha -\gamma \sqrt{k}\beta (1+\xi )=0,\end{array}$$

(27)

$$\begin{array}{c}-{\omega }_0\beta +\gamma \alpha \left(\frac{k-{\beta }^2}{\sqrt{k}}\right)(1+\xi )-{\eta }_{x}k\beta \left(1-\frac{{\beta }^2}{k}\right)-{\eta }_z(2{\beta }^2-1)\beta =0\end{array}$$

(28)

Solving the system of equations leads to

$$\begin{array}{c}\alpha =\frac{\gamma (1+\xi )}{2\omega }\sqrt{1-{\mu }_x^2}\mathrm{and}\beta =\sqrt{\frac{1}{2}\left(1-{\mu }_x\right)},\end{array}$$

(29)

where

$$\begin{array}{c}{\mu }_x={\left[{\tilde{\omega }}_{zx}^2\left(f_{\xi x}-1\right)+1\right]}^{-1}.\end{array}$$

(30)

The critical coupling ${\gamma }_{\xi x}^c$ in Equation (17) for the existence of superradiant behavior is determined by ${\mu }_x=1$ because for ${\mu }_x>1$, $\alpha$ and $\beta$ in Equation (29) become complex. Writing the resulting Hamiltonian in terms of Equation (29), we have:

$$\begin{array}{c}\widehat{H}=\omega {\widehat{c}}^{†}\widehat{c}+{\omega }_A{\widehat{d}}^{†}\widehat{d}+{\omega }_B{({\widehat{d}}^{†}+\widehat{d})}^2+{\omega }_C{({\widehat{d}}^{†}-\widehat{d})}^2+{\omega }_D({\widehat{c}}^{†}+\widehat{c})({\widehat{d}}^{†}+\widehat{d})\\ +{\omega }_E\left[(c{d}^{†}+{\widehat{c}}^{†}\widehat{d})+\xi ({\widehat{c}}^{†}{\widehat{d}}^{†}+\widehat{c}\widehat{d})\right]+{\omega }_F\end{array}$$

(31)

where we have defined a set of reduced variables

$$\begin{array}{c}{\omega }_A={\displaystyle \frac{{\omega }_0}{2}}\left({\displaystyle \frac{1}{{\mu }_x}}-{\displaystyle \frac{{\eta }_z}{{\omega }_0}}\right)(1+{\mu }_x),\\ {\omega }_B={\displaystyle \frac{{\omega }_0}{2}}{\displaystyle \frac{1-{\mu }_x}{4}}\left[{\displaystyle \frac{1}{1+{\mu }_x}}{\displaystyle \frac{3+{\mu }_x}{{\mu }_x}}+{\displaystyle \frac{{\eta }_z}{{\omega }_0}}{\displaystyle \frac{1+3{\mu }_x}{1+{\mu }_x}}+{\displaystyle \frac{{\eta }_x}{{\omega }_0}}{\displaystyle \frac{4{\mu }_x^2}{1-{\mu }_x^2}}\right],\\ {\omega }_C=-\frac{{\omega }_0}{8}\frac{{\eta }_y}{{\omega }_0}(1+{\mu }_x),{\omega }_D=-\frac{\sqrt{2}}{4}{f}_{\xi x}^{1/2}\sqrt{\omega {\omega }_0}{\tilde{\omega }}_{zx}\frac{1-{\mu }_x}{\sqrt{1+{\mu }_x}},{\omega }_E=\gamma \sqrt{\frac{1}{2}(1+{\mu }_x)},\\ {\omega }_F=-{\displaystyle \frac{{\omega }_0}{2}}\left\{\left[{\displaystyle \frac{1-{\mu }_x^2}{2{\mu }_x}}+\left({\mu }_x-{\displaystyle \frac{{\eta }_z}{2{\omega }_0}}\right)\right]2j+{\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{1}{{\mu }_x}}-{\displaystyle \frac{{\eta }_z}{{\omega }_0}}\right](1-{\mu }_x)\right\}.\end{array}$$

Just as we did in the deformed normal phases, we will proceed to diagonalize the Hamiltonian in Equation (31) by transforming it into quadratures and then applying rotations that decoupled them. This transformation reads

$$\begin{array}{c}\widehat{c}=\sqrt{\frac{\omega }{2}}\left(\widehat{x}+\frac{i}{\omega }{\widehat{p}}_x\right),\widehat{d}=\sqrt{\frac{{\omega }_A}{2}}\left(\widehat{y}+\frac{i}{{\omega }_A}{\widehat{p}}_y\right),\end{array}$$

(32)

so the Hamiltonian is expressed as

$$\begin{array}{c}\widehat{H}=\frac{1}{2}\left({\omega }^2{\widehat{x}}^2+{\omega }_A^2{\chi }_{+}{\widehat{y}}^2+{\widehat{p}}_x^2+{\chi }_{-}{\widehat{p}}_y^2+\sqrt{\omega {\omega }_A}{\kappa }_{+}\widehat{x}\widehat{y}+\frac{{\kappa }_{-}}{\sqrt{\omega {\omega }_A}}{\widehat{p}}_x{\widehat{p}}_y-{ϵ}_1\right).\end{array}$$

(33)

where

$$\begin{array}{c}{ϵ}_1={\omega }_0\left[{\displaystyle \frac{1-{\mu }_x^2}{2{\mu }_x}}+\left({\mu }_x-{\displaystyle \frac{{\eta }_z}{2{\omega }_0}}\right)\right]2j+{\omega }_0\left({\displaystyle \frac{1}{{\mu }_x}}-{\displaystyle \frac{{\eta }_z}{{\omega }_0}}\right)+\omega \end{array}$$

(34)

Next, we consider the following variable changes

$$\begin{array}{c}{\chi }_{+}=\left(1+4\frac{{\omega }_B}{{\omega }_A}\right),{\kappa }_{+}=4{\omega }_D+\sqrt{2}{f}_{\xi x}^{1/2}\sqrt{\omega {\omega }_0}{\tilde{\omega }}_{zx}\sqrt{1+{\mu }_x},\\ {\chi }_{-}=\left(1-4\frac{{\omega }_C}{{\omega }_A}\right),{\kappa }_{-}=\sqrt{2}{g}_{\xi y}^{1/2}\sqrt{\omega {\omega }_0}{\tilde{\omega }}_{zy}\sqrt{1+{\mu }_x}.\end{array}$$

Again, we apply two rotations, ${\varphi }_1$ for x and y, and ${\varphi }_2$ for $p_x$ and $p_y$,

$$\begin{array}{c}\left(\begin{array}{c}\widehat{x}\\ \widehat{y}\end{array}\right)\left(\begin{array}{cc}cos{\varphi }_1& sin{\varphi }_1\\ -sin{\varphi }_1& cos{\varphi }_1\end{array}\right)\left(\begin{array}{c}{\widehat{q}}_1\\ {\widehat{q}}_2\end{array}\right),\left(\begin{array}{c}{\widehat{p}}_x\\ {\widehat{p}}_y\end{array}\right)\left(\begin{array}{cc}cos{\varphi }_2& sin{\varphi }_2\\ -sin{\varphi }_2& cos{\varphi }_2\end{array}\right)\left(\begin{array}{c}{\widehat{p}}_1\\ {\widehat{p}}_2\end{array}\right).\end{array}$$

(35)

We obtain the following conditions for the angles ${\varphi }_1$ and ${\varphi }_2$ from eliminating the crossed terms ${\widehat{q}}_1{\widehat{q}}_2$ and ${\widehat{p}}_1{\widehat{p}}_2$,

$$\begin{array}{c}tan2{\varphi }_1=\frac{\sqrt{\omega {\omega }_A}{\kappa }_{+}}{{\omega }_A^2{\chi }_{+}-{\omega }^2}ytan2{\varphi }_2=\frac{{\kappa }_{-}}{\sqrt{\omega {\omega }_A}({\chi }_{-}-1)}.\end{array}$$

(36)

The Hamiltonian reduces to:

$$\begin{array}{c}\widehat{H}=\frac{1}{2}\left[{\left({ϵ}_{1-}^{\mathrm{S}}{\widehat{q}}_1\right)}^2+{\left({ϵ}_{1+}^{\mathrm{S}}{\widehat{q}}_2\right)}^2+{\left({ϵ}_{2-}^{\mathrm{S}}{\widehat{p}}_1\right)}^2+{\left({ϵ}_{2+}^{\mathrm{S}}{\widehat{p}}_2\right)}^2-{ϵ}_1\right],\end{array}$$

(37)

with energies given by

$$\begin{array}{c}{ϵ}_{1\pm }^{\mathrm{S}}=\sqrt{\frac{1}{2}\left(({\omega }^2+{\omega }_A^2{\chi }_{+})\pm \sqrt{{({\omega }_A^2{\chi }_{+}-{\omega }^2)}^2+\omega {\omega }_A{\kappa }_{+}^2}\right)},\end{array}$$

(38)

$$\begin{array}{c}{ϵ}_{2\pm }^{\mathrm{S}}=\sqrt{\frac{1}{2}\left((1+{\chi }_{-})\pm \sqrt{{({\chi }_{-}-1)}^2+\frac{{\kappa }_{-}^2}{\omega {\omega }_A}}\right)}.\end{array}$$

(39)

We decouple the Hamiltonian using quadratures:

$$\begin{array}{c}{\widehat{q}}_1=\frac{({\widehat{a}}_1^{†}+{\widehat{a}}_1)}{\sqrt{2{\omega }_{{\mathrm{S}}_{-}}}},{\widehat{p}}_1=i\sqrt{\frac{{\omega }_{{\mathrm{S}}_{-}}}{2}}({\widehat{a}}_1^{†}-{\widehat{a}}_1),\\ {\widehat{q}}_2=\frac{({\widehat{a}}_2^{†}+{\widehat{a}}_2)}{\sqrt{2{\omega }_{{\mathrm{S}}_{+}}}},{\widehat{p}}_2=i\sqrt{\frac{{\omega }_{{\mathrm{S}}_{+}}}{2}}({\widehat{a}}_2^{†}-{\widehat{a}}_2),\end{array}$$

where, as in the deformed normal case, ${\omega }_{{\mathrm{S}}_{-}}={ϵ}_{1-}^{\mathrm{S}}/{ϵ}_{2-}^{\mathrm{S}}$, ${\omega }_{{\mathrm{S}}_{+}}={ϵ}_{1+}^{\mathrm{S}}/{ϵ}_{2+}^{\mathrm{S}}$. Finally, the Hamiltonian reads in second quantized form as

$$\begin{array}{c}\widehat{H}={ϵ}_0^{\mathrm{S}}+{ϵ}_{-}^{\mathrm{S}}{\widehat{a}}_1^{†}{\widehat{a}}_1+{ϵ}_{+}^{\mathrm{S}}{\widehat{a}}_2^{†}{\widehat{a}}_2\end{array}$$

(40)

with

$$\begin{array}{c}{ϵ}_0^{\mathrm{S}}=\frac{1}{2}\left({ϵ}_{-}^{\mathrm{S}}+{ϵ}_{+}^{\mathrm{S}}-{ϵ}_1\right),\end{array}$$

(41)

with the phase ${ϵ}_{-}^{\mathrm{S}}$ and amplitude ${ϵ}_{+}^{\mathrm{S}}$ modes being

$$\begin{array}{c}{ϵ}_{-}^{\mathrm{S}}={ϵ}_{2-}^{\mathrm{S}}{ϵ}_{1-}^{\mathrm{S}},{ϵ}_{+}^{\mathrm{S}}={ϵ}_{2+}^{\mathrm{S}}{ϵ}_{1+}^{\mathrm{S}}.\end{array}$$

(42)

3.3. Deformed Phase

Here, we comment on the deformed phase in the Dicke limit ($\xi =1$) that suppresses the superradiant-y phase [117,122]. It exists independently of the spin–boson coupling $\gamma$, while $\Delta {\eta }_{zy}\le {\omega }_0$, and it is characterized by the zero expectation value of the photon number but is not a normal phase. This can be recognized in the energy surface by a $\pi$ rotation of the extreme points. If we take $\xi =1$ in Equations (19) and (38), $g_{\xi y}\to 0$, we obtain

$$\begin{array}{c}{ϵ}_{2\pm }^N=\sqrt{1-(1\mp 1){\displaystyle \frac{\Delta {\eta }_{zy}}{2{\omega }_0}}},{ϵ}_{2\pm }^S=\sqrt{{\displaystyle \frac{\left(1\mp 1\right)\left(1-\frac{\Delta {\eta }_{zy}{\mu }_x}{{\omega }_0}\right)+\left(1\pm 1\right)\left(1-\frac{{\eta }_z{\mu }_x}{{\omega }_0}\right)}{2\left(1-\frac{{\eta }_z{\mu }_x}{{\omega }_0}\right)}}}.\end{array}$$

(43)

Because ${ϵ}_{\pm }^{\mathrm{S},\mathrm{N}}={ϵ}_{2\pm }^{\mathrm{S},\mathrm{N}}{ϵ}_{1\pm }^{\mathrm{S},\mathrm{N}}$, the phase mode becomes undefined when $\Delta {\eta }_{zy}\ge {\omega }_0$, or $\Delta {\eta }_{zy}{\mu }_x\ge {\omega }_0$ in the normal and superradiant phases, respectively, and the amplitude mode turns to ${ϵ}_{-}^{\mathrm{N},\mathrm{S}}={ϵ}_{1-}^{\mathrm{N},\mathrm{S}}$.

4. Role of Matter–Matter Interactions

Next, we explore the behavior of the phase and amplitude modes as a function of the matter–matter interaction parameters for the Dicke ($\xi =1$), TC ($\xi =0$) and an anisotropic case, particularly $\xi =0.5$ as an example.

4.1. Absence of Matter–Matter Interactions

As shown in Figure 1, first, we revisit the polariton modes without matter interactions. In Figure 1(a1), we show the polariton modes in the Tavis–Cummings limit. There, the lower curve corresponds to the phase mode ${ϵ}_{-}^N$, and the upper one to the amplitude mode ${ϵ}_{+}^N$. In the normal phase, the parametric evolution of both modes increases proportionally to the light–matter coupling. Their energy gap is just the Rabi splitting $2\gamma$. This agrees with the overall behavior of the energy surface, which adopts a spherical well shape, a reflection of the conserved $U\left(1\right)$ symmetry [see Figure 1(a2)]. The dashed, vertical line marks the critical coupling ${\gamma }_{0x}^{\left(c\right)}={\gamma }_{0+}^{\left(c\right)}$, separating the normal from the superradiant phases. There, the energy surface experiences a sudden change, taking the shape of a Mexican hat (see [Figure 1(a3,a4)]. Because of the onset of degeneration, we have spontaneous symmetry breaking. Once we enter the superradiant phase, the phase mode ${ϵ}_{-}^S$ becomes a Goldstone mode with zero energy, depicted on the energy surface as the minimum energy ring of the Mexican hat [22,23,88,106,170,171]. Instead, the amplitude mode ${ϵ}_{+}^S$ depends on the brim of the Mexican hat. As it grows, it requires increasingly higher energy costs than the phase modes, which remain constant. Note that in this case, ${\gamma }_{0x}={\gamma }_{0y}={\gamma }_{0-}$, so there is no presence of the superradiant-$(-)$ effects. In polariton terms, the two modes represent a maximal light–matter hybridization in resonance ($\omega ={\omega }_0$). The amplitude mode, being massive, can be identified as an upper-polariton branch that contains the most matter content, while the phase mode, being massless, becomes the lower-polariton branch with the most photonic content. Here, entering the superradiant phase implies a vanishment of the lower-polariton energy.

The energy spectra on Figure 1(b1,c1) correspond to the standard Dicke model $(\xi =1)$ and the anisotropic case when $\xi =0.5$, respectively. Both cases exhibit quite similar behavior. In the normal phase, the behavior resembles the TC limit because the system is in the strong-interacting regime. At the critical coupling ${\gamma }_{\xi x}^c={\gamma }_{\xi +}^c$, however, the phase mode tends to zero, and the amplitude mode to a fixed value given by

$$\begin{array}{c}{ϵ}_{+}^{\mathrm{N},c}=\sqrt{{\displaystyle \frac{1}{2}}\left({\omega }^2+{\omega }_0^2{\omega }_{zx}^2\right)\left[(1+{\tilde{\omega }}_{zy}^2)+\sqrt{{(1-{\tilde{\omega }}_{zy}^2)}^2+4{\tilde{\omega }}_{zx}^2{\displaystyle \frac{{(1-\xi )}^2}{{(1+\xi )}^2}}}\right]},\end{array}$$

(44)

in the normal phase, and by

$$\begin{array}{c}{ϵ}_{+}^{\mathrm{S},c}=\sqrt{{\displaystyle \frac{1}{2}}\left[{\omega }^2+\left(1-{\displaystyle \frac{{\eta }_z}{{\omega }_0}}\right){\omega }_0^2{\tilde{\omega }}_{zx}^2\right]\left[\left(1+{\displaystyle \frac{{\tilde{\omega }}_{zy}^2}{1-\frac{{\eta }_z}{{\omega }_0}}}\right)+\sqrt{{\left(1-{\displaystyle \frac{{\tilde{\omega }}_{zy}^2}{1-\frac{{\eta }_z}{{\omega }_0}}}\right)}^2+{\displaystyle \frac{4{\tilde{\omega }}_{zx}^2}{1-\frac{{\eta }_z}{{\omega }_0}}}{\displaystyle \frac{{(1-\xi )}^2}{{(1+\xi )}^2}}}\right]}\end{array}$$

(45)

in the superradiant phase. Discontinuities may arise for ${\eta }_z\ne 0$ at the critical coupling. In the absence of interactions it becomes

$$\begin{array}{c}{ϵ}_{+}^c=\sqrt{\left({\omega }^2+{\omega }_0^2\right)\left(1+\left|\frac{1-\xi }{1+\xi }\right|\right)},\end{array}$$

(46)

and for the standard Dicke model, it is ${ϵ}_{+}^c={({\omega }^2+{\omega }_0^2)}^{1/2}$ [80]. When one passes the QPT, the phase mode gains a finite energy, becoming the rotonic mode [23,88,170], which converges to the finite value

$$\begin{array}{c}{ϵ}_{-}^{\gamma \gg 1}=\omega \sqrt{1+\left|\frac{1-\xi }{1+\xi }\right|},\end{array}$$

(47)

for larger couplings. This value comes from the energy cost to pass between the two degenerate minima (whose phases are 0 and $\pi$) observed in Figure 1(b3,b4). The amplitude mode behaves quite similarly to the TC since the center of the energy surface presents an unstable point and tends quadratically as a function of $\gamma$, ${ϵ}_{+}^{\gamma \gg 1}={(1+\xi )}^2{\gamma }^2/\omega$ in the absence of matter interactions, which has the value $4{\gamma }^2/\omega$ in the standard Dicke model [80].

It has been shown that increasing the anisotropy decreases the amplitude mode gap until one obtains the Goldstone mode. This can be seen from Equation (47), where the limit $\xi \to 0$ leads to $\omega$ returning to the non-interacting value of photon field for the isotropic model [80]. This critical feature has been proposed as an experimental signature for detecting the Goldstone mode [18,22]. We observe this feature in Figure 1(c1), and it can be explained in terms of the change in the energy surface that leads to a new set of extremal points related to the superradiant-$(-)$ phase that reduces the energy of the phase mode [122]. In this case, the presence of the counter-rotating terms that do not conserve the excitation number is reflected in the finite energy of the phase mode as a lower-polariton mode, where light gains “mass” close to the QPT and then becomes a massless but highly correlated state. In polariton terms, the lower-polariton energy becomes gapped. We notice, however, that the criticality of these cases is not affected by the superradiant-$(-)$ extreme points because, in the absence of interactions, they cannot be shifted as it happens below. For the Dicke model ${\gamma }_{1-}^c\to \infty$, and for the anisotropic case with $\xi =0.5$, the critical coupling ${\gamma }_{0.5-}^c=3{\gamma }_{0.5+}^c$ is deep in the superradiant phase.

4.2. Presence of Matter–Matter Interactions

Next, we explore the effect of matter interactions over the polariton modes. In general, the z-interactions alone have the role of shifting both the critical coupling and the energy spectrum. To simplify the analysis, we will focus on the role of the x and y interactions via the relevant parameters $\Delta {\eta }_{zx}$ and $\Delta {\eta }_{zy}$, respectively.

In Figure 2, we show a representative value of the interactions given by ${\eta }_y=0.9{\omega }_0$ ($\Delta {\eta }_{zy}=-0.9{\omega }_0$), keeping ${\eta }_x=0$. Because of the interactions, in the normal phase, the energy surface widens along the u direction for all $\xi$, so we call it the deformed normal phase [122] [See Figure 2(a2–c2)]. In this case, the deformation does not substantially affect the phase and amplitude modes in the normal phase, except for an energy gap of the amplitude mode for zero spin–boson coupling given by

$$\begin{array}{c}{ϵ}_{+}^N={\omega }_0{\tilde{\omega }}_{zx}{\tilde{\omega }}_{zy}.\end{array}$$

(48)

and representing the cost to climb the well in the presence of interactions. Next, we cross the superradiant QPT. For the three cases under study ($\xi =0,0.5,1$), the amplitude mode is not significantly changed by the interactions compared to the case without. However, this is not the case for the phase mode in the TC limit ($\xi =0$). As the system enters the superradiant phase, the excitation number is not conserved anymore because of the deformation, so we do not have the Mexican hat potential anymore. As shown in Figure 2(a1), the phase mode gains energy, becoming rotonic. As we increase coupling, the onset of new extreme points at higher energies contributes to somehow restoring the original symmetry of the Mexican hat potential in the absence of interactions and reducing the mode’s energy. Thus, as the spin–boson coupling increases, the phase mode slowly tends to a Goldstone-like phase mode. As a result, interactions produce a “variable” mass roton mode. In both the anisotropic and Dicke limit cases, the results are similar to the case without interactions. After transitioning to the superradiant phase, two stable degenerate minimum points form, so there is no mechanism capable of decreasing the energy of the phase mode and eventually converging to a general energy value given by Equation (47).

The situation is different for x interactions as shown in Figure 3 for a representative value ${\eta }_x=0.9$ ($\Delta {\eta }_{zx}=-0.9$). We observe similar dynamics for the upper- and lower-polariton modes to the previous case. However, new phenomena emerge from the interplay between the anisotropy and matter interactions, i.e., the ability of the latter to shift the position of the superradiant-$(-)$ phenomena. Again, the normal phase becomes deformed for every $\xi$ but now with a different orientation, alongside the v direction corresponding to the change from y to x interactions [122] [See Figure 3(a2–c2)]. While the overall behavior of the amplitude mode remains the same, in the TC limit [Figure 3(a1)], we observe that the phase mode vanishes in the normal phase before the onset of the critical coupling. This is not a Goldstone mode but a suppression due to matter–matter interactions. As it can be seen from Equation (19), two conditions make the phase mode vanish, either $f_{\xi x}=1$ or $g_{\xi y}=1$. As it can be shown in Figure 3(a1), the relative position of the critical coupling ${\gamma }_{0x}^c$ with respect to ${\gamma }_{0y}^c$ is shifted thanks to ${\eta }_x\ne 0$, creating a coupling regime where the phase mode cannot exist, although the amplitude mode does. This can be seen from the relation between the two critical couplings

$$\begin{array}{c}\Delta \gamma ={\gamma }_{\xi x}^c-{\gamma }_{\xi y}^c={\gamma }_{\xi x}^c\left(1-\frac{{\tilde{\omega }}_{zy}}{{\tilde{\omega }}_{zx}}\frac{1+\xi }{1-\xi }\right).\end{array}$$

(49)

and it occurs when $\Delta \gamma <0$. This unique effect results from the influence of the matter interactions over the competence between the superradiant-+ and superradiant-y phases resulting from the anisotropy. This behavior does not occur in the Dicke and anisotropic cases [See Figure 3(a1,c1)], provided the onset of the superradiant-y behavior occurs deep in the superradiant-x one. However, for the anisotropic case, the interplay between the anisotropy and the matter interactions might shift the critical coupling ${\gamma }_{\xi y}^c$ to suppress the mode. Beyond that, for the Dicke and anisotropic cases, the behavior of the polariton modes remains qualitatively similar to previous cases, except for the $\pi$ rotation of the energy surfaces as one crosses the superradiant QPT as it can be seen in Figure 3(b2,b3), as well as in Figure 3(c2,c3).

5. Geometric Phase in Presence of Matter Interactions

Next, we investigate the influence of matter interactions over the Berry or geometric phase induced by the bosonic field and the collective pseudospin over both polariton modes. The geometric phase is a topological quantum phase that the eigenstates of a Hamiltonian acquire when one goes adiabatically over a close path in the system’s parameter space, accounting for the Hilbert space geometry [129,130]. The geometric phase has already been studied for the Dicke model [137,138], extensions [140,141], and under the presence of matter interactions in the z-direction [142]. Mainly because it becomes singular at criticality, it has been shown to serve as a signature for the superradiant QPT. Also, matter interactions could modify its scaling behavior [142].

First, for comparison, we obtain the approximated ground-state energy without excitations. According to Equation (3), the scaled energy reads

$$\begin{array}{c}{ϵ}_0=E_0(\alpha ,\beta )/2j=\omega {\alpha }^2+{\omega }_0{\beta }^2-\frac{{\omega }_0}{2}-\gamma \sqrt{k}2\alpha \beta (1+\xi )+{\eta }_{x}k{\beta }^2+{\eta }_z{\left({\beta }^2-\frac{1}{2}\right)}^2.\end{array}$$

(50)

It agrees with that calculated using coherent states as trial states [122]. Then, by substituting Equation (29) into Equation (50), we obtain the ground-state energy for both the normal and superradiant phases

$$\begin{array}{c}{ϵ}_0=\langle {\Psi }_0|2j{H}_0|{\Psi }_0\rangle =-\frac{j{\omega }_0}{2}\left\{\begin{array}{cc}1-\frac{{\eta }_z}{2{\omega }_0}& \gamma <{\gamma }_{\xi x}^c\\ \frac{1}{2}\left({\mu }_x+\frac{1}{{\mu }_x}\right)-\frac{{\eta }_z}{2{\omega }_0}& \gamma >{\gamma }_{\xi x}^c\end{array}\right.\end{array}$$

(51)

where $|{\Psi }_0\rangle$ is the ground-state of the system, and ${\mu }_x(\Delta {\eta }_{zx},\xi )={\left(\Delta {\eta }_{zx}/{\omega }_0+f_{\xi +}\right)}^{-1}$, with $f_{\xi +}={\gamma }^2/{\gamma }_{\xi +}^2$. The ground-state energy signals the superradiant QPT modified by the anisotropy and matter interactions at ${\gamma }_{\xi x}^c$ (${\mu }_x=1$).

To calculate the geometric phase, we follow the standard procedure [138,139,172,173,174,175], where one considers a unitary transformation producing a cyclic trajectory over the parameter space. It may be an arbitrary adiabatic circulation generated by the photon number [140] or over the Bloch sphere [175]. First, let us consider an adiabatic circulation generated by the photon number

$$\begin{array}{c}\widehat{U}\left({\varphi }_n\right)=exp\left(-i{\varphi }_n\widehat{n}\right)=exp\left(-i{\varphi }_n{\widehat{a}}^{†}\widehat{a}\right).\end{array}$$

(52)

It provides an additional phase to the bosonic annihilation (creation) operator $\widehat{a}\to \widehat{a}{e}^{i{\varphi }_n}$ (${\widehat{a}}^{†}\to {\widehat{a}}^{†}{e}^{-i{\varphi }_n}$), effectively displacing them. Next, we introduce a time-dependent unitary transformation, where the parameter ${\varphi }_n\left(t\right)={\omega }_{n}t$ varies adiabatically in the interval $[0,2\pi )$ with angular frequency ${\omega }_n$. The geometric phase is calculated in terms of the circuit integral of the Berry connection

$$\begin{array}{c}{\Gamma }_n=i{\int }_0^{2\pi }d{\varphi }_n\langle \Psi \left({\varphi }_n\right)|\frac{d}{d{\varphi }_n}|\Psi \left({\varphi }_n\right)\rangle =2\pi \langle {\Psi }_0|\widehat{n}|{\Psi }_0\rangle ,\end{array}$$

(53)

where $\langle {\Psi }_0|\widehat{n}|{\Psi }_0{\rangle =\left|\alpha \right|}^2$ is just the expectation value of the photon number within the HPA. Hence, the geometric phase reads

$$\begin{array}{c}{\Gamma }_n=j\left\{\begin{array}{cc}0& \gamma <{\gamma }_{\xi x}^c\\ \frac{\pi }{2}{\left(\frac{\gamma }{{\gamma }_{\xi x}^c}\right)}^2\frac{{\omega }_{zx}^2}{{\omega }_0}(1-{\mu }_x^2)& \gamma >{\gamma }_{\xi x}^c\end{array}\right.\end{array}$$

(54)

For $\xi =1$, we recover the results of the standard [138] and modified Dicke models [142]. Similarly, if now we consider a circulation generated by the excited pseudospin population ${\widehat{J}}_z-j\mathbb{I}$

$$\begin{array}{c}\widehat{V}\left({\varphi }_m\right)=exp\left[-i{\varphi }_m\left({\widehat{J}}_z-j\mathbb{I}\right)\right]=exp\left(-i{\varphi }_m{\widehat{b}}^{†}\widehat{b}\right).\end{array}$$

(55)

Like in the boson case, this shifts the bosonic annihilation (creation) operator $\widehat{b}\to \widehat{b}{e}^{i{\varphi }_m}$ (${\widehat{b}}^{†}\to {\widehat{b}}^{†}{e}^{-i{\varphi }_m}$), effectively displacing them. Changing adiabatically in time allows us to calculate the geometric phase as

$$\begin{array}{c}{\Gamma }_m=i{\int }_0^{2\pi }d{\varphi }_m\langle \Psi \left({\varphi }_m\right)|\frac{d}{d{\varphi }_m}|\Psi \left({\varphi }_m\right)\rangle =2\pi \langle {\Psi }_0|{\widehat{J}}_z-j\mathbb{I}|{\Psi }_0\rangle ,\end{array}$$

(56)

where $\langle {\Psi }_0|{\widehat{J}}_z-j\mathbb{I}|{\Psi }_0{\rangle =\left|\beta \right|}^2$, so

$$\begin{array}{c}{\Gamma }_m=j\left\{\begin{array}{cc}0& \gamma <{\gamma }_{\xi x}^c\\ \pi \left(1-{\mu }_x\right)& \gamma >{\gamma }_{\xi x}^c\end{array}\right.\end{array}$$

(57)

We notice that, for these contours, both the ${\Gamma }_n$ and ${\Gamma }_m$ geometric phases depend only on the ratio $f_{\xi x}$ as in Equation (30), where ${\mu }_x^{-1}={\tilde{\omega }}_{zx}/{\omega }_0\left(f_{\xi x}-1\right)+1$. As a result, the behavior is qualitatively independent from $\xi$ and ${\eta }_y$. That it does not depend on ${\eta }_y$ is because the extreme points in the energy surface belonging to the superradiant-y occur higher in energy, far away from the ground-state energy. Instead, the independence from $\xi$ is a feature that captures the universal behavior of the superradiant QPT. In Figure 4, we show the results for both geometric phases ${\Gamma }_n$ and ${\Gamma }_m$ and compare them to the scaled ground-state energy. We observe that the geometric phase identifies the position of the superradiant QPT as expected. The geometric phases are sensitive to the sign of ${\eta }_x$, increasing with ${\eta }_x$. Furthermore, we observe that the geometric phase also helps to identify the first-order phase transition occurring at $\Delta {\eta }_{zx}={\omega }_0$, where the geometric phase becomes zero independently of $\gamma$ in both cases and then turns negative.

6. Discussion and Conclusions

We derived analytic expressions for low-energy excitations of the anisotropic Dicke model in the presence of matter interactions by employing the Holstein–Primakoff transformation and neglecting higher-order terms in the thermodynamic limit. We explored the modes in the normal (ordered) phase, deformed by matter interactions, and the superradiant (ordered) phase from a unified point of view, whose interpretation is supported by the analysis of energy surfaces obtained using the corresponding classical Dicke model via coherent states [122].

The anisotropy in the Dicke model that, in this work, tunes the Hamiltonian between the Tavis–Cummings to the standard Dicke models produces a new set of extreme points and quantum phases, the superradiant-$(+)$ and $(-)$ phases, where superradiance is determined by the $(+)$ phase, given that it lies lower in energy. Matter interactions serve two purposes: to shift the values of critical couplings and ground-state energy and to produce new extreme points in the energy surface, reflecting the onset of new quantum phases. Particularly, interactions in the x and y directions directly modify the criticality of the $(+)$ and $(-)$ phases, respectively, and shift their relative position in terms of the spin–boson coupling.

Beyond the shifting of critical couplings and ground-state energy associated to matter–matter interactions, we found unique effects from combining them with the anisotropy. In most cases, amplitude, massive, or upper-polariton mode behave qualitatively the same, just changing energetically to reflect the deformations of the energy surface, so interactions increase the gap to the phase mode, creating an energy cost given by ${\omega }_0{\tilde{\omega }}_{zx}{\tilde{\omega }}_{zy}$. Novel behavior appears mostly over the ground-state mode, i.e., the phase or lower-polariton mode. For the Dicke model, it becomes a roton mode that, for larger couplings, tends to the non-interacting photon energy $\omega$ and decreases in energy with the anisotropy. Without interactions, it eventually transforms into the Goldstone mode. Instead, the TC limit exhibits a Goldstone mode in the absence of matter–matter interactions that becomes rotonic as the latter breaks the $U\left(2\right)$ symmetry, leaving conserved only the $Z_2$ one.

A major result is that tuning the combination of matter interactions in the z and y directions can reduce the rotonic mode energy in the TC Hamiltonian, producing a “variable” mass roton mode, or a lower polariton that starts losing matter content with interactions. This is understood in terms of the energy surface, as the interactions’ influence on higher energies becomes less relevant to the ground state, thus restoring to some extent the Goldstone-like behavior. Second, we discover that the Goldstone mode can be suppressed by tuning the interactions in the z and x directions by shifting the relative position between the critical couplings of the superradian-x and -y phases. This suppression does not affect the amplitude mode.

We anticipate these features could be detected and employed as an experimental signature to observe the Goldstone mode in novel setups, such as spin–magnon systems, where the interactions are tunable [109,110,111], as well as in others within the broad range of setups where the Dicke model can be realized, from quantum optics and atomic physics to condensed matter.

We also calculated the geometric phase for both boson and pseudospin contours, which equals the expectation value of the photon number and excited pseudospin in the phase mode within this approximation. We found that the role of the z and x interactions is to shift the phase. Still, it is independent of the y interactions and the anisotropy parameter $\xi$, as it is a low-energy approach capable only of capturing the universal behavior of the superradiant-x transition. Furthermore, the geometric phase becomes singular at the superradiant QPT, serving as a signature of criticality, where we recover previous results from refs. [119,142]. It also detects the first-order QPT that results from changing matter–matter interactions by changing sign.

We conjecture that the predominance of the superradiant-x phase will be inverted to the superradiant-y phase in the case of the anti-Tavis–Cummings model, where now the rotating terms become smaller than the counter-rotating ones [18,92,176]. Our model can be extended to explain this case, but this goes beyond the scope of this work. Finally, we note that establishing a quantitative comparison of the polariton modes with finite-size quantum exact numerical solutions as a function of energy is desirable [150]. Nevertheless, this analysis and exploring the effect of the anisotropy and matter interactions over the polariton modes beyond the mean field are left for future work [82,177].