Unified Superradiant phase transitions

TL;DR

This paper demonstrates that superradiant phase transitions in Dicke and classical oscillator models are fundamentally the same, providing a unified theoretical framework and analyzing their structure and symmetry properties.

Contribution

It offers a unified proof that superradiant phase transitions are of the same type across different models, with exact mean-field analysis and symmetry considerations.

Findings

Superradiant transitions are of the same type in Dicke and classical oscillator models.

Mean-field approximation is exact for these transitions.

General features of phase structure are uncovered across models.

Abstract

We prove, by means of a unified treatment, that the superradiant phase transitions of Dicke and classical oscillator limits of simple light-matter models are indeed of the same type. We show that the mean-field approximation is exact in both cases, and compute the structure and location of the transitions in parameter space. We extend this study to a fuller range of models, paying special attention to symmetry considerations. We uncover general features of the phase structure in the space of parameters of these models.

Figures

Equations

$$H=\omega a^{†}a+\sum _{i=1}^N\frac{g}{\sqrt{N}}{\sigma }_{ix}(a+a^{†})+\sum _i\Delta {\sigma }_{iz}.$$

$$H=\Delta \sum _{i=1}^N[\frac{\omega C}{N\Delta }\frac{a^{†}a}{C}+\frac{g\sqrt{C}}{\sqrt{N}\Delta }{\sigma }_{ix}\frac{(a+a^{†})}{\sqrt{C}}+{\sigma }_{iz}].$$

$$\frac{⟨a^{†}a⟩}{N\Delta /\omega }\approx -\frac{\gamma }{2}{\partial }_{\gamma }\varphi (u_{\min })=u_{\min }^2.$$

$$H$$

$$H=\omega a^{†}a+\sum _{i=1}^N[{\sigma }_{ix}[\frac{g}{\sqrt{N}}(a+a^{†})+\frac{g^{\prime }}{N}(a^2+(a^{†}{)}^2)]+\Delta {\sigma }_{iz}].$$

$$H=\omega a^{†}a+\sum _{j=1,2}[g_j{\sigma }_{jx}(a+a^{†})+{\Delta }_j{\sigma }_{jz}]+\sum _{\alpha =x,y,x}{J}^{(\alpha )}{\sigma }_{1\alpha }{\sigma }_{2\alpha }.$$

$$H=\sum _{\nu =1}^M{\omega }_{\nu }{a}_{\nu }^{†}{a}_{\nu }+\sum _{i=1}^N\sum _{\nu =1}^M\frac{g_{\nu }}{\sqrt{N}}{\sigma }_{ix}(a_{\nu }+a_{\nu }^{†})+\sum _{i=1}^N{\Delta }_i{\sigma }_{iz}.$$

$$\frac{⟨a^{†}a⟩}{C}=\frac{\partial f}{\partial \Omega },\frac{1}{N}⟨\sum _{i=1}^N{\sigma }_{ix}\frac{a+a^{†}}{\sqrt{C}}⟩=\frac{\partial f}{\partial \gamma }.$$

$$Z=\mathrm{Tr}\exp (-\beta H)=\int {\mathrm{Tr}}_{sp}⟨\alpha |e^{-\beta H(a^{†},a)}|\alpha ⟩\frac{d^2\alpha }{\pi },$$

$$\tilde{Z}\le Z\le \exp \{\beta \omega \}\tilde{Z},$$

$$\tilde{Z}=\frac{N\Delta }{\omega }\int {\mathrm{Tr}}_{spin}\exp \{-\beta N\Delta h({\alpha }^{\prime },{\alpha }^{\prime \ast })\}\frac{d^2{\alpha }^{\prime }}{\pi },$$

$$\frac{-\ln \tilde{Z}}{\beta N\Delta }\ge \frac{-\ln Z}{\beta N\Delta }\ge -\frac{\omega }{N\Delta }-\frac{\ln \tilde{Z}}{\beta N\Delta }.$$

$$H=\sum _{\nu }\left[{\omega }_{\nu }{a}_{\nu }^{†}{a}_{\nu }+a_{\nu }^{†}B+a_{\nu }^{†}{B}^{†}\right]+\sum _i\Delta {\sigma }_{iz},$$

$$\tilde{Z}\le Z\le \exp (\beta \sum _{\nu }{\omega }_{\nu })\tilde{Z}.$$

$$H=\omega a^{†}a+A+a^{†}B+a{B}^{†}+a^{{†}^2}D+a^2{D}^{†}.$$

$$|\alpha ,n⟩=P_n|\alpha ⟩,$$

$$⟨\alpha ,n|\beta ,n⟩=K_n(\alpha ,\beta )=\exp [-\frac{1}{2}(|\alpha {|}^2+|\beta {|}^2)]\sum _{m=0}^n({\alpha }^{\ast }\beta {)}^m/m!,$$

$$⟨\alpha ,n|a|\alpha ,n⟩=\alpha K_{n-1}(\alpha ,\alpha ),⟨\alpha ,n|a^{†}|\alpha ,n⟩={\alpha }^{\ast }{K}_{n-1}(\alpha ,\alpha ),$$

$$⟨\alpha ,n|a^2|\alpha ,n⟩={\alpha }^2{K}_{n-2}(\alpha ,\alpha ),⟨\alpha ,n|a^{{†}^2}|\alpha ,n⟩={\alpha }^{{\ast }^2}{K}_{n-2}(\alpha ,\alpha ).$$

$$Z_n=\mathrm{Tr}\left[P_n{e}^{-\beta H}{P}_n\right]\ge {\mathrm{Tr}}_{\mathfrak{H}}\int \frac{d^2\alpha }{\pi }{K}_n(\alpha ,\alpha )\times \exp \{-|\alpha {|}^2{K}_{n-1}(\alpha ,\alpha )/K_n(\alpha ,\alpha )+A$$

$$+[B{\alpha }^{\ast }{K}_{n-1}(\alpha ,\alpha )/K_n(\alpha ,\alpha )+H.C.]+[D{\left({\alpha }^{\ast }\right)}^2{K}_{n-2}(\alpha ,\alpha )/K_n(\alpha ,\alpha )+H.C.]\}$$

$$Z_n(ϵ)\le {\pi }^{-1}\int {\mathrm{d}}^2\alpha K_n(\alpha ,\alpha ){\mathrm{Tr}}_{\mathfrak{H}}\exp \{-\beta H_n(\alpha ,ϵ)\},$$

$$H_n=P_{n}H{P}_n=\frac{1}{\pi }\int \mathrm{d}\alpha H_n(\alpha )|\alpha ,n⟩⟨\alpha ,n|,$$

$$H_n(\alpha ,ϵ)=H_n(\alpha )e^{-ϵ|\alpha {|}^2}.$$

$$H_n(\alpha )=\omega (|\alpha {|}^2-1)+A+B{\alpha }^{\ast }+B^{†}\alpha +D{\alpha }^{{\ast }^2}+D^{†}{\alpha }^2.$$

$$H=\omega a^{†}a+\sum _{i=1}^N\left[\frac{g_i}{\sqrt{N}}(a^{†}{\sigma }_i^{-}+a{\sigma }_i^{+})+\frac{g_i{\lambda }_i}{\sqrt{N}}(a{\sigma }_i^{-}+a^{†}{\sigma }_i^{+})+{\Delta }_i{\sigma }_z\right].$$

$${\varphi }_1(x,y)=\beta \omega (x^2+y^2)-\sum _{i=1}^N\ln 2\cosh [\beta \sqrt{{\Delta }_i^2+\frac{g_i^2(1+{\lambda }_i{)}^2{x}^2}{N}+\frac{g_i^2(1-{\lambda }_i{)}^2{y}^2}{N}}].$$

$$\varphi (u,v)=u^2+v^2-\sum _{i=1}^N\frac{1}{\beta N\Delta }\ln 2\cosh [\beta \Delta \sqrt{{\delta }_i^2+{\gamma }_i^2(1+{\lambda }_i{)}^2{u}^2+{\gamma }_i^2(1-{\lambda }_i{)}^2{v}^2}],$$

$$\frac{\partial \varphi }{\partial u}=2u-\frac{1}{N}\sum _{i=1}^N\tanh [\beta \Delta \sqrt{{\delta }_i^2+{\gamma }_i^2(1+{\lambda }_i{)}^2{u}^2+{\gamma }_i^2(1-{\lambda }_i{)}^2{v}^2}]\frac{{\gamma }_i^2(1+{\lambda }_i{)}^{2}u}{\sqrt{{\delta }_i^2+{\gamma }_i^2(1+{\lambda }_i{)}^2{u}^2+{\gamma }_i^2(1-{\lambda }_i{)}^2{v}^2}}=0,$$

$$\frac{\partial \varphi }{\partial v}=2v-\frac{1}{N}\sum _{i=1}^N\tanh [\beta \Delta \sqrt{{\delta }_i^2+{\gamma }_i^2(1+{\lambda }_i{)}^2{u}^2+{\gamma }_i^2(1-{\lambda }_i{)}^2{v}^2}]\frac{{\gamma }_i^2(1-{\lambda }_i{)}^{2}v}{\sqrt{{\delta }_i^2+{\gamma }_i^2(1+{\lambda }_i{)}^2{u}^2+{\gamma }_i^2(1-{\lambda }_i{)}^2{v}^2}}=0.$$

$$\frac{1}{N}\sum _{i=1}^N\tanh [\beta \Delta {\delta }_i]\frac{{\gamma }_i^2(1-{\lambda }_i{)}^2}{{\delta }_i}=\sum _{i=1}^N\tanh [\beta {\Delta }_i]\frac{g_i^2(1-{\lambda }_i{)}^2}{N\omega {\Delta }_i}\ge 2.$$

$$\frac{1}{N}\sum _{i=1}^N\tanh [\beta \Delta {\delta }_i]\frac{{\gamma }_i^2(1+{\lambda }_i{)}^2}{{\delta }_i}=\sum _{i=1}^N\tanh [\beta {\Delta }_i]\frac{g_i^2(1+{\lambda }_i{)}^2}{N\omega {\Delta }_i}\ge 2,$$

$$\sum _{i=1}^N\tanh [\beta {\Delta }_i]\frac{2g_i^2}{N\omega {\Delta }_i}>1.$$

Full text

Introduction.– Since the groundbreaking result of Dicke, in which he found a superradiant state for a localized radiating gas PhysRev.93.99 , the very concept of superradiance has been actively investigated. In a superradiant state, the average photon number is macroscopic, in that the ratio $\langle a^{\dagger}a\rangle/N$ is finite in the limit $N\to\infty$, where $N$ is the number of intervening atoms. The existence of superradiant phase transitions (SPT) was rigourously shown hepp ; wang in the Tavis–Cummings model PhysRev.170.379 , which describes coupling of $N$ spins to a single bosonic mode with dipolar coupling, but restricts it the Rotating Wave Approximation (RWA). Those early works also considered the extension to the multimode case, beyond the RWA and for a multitude of analogous models PhysRevA.8.1440 ; PhysRevA.8.2517 ; zhang ; PhysRevA.95.053854 ; guo ; li ; li1 ; GJC18 .

One more recent development is the realization that a similar phase transition is also present in a different limit, which in the Rabi model case corresponds to the photon frequency tending to zero compared to qubit frequency Bakemeier ; PhysRevA.87.013826 ; hwang1 ; hwang2 ; liu1 ; lv . This limit is understood as a classical oscillator Bakemeier . They are realized at zero temperature larson , thus falling prima facie under the purview of Quantum Phase Transitions (QPT) sachdev_2011 (see however the discussion in larson to this point). Following the theoretical developments, there has also been keen experimental interest in realizing at least some of these transitions Baumann:2010aa ; PhysRevLett.118.073001 .

In this letter we put forward a unifying treatment of SPTs in a wide class of models br ; chen ; jorge ; ion ; jin ; pj0 ; guxin ; jorge1 ; KMLSN19 , including both the thermodynamical and the classical oscillator limit. In all these models there is a competition between the free hamiltonian for photons and spins, which tends to the normal phase, and an interaction that, when increased, produces the superradiant phase. Furthermore, in all these models the mean field approximation for the bosons PhysRevA.8.2517 will be exact in the relevant limit. The physical phases are then studied by analysing the extrema of a Landau potential. The order parameter at hand is the macroscopic boson expectation number. The symmetry or lack thereof of the underlying Hamiltonian will be apparent in the Landau potential, and the presence of second and first order superradiant phase transitions will be directly related to it.

Unified analysis.– We present a unified treatment of the thermodynamic limit of the Dicke model and the classical oscillator model of the quantum Rabi model and SPTs thereof. Thus, we first consider the family of models ($\hbar=1$)

[TABLE]

As usual, $a$ and $a^{\dagger}$ denote annihilation and creation operators for a bosonic mode (the “photon”), while the Pauli matrices also carry a spin index. The parameters of this family are $\omega$ (photon energy), $\Delta$ (spin energy), $g$ (coupling strength) and $N$ (number of spins), together possibly with inverse temperature $\beta=1/k_{B}T$. We will henceforward denote this model as the Dicke–Quantum Rabi (DQR) model.

In the DQR model in the thermodynamic limit the order parameter is the ratio $\langle a^{\dagger}a\rangle/N$, which, if finite in the limit $N\to\infty$ for either the thermal or the ground state is the diagnostic of the superradiant phase.

We now endeavour to identify regions of parameter space that can be understood as the infinite level set of a function $C$ on parameter space and for which the ratio $\langle a^{\dagger}a\rangle/C$ is an order parameter, as first hinted at in larson . We thus rewrite the Hamiltonian in the form

[TABLE]

If indeed a transition is to be found, all three terms in the expression above must appear, since otherwise there will be no competition that accounts for the change of phase. Furthermore, they are to be examined under the hypothesis that $\langle a^{\dagger}a\rangle/C$ be an order parameter. This amounts to requiring that in the limit $C\to\infty$ the effective couplings $\Omega=\omega C/N\Delta$ and $\gamma=g\sqrt{C}/\Delta\sqrt{N}$ tend to a finite value when $C\to\infty$.

Thus the region of parameter space of interest is given by the functions $C=\Omega\left(N\Delta/\omega\right)\to\infty$, with $\Omega$ a finite nonzero quantity, while at the same time $g^{2}/\omega\Delta$ is kept finite. One way of reaching this region is to take $N\to\infty$, while keeping finite all other parameters; this corresponds to the usual thermodynamical limit. An alternative, the classical oscillator limit, is achieved by considering $\Delta/\omega\to\infty$, as well as $g/\omega=(\gamma/\sqrt{\Omega})\sqrt{\Delta/\omega}\to\infty$, while $\gamma$ and $\Omega$ are finite, together with $N$. This analysis can be extended to the finite temperature case sl .

Further to the thermodynamical and the classical oscillator limits both being included in the infinite level set of the combination of parameters $N\Delta/\omega$, they are also determined there by the mean field approximation. Namely in that the mean field partition function $\tilde{Z}$ and the full partition function $Z$ will coincide in that limit. In a more concrete way, the mean field reduced free energy $\tilde{f}=-\ln\tilde{Z}/\beta N\Delta$ and the full reduced free energy $f=-\ln Z/\beta N\Delta$ are equal up to $O(\omega/N\Delta)$ sl .

Following this result, the analysis of phases in terms of the mean field approximation is exact in the relevant limit $\omega/N\Delta\to 0$ with finite $g^{2}/\omega\Delta$. For the case at hand, namely for hamiltonian (1), the mean field partition function is computed to be $\tilde{Z}=\sqrt{N\Delta/\pi\beta\omega^{2}}\int_{-\infty}^{+\infty}\mathrm{d}u\,\exp\left[-\beta N\Delta\phi(u)\right]\,,$ where $\phi(u)=u^{2}-\ln\left[2\cosh\left(\beta\Delta\sqrt{1+4\gamma^{2}u^{2}}\right)\right]/\beta\Delta$ plays the role of a Landau potential. This last statement is the consequence the integral being well approximated by $\exp\left[-\beta N\Delta\phi(u_{\mathrm{min}})\right]$ in the limit $\beta N\Delta\to\infty$, with $u_{\mathrm{min}}$ the location of the global minimum of $\phi(u)$. In other words, $f=-\ln\tilde{Z}/\beta N\Delta\approx\phi(u_{\mathrm{min}})$. Notice that $\beta N\Delta=\beta\omega(N\Delta/\omega)$, so the limit we were considering from the outset implies the validity of the further approximation of the integral by Laplace’s method, as long as $\beta\omega\neq 0$.

The superradiant or normal character of each relevant phase is determined by the order parameter $\langle a^{\dagger}a\rangle/(N\Delta/\omega)$, which is well approximated by the mean field prediction $\omega\partial_{\omega}f=\omega\partial_{\omega}\phi(u_{\mathrm{min}})$. Since $\omega\partial_{\omega}\gamma=-\gamma/2$ and $\gamma\partial_{\gamma}\phi(u)=u\phi^{\prime}(u)-2u^{2}$, we conclude

[TABLE]

Thus, the system is in a normal phase if the global minimum of the Landau potential $\phi(u)$ is zero, and in a superradiant phase if the origin is no longer the global minimum. The critical line separating the two regions is $\tanh(\beta\Delta)-\frac{1}{2\gamma_{c}^{2}}=0$. For couplings smaller than the critical coupling $\gamma_{c}(\beta\Delta)$ determined by this equation the system is in its normal phase, and for $\gamma>\gamma_{c}$ in a superradiant one. Close to criticality $u^{2}_{\mathrm{min}}$ is proportional to $\gamma-\gamma_{c}$, and the transition is second order. Notice that in this unified treatment we have recovered both the $N\to\infty$ Dicke SPT PhysRevA.8.1440 and the Quantum Rabi classical oscillator SPT hwang1 .

Symmetry.– The fact that Hamiltonian (1) presents a $\mathbbm{Z}_{2}$ symmetry is reflected in $\phi(u)$ being an even function. This, together with continuity and $u^{2}$ being the dominant term of $\phi(u)$ for large $|u|$, fully determines that the phase transition is second order. To better establish the connection between the $\mathbbm{Z}_{2}$ symmetry and the Landau potential being an even function, we consider a single mode inhomogenous anisotropic extension of the DQR model,

[TABLE]

We term this model inhomogeneous because the dipole coupling constants $g_{i}$ are different for different qubits, and similarly the spin energies $\Delta_{i}$ and the anisotropy parameters $\lambda_{i}$. These are anisotropies in that rotating and counterrotating terms are differently presented. In the general case ($\lambda_{i}\neq 0$ at least for one value of $i$), this model has a $\mathbbm{Z}_{2}$ symmetry with non trivial group element $\Pi=\exp[i\pi(a^{\dagger}a+\sum_{i}^{N}(1+\sigma_{iz})/2)]$. Since, for coherent states $|\alpha\rangle$, $\langle\alpha|\Pi H\Pi|\alpha\rangle=U\langle-\alpha|H|-\alpha\rangle U^{\dagger}$, with $U=U^{\dagger}=\prod_{i=1}^{N}\sigma_{i}^{z}$ unitary, the mean field Landau potential is even in $\alpha$.

Let $\Delta$ be a positive generalized mean of the set $\{\Delta_{i}\}$. Then the mean field approximation is asymptotically exact as $N\Delta/\omega\to\infty$, as before, and again the phase diagram is determined by the mean field Landau potential. As shown in sl there are indeed in this limit SPT, all of them second order and of mean field type, with $\mathbbm{Z}_{2}$ symmetry breaking. For finite $N\Delta/\omega\gg 1$, the mean field approach is still approximately valid, and the $\mathbb{Z}_{2}$ symmetry should be preserved. So the ground state will take the approximate form of a “Schrödinger cat”-like state for finite $N$ sl . Notice that this model includes as a special case the DQR Hamiltonian (1). A very different symmetry is present if $\lambda_{i}=0$ for all $i$, i.e. for the inhomogeneous Tavis–Cummings model. Namely a $U(1)$ symmetry generated by the conserved quantity $a^{\dagger}a+\sum_{i}^{N}(1+\sigma_{iz})/2$. In this case the Landau potential depends on $\alpha$ only through $|\alpha|^{2}$, reflecting the underlying continuous symmetry, and the spontaneous symmetry breaking of the second order SPT will be accompanied by the presence of a Goldstone mode sl .

To stress the role of symmetry, consider the explicit breaking of the $\mathbbm{Z}_{2}$ by the addition of a static bias $\epsilon_{i}\sigma_{ix}$ in the previous Hamiltonian. In this case there is no SPT, as expected, but rather a continuous increase of the order parameter with the couplings sl .

1st order SPT without symmetry.–

A second order phase transition controlled by an order parameter requires a symmetry to be broken. Not so however for a first order transition to take place. As an example of this assertion in the context at hand, consider the QRM Hamiltonian extended, homogeneously, with a two-photon coupling term ying ; chen1 ,

[TABLE]

The stability of the system requires $g^{\prime}/\omega<1/2$ pj1 ; PhysRevA.92.033817 . In the generic case $g,g^{\prime}\neq 0$ this Hamiltonian has no symmetry. Following the by now well established routine, we obtain the reduced free energy as the global minimum of a potential $\phi(u)=u^{2}-\frac{1}{\beta\Delta}\ln\left[2\cosh\left(\beta\Delta\sqrt{1+\left[2\gamma u+2\gamma^{\prime}u^{2}\right]^{2}}\right)\right]\,,$ where $\gamma=\frac{g}{\sqrt{\Delta\omega}}$ and $\gamma^{\prime}=g^{\prime}/\omega$ sl . The inexistence of the symmetry is reflected in the Landau potential not being even for the generic case. Nonetheless, $u=0$ is always an extreme, and indeed the global minimum when $\gamma$ and $\gamma^{\prime}$ are close to [math]. On further analysis, there is a critical region in parameter space that separates normal and superradiant phases, determined by the existence of a nonzero solution of $\phi(u)=\phi(0)$ and $\phi^{\prime}(u)=0$ simultaneously. When $\beta\Delta\to\infty$, this set of equations is algebraic and allows the determination of the explicit critical line $2\gamma^{2}+4\gamma^{\prime 2}=1$, coinciding with the single qubit QPT case ying . We present the phase diagram for the general case in Fig. 1(a). Given the lack of symmetry, the transition has to be first order, and this is indeed checked for $\beta\Delta\to\infty$, since just above the critical line the rescaled photon number $u(\gamma_{c},\gamma^{\prime}_{c})_{\min}^{2}=\left(2\gamma^{\prime}_{c}/\gamma_{c}\right)^{2}\neq 0$, as can be seen from Fig. 1(b).

General structure of the phase diagram.–

We now present a model with a discrete symmetry for which both continuous and discontinuous SPTs exist, and use it to illustrate general properties of the phase diagrams of the family of models under consideration. In particular, we now study two spins which, on top of being dipole coupled to a photon mode, present XYZ spin-spin interactions pj2 , with $\mathbbm{Z}_{2}$ symmetry,

[TABLE]

We study its phase diagram in the classical oscillator limit $\Delta/\omega\to\infty$. Defining the reduced free energy $f$ from the partition function $Z=\exp\left(-\beta\Delta f\right)$, under the by now standard assumptions it will be determined as $\phi(u_{\mathrm{min}})$ at the global minimum of a Landau potential $\phi(u)=u^{2}+\lambda(u)$, where $\lambda(u)$ is the smallest eigenvalue of the two spin operator $h(u)=\sum_{j=1}^{2}\left[2\gamma_{j}u\sigma_{jx}+\delta_{j}\sigma_{jz}\right]+\sum_{\alpha=x,y,z}\epsilon_{\alpha}\sigma_{1\alpha}\sigma_{2\alpha}\,,$ where $\Delta$ is some positive generalised mean of $\{\Delta_{j}\}$, $\gamma_{j}=g_{j}/\sqrt{\Delta\omega}$, $\delta_{j}=\Delta_{j}/\Delta$, $\epsilon_{\alpha}=J^{(\alpha)}/\Delta$, and the order parameter $\omega\langle a^{\dagger}a\rangle/\Delta$ is computed to be $u_{\mathrm{\min}}^{2}$.

The phase space of interest, thus, is seven dimensional, and we look at two dimensional sections controlled by the dipole couplings, organised in a vector $\vec{\gamma}$. In a generic section there will be both continuous and discontinuous phase transitions. The origin $\vec{\gamma}=0$ is always in normal phase, since the effective hamiltonian $h(u)$ is then independent of $u$. Furthermore, since the Landau function presents as $\phi(u)=u^{2}+f(u\vec{\gamma})$, with $f$ a function of several (here two) variables, we see that at critical points $u_{*}$ such that $\phi^{\prime}(u_{*})=0$ one has $\vec{\gamma}\cdot\nabla_{\gamma}\phi(u_{*})=-2u_{*}^{2}$. Thus, for non-zero critical points, negative. Additionally, as $h(u)$ is a bounded operator, one sees that for large $\|\vec{\gamma}\|$ the minimum of $\phi$ can be approximately located at $|\gamma_{1}|+|\gamma_{2}|$, with the prediction that the system will be in its superradiant phase. It follows that there will be a region of normal phase around the origin, and that the radial component of the normal to its boundary, i.e. to the critical line, will not be zero.

For fixed $\vec{\delta}$, $\{\epsilon_{\alpha}\}$, one can identify whether the transition is first or second order in a given direction in the $\vec{\gamma}$ plane by several methods. First, the second derivative of $\phi(u)$ at zero can be obtained by a second order perturbative computation. If $\phi^{\prime\prime}(0)$ does not change sign as $\|\vec{\gamma}\|$ grows in a fixed direction, the transition in that direction will be necessarily discontinuous. Alternatively, since the origin in $u$ space cannot be the (unique) global minimum if there exists a point $u_{s}\neq 0$ such that $\phi(0)=\phi(u_{s})$, this equation is analysed. When this equation is inserted into the secular equation for $h(u)$ a cubic for the variable $u_{s}^{2}$ is obtained. The study of its discriminant provides us with criteria for the existence and kind of solution, and thus for the location and character of the transition. In fact, the quantity $\phi^{\prime\prime}(0)$ will determine the value at the origin of the cubic, thus connecting both methods.

Notice that in any case there are symmetries in phase space on which we have relied in the discussion above. Thus, the effective spin hamiltonian $h(u)$ is isospectral with $h(-u)$ as a consequence of the $\mathbb{Z}_{2}$ symmetry of the initial Hamiltonian. This leads to the symmetry of phase space with respect to $\vec{\gamma}\to-\vec{\gamma}$. There are additional symmetries of phase space that can be identified in a similar manner. Thus, for instance, under the transformation $\vec{\delta}\to-\vec{\delta}$.

We illustrate this general analysis with the portrayal of one $\vec{\gamma}$ plane section in Fig. 2.

We see in this manner that in the general class of spin models dipole coupled to a bosonic mode, if there is a symmetry, there will be a region of normal phase in the neighbourhood of the origin of the space of dipolar couplings, and that generally one flows radially in this space to a superradiant phase. The boundary will be composed generically of patches of first order and second order transition lines.

The multimode DQR model.- In the single mode analysis above we observe that the possiblity exists that by progressively increasing one dipole coupling the system can go from superradiant to normal to superradiant again. Thus, by using multiqubit systems we can have sequences of transitions, whose experimental realizations would be of great interest.

Even though the role of the bosonic mode and the spins are radically different, the fact above suggests that an extension to several bosonic modes (multi-mode in what follows) might provide us with rich phenomenology, even more amenable to experimental realization. Indeed multimode models with dipole couplings can be analyzed with the techniques presented here sl , and will have superradiant phases.

In particular, let us examine the following Hamiltonian:

[TABLE]

In this case the dipole couplings are homogeneous for the qubits. Even if that were not the case, it is easy to observe that a $\mathbbm{Z}_{2}$ exists, with generator $\Pi=\exp\left\{i\pi\left[\sum_{\nu}a_{\nu}^{\dagger}a_{\nu}+(1/2)\sum_{i}\left(1+\sigma_{iz}\right)\right]\right\}$. Furthermore, the system is in normal phase if the vector of dipolar couplings is zero. Following the steps of the previous analyses, we identify a set of order parameters $\omega_{\nu}\langle a_{\nu}^{\dagger}a_{\nu}\rangle/N\Delta$, with $\Delta$ some average of the spin frequencies $\Delta_{i}=\Delta\delta_{i}$. Similarly we arrive at a Landau function of $M$ variables $\{u_{\nu}\}_{\nu=1}^{M}$, whose global minimum determines the phase. In fact, these variable appear in the Landau function only through the combinations $\vec{u}^{2}$ and $\vec{\gamma}\cdot\vec{u}$, where we use vector notation for the corresponding $M$ dimensional objects, in the form $\phi(\vec{u})=\vec{u}^{2}-s\left[\left(\vec{\gamma}\cdot\vec{u}\right)^{2}\right]$ for some function $s$. The critical points of the Landau function, other than the origin, are therefore parallel to the vector $\vec{\gamma}$. Analogously to the analysis above, we see that for critical points other than the origin $\vec{\gamma}\cdot\nabla_{\gamma}\phi=-2\vec{u}^{2}$ is negative, whence the radial flow in coupling space is towards a superradiant phase. The coupling line is determined by $1/\vec{\gamma}_{c}^{2}=2\sum_{i}\tanh(\beta\Delta_{i})/N|\delta_{i}|$. The transition is continuous, with the order parameters behaving radially with mean field exponent. Since criticality is controlled by the length of the effective coupling vector $\vec{\gamma}$, one can achieve superradiance in a cooperative way. We term this phenomenon assisted SPT.

Conclusion and perspective.– We put forward the analytic study of SPT in a unified manner for both the thermodynamic and the classical oscillator limits of a wealth of boson-spin models with dipolar coupling. We show that in both cases the mean field approximation is exact, as was expected in the literature. Since we study the transitions in terms of a Landau potential, the analysis of symmetry comes to the fore, and we highlight its relevance in a general family of models. We demonstrate that general features of phase space in terms of dipolar couplings can be recovered from this approach, and portray anisotropies and first order/continuous boundaries. We expect that this general framework will be a good setting for new experimental proposals, given the wide range of phenomenology now available, including assisted SPT and the identification of first order transitions.

Acknowledgements.– We are thankful to Ricardo Puebla for helpful discussions. This work was supported by the the program of China Scholarship Council (No. 201707230010), National Natural Science Foundation of China (11704320, 11604240), Natural Science Foundation of Hunan Province, China ( 2018JJ3482), the National Basic Research Program of China (2015CB921103), the Program for Changjiang Scholars and Innovative Research Team in University (No. IRT13093), Spanish MINECO/FEDER FIS2015-69983-P, Basque Government IT986-16, QMiCS (820505) and OpenSuperQ (820363) of the EU Flagship on Quantum Technologies.

Supplementary Material: Superradiant phase transitions: exact unified presentation

S1 Magnitude of parameters for SPT at finite temperature

In the discussion of Eq. (2) we have stressed the need for two parameters to be finite. Indeed, the finiteness of $\Omega$ and $\gamma$ has been predicated on our desire to identify a SPT at zero temperature. Nonetheless, even at finite temperature we would require their presence, as will be warranted by considering the partition function $Z=\mathrm{Tr}\left[\exp(-\beta H)\right]$ for the Hamiltonian in (2). The reduced free energy $f$, defined as $f=F/(N\Delta)=-(\ln Z)/(\beta N\Delta)$, should be finite for SPT to be possible in all the (generalised) thermodynamic limits we examine. Denoting, as in the main text, the purportedly nonzero coefficients $\omega C/(N\Delta)=\Omega$ and $g\sqrt{C}/(\Delta\sqrt{N})=\gamma$, we find that in a thermal state,

[TABLE]

Since SPT takes place if the first quantity goes from zero in a region to non-zero elsewhere, the reduced free energy and its derivatives have to be finite, whence we require finiteness of $\Omega$ and $\gamma$ also in a thermal state.

S2 Validity of mean field approach on SPT

Here we study the validity of the mean field approximation in

[TABLE]

where $H(a^{\dagger},a)=\omega a^{\dagger}a+\sum_{i=1}^{N}\left[\frac{g}{\sqrt{N}}\sigma_{ix}(a+a^{\dagger})+\Delta\sigma_{iz}\right]$.

Hepp and Lieb PhysRevA.8.25170 obtained the upper and lower bound

[TABLE]

where $\tilde{Z}=\int\mathrm{Tr}_{spin}\exp\{-\beta\left[\omega|\alpha|^{2}+\sum_{i=1}^{N}\frac{g}{\sqrt{N}}\sigma_{ix}(\alpha+\alpha^{*})+\sum_{i=1}^{N}\Delta\sigma_{iz}\right]\}\frac{d^{2}\alpha}{\pi}$ is the mean field approximation of $Z$. Here we proceed to prove the asymptotic validity of the mean field approximation when ${N\Delta}/{\omega}$ tends to infinity, while maintaining ${g}/{\sqrt{\Delta\omega}}$ finite.

After some rearrangements and changing variables to $\alpha=\sqrt{{N\Delta}/{\omega}}\alpha^{\prime}$, we obtain

[TABLE]

where the effective spin hamiltonian $h(\alpha^{\prime},\alpha^{*\prime})$ is defined as $|\alpha^{\prime}|^{2}+\sum_{i=1}^{N}({g}/{N\sqrt{\omega\Delta}})\sigma_{ix}(\alpha^{\prime}+\alpha^{\prime*})+{\sum_{i=1}^{N}\sigma_{iz}}/{N}$. After tracing out the qubit part, $\tilde{Z}=({N\Delta}/{\omega})\int\exp\{-\beta N\Delta\phi_{1}(\alpha^{\prime},\alpha^{*\prime})\}{d^{2}\alpha^{\prime}}/{\pi}$, with all coefficients of $\alpha^{\prime}$ and $\alpha^{*\prime}$ in $\phi_{1}(\alpha^{\prime},\alpha^{*\prime})$ finite in the limit of interest.

The bounds in Eq. (S3) can be rewritten to

[TABLE]

Thus, if indeed $f=-{\ln\tilde{Z}}/{\beta N\Delta}$ is finite in the limit ${N\Delta}/{\omega}\to\infty$, the reduced free energy ${F}/{N\Delta}={-\ln Z}/{\beta N\Delta}$ is asymptotically determined by $f$. For all finite values of $\beta\omega$, it is the case that $f$ is asymptotically given by the value of $\phi_{1}$ at its minimum, whence the result follows.

This analysis can be extended to multimode models with the form

[TABLE]

where B is a atomic operator, i.e. an operator acting on a finite dimensional Hilbert space. The bounds on their corresponding partition function are PhysRevA.8.25170

[TABLE]

We again observe that the mean field approximation is asymptotically exact as $\sum_{\mu}\omega_{\mu}/N\Delta\to 0$.

Next we extend the study on bounds PhysRevA.8.25170 to the new case of interaction terms with two photons,

[TABLE]

$B$ and $D$ are again operators acting on a finite dimensional space $\mathfrak{H}$.

Following similar procedures to those presented in PhysRevA.8.25170 , we define the cutoff coherent state as

[TABLE]

where $P_{n}$ is the projector onto the states with $n$ photons, so that $P_{n}\to I$ strongly. It is useful to define $K_{n}$ and identify some formulae in which it appears, as follows;

[TABLE]

The standard convexity inequality $\langle\psi|e^{X}|\psi\rangle\geq\langle\psi|\psi\rangle\exp(\langle\psi|X|\psi\rangle/\langle\psi|\psi\rangle)$, or its more general Schwartz inequality form $P_{n}\exp(X)P_{n}\geq\exp(P_{n}XP_{n})$ Davis:1957wj provide us with the corresponding bound with cutoffs

[TABLE]

Then we take the limit $n\to\infty$ in both sides to obtain $Z\geq\tilde{Z}$, where $\tilde{Z}$ is the partition function in mean field approximation.

On the other hand PhysRevA.8.25170 ; lieb ,

[TABLE]

where

[TABLE]

Here

[TABLE]

We then take the limit $n\to\infty$ and $\epsilon\to 0$ to obtain $Z\leq e^{\beta\omega}\tilde{Z}$. Implicit in this analysis is the requirement that $\omega^{2}/4>\|D^{\dagger}D\|$, as required for stability.

In conclusion, the bounds $\tilde{Z}\leq Z\leq e^{\beta\omega}\tilde{Z}$ do hold for the (stable) two-photon interaction case, and the reduced free energy is given asymptotically by the mean field value.

S3 2nd order SPT with symmetry: anisotropic inhomogeneous DQR model

Here we consider the anistropic inhomogeneous Dicke-quantum Rabi model,

[TABLE]

As before, we consider the limit $N{\Delta}/{\omega}\to\infty$, in which case the mean field approximation of the partition function $\tilde{Z}=\int({\mathrm{d}^{2}\alpha}/{\pi})\,\mathrm{Tr}_{spin}\left[e^{-\beta\langle\alpha|H(a^{\dagger}\to\alpha^{*},a\to\alpha)|\alpha\rangle}\right]$ is asymptotically exact. Tracing out the qubit part, we obtain $\tilde{Z}=\int\exp[-\phi_{1}(x,y)]{dxdy}/{\pi}$ with

[TABLE]

Let us change variables $u=\sqrt{{\omega}/{N\Delta}}x$ and $v=\sqrt{{\omega}/{N\Delta}}y$, where $\Delta$ is again a positive generalised mean of $\{\Delta_{i}\}$. Then the reduced free energy $f=-\ln\tilde{Z}/(\beta N\Delta)$ is well approximated by the value at its global minimum of the Landau potential

[TABLE]

where $\delta_{i}={\Delta_{i}}/{\Delta}$ and $\gamma_{i}={g_{i}}/{\sqrt{\Delta\omega}}$.

Extrema of this function satisfy

[TABLE]

Clearly the origin $(u,v)=(0,0)$ is always an extremum, and the global minimum for $\gamma_{i}$ close to [math]. If $\lambda_{i}\neq 0$, Eqs. (S20) and (S21) cannot be simultaneously satisfied with both $u$ and $v$ different from 0. Let us first assume $u=0$. In this case, $v\neq 0$ solutions to Eq. (S21) exist if

[TABLE]

The transition is continuous.

Passing now to the case $v=0$, we see that a $u\neq 0$ solution to Eq. (S20) exist if

[TABLE]

again leading to a continuous transition.

To illustrate this assertion, let us consider the homogeneous case in the limit $\beta\Delta\to\infty$. Criticality is achieved at $\gamma^{2}(1+\lambda)^{2}=2$ for $u\neq 0$ solutions and at $\gamma^{2}(1-\lambda)^{2}=2$ for $v\neq 0$. Close to criticality in the first case, $u^{2}_{\mathrm{min}}\approx 2(|\gamma(1+\lambda)|-\sqrt{2})$.

As particular cases of this Hamiltonian Eq. (S20), notice that when all $\lambda_{i}=1$, then the inhomogeneous Dicke model is recovered. The SPT presents itself as $u$ type, with the condition for superradiance

[TABLE]

The phase diagram for homogeneous Dicke model is shown in Fig. S1. Notice now that $N=1$ with $\lambda=1$ corresponds to the Rabi model.

For these special situations, and for generic values of the parameters as well, the Hamiltonian presents a $\mathbb{Z}_{2}$ symmetry, with generator $\Pi=\exp\left\{i\pi\left[a^{\dagger}a+\sum_{i=1}^{N}(1+\sigma_{iz})/2\right]\right\}$. This is reflected by the two $u\leftrightarrow-u$ and $v\leftrightarrow-v$ symmetries of $\phi(u,v)$. In fact it is a function of $u^{2}$ and $v^{2}$. It follows that the Hessian matrix at the origin is diagonal for all values of the parameters. The critical lines identified above correspond to a first sign change in the Hessian determinant at the origin, that is, when the origin in $(u,v)$ space goes from being a minimum to being a saddle point. We see here the usual behaviour of continuous transition Landau functions for $\mathbb{Z}_{2}$ symmetry.

The symmetry is enhanced if all $\lambda_{i}$ are zero, in which case the Hamiltonian Eq. (S20) corresponds to the inhomogeneous Tavis–Cummings model. Eq. (S19) reduces to

[TABLE]

This function presents rotational symmetry in the $(u,v)$ plane, which is a reflection of the $U(1)$ symmetry of the inhomogeneous Tavis–Cummings model. The Hessian matrix at the origin is proportional to the identity, with proportionality constant $2-(1/N)\sum_{i}\gamma_{i}^{2}\tanh(\beta\Delta_{i})/\delta_{i}$. Criticality is determined by this quantity being zero, or, equivalently, by the Hessian determinant at origin being zero. That is, the system is in a superradiant phase if

[TABLE]

The Landau function goes from having a global minimum at the origin to a Mexican hat like shape. Using polar coordinates, we can interpret the symmetry breaking as forcing the choice of a particular direction $\theta$, with $\theta$ the polar angle. All these directions are equivalent, so we expect the gapless excitations.

Let us now study ground states, of particular relevance in the Quantum SPT case. Namely, once the minimum point has been determined, we can produce a variational estimate also for the ground state wavefunction. For concreteness, take the case $\lambda_{i}=1$, which is the inhomogeneous DQR Hamiltonian. The effective mean field Hamiltonian reads

[TABLE]

which is a sum of independent spin terms, readily diagonalizable. The ground state in spin space of $H\left(u_{\mathrm{min}}\right)$ is

[TABLE]

where the angles are defined by

[TABLE]

The average spin energy in the ground state, in units of $\Delta$, is therefore

[TABLE]

When the symmetry is broken, the ground state is determined by a non-zero $u_{\mathrm{min}}$, resulting in the total ground state (both bosonic and spin spaces) $|gs\rangle=|\alpha_{\mathrm{min}}\rangle\otimes|\psi_{sp}\rangle$.

Outside of the generalised thermodynamic limit the mean field approach still provides us with a good approximation. The $\mathbb{Z}_{2}$ symmetry is however not broken, so the ground state, which presents odd parity, takes the approximate form

[TABLE]

where $|\psi_{sp\pm}\rangle$ is the $|\pm\rangle$ eigenstate of $\prod_{i=1}^{N}\sigma_{iz}$. Given the effectively zero overlap $\langle\alpha_{\mathrm{min}}|-\alpha_{\mathrm{min}}\rangle$, this state is normalized to be a “Schrödinger cat”-like state for finite $N$. In single qubit Rabi case, $u_{\mathrm{min}}=\frac{1}{\sqrt{2}}\sqrt{\frac{\gamma^{2}}{\gamma_{c}^{2}}-\frac{\gamma_{c}^{2}}{\gamma^{2}}}\,$, $\gamma_{c}=\sqrt{1/2}$, so

[TABLE]

Knowing the explicit ground state (in the variational approximation), we can portray many of its characteristics. For illustration, we consider the photon number distribution of the ground state of the quantum Rabi model ($N$=1) with $\Delta/\omega\gg 1$. Now, the probability of measuring $n$ photons in this state while having spin up, $P_{ne}$, will be zero for even $n$, and, similarly, the probability of odd photons with spin down will be zero. This results in

[TABLE]

We portray these expressions, together with those obtained by direct diagonalization, in Fig. S2.

S4 No symmetry

We have stressed in the analysis above that the $\mathbbm{Z}_{2}$ symmetry has guaranteed that the relevant $\phi(u)$ function is even. This has allowed us to establish the existence of a SPT and its second order character. Now we break this symmetry by adding a static bias. More concretely, we consider the Hamiltonian

[TABLE]

where, as before, $C$ will be the parameter or combination of parameters that give a macroscopic limit, while $\Omega$, $\gamma_{i}$, $\delta_{i}=\Delta_{i}/\Delta$ and $\epsilon_{i}=e_{i}/\Delta$ are finite, with $\Delta$ being some average of $\Delta_{i}$. Choosing $C=N\Delta/\omega$ we set $\Omega=1$, while $\gamma_{i}=\frac{g_{i}}{\sqrt{\Delta_{i}\omega}}$.

This model presents no parity symmetry. Nonetheless, the analysis above will be applicable in that in both the thermodynamical and the classical oscillator limit the mean field approximation of the partition function will be adequate. Following indeed the steps presented above we obtain the corresponding $\phi(u)$ function,

[TABLE]

Clearly this is no longer an even function of $u$ if not all biases $\epsilon_{i}$ are zero, whence it follows that $u=0$ is not generically an extremum. For definiteness, let us concentrate on the homogenous case with finite $N$, in the classical oscillator limit. The Landau potential $\phi(u)$ is

[TABLE]

Its derivative is

[TABLE]

One readily sees that the origin is never an extremum for generic $\epsilon$ and $\gamma$: The system always presents a macroscopic photon number.

S5 Dicke model with one and two-photon terms

A relevant question is whether models with dipolar coupling to two photons can also present superradiance. As we shall see, in general the parity symmetry we have so much relied upon will not be available, so any transition, if present, cannot be continuous. In order to assess this point, we will, as before, rewrite the Hamiltonian Eq. (5) in terms of a set of finite parameters $\Omega$, $\gamma$ and $\gamma^{\prime}$ and a macroscopicity parameter $C=N\Delta/\omega$. That is,

[TABLE]

where $\gamma$=$g/\sqrt{\Delta\omega}$, $\gamma^{\prime}=g^{\prime}/\omega$. Following the by now well established routine, we obtain the reduced free energy $F/(N\Delta)$ as the minimum of a potential, now in two variables,

[TABLE]

The stability of the system requires that the minima be located at a finite point in $(u,v)$ space. Using polar coordinates, the dominant large radius behaviour of the potential is $r^{2}\left(1-2\gamma^{\prime}\cos(2\theta)\right)$. Thus, stability requires $|\gamma^{\prime}|<1/2$, or, in the original parameters, $|g^{\prime}|<\omega/2$. This result has been obtained in different forms in the literature, and the limiting value corresponds to the spectral collapse due to the two-photon term.

This model has no symmetry for nonzero $\gamma^{\prime}$ and $\gamma$. Nonetheless, $\phi(u,v)$ is an even function of $v$ for all values of the parameters. As is only to be expected, there are some regions in parameter space for which a symmetry does exist. Namely, if $\gamma^{\prime}=0$ we recover the Dicke case with $\mathbbm{Z}_{2}$ symmetry, and this fact is reflected in $\phi$ being an even function of $u$ as well. If $\gamma=0$ we have a $\mathbbm{Z}_{4}$ symmetry arising from $u\leftrightarrow v$ exchange symmetry and parity. Now we search for extrema of $\phi(u,v)$ on computing the its derivative

[TABLE]

From Eq.(S42), one immediately sees that $\partial_{v}\phi(u,0)$ is identically zero. In the alternative, one obtains a line in the $(u,v)$ plane. On computing the first equation on this line, one obtains a linear equation for $u$. On substituting back into the line to determine critical points, one sees that there are no real solutions for $v$, as long as the stability condition is satisfied. It follows that all the extrema are located on the $v=0$ line. We thus have to analyze the extrema of

[TABLE]

Clearly $u=0$ is the global minimum when $\gamma$ and $\gamma^{\prime}\to 0$.

In the limit $\beta\Delta\to\infty$ the crucial equations become algebraic. In particular, in this limit a nonzero solution to $\phi(u)=\phi(0)$ will exist when $2\gamma^{2}+4\gamma^{\prime 2}\geq 1$, and in fact $\phi^{\prime}(u)=0$ at criticality, $2\gamma^{2}+4\gamma^{\prime 2}\geq 1$. We thus establish the existence of a first order superradiant transition with the aforementioned critical line.

We illustrate our results in the main text by diagonalizing the Hamiltonian with $N=1$ and ${\Delta}/{\omega}=200$. We portray the expected first order quantum phase transition in Fig. S3(a) and S3(b). The crossing between ground and first excited state is necessarily avoided with finite parameters. Nonetheless, the change from zero to macroscopic photon number is extremely sharp, as expected from the analytic results for the limiting case of $\Delta/\omega\to\infty$.

S6 General structure of the phase diagram

The single mode models we are considering have the general structure

[TABLE]

with $A$, $B$ and $D$ operators that act on a finite dimensional Hilbert space $\mathfrak{H}$. These operators depend continuously on the (reduced) parameters that define the phase space of interest. We have shown above that, as long as $\|D^{\dagger}D\|<\omega^{2}/4$, the mean field approximation will be valid for (generalised) thermodynamic limits, and that the phase structure will be determined by a Landau potential of the form

[TABLE]

or possibly reduced to one variable $u$. The symmetries of $\phi(u,v)$ will be inherited from those of the Hamiltonian. This function is continuous and differentiable in the situations we study. Furthermore, it is also continuous and even differentiable with respect to the parameters of the model. In all cases of interest, the dominant term for large $u^{2}+v^{2}$ will be quadratic in these variables. Additionally, in most cases the origin will be an extremum of the Landau potential. Not so, let it be noted, when a bias is introduced. Given these two considerations, namely that for large order parameter the Landau potential grows quadratically, and that the origin is an extremum for all values of the parameters, it follows that there are two alternatives: either the origin is the unique global minimum or there exists at least one non zero point $(u_{*},v_{*})$ such that $\phi(u_{*},v_{*})=\phi(0,0)$. This special point (or points) will depend on the reduced parameters.

If indeed the system is in the superradiant region of phase space, the potential will have a global minimum away from the origin. The location of the minimum will also depend on the parameters, non analytically at the phase transition.

In the cases of interest, there is a set of dipolar couplings $\vec{\gamma}$, such that the potential can be written as

[TABLE]

with $\mu$ a function of several variables, $\mu(\vec{x})$. Notice that this indeed applicable to the Dicke model with both one and two photon terms presented in Eq. (S39) of Section S5, with $\vec{\gamma}=(\gamma,\sqrt{\gamma^{\prime}})$. We readily write $\partial_{u}\phi=2u+\vec{\gamma}\cdot\nabla\mu$, where we denote with $\nabla$ the gradient with respect to the variables of function $\mu$. Since $\phi$ depends on the dipolar couplings, we can also consider its gradient in dipolar coupling space, $\nabla_{\gamma}$. Clearly we have $\nabla_{\gamma}\phi=u\nabla\mu$. Combining these two results, we conclude

[TABLE]

We see that if there is an extremum of $\phi$ away from the origin, $\vec{\gamma}\cdot\nabla_{\gamma}\phi$ will be negative at that point. Hence, the radial derivative in dipolar coupling space is negative at that point, and if the nonzero extremum is a minimum, an increase in coupling it will deepen. Thus, once in a superradiant phase, a radial increase in the dipolar couplings will move us further into the superradiant phase.

S7 Two-qubit Rabi model with Heisenberg interaction

As we have stated, in the family of models under consideration going radially in dipolar coupling space one moves deeper into the superradiant phase, once that superradiance has been achieved. It is not the case, however, that there is isotropy in dipolar couplings. To illustrate this point, we have mentioned in the main text, Eq. (6), the following model that we restate here,

[TABLE]

This model presents generically a $\mathbb{Z}_{2}$ symmetry, with operator $\Pi=\exp\left\{i\pi\left[a^{\dagger}a+\sum_{i=1}^{2}\left(1+\sigma_{iz}\right)/2\right]\right\}$.

Let us define $\Delta$ as some average of the spin energies $\Delta_{i}$, be it $\Delta=\sum_{i}\Delta_{i}/N$ if not zero, or $\Delta_{i}=\left(\sum_{i}\Delta_{i}^{2}/N\right)^{1/2}$, or some other generalised mean. We now define $\gamma_{j}=g_{j}/\sqrt{\Delta\omega}$, $\delta_{j}=\Delta_{j}/\Delta$ and $\epsilon_{\alpha}=J^{(\alpha)}/\Delta$. After following the same analysis as before for other models, we arrrive at an effective mean field spin hamiltonian

[TABLE]

We study this model in the limit $\beta\Delta\to\infty$. In this case the phase structure is controlled by the Landau potential $\phi(u)=u^{2}+\lambda(u)$, with $\lambda(u)$ the smallest eigenvalue of $h(u)$.

In this case, as we shall see, there is generically an anisotropy in the space of dipolar couplings due to the presence of the spin interactions and different spin energies, together with inhomogeneous dipolar couplings. In fact, in the limiting case $\epsilon_{\alpha}\to 0$, which we have studied earlier, it is trivial to determine that $\lambda(u)=-\sum_{j}\sqrt{\delta_{j}^{2}+4\gamma_{j}^{2}u^{2}}$. This provides us with the critical line

[TABLE]

This critical line is an ellipse, more or less elongated depending on the inhomogeneity of $\delta_{j}$. It is important to notice that the actual shape of the critical line can be computed perturbatively in this case. To understand this statement, consider the case in which the spin energies are the same and positive, $\delta_{1}=\delta_{2}=\delta>0$. The smallest eigenvalue of $h(0)$ will be $-2\delta$. Computing $\lambda(u)$ perturbatively to second order, the corresponding eigenvalue reads $-2\delta-2(\gamma_{1}^{2}+\gamma_{2}^{2})u^{2}/\delta$. This comes about because in these cases the transition is second order in all directions.

In this general model, however, the possibility exists of first order transitions. To show that this is indeed the case, we consider the homogeneous model, in which $\gamma_{1}=\gamma_{2}=\gamma$ and $\delta_{1}=\delta_{2}=\delta$. In this case we rewrite the effective spin hamiltonian as

[TABLE]

where $\mathbf{S}=(\bm{\sigma}_{1}+\bm{\sigma}_{2})/2$. In the singlet subspace the eigenvalue is $-\sum_{\alpha}\epsilon_{\alpha}$. This eigenvalue is independent of $u$. Thus, consider the situation in which there is a point $u_{*}$ at which the smallest eigenvalue of $h(u)$ in the triplet space crosses the constant eigenvalue in the singlet subspace. The Landau potential, for $u<u_{*}$, is of the form $u^{2}-\sum_{\alpha}\epsilon_{\alpha}$, while for $u>u_{*}$ there is an additional $u$ dependence, which can introduce a new minimum at finite distance from the origin. Consider, for example, a situation in which $\epsilon_{y}=0$, $\epsilon_{x}>0$ and $\epsilon_{z}>\sqrt{\epsilon_{x}^{2}+4\delta^{2}}-\epsilon_{x}>0$. In this case the smallest eigenvalue of $h(u)$ close to $u=0$ is the singlet space one. As $u$ grows, the dominant term of $h(u)$ will be the first one, and the smallest eigenvalue of $h(u)$ for large $u$ will be $-4\gamma u$. Clearly at some point the smallest eigenvalue of the triplet subspace will cross the constant singlet eigenvalue. Furthermore, from that point on the derivative of the Landau potential will behave as $u^{2}-4\gamma u$ plus correction terms, thus having a minimum with smaller value than that at the origin. This means that the transition can in this case be a first order one.

This analysis suggests rewriting the effective spin hamiltonian as

[TABLE]

In this form, the first line acts only on the triplet space, the second line is proportional to the identity, and the third line connects the singlet and the triplet spaces. At large $u>0$, the smallest eigenvalue will be $-\left(|\gamma_{1}|+|\gamma_{2}|\right)u$. One sees that the transition can be first order if this behaviour sets in before being overwhelmed by the $u^{2}$ term and the behaviour close to the origin of the smallest eigenvalue is very weakly dependent on $u$.

We have thus shown, qualitatively, that anisotropy is present in dipole coupling space, and that both first and second order transitions are feasible for some range of parameters.

Let us consider now an specific case for which a complete analysis can be easily provided; namely, identical spins ($\gamma_{1}=\gamma_{2}=\gamma$, $\delta_{1}=\delta_{2}=\delta$) with isotropic spin–spin coupling ($\epsilon_{x}=\epsilon_{y}=\epsilon_{z}=\epsilon$). Under the condition $\epsilon>|\delta|/2$ the Landau potential reads

[TABLE]

Differentiating the second expression, we see that its minimum (other than zero) would be located at $|u_{\mathrm{m}}|=2|\gamma|\sqrt{1-\delta^{2}/16\gamma^{2}}$. This is the minimum, corresponding to a superradiant phase, if $16\gamma^{4}>4\epsilon^{2}>\delta^{2}$. Once $\epsilon$ and $\delta$ are fixed (and the condition $\epsilon>|\delta|/2$ holds), the critical point for a first order transition is given by $\gamma_{c}^{2}=\epsilon/2$.

On the other hand, if $\epsilon<|\delta|/2$ the smallest eigenvalue is always $\epsilon-2\sqrt{\delta^{2}+4\gamma^{2}u^{2}}$, for all values of $u$. There will be a continuous transition with critical dipolar coupling $\gamma_{c}^{2}=|\delta|/4$. The first and second order critical surfaces in the three dimensional space of $\gamma$, $\epsilon$ and $\delta$ have the common line $\gamma^{2}=\epsilon/2=|\delta|/4$.

In order to go forward and provide a more quantitative assessment, it is convenient to remember that the Landau potential will not have a unique global minimum at the origin $u=0$ if there exists a finite, nonzero value of $u$ such that $\phi(u)=\phi(0)$ holds. In the case at hand, this corresponds to the equation $\lambda(u)=\lambda(0)-u^{2}$. Since the determination of the smallest eigenvalue for generic parameters and order parameter $u$ is rather involved as a first step, and the computation of the derivative of the Landau potential hinges on this first computation, we sidestep the process in some cases by inserting the equation $\lambda(u)=\lambda(0)-u^{2}$ in the secular equation $\mathrm{Det}\left[h(u)-\lambda(u)\mathbbm{1}\right]$. This technique provides us with an analytical tool that can be also implemented numerically. We portray a number of cases in Figures S4, S5 and S6. To be more explicit, the substitution we have mentioned leads us to a cubic equation in the variable $u^{2}$. We are looking for real positive roots of the corresponding cubic polynomial. In particular we are looking for the existence of a real positive double root in order to identify criticality. This entails the discriminant of the cubic being zero. One also must look at the value of the polynomial at $u^{2}=0$, to assess whether the number of positive zeros (with multiplicities) is even or odd. Given some parameters for the spin energies and spin–spin couplings, this provides us with analytic equations for the dipolar couplings, and the result can be checked by numerical diagonalization.

This analysis provides us with a specific check on the general phase structure on the space of dipolar couplings presented above in Section S6. Indeed, as shown in the main text, the region close to the origin of dipolar couplings is in normal phase, and one goes radially into a superradiant one in the cases portrayed there.

S8 Multi-mode Dicke-QR model

According to our analysis, the partition function can be written as

[TABLE]

After tracing out the qubit and $y$ part, we obtain

[TABLE]

Making the transformation $x_{\nu}\rightarrow\sqrt{\frac{N\Delta}{\omega_{\nu}}}u_{\nu}$, we obtain $Z=\mathbf{N^{\prime}}\int\exp\{-\beta N\Delta\phi(u_{\nu})\}du_{\nu}$, with

[TABLE]

where $\gamma_{\nu}=\frac{g_{\nu}}{\sqrt{\omega_{\nu}\Delta}}$, $\delta_{i}=\Delta_{i}/\Delta$, $\Delta$ is the average of $\Delta_{i}$. For finite temperature $\beta\gg\frac{1}{N\Delta}$, $\beta N\Delta$ is infinite, so according to Laplace’s method, the reduced free energy $\frac{F}{N\Delta}$ is determined by the global minimum of $\phi(u_{m})$, $m=1,2,\ldots,M$, which should satisfy

[TABLE]

so $u_{\nu}=\frac{\gamma_{\nu}u_{1}}{\gamma_{1}}$. Substituting this into Eq. (S56) for $m=1$, and defining $\gamma^{2}=\sum_{\nu}\gamma_{\nu}^{2}$, we obtain

[TABLE]

One solution is $u_{1}=0$ and hence $u_{\nu}=0$, which is clearly global minimum when $g\sim 0$. $\gamma$-dependent solution exists for

[TABLE]

At the critical $\gamma_{c}$, $u_{\nu}(\gamma_{c})=0$, so $f(u_{\nu}(\gamma_{c}))=f(0)$, but as $\gamma_{\nu}$ increases, all nonzero $\phi(u_{\nu})$ are decreasing, so the global minimum should be this $\gamma$-dependent extreme. For identical qubit $\Delta_{i}=\Delta$, it is easy to obtain the critical $\gamma_{c}^{2}=1/2$ and $\frac{\omega_{\nu}a^{\dagger}a}{N\Delta}=u_{\nu}^{2}=\frac{\gamma_{\nu}^{2}}{4}(\frac{1}{\gamma_{c}^{4}}-\frac{1}{\gamma^{4}})$ for SRP in classical oscillator or zero temperature case ($\beta\Delta\to\infty$) from Eq. (S57).