# Cavity quantum electrodynamics in the non-perturbative regime

Daniele De Bernardis, Tuomas Jaako, and Peter Rabl  Vienna Center for Quantum
Science and Technology, Atominstitut, TU Wien, 1040 Vienna, Austria

###### Abstract

We study a generic cavity-QED system where a set of (artificial) two-level
dipoles is coupled to the electric field of a single-mode $LC$ resonator. This
setup is used to derive a minimal quantum mechanical model for cavity QED,
which accounts for both dipole-field and direct dipole-dipole interactions.
The model is applicable for arbitrary coupling strengths and allows us to
extend the usual Dicke model into the non-perturbative regime of QED, where
the dipole-field interaction can be associated with an effective finestructure
constant of order unity. In this regime, we identify three distinct classes of
normal, superradiant and subradiant vacuum states and discuss their
characteristic properties and the transitions between them. Our findings
reconcile many of the previous, often contradictory predictions in this field
and establish a common theoretical framework to describe ultrastrong coupling
phenomena in a diverse range of cavity-QED platforms.

##  I Introduction

Quantum electrodynamics (QED) is the fundamental theory of charges and
electromagnetic fields, which in its low-energy limit describes the physics of
photons interacting with atoms, molecules and solids. Cavity QED HarocheCQED
is a minimal framework within which such light-matter interactions are studied
at the quantum level in terms of two-level emitters coupled to a single
radiation mode. A hallmark of cavity QED is the strong coupling between single
atoms and single photons, which has been the subject of many theoretical and
experimental works in this field. These strong interactions between excited
atomic and photonic states are, however, still perturbative in the sense that
the coupling is much smaller than the absolute atomic or photonic energy
scales involved. Indeed, it follows from basic geometric considerations that
the coupling strength $g$ between an elementary electric dipole and a cavity
mode of frequency $\omega_{c}$ is limited to HarocheCQED ; Devoret2007 ;
Comment

$\frac{g}{\omega_{c}}\lesssim\sqrt{2\pi\alpha_{\rm fs}},$ (1)

where $\alpha_{\rm fs}\simeq 1/137$ is the finestructure constant. As a
consequence, the vacuum of (cavity) QED is to a good approximation represented
by the trivial state with all atoms in their ground state and no photons. This
is in stark contrast to the theory of quantum chromodynamics, where much more
complex ground states of strongly interacting quarks and gluons arise.

The interest in the physics of light-matter interactions beyond this ‘weak-
coupling’ regime dates back to the early days of cavity QED and is
traditionally closely connected to the Dicke model Dicke1954 ; Brandes2005 ;
Garraway2011  describing the coupling of $N$ two-level atoms to a single
optical mode. For a sufficiently strong _collective_ coupling, $G=g\sqrt{N}$,
the ground state of this model undergoes a quantum phase transition from the
normal vacuum into a so-called superradiant phase, where the atoms
spontaneously polarize and the field acquires a non-vanishing expectation
value Hepp1973 ; Wang1973 ; Hioe1973 . Over the past decades the existence of
such a cavity-induced instability has been subject of many controversial
debates. Most notably, it has been argued Rzazewski1975  that the superradiant
phase does not occur in more realistic models when the usually neglected
diamagnetic “$A^{2}$-term” is correctly taken into account. This famous no-go-
theorem has both been confirmed as well as rejected by many subsequent studies
of various cavity QED Kudenko1975 ; Emelanov1976 ; Knight1978 ; Bialynicki1979
; Bowden1979 ; Yamanoi1979 ; Keeling2007 ; Hagenmuller2010 ; Todorov2012 ;
VukicsPRA2012 ; DeLiberato2014 ; Vukics2014 ; Bamba2014 ; Vukics2015 ;
Tufarelli2015 ; Griesser2016  and analogous circuit QED CiutiNatComm2010 ;
NatafPRL2010 ; ViehmannPRL2011 ; Leib2012 ; Jaako2016 ; Bamba2016 ; Bamba2017
setups, but despite its fundamental relevance, this matter is still not fully
resolved.

More recently, the development of various solid-state cavity QED platforms has
led to a growing number of experimental activities related to what is now
quite generally called the ultrastrong coupling (USC) regime Ciuti2005  of
light-matter interactions. By using, for example, organic materials
Schwartz2011 ; KenaCohen2013 ; Mazzeo2014 , intersubband transitions
Todorov2010 ; Geiser2012 ; Benz2013 ; Dietze2013 ; Vasanelli2016 , or 2D
electron gases Scalari2012 ; Maissen2014 ; Zhang2016 ; Keller2017 , the
collective coupling of such dense dipolar ensembles to optical or THz modes
can reach a considerable fraction of the bare photon frequency. In parallel,
it has been demonstrated in the context of circuit QED Wallraff2004 ;
Blais2004  that artificial atoms, like superconducting qubits YouNature2011
or quantum dots Mi2017 ; Stockklauser2017 ; Bruhat2016 ; Cottet2017 , can be
coupled very efficiently to microwave resonators, in which case the USC regime
becomes accessible even at the single-qubit level Forndiaz2017 ; Yoshihara2017
; Niemczyk2010 ; Forndiaz2010 ; Baust2016 ; Chen2017 ; Bosman2017 . In light
of these experimental developments and potential applications ranging from
USC-assisted chemical reactions Hutchison2012 ; Wang2014 ; Galego2015 ;
Herrera2016 ; Flick2016  to ultra-fast superconducting quantum information
processing schemes Romero2012 ; Armata2017 , a refined understanding of the
basic principles of USC cavity QED on the single-, few- and many-particle
level becomes of uttermost importance.

In this work we analyze a generic cavity QED setup where multiple two-level
dipoles are coupled to a single electromagnetic mode of a lumped-element $LC$
resonator. In this setup the limit on the interaction strength stated above
can be overcome by coupling (artificial) dipoles to the electric field of a
tailored circuit mode with an impedance much higher than that of free space
Devoret2007 . In view of Eq. (1), one can then associate with this system an
effective finestructure constant of order unity, meaning that already for a
single dipole a non-perturbative treatment of electromagnetic interactions
must be taken into account. The purpose of this study is, first of all, to
derive a minimal consistent model for cavity QED, which is applicable in this
non-perturbative regime Comment2 , and second, to evaluate and describe the
resulting vacuum states under various conditions. In contrast to most previous
studies on this subject, we here focus explicitly on the long-wavelength and
low-frequency regime to avoid many of the complications related to the
quantization of the full electromagnetic field PhotonsAndAtoms ; ScheelBook .
This approach still captures correctly the relevant low-energy physics and
allows us to rigorously separate the collective coupling to a single dynamical
field mode from all direct dipole-dipole interactions. Thereby, most of the
ambiguities about the existence or non-existence of superradiant instabilities
can be resolved and explained in terms of basic electrostatic considerations.
Our analysis also addresses several other subtle issues, like the breakdown of
gauge invariance, which results in a unique extension of the Dicke model into
the USC regime.

From the analysis of the ground states of this model we identify three
distinct classes of normal, superradiant and subradiant vacuum states, which
arise from the competition between direct dipole-dipole and cavity-mediated
interactions. Our study first of all shows that a superradiant phase
transition (in the conventional sense) can exist for very specific geometries,
but must be understood as a ferroelectric instability Keeling2007 ;
Griesser2016 ; Emelanov1976 , which is essentially unaffected by the coupling
to the resonator mode. Nevertheless, this transition is still associated with
a characteristic kink in the vacuum fluctuations of the gauge-invariant
voltage and flux degrees of freedom. In the non-perturbative regime
significant corrections from this classical picture arise due to a
hybridization of individual dipoles and photons. Most importantly, in this
regime the cavity induces a collective anti-ferroelectric interaction, which
favors subradiant ground states where the dipoles tend to anti-align and
decouple from the field mode Jaako2016 . In this regime also a new transition
between superradiant and subradiant ground states becomes possible. These
preliminary findings already show that for very strong interactions the
physics of cavity QED can differ significantly from the usual picture conveyed
by discussions of Dicke or Hopfield-type Hopfield  models and that many
surprising aspects of USC physics are still unexplored.

The remainder of the paper is structured as follows. After introducing in Sec.
II the setup and the quantities of interest, we first discuss in Sec. III the
polaritonic eigenmodes and instabilities of classical systems of dipoles in a
cavity. In Sec. IV we then derive a minimal quantum mechanical model for this
system, which after some further simplifications is used in Sec. V to
investigate the different ground states of cavity QED. We conclude our work in
Sec. VII by connecting the findings of this work to different experimental
platforms.

##  II Cavity QED: A toy model

Figure 1: (a) Sketch of the cavity QED setup considered in this work. (b)
Different effective potentials $V(\xi)$ for the dipole variable $\xi$ are used
to model either harmonic or two-level dipoles of frequency $\omega_{0}$. (c)
Illustration of the two different contribution to the total charge
$Q=Q_{U}+Q_{\rm in}$ on the upper capacitor plate. See text for more details.

We consider a setup as shown in Fig. 1(a), where $N$ dipoles are coupled to
the electric field of a lumped-element $LC$ resonator. The resonator has a
bare resonance frequency $\omega_{c}=\sqrt{1/LC}$, where $C$ is the
capacitance and $L$ the inductance of the circuit. This frequency is far
separated from all higher order electromagnetic resonances such that the $LC$
resonator is well-described by a single harmonic oscillator mode. The dipoles
are assumed to be fixed at positions $\vec{r}_{i}=(x_{i},y_{i},z_{i})$ and
formed by a pair of charges $+q$ and $-q$, which are displaced by an amount
$\xi_{i}$ in the direction perpendicular to the plates. The dynamics of each
dipole is modelled as an (effective) particle of mass $m$ moving in a
potential $V(\xi)$, as indicated in Fig. 1(b). This allows us to treat both
harmonically bound dipoles as well as two-level systems by changing from a
quadratic to a double-well potential. For all the following derivations it is
assumed that the dipole approximation is valid and that radiative effects as
well as magnetic interactions can be neglected.

The dynamics of the $LC$ resonator is governed by the circuit relations
$\dot{\Phi}=U$ and $\dot{Q}=-\Phi/L$, where $U$ is the voltage drop across the
capacitor, $Q$ is the total charge on the upper capacitor plate and $\Phi$ is
the magnetic flux through the inductor. In the following we write
$Q=Q_{U}+Q_{\rm in}$, where $Q_{U}=CU$ is the charge in the absence of the
dipoles and $Q_{\rm in}=\int_{A}dxdy\,\sigma_{\rm in}(x,y)$ is the total
charge induced by the dipoles when $U=0$ [cf. Fig. 1(c)]. The induced surface
charge density, $\sigma_{\rm in}(x,y)$, depends on the exact distribution of
dipoles and can vary strongly across the capacitor plate of total area $A$.
With these definitions we obtain the equation of motion for the flux variable
$\Phi$,

$C\ddot{\Phi}+\frac{\Phi}{L}=-\dot{Q}_{\rm
in}\simeq\frac{q}{d}\sum_{i}\dot{\xi}_{i}.$ (2)

In the last step we have used the fact that sufficiently far away from the
edges of the capacitor the total surface charge induced by a single dipole is
$-q\xi_{i}/d$, where $d$ is the distance between the plates (see App. A). This
approximation is not essential, but results in a convenient homogeneous
dipole-resonator interaction.

Based on the assumptions stated above, the equations of motion for the dipole
variables $\xi_{i}$ are
$m\ddot{\xi}_{i}+V_{i}^{\prime}(\xi_{i})=qE(\vec{r}_{i})$, where
$E(\vec{r}_{i})$ is the total electric field at the position of the $i$-th
dipole. We decompose this field into two parts,

$qE(\vec{r}_{i})=-\frac{q}{d}\dot{\Phi}-m\omega^{2}_{\rm
p}\sum_{i,j}\mathcal{D}_{ij}\xi_{j},$ (3)

where we introduced the plasma frequency,

$\omega_{\rm p}^{2}=\frac{q^{2}}{\varepsilon_{0}mr_{0}^{3}},$ (4)

as the characteristic frequency scale related to the interaction between two
neighboring dipoles separated by a distance $r_{0}$. The dimensionless
coupling parameters $\mathcal{D}_{ij}\sim\mathcal{O}(1)$ account for the exact
spatial dependence of dipole-dipole interactions. In free space we would
simply obtain

$\mathcal{D}_{ij}=\frac{r_{0}^{3}}{4\pi}\frac{|\vec{r}_{ij}|^{2}-3(\vec{r}_{ij}\cdot\vec{e}_{z})^{2}}{|\vec{r}_{ij}|^{5}},$
(5)

where $\vec{r}_{ij}=\vec{r}_{i}-\vec{r}_{j}$, but the capacitor plates can
strongly modify this dependence due to the presence of additional image
charges Perram1996 . The numerical evaluation of the $\mathcal{D}_{ij}$ in
this confined geometry is detailed in App. A. Note that each dipole also
interacts with its own image charges and $\mathcal{D}_{ii}\neq 0$. In the
following we absorb this self-interaction into a redefinition of the
potential, i.e., $V_{i}(\xi_{i})+m\omega_{\rm
p}^{2}\mathcal{D}_{ii}\xi_{i}^{2}/2\rightarrow V^{\prime}_{i}(\xi_{i})\simeq
V(\xi_{i})$, which, for the sake of simplicity, is assumed to be approximately
the same for all dipoles.

In summary, we obtain the following equations of motion for the dynamical
variables $\xi_{i}$,

$m\ddot{\xi}_{i}+V^{\prime}(\xi_{i})+m\omega_{\rm p}^{2}\sum_{j\neq
i}\mathcal{D}_{ij}\xi_{j}=-\frac{q}{d}\dot{\Phi}.$ (6)

Together with Eq. (2), this set of equations specifies a minimal model for a
cavity QED system consisting of multiple electric dipoles coupled to a single
electromagnetic mode.

##  III Polaritons, instabilities and geometry

Before we proceed with the quantization of our model, it is instructive to
consider first a few basic properties of this system in the limit of a large
number of harmonically bound dipoles, i.e., $V(\xi)=m\omega_{0}^{2}\xi^{2}/2$.
For a sufficiently homogeneous system the cavity will couple primarily to the
collective variable $\mathcal{Z}=\sum_{i}\xi_{i}/\sqrt{N}$, where all dipoles
oscillate in phase. By ignoring for now the weak admixing of other excitation
modes due to dipole-dipole interactions, we arrive at a reduced set of two
coupled oscillator equations

$\displaystyle\ddot{\mathcal{Z}}+\left(\omega_{0}^{2}+\eta\omega_{\rm
p}^{2}\right)\mathcal{Z}$ $\displaystyle=$
$\displaystyle-\frac{q}{dm}\sqrt{N}\dot{\Phi},$ (7)
$\displaystyle\ddot{\Phi}+\omega_{c}^{2}\Phi$ $\displaystyle=$
$\displaystyle\frac{qd}{V\varepsilon_{0}}\sqrt{N}\dot{\mathcal{Z}}.$ (8)

Here we have assumed a parallel plate capacitor with volume $V=Ad$ and
capacitance $C=\varepsilon_{0}A/d$ and introduced the dimensionless parameter

$\eta=\frac{1}{N}\sum_{i\neq j}\mathcal{D}_{ij}.$ (9)

This geometrical constant captures the average influence of dipole-dipole
interactions in a homogeneously polarized sample and is closely related (but
not identical) to the usual depolarization factor of dielectric bodies
LandauECM . Its value depends on the lattice configuration, the shape of the
dipole ensemble, and the metallic boundaries, but for a fixed minimal distance
$r_{0}$, it does not scale with the number of dipoles.

###  III.1 Polaritons and instabilities

From Eqs. (7) and (8) we readily obtain two polaritonic eigenmodes with
frequencies (see App. B)

$\Omega_{\pm}^{2}=\frac{\omega_{d}^{2}+\omega_{c}^{2}+\nu\omega_{\rm
p}^{2}\pm\sqrt{(\omega_{d}^{2}+\omega_{c}^{2}+\nu\omega_{\rm
p}^{2})^{2}-4\omega_{d}^{2}\omega_{c}^{2}}}{2},$ (10)

where $\omega_{d}=\sqrt{\omega_{0}^{2}+\eta\omega_{\rm p}^{2}}$ denotes the
bare oscillation frequency of the interacting ensemble of dipoles. In Eq. (10)
we have used the identity
$q^{2}N/(mCd^{2})=q^{2}N/(\varepsilon_{0}Vm)=\nu\omega_{\rm p}^{2}$ to express
the collective dipole-field coupling in terms of the plasma frequency and the
filling factor $\nu=Nr_{0}^{3}/V$.

Figure 2: The spectrum of the two bright polariton branches is plotted as a
function of $\omega_{\rm p}$ and for $\omega_{c}=\omega_{0}$. In (a) a
positive value of $\eta\approx 0.3$ and in (b) a negative value of
$\eta\approx-0.2$ has been assumed. In both plots the orange (lower) and the
dark blue (upper) lines represent the spectrum obtained from Eq. (10), while
the shaded area indicates the range of frequencies of all other dark polariton
modes obtained from the numerical solution of the full eigenvalue problem (see
App. B). (c) Sketch of an ensemble of $N=3\times N_{x}^{2}$ dipoles, which are
arranged in three layers on a square lattice with spacing $r_{0}$ and placed
between two capacitor plates. For this configuration the resulting value of
$\eta$ is plotted in (d) for varying $d\geq h$ and different $N_{x}$. The
values of $\eta$ and the full coupling matrix $\mathcal{D}_{ij}$ used in the
calculations of the polariton spectra in (a) and (b) have been obtained for
the case $N_{x}=10$ and the values of $h/d=\nu\approx 0.2$ and $h/d=\nu\approx
0.9$, respectively.

Figure 2 shows examples of polaritonic spectra plotted as a function of
increasing plasma frequency, i.e., increasing density of dipoles, and for two
different values of $\eta$. For resonant interactions,
$\omega_{0}\approx\omega_{c}$, and for small values of $\omega_{\rm p}$, we
observe in both cases the expected normal mode splitting,
$\Delta\Omega=\Omega_{+}-\Omega_{-}\approx\sqrt{\nu}\omega_{\rm p}$. For
larger $\omega_{\rm p}$ and $\eta>0$ the lower branch approaches a finite
value $\Omega_{-}\approx\omega_{c}\sqrt{1-\nu/(\nu+\eta)}$ and remains stable
for all parameters. This behavior is well-known from the study of various
solid-state cavity QED systems Ciuti2005 ; Todorov2010 ; Todorov2012 ;
KenaCohen2013 ; Geiser2012 ; Maissen2014 ; Zhang2016 ; Hagenmuller2010 , where
the regime $\Delta\Omega\sim\omega_{c}$ is experimentally accessible. In these
systems, the observed deviation from a linearly increasing mode splitting is
usually derived from the Hopfield model Hopfield , where the $A^{2}$-term is
taken into account. In the opposite case, $\eta<0$, i.e., when dipole-dipole
interactions are on average attractive, there exists a critical density or
critical plasma frequency, $\omega_{\rm p}^{c}=\omega_{0}/\sqrt{\eta}$, at
which the frequency of the lower polariton mode vanishes. Beyond this point
the eigenfrequency $\Omega_{-}$ is imaginary, which means that any excitation
of this mode will be exponentially amplified. Therefore, the linear system
becomes unstable and a more accurate description of the dipoles must be taken
into account. As shown in more detail in App. B, the critical density at which
this instability occurs is determined by a vanishing frequency of the
interacting dipole ensemble, i.e. $\omega_{d}=0$, and is independent of the
cavity frequency. An experimental signature consistent with such an
instability has recently been reported for a 2D hole gas coupled to a THz
resonator Keller2017 .

In typical cavity QED experiments the excitation spectrum is inferred from the
cavity output field and therefore only the ‘bright’ polariton modes, which are
described by Eq. (10) and involve a large photonic component, are observable.
However, there also exist $N-1$ unobservable, i.e., ‘dark’ excitation modes of
the dipole ensemble, which due to their spatial profile are almost decoupled
from the cavity field (see App. B). In the examples shown in Fig. 2(a) and (b)
the frequency range of these additional modes is indicated by the shaded area.
We see that even for $\eta>0$ some of these modes become unstable at high
enough densities. Thus, the stability of the experimentally observable bright
polarition modes does not necessarily imply the linear stability of the system
as a whole.

###  III.2 Geometrical considerations

The shape-dependence of macroscopic thermodynamic properties is a peculiarity
of systems interacting via long-range dipole-dipole interactions and is well-
known from the study of magnetic or ferroelectric systems. In the limit $N\gg
1$ approximate expressions for $\eta$ can be derived, for example, by treating
the dipoles as a continuous medium with polarization density $P(\vec{r})$ and
solving for the macroscopic electric field $E_{M}(\vec{r})$. The local field
$E(\vec{r})$ can then be obtained from the relation
$E(\vec{r})=E_{M}(\vec{r})+E_{\rm near}(\vec{r})-E_{P}(\vec{r})$, where
$E_{\rm near}$ is the exact field and $E_{P}=-P/(3\epsilon_{0})$ the average
field from neighboring dipoles inside a small Lorentz sphere centered around
$\vec{r}$ LandauECM . In free space and for dipoles placed on a regular cubic
lattice, where $E_{\rm near}\approx 0$, one obtains $\eta\approx 2/3$ for a
disc-shaped ensemble, $\eta\approx 0$ for a spherical ensemble and
$\eta\approx-1/3$ for an elongated, cigar-like configuration LandauECM .
Importantly, these values are modified in the presence of the capacitor
Takae2013 ; Perram1996 , as illustrated in Fig. 2(c) and (d) for the case of a
flat layer of dipoles placed between two metallic plates. For this geometry we
obtain

$\eta\approx\frac{2}{3}-\frac{h}{d},$ (11)

where $h$ is the thickness of the dipole layer. Therefore, the presence of the
metallic boundaries can have a substantial effect and bring the system from a
stable to an unstable configuration as the distance between the plates is
decreased Bratkovsky2008 ; Levanyuk2016 . For a single layer of dipoles placed
on a triangular lattice and $d\gg r_{0}$ we obtain $\eta\approx 0.88$
Nijboer1958 ; Dell'Anna2016 , while the minimal possible value of $\eta=-{\rm
Zeta}(3)/\pi\approx-0.38$ is obtained for a line of dipoles placed on top of
each other.

In view of Eq. (11) it is important to keep in mind that our definition of the
potential $V(\xi)$ includes, apart from the external confining potential, also
the energy that it takes to separate the charges $+q$ and $-q$. When applying
the current analysis to the case of a free electron gas, the limit
$\omega_{0}^{2}\rightarrow\omega_{\rm p}^{2}/3$ must be taken to retain this
local field contribution. In this limit we recover the usual plasma
oscillations, $\omega_{d}\approx\omega_{\rm p}$, for $\nu\rightarrow 0$ and
$\omega_{d}\rightarrow 0^{+}$ for $\nu\rightarrow 1$.

###  III.3 Discussion

From the basics properties of polaritonic systems discussed in this section we
can already make the following important observations. (i) Both the collective
dipole-field coupling, $\Delta\Omega\sim\omega_{\rm p}$, as well as the
strength of direct dipole-dipole interactions, $\sim\omega_{\rm p}^{2}$, scale
with the density and cannot be treated as independent effects. In particular,
in the USC regime, where $\omega_{\rm p}\sim\omega_{c},\omega_{0}$, the effect
of dipole-dipole interactions plays a dominant role and must be fully taken
into account. (ii) A cavity QED system of dipoles coupled to a single
electromagnetic mode can exhibit an instability. This instability is induced
by dipole-dipole interactions and therefore depends on details like the shape
of the ensemble or the lattice configuration. This explains, why many no-go-
and counter-no-go-theorems for superradiant phase transitions, which either
completely omit dipole-dipole interactions or do not treat them in all detail,
come to very different conclusions. (iii) Most importantly, if an instability
exists, it is _solely_ induced by dipole-dipole interactions and not
influenced by the frequency or other properties of the resonator mode. This
observation contradicts the usual picture conveyed by discussions of the Dicke
model, where the transition into the superradiant phase is commonly
misinterpreted as being induced by the coupling to a dynamical field mode. Of
course, adding the metallic plates in the first place can still substantially
modify the properties of the confined system of dipoles compared to its
counterpart in free space.

##  IV Cavity QED Hamiltonian

Our goal is now to derive a minimal quantum mechanical model for the cavity
QED system described in Sec. II, which is also applicable for highly non-
linear dipolar systems and for arbitrary coupling strengths. As a starting
point for this derivation we consider the Lagrangian
$\mathcal{L}\equiv\mathcal{L}(\Phi,\dot{\Phi},\\{\xi_{i}\\},\\{\dot{\xi}_{i}\\})$
of the form

$\begin{split}\mathcal{L}=&C\frac{\dot{\Phi}^{2}}{2}-\frac{\Phi^{2}}{2L}+\dot{\Phi}Q_{\rm
in}\\\
&+\sum_{i}\left[\frac{m}{2}\dot{\xi}_{i}^{2}-V(\xi_{i})\right]-\frac{m\omega_{\rm
p}^{2}}{2}\sum_{i,j}\mathcal{D}_{i\neq j}\xi_{i}\xi_{j},\end{split}$ (12)

from which the equations of motion (2) and (6) can be derived. For this
Lagrangian, the resulting canonical momenta are

$\Pi=\frac{\partial\mathcal{L}}{\partial\dot{\Phi}}=C\dot{\Phi}+Q_{\rm
in},\qquad\pi_{i}=\frac{\partial\mathcal{L}}{\partial\dot{\xi}_{i}}=m\dot{\xi}_{i}=p_{i},$
(13)

and correspond to the total charge $Q$ on the capacitor plate and the kinetic
momenta, respectively. By following the usual quantization procedure we obtain
the Hamilton operator

$H=\frac{(Q-Q_{\rm
in})^{2}}{2C}+\frac{\Phi^{2}}{2L}+\sum_{i}H_{d}^{i}+\frac{m\omega_{\rm
p}^{2}}{2}\sum_{i,j}\mathcal{D}_{i\neq j}\xi_{i}\xi_{j},$ (14)

where $H_{d}^{i}=p_{i}^{2}/(2m)+V(\xi_{i})$ and $\Phi$, $Q$, $\xi_{i}$ and
$p_{i}$ are now operators obeying the commutation relations
$[\Phi,Q]=[\xi_{i},p_{j}]=i\hbar\delta_{ij}$. Using $Q_{\rm
in}\simeq-q\sum_{i}\xi_{i}/d$, Hamiltonian (14) can further be expanded in
terms of the operators $\xi_{i}$,

$\begin{split}H=&\frac{Q^{2}}{2C}+\frac{\Phi^{2}}{2L}+\frac{q}{Cd}Q\sum_{i}\xi_{i}+\frac{q^{2}}{2Cd^{2}}\sum_{i,j}\xi_{i}\xi_{j}\\\
&+\sum_{i}\left[\frac{p_{i}^{2}}{2m}+V(\xi_{i})\right]+\frac{m\omega_{\rm
p}^{2}}{2}\sum_{i\neq j}\mathcal{D}_{ij}\xi_{i}\xi_{j}.\end{split}$ (15)

This result represents the full Hamilton operator for the model cavity QED
system considered in this work.

Equation (15) shows that apart from the expected collective coupling of all
dipoles to the ‘charge’ of the $LC$ resonator, there are two additional
dipole-dipole interaction terms $\sim\xi_{i}\xi_{j}$. Since, by construction
of our model, the term $\sim\mathcal{D}_{ij}$ already accounts for all direct
interactions between the dipoles, the additional presence of the last term in
the first line of Eq. (15) is very counterintuitive. However, as can be seen
from Eq. (14), this term simply arises from expressing the electrostatic
energy contribution, $CU^{2}/2$, in terms of the canonical charge $Q$. It is
thus merely an artifact of our choice of variables and should not be
associated with a physical interaction. The inclusion of this term is
nevertheless crucial to recover the correct equations of motion (6) from the
relation $m\ddot{\xi}_{i}=i/\hbar[H,p_{i}]$. This subtle difference between
apparent interaction terms in the Hamiltonian and real physical couplings is a
common source of confusion in the interpretation of cavity QED models.

###  IV.1 The $P^{2}$-term

In our model, the distribution of point-like dipoles corresponds to a
polarization density

$\vec{P}(\vec{r})=q\vec{e}_{z}\sum_{i}\xi_{i}\delta(\vec{r}_{i}-\vec{r}).$
(16)

For a sufficiently dense and homogeneous ensemble, where
$\vec{P}(\vec{r})\simeq\vec{e}_{z}P$, we can identify $Q_{\rm
in}\epsilon_{0}/(Cd)\simeq-P$ with the polarization density and
$Q\epsilon_{0}/(Cd)\simeq\epsilon_{0}U/d-P=-D$ with the displacement field
$\vec{D}(\vec{r})\simeq\vec{e}_{z}D$. With these identifications, Hamiltonian
(15) can be directly related to the Hopfield model expressed in the electric
dipole gauge Todorov2012 ; Bamba2014 ; PhotonsAndAtoms ; ScheelBook ,

$\begin{split}H_{\rm HM}=&H_{\rm matter}+\int
d^{3}r\,\frac{[\vec{D}(\vec{r})-\vec{P}(\vec{r})]^{2}}{2\epsilon_{0}}+\frac{\vec{B}^{2}(\vec{r})}{2\mu_{0}}\\\
=&H_{\rm matter}+\int
d^{3}r\,\frac{\vec{D}^{2}(\vec{r})}{2\epsilon_{0}}+\frac{\vec{B}^{2}(\vec{r})}{2\mu_{0}}\\\
&-\frac{1}{\epsilon_{0}}\int
d^{3}r\,\vec{D}(\vec{r})\cdot\vec{P}(\vec{r})+\frac{1}{2\epsilon_{0}}\int
d^{3}r\vec{P}^{2}(\vec{r}).\end{split}$ (17)

Here $\vec{B}(\vec{r})$ is the magnetic field and $H_{\rm
matter}=\sum_{i}H_{d}^{i}$ is the Hamiltonian for the matter part. Therefore,
the above-discussed $Q_{\rm in}^{2}$-contribution plays an equivalent role as
the polarization self-interaction or “$P^{2}$-term”, which appears in the
description of macroscopic polarizable media Todorov2010 ; Todorov2012 . It
should be emphasized though that for a discrete polarization density as in Eq.
(16), this polarization self-interaction term results in purely local
interactions PhotonsAndAtoms ; VukicsPRA2012 ; Vukics2014

$\int d^{3}r\,\vec{P}^{2}(\vec{r})\,\,\longrightarrow\,\,\sum_{i}\xi_{i}^{2}.$
(18)

The apparent discrepancy between such a local $P^{2}$-term and the non-local
coupling derived in Eq. (15) can be resolved by taking into account that
Hamiltonian $H_{\rm HM}$ still contains the coupling of the dipoles to _all_
electromagnetic modes. As illustrated in App. C for a basic geometry, the
coupling to these other high-frequency modes introduces effective
interactions, which restore the correct non-local $P^{2}$\- and direct dipole-
dipole interaction terms. In other words, starting from the full model $H_{\rm
HM}$ in the electric dipole gauge, a single-mode approximation is—independent
of the frequency separation—not permitted and various approximate treatments
of the $P^{2}$-term lead to very different physical predictions Todorov2012 ;
VukicsPRA2012 ; Vukics2014 ; Bamba2014 . Our derivation avoids such
complications by including the correct dipole-field and dipole-dipole
interactions before passing to a quantum description.

###  IV.2 Two-level-approximation

Of primary interest in the field of cavity QED is the study of nonlinear
quantum phenomena, which arise from the coupling of the harmonic field mode to
nonlinear matter, in the simplest case represented by two-level dipoles. In
our model we can describe this scenario by considering for each dipole a
double-well potential with eigenstates $|\psi_{n}\rangle$ of energy
$\hbar\omega_{n}$ [cf. Fig. 1(b)]. For an appropriate choice of parameters the
two lowest tunnel-coupled states $\ket{\downarrow}\equiv|\psi_{1}\rangle$ and
$\ket{\uparrow}\equiv|\psi_{2}\rangle$ are energetically well-separated from
all higher excited states and the dynamics of the dipoles can be restricted to
this two-level subspace. Under such conditions we can approximate

$H_{d}^{i}\approx\frac{\hbar\omega_{0}}{2}\sigma_{z}^{i},\qquad\xi_{i}\approx\frac{\xi_{0}}{2}\sigma_{x}^{i},$
(19)

where the $\sigma_{k}$ are the usual Pauli operators,
$\omega_{0}=\omega_{2}-\omega_{1}$ is the transition frequency between the two
lowest states and $\xi_{0}=2\langle\downarrow\\!|\xi_{i}|\\!\uparrow\rangle$
is the separation between the wells. According to Eq. (19), the definition of
$\omega_{0}$ does not include a small renormalization of the potential from
the additional term $\sim\xi_{i}^{2}$ in Eq. (15). This approximation is
justified when $\nu m\omega_{\rm
p}^{2}\xi_{0}^{2}/N\ll\hbar|\omega_{3}-\omega_{2}|$, which can be achieved for
a sufficiently nonlinear potential. Note, however, that for weakly nonlinear
systems, for example, superconducting transmon qubits, this renormalization
term is highly relevant and constraints the resulting coupling constant to
$g<\sqrt{\omega_{c}\omega_{0}}$ Jaako2016 ; Bosman2017 .

Within the validity of the two-level approximation and by expressing the
resonator variables in terms of annihilation and creation operators,
$\Phi=-i\sqrt{\hbar/(2C\omega_{c})}(a^{\dagger}-a)$ and $Q=\sqrt{\hbar
C\omega_{c}/2}(a^{\dagger}+a)$, we finally obtain the cavity QED Hamiltonian

$\begin{split}H_{\rm cQED}=&\,\,\hbar\omega_{c}a^{\dagger}a+\frac{\hbar
g}{2}(a+a^{\dagger})\sum_{i}\sigma_{x}^{i}+\frac{\hbar\omega_{0}}{2}\sum_{i}\sigma^{i}_{z}\\\
&+\frac{\hbar
g^{2}}{4\omega_{c}}\sum_{i,j}\left(1+\frac{N}{\nu}\mathcal{D}_{ij}\right)\sigma_{x}^{i}\sigma_{x}^{j}.\end{split}$
(20)

In this expression we have adopted a notation more familiar in the field of
quantum optics and introduced the single-dipole coupling constant

$g=\frac{q\xi_{0}}{Cd\hbar}\sqrt{\frac{\hbar
C\omega_{c}}{2}}=\sqrt{\omega_{c}\frac{m\omega_{\rm
p}^{2}\xi_{0}^{2}}{2\hbar}\frac{\nu}{N}}.$ (21)

For weak couplings, $g\rightarrow 0$, the second line in Eq. (20) can be
neglected and $H_{\rm cQED}$ reduces to the standard Dicke model with a
collective coupling constant $G=g\sqrt{N}$. When this coupling becomes
comparable to $\omega_{c}$, the Dicke model is no longer valid and the effect
of dipole-dipole interactions and the $P^{2}$-term must be taken into account.
Note that both contributions scale as $\sim G^{2}/\omega_{c}$, as will become
more apparent in the discussion below.

###  IV.3 Coupling parameter

In Eq. (21) we have related the coupling constant $g$ to the cavity frequency
$\omega_{c}$ and the plasma frequency $\omega_{\rm p}$ such that for a
harmonic dipole, where $\xi_{0}/2=\xi_{\rm HO}=\sqrt{\hbar/(2m\omega_{0})}$,
the results of Sec. III are recovered, i.e.,
$\Delta\Omega=G\approx\sqrt{\nu}\omega_{p}$, when
$\omega_{0}\approx\omega_{c}$. In the few dipole, quantum regime the key
quantity of interest is the ratio $g/\omega_{c}=\sqrt{2\pi\alpha}$, which can
be expressed in terms of the dimensionless parameter (see also Ref.
Devoret2007 )

$\alpha=\alpha_{\rm
fs}\left(\frac{\xi_{0}}{d}\right)^{2}\left(\frac{q}{e}\right)^{2}\frac{Z}{Z_{0}}.$
(22)

Here $e\simeq 1.6\times 10^{-19}$ C is the elementary charge, $Z=\sqrt{L/C}$
the circuit impedance, and $Z_{0}=\mu_{0}/\varepsilon_{0}\approx 377\,\Omega$
the impedance of free space. For an electromagnetic mode with $Z\approx Z_{0}$
and elementary dipoles of charge $q=e$ the maximal value of $\alpha$ is set by
the finestructure constant $\alpha_{\rm fs}$, which is reached when the size
of a dipole is comparable to the size of the cavity, $\sim d$. This
illustrates the natural bound on the coupling parameter stated in Eq. (1),
which can also be obtained for an optical mode confined to a volume
$V\approx\xi_{0}^{3}$, electric transitions between Landau levels
Hagenmuller2010 , etc. Eq. (22) shows that this bound can be reached or even
overcome by using artificial atoms like superconducting qubits or quantum dots
coupled to tailored circuit resonances with $Z\gg Z_{0}$ Bosman2017 ;
Stockklauser2017 . In view of Eq. (1), one can then reinterpret such
artificial setups as regular cavity QED systems with an effective
finestructure constant $\alpha\sim\mathcal{O}(1)$. This analogy establishes an
interesting connection to the underlying theory of QED and motivates the study
of cavity QED systems in the regime $\alpha\gtrsim 1$ ($g/\omega_{c}\gtrsim
2.5$), where the electromagnetic interaction can no longer be considered as
weak.

Another interesting and in practice useful observation is that the coupling
strength is bounded by

$\frac{g}{\omega_{c}}\leq\frac{q}{2Q_{0}},$ (23)

where $Q_{0}=\sqrt{\hbar C\omega_{c}/2}$ is the magnitude of the zero-point
charge fluctuations. The non-perturbative regime is thus equivalent to the
condition that the charge induced by single dipole exceeds the quantum
fluctuations of the charge on the capacitor plate. Note that similar bounds
can also be obtained for other cavity QED implementations. For example, for a
flux qubit coupled inductively to a microwave cavity the coupling is bounded
by the ratio $g/\omega_{c}\leq\Phi_{q}/(2\Phi_{0})$, where $\Phi_{q}$ is the
flux of a qubit state and $\Phi_{0}=\sqrt{\hbar/(2C\omega_{c})}$ the magnitude
of the zero-point flux fluctuations of the cavity.

###  IV.4 Gauge non-invariance

As discussed above, Hamiltonian $H_{\rm cQED}$ represents a cavity QED model
in the electric dipole gauge, which is derived from the Lagrangian
$\mathcal{L}$ in Eq. (12). The quantization of the electromagnetic field in
free space is commonly performed in the Coulomb gauge, where the so-called
minimal coupling Hamiltonian emerges as the fundamental model for light-matter
interactions PhotonsAndAtoms . In the current setup, the Coulomb gauge is
represented by the Lagrangian,

$\mathcal{L}_{C}=\mathcal{L}-\frac{d}{dt}(\Phi Q_{\rm in}),$ (24)

which is related to $\mathcal{L}$ by a canonical transformation. In this gauge
the canonical momenta are

$\Pi=\frac{\partial\mathcal{L}}{\partial\dot{\Phi}}=C\dot{\Phi}\equiv
Q_{U},\qquad\pi_{i}=\frac{\partial\mathcal{L}}{\partial\dot{\xi}_{i}}=p_{i}+\frac{q}{d}\Phi.$
(25)

The canonical charge is now proportional to the voltage across the capacitor,
while the canonical momenta of the dipoles contain an additional magnetic
component. The resulting Hamilton operator reads

$\begin{split}H_{C}=&\frac{Q_{U}^{2}}{2C}+\frac{\Phi^{2}}{2L}+\sum_{i}\left[\frac{(p_{i}-\frac{q}{d}\Phi)^{2}}{2m}+V(\xi_{i})\right]\\\
&+\frac{m\omega_{p}^{2}}{2}\sum_{i\neq
j}\mathcal{D}_{ij}\xi_{i}\xi_{j},\end{split}$ (26)

and by identifying $\Phi$ with the magnetic vector potential $\vec{A}$, it can
be directly mapped on the minimal coupling Hamiltonian of QED. In this
representation there are no spurious dipole-dipole interactions, but when
expanding the kinetic energy term, we obtain an additional contribution
$\sim\Phi^{2}$. This is the analogue of the diamagnetic $A^{2}$-term and leads
to a positive frequency renormalization of the cavity mode. Although
Hamiltonian (14) and (26) have a different structure, the canonical
transformation in Eq. (24) ensures that both Hamiltonians represent the same
physical system. Indeed, they are related by the unitary transformation
$H_{\rm C}=\mathcal{U}H\mathcal{U}^{\dagger}$, where

$\mathcal{U}=e^{iq\sum_{i}\xi_{i}\Phi/(d\hbar)}=e^{-i\Phi Q_{\rm in}/\hbar},$
(27)

and for harmonic dipoles it can be explicitly shown that both Hamiltonians
reproduce the same spectra Todorov2012 .

This gauge equivalence, however, is only guaranteed when the full Hilbert
spaces in each representation are considered. By applying to $H_{C}$ the same
two-level approximation as in Sec. IV.2, we obtain an alternative cavity QED
Hamiltonian

$\begin{split}H_{\rm
cQED}^{\prime}=\hbar\omega_{c}a^{{\dagger}}a+\frac{\hbar\omega_{0}}{2}\sum_{i}\sigma_{z}^{i}-i\frac{\hbar
g_{C}}{2}(a^{{\dagger}}-a)\sum_{i}\sigma_{y}^{i}\\\
-\bar{\xi}^{\,2}\frac{\hbar
g_{C}^{2}N}{4\omega_{0}}(a^{\dagger}-a)^{2}+\frac{\hbar
g^{2}}{4\omega_{c}}\frac{N}{\nu}\sum_{i,j}\mathcal{D}_{ij}\sigma_{x}^{i}\sigma_{x}^{j},\end{split}$
(28)

where $g_{C}=g\omega_{0}/\omega_{c}$ and $\bar{\xi}=2\xi_{\rm HO}/\xi_{0}\geq
1$. The last inequality follows from the Thomas-Reiche-Kuhn sum rule
Lipparini2008 .

By comparing Eq. (20) and Eq. (28) it can be readily shown, for example, by
setting $\mathcal{D}_{ij}=0$ and $\omega_{0}=\omega_{c}$, that after the two-
level approximation the unitary equivalence is lost, i.e., $H_{\rm
cQED}\neq\mathcal{U}H_{\rm cQED}^{\prime}\mathcal{U}^{\dagger}$. While the
difference is negligible for weakly coupled systems, the obvious question
arises: Which is the appropriate model for cavity QED systems in the USC
regime? The answer is suggested by the following general relation between the
matrix elements of the position and the momentum operator,

$\langle\psi_{n}|p|\psi_{m}\rangle=im(\omega_{n}-\omega_{m})\langle\psi_{n}|\xi|\psi_{m}\rangle.$
(29)

This relation shows that for the momentum operator, the coupling to
energetically higher states increases with the energy gap. Therefore, for the
$(\Phi\sum_{i}p_{i})$-type coupling transitions to states out of the two-level
subspace are not systematically suppressed by a large energy denominator and
in the Coulomb gauge a two-level approximation is in general not permitted.
This basic argument can be verified numerically for explicit examples, which
will be detailed elsewhere GaugeNonInvariance .

Note that this gauge non-invariance does not contradict any of the previous
models for QED systems with atoms, molecules or intersubband transitions. In
these systems $\bar{\xi}\approx 1$ and the single-dipole coupling $g$ is very
weak such that the equivalence between the dipole and the Coulomb gauge still
holds. However, once the regime $g/\omega_{c}\sim 1$ is reached, the effective
cavity QED Hamiltonians derived in different gauges do no longer agree and
lead to qualitatively different predictions. From the analysis presented in
this section we conclude that $H_{\rm cQED}$ given in Eq. (20) represents
indeed the correct effective model for two-level dipoles coupled to a single
cavity mode, which is valid both in the weak and USC regime.

##  V The vacua of cavity QED

Due to the presence of both short- and long-range dipole-dipole interactions,
the properties of Hamiltonian $H_{\rm cQED}$ can be very complex and will also
depend in detail on the specific configuration of dipoles. For a qualitative
discussion of the possible ground states of cavity QED it is thus preferential
to proceed with a further simplification and replace the actual dipole-dipole
interactions by the corresponding all-to-all coupling,

$\frac{N}{4}\sum_{i,j}\mathcal{D}_{ij}\sigma_{x}^{i}\sigma_{x}^{j}\,\rightarrow\,\eta
S_{x}^{2}.$ (30)

Here $\eta$ is the dimensionless configuration parameter already defined in
Eq. (9) and we have introduced the collective angular momentum operators
$S_{k}=\sum_{i}\sigma_{k}^{i}/2$. This substitution maps the full Hamiltonian
$H_{\rm cQED}$ onto the extended Dicke model $(\hbar=1)$,

$H_{\rm
EDM}=\omega_{c}a^{\dagger}a+\omega_{0}S_{z}+g(a^{\dagger}+a)S_{x}+\frac{g^{2}}{\omega_{c}}\left(1+\varepsilon\right)S_{x}^{2}.$
(31)

In this model, dipole-dipole interactions are treated in an averaged way and
can be described by a single parameter $\varepsilon=\eta/\nu$. The usual Dicke
model is recovered as a specific instance of strong ferroelectric couplings,
i.e. $\varepsilon=-1$, while the cases of non-interacting ($\varepsilon=0$) or
repulsive ($\varepsilon>0$) dipoles appear, for example, in the description of
intersubband transitions Todorov2014  or certain circuit QED settings
Jaako2016 ; Armata2017 . Therefore, $H_{\rm EDM}$ interpolates between and
extends various other collective cavity-QED Hamiltonians and shows that each
of these models can be associated with a different arrangement of dipoles. We
emphasize though that the replacement in Eq. (30) is not a systematic
approximation and we will discuss some important differences between
collective spin models like $H_{\rm EDM}$ and the full Hamiltonian $H_{\rm
cQED}$ in Sec. VI.3 below.

Figure 3: Ground-state phase diagram of the extend Dicke model $H_{\rm EDM}$
as a function of the effective finestructure constant
$\alpha=g^{2}/(2\pi\omega_{c}^{2})$ (horizontal axis) and the average dipole-
dipole interaction strength $\varepsilon=\eta/\nu$ (vertical axis). For this
plot $\omega_{0}/\omega_{c}=1$ and $N=8$. The red dotted line indicates the
value of the critical coupling strength given in Eq. (34) and the other phase
boundaries are defined in the text. For each phase, the insets illustrate the
reduced state of the dipoles, $\rho_{d}$, in terms of a Bloch-sphere
representation. The color shows the value of the Q-function
$Q(\vec{n})=\langle\vec{n}|\rho_{d}|\vec{n}\rangle\in[0,1]$, where $\vec{n}$
is a unit vector and $|\vec{n}\rangle$ the corresponding coherent spin state.
Note that for a better visibility, the three insets have been plotted with
different colorscales.

Figure 3 shows a diagram of the ground states of $H_{\rm EDM}$ for different
parameters $\alpha$ and $\varepsilon$, which separates into three distinct
regimes. For weak couplings the system is in a _normal_ phase, where $\langle
a^{\dagger}a\rangle\approx 0$ and $\langle S_{z}\rangle\approx-N/2$. For
increasing $\alpha$ and $\varepsilon<0$ this phase becomes unstable and the
system undergoes a transition into a _superradiant_ phase. This phase breaks
the $Z_{2}$ symmetry of $H_{\rm EDM}$ and is characterized by a finite
expectation value $\langle a\rangle\neq 0$ and a finite polarization $\langle
S_{x}\rangle\neq 0$. In the opposite case, $\varepsilon>0$, there is a smooth
crossover into a _subradiant_ phase. This symmetry-preserving phase is
characterized by an anti-aligned spin configuration with vanishing
polarization, $\langle S_{z}\rangle\approx\langle S_{x}\rangle\approx 0$,
which decouples from the field and therefore $\langle
a^{\dagger}a\rangle\approx 0$. For $\alpha>1$ the superradiant and subradiant
phase merge and an additional sharp transition between these two phases
appears.

###  V.1 “Normal phase”

In the limit $g\rightarrow 0$ the ground state of a cavity QED system is the
normal vacuum state with $\langle a^{\dagger}a\rangle=0$ and $\langle
S_{z}\rangle=-N/2$. For finite $g$, corrections to this state can be taken
into account by a Holstein-Primakoff approximation HolsteinPrimakoff , where
the spins are replaced by harmonic oscillators, i.e., $S_{z}\rightarrow
b^{\dagger}b-N/2$, $S_{x}\rightarrow\sqrt{N}(b+b^{\dagger})/2$ and
$[b,b^{\dagger}]=1$. Under this approximation we obtain the quadratic
Hamiltonian

$H_{\rm
HP}=\omega_{c}a^{\dagger}a+\omega_{0}b^{\dagger}b+\frac{G}{2}(a+a^{\dagger})(b+b^{\dagger})+\frac{D}{4}(b+b^{\dagger})^{2},$
(32)

where $D=\left(1+\varepsilon\right)G^{2}/\omega_{c}$. $H_{\rm HP}$ can be
diagonalized by a Bogoliubov transformation and written in terms of a new set
of eigenmode operators $d_{\pm}$ as $H_{\rm
HP}=\Omega_{+}d^{\dagger}_{+}d_{+}+\Omega_{-}d_{-}^{\dagger}d_{-}$. By
identifying $\sqrt{\nu}\omega_{\rm p}\leftrightarrow
G\sqrt{\omega_{0}/\omega_{c}}$, the eigenfrequencies $\Omega_{\pm}$ are the
same as already obtained for the classical system in Eq. (10). This shows that
the vacuum state $|G\rangle$ in the normal phase is simply the ground state of
the two bright polariton modes described in Sec. III. However, one should keep
in mind that $H_{\rm HP}$ does not account for other dark polariton modes,
which in the presence of dipole-dipole interaction can lead to important
corrections and additional instabilities in the USC regime.

The presence of excitation number non-conserving terms $\sim
b^{\dagger}a^{\dagger}$ and $\sim(b^{\dagger})^{2}$ in $H_{\rm HP}$ implies
that the ground state in the normal phase still exhibits many nontrivial
properties when expressed in terms of the original field and matter modes
Ciuti2006 ; Fedortchenko2016 . A quantity of interest for the discussion below
is the ground state ‘photon number’ $\langle a^{\dagger}a\rangle$, which for
moderate couplings is approximately given by

$\langle
a^{\dagger}a\rangle\simeq\frac{Ng^{2}\omega_{0}}{4(\omega_{c}+\omega_{0})^{2}(\omega_{0}+\varepsilon
Ng^{2}/\omega_{c})}.$ (33)

Many other ground-state properties of light-matter systems in the linearized
regime have been extensively studied in the literature and will not be further
elaborated here.

###  V.2 “Superradiant phase”

For increasing coupling $g$ and $\varepsilon<0$ the normal phase eventually
becomes unstable and for $N\gg 1$ a second order phase transition into a
superradiant phase occurs. This superradiant phase exists for

$g\geq g_{c}=\sqrt{\frac{\omega_{c}\omega_{0}}{-\varepsilon N}},$ (34)

and is characterized by a finite polarization of the spins, $\langle
S_{x}\rangle$, and a finite expectation value of the field mode $\langle
a\rangle$. For $g$ close to $g_{c}$ we obtain Dimer2007

$\langle
a\rangle\simeq\pm\dfrac{Ng}{2\omega_{c}}\sqrt{1-\left(\dfrac{g_{c}}{g}\right)^{4}},\quad\langle
S_{x}\rangle\simeq\mp\dfrac{N}{2}\sqrt{1-\left(\dfrac{g_{c}}{g}\right)^{4}},$
(35)

and $\langle a\rangle\simeq\pm gN/(2\omega_{c})$ and $\langle
S_{x}\rangle\simeq\mp N/2$ for very large couplings. As shown in Fig. 4(a),
the transition into the superradiant phase is indicated by a sharp peak in the
fluctuations of the polarization, $\Delta S_{x}^{2}=\langle
S_{x}^{2}\rangle-\langle S_{x}\rangle^{2}$, and the field, $\Delta
a^{2}=\langle a^{\dagger}a\rangle-|\langle a\rangle|^{2}$, and a continuous
increase of the order parameter, $\langle
a\rangle\sim(g-g_{c})^{\frac{1}{2}}$. Note that in all our numerical
simulations we have added a small symmetry-breaking bias field, which is
necessary to deterministically pick one of the two degenerate ground states in
the symmetry-broken regime. In Fig. 3 the maximum of $\Delta S_{x}^{2}$ is
used to mark the boundary between the normal and the superradiant phase, which
even for moderate $N$ agrees reasonably well with the value of $g_{c}$
obtained from the divergence of $\langle a^{\dagger}a\rangle$ in the Holstein-
Primakoff approximation. Thus, the described transition is identical to the
conventional superradiant phase transitions discussed within the framework of
the Dicke model, but generalized to arbitrary negative values of the
interaction parameter $\varepsilon$.

Figure 4: Superradiant phase transition. (a) Dependence of the mean value
$\langle a\rangle$ and the spin- and field fluctuations across the
superradiant phase transition point. (b) Illustration of the two possible
superradiant ground states in terms of polarized dipoles and the corresponding
induced charges. (c) Comparison of the spin fluctuations $\Delta S_{x}$
evaluated with the extended Dicke model and the Lipkin-Meshkov-Glick (LMG)
Hamiltonian for two different numbers of dipoles. The value of $\alpha_{c}$,
where the fluctuations reach there maximum, as well as the width of the
fluctuations at half of the maximum, $\delta\alpha$ are plotted in (d) for
varying $N$. (e) Plot of the fluctuations of the voltage operator $U$ (solid
lines) and flux operator $\Phi$ (dashed lines) for different numbers of
dipoles. The horizontal dashed line marks the approximate analytic result
given in Eq. (39). In all plots $\omega_{0}=\omega_{c}$ and a value of
$\varepsilon=-0.1$ have been assumed and the vertical dotted lines indicate
the analytic phase transition point given in Eq. (34). In all numerical
simulations a symmetry-breaking bias field, $H_{\rm bias}=\lambda S_{x}$,
where $\lambda/\omega_{c}=10^{-3}$, has been added to the bare Hamiltonians
$H_{\rm EDM}$ and $H_{\rm LMG}$.

In Eq. (34) the phase boundary between the normal and the superradiant phase
is expressed as usual in terms of the dipole-field coupling and the cavity
frequency. This form can be very deceiving for identifying the physical origin
of this phase transition. By reexpressing $g$ instead in terms of the original
system parameters, Eq. (34) can be rewritten as

$\hbar\omega_{d}=\hbar\omega_{0}+\dfrac{\eta m\omega_{\rm
p}^{2}\xi_{0}^{2}}{2}\leq 0,$ (36)

and all cavity-related parameters disappear. While the same is true for the
phase transition point of the original Dicke model Rzazewski1975 ; Keeling2007
; Vukics2015 , Eq. (36) shows that the cancellation of $\omega_{c}$ is not
simply a coincidence. By setting $\eta=0$, the dipole-field coupling can still
be arbitrary strong, but no instability occurs. This confirms our observation
from above, namely that the superradiant instability is in essence a
ferroelectric instability and not related to the coupling to the dynamical
field mode.

We can further elaborate this point by looking more closely at the physical
properties of the superradiant phase. In the current setup a finite
expectation value $\langle a\rangle\in\mathbbm{R}$ corresponds to a finite
expectation value of the total charge $\langle Q\rangle\sim\langle
a+a^{\dagger}\rangle$, which includes charges induced by the dipoles.
Therefore, as illustrated in Fig. 4(b), the superradiant ground state simply
corresponds to a state of polarized dipoles and the corresponding induced
image charge on the capacitor plate. Since the total charge $Q$ is not
directly accessible, the more relevant resonator variables to consider are the
magnetic flux $\Phi$ and the voltage drop $U$. The latter can be expressed as

$U=U_{0}\left[a+a^{\dagger}+\frac{2g}{\omega_{c}}S_{x}\right],$ (37)

where $U_{0}=Q_{0}/C$. Importantly, the expectation values of the flux and
voltage operators are unaffected by the phase transition and we have
$\langle\Phi\rangle=\langle U\rangle=0$ in the normal as well as in the
superradiant phase. This can be seen directly from Eq. (35), or more generally
from the fact that for any stationary state $\langle\dot{\Phi}\rangle=\langle
U\rangle=0$ and $L\langle\dot{Q}\rangle=-\langle\Phi\rangle=0$. Therefore,
although the superradiant phase is conventionally characterized by a finite
‘field’ expectation value $\langle a\rangle\neq 0$, the transition affects the
displacement field, $D\sim Q$, and not the electric field, $E\sim U$
Keeling2007 . On the mean-field level the actual physical properties of the
cavity do not change when transitioning between the normal and the
superradiant phase. This example illustrates that the imprecise notion of a
‘photon’ annihilation operator $a$ can be very misleading, since depending on
the choice of gauge and the setup under consideration this operator can
represent very different physical quantities.

Given this interpretation in terms a ferroelectric phase transition, where the
‘radiation mode’ does not play a role, it would seem natural to abandon
notions of a superradiant transitions and phases all together. It should be
kept in mind though that this analogy only concerns average quantities and
strictly holds only in the limit of $N\gg 1$ and $\alpha\ll 1$. When finite-
size and strong-coupling effects are taken into account, the presence of the
electromagnetic mode can substantially influence the transition as well as all
thermodynamic properties, where excitations on top of the ground state must be
taken into account. As an example, we compare in Fig. 4(c) and (d), the
predictions from the extended Dicke model with the predictions from the
corresponding Lipkin-Meshkov-Glick Hamiltonian LMG

$H_{\rm LMG}=\omega_{0}S_{z}+\varepsilon\frac{g^{2}}{\omega_{c}}S_{x}^{2}.$
(38)

For $\varepsilon<0$ this Hamiltonian represents a model for ferroelectricity
with infinite-range interactions and can be obtained from $H_{\rm EDM}$ by
taking the limit $\nu\rightarrow 0$, but keeping $g^{2}/\nu$ fixed. We see
that for values of $|\varepsilon|N\lesssim 1$, i.e., when the transition
already happens at rather large values of $g_{c}/\omega_{c}$, the range of
fluctuations of $\Delta S_{x}^{2}$ as well as the transition point itself are
still considerably different. Only for larger numbers, $|\varepsilon|N\gg 1$,
the two models start to agree better. Overall we find that the coupling to the
cavity mode suppresses fluctuations and generates a sharper transition even
for small $N$. This can be in part explained by a dressing of the dipoles with
photons, as explained further below.

Finally, a unique signature of a superradiant transition can be obtained by
looking at cavity observables, which do not have a counterpart in
ferroelectric models. As an example, we plot in Fig. 4(e) the behavior of the
voltage fluctuations, $\langle U^{2}\rangle$, across the transition point.
While the gauge non-invariant photon number $\langle a^{\dagger}a\rangle$
diverges at the transition point, the voltage fluctuations remain finite and
show a characteristic kink. For $N\gg 1$ the position of this kink coincides
with the classical transition point $g_{c}$ and the maximal value of the
fluctuations scales approximately as

$\left.\frac{\langle
U^{2}\rangle}{U^{2}_{0}}\right|_{g=g_{c}}\approx\sqrt{1+\dfrac{1}{|\varepsilon|}\left(\dfrac{\omega_{0}}{\omega_{c}}\right)^{2}}.$
(39)

Interestingly, this maximum does neither scale with $N$ nor the coupling
parameter $\alpha$, but the kink vanishes for an interaction dominated system,
$\varepsilon\rightarrow\infty$. It thus represents a quantum mechanical
signature of a superradiant phase transition, which involves the dynamical
cavity mode. Note that the maximum of $\langle U^{2}\rangle$ is accompanied a
corresponding minimum of the flux fluctuations,
$\langle\Phi^{2}\rangle_{g=g_{c}}\simeq\hbar^{2}/(4C^{2}\langle
U^{2}\rangle)$, as expected for a minimum uncertainty squeezed state.

###  V.3 “Subradiant phase”

For repulsive dipole-dipole interactions, i.e. $\eta>0$, the linearized
Hamiltonian $H_{\rm HP}$ predicts that the normal phase remains stable for
arbitrary interaction strengths. For this reason the parameter regime
$\eta\geq 0$ and $g\gtrsim\omega_{c}$ has received little attention in the
discussion of USC cavity QED so far. However, although there is indeed no
sharp phase transition, Fig. 5(a) clearly shows that the properties of the
ground state change significantly when we go beyond the validity of the
Holstein-Primakoff approximation into the non-perturbative regime
$\alpha\gtrsim 1$. In stark contrast to the superradiant phase, the ground
state photon number in this regime decreases with increasing coupling strength
and approaches zero for very large couplings. This behavior has recently been
described in the context of circuit QED Jaako2016  and explained in terms of
an anti-ferroelectric alignment of the dipoles, which then decouple from the
cavity mode. For $N$ even, the resulting ground state is approximately of the
form

$|G\rangle\simeq|0\rangle\otimes|D_{0}\rangle,$ (40)

where $|D_{0}\rangle=|s=N/2,m_{x}=0\rangle$ is the fully symmetric Dicke state
with vanishing projection along $S_{x}$ (see Fig. 3), i.e.,
$S_{x}|D_{0}\rangle=0$. For an odd number of dipoles a perfect anti-alignment
is not possible and, as shown in Fig. 5(a), the dipoles and the cavity remain
coupled.

Figure 5: Subradiant phase. (a) Plot of the ground state photon number as a
function of the coupling parameter $\alpha$ for an even and an odd number of
dipoles. (b) Dependence of the voltage (solid lines) and flux (dashed lines)
fluctuations on the coupling strength $g$. In both plots the parameters
$\varepsilon=0.05$ and $\omega_{c}=\omega_{0}=1$ have been assumed and the
dotted lines shows the corresponding results for $\langle a^{\dagger}a\rangle$
and $\langle U^{2}\rangle$ obtained from the ground state of the Holstein-
Primakoff (HP) Hamiltonian $H_{\rm HP}$ for $N=10$. (c) Plot of the residual
single-spin entropy $\Delta\mathcal{S}$ of the ground state of $H_{\rm EDM}$
for $\omega_{0}=\omega_{c}$ and $N=4$. (d) The entanglement entropies for a
single dipole and for all dipoles are plotted for a fixed $g/\omega_{c}=2$ and
otherwise the same parameters as in (c). In all numerical simulations a
symmetry-breaking bias field, $H_{\rm bias}=\lambda S_{x}$, where
$\lambda/\omega_{c}=10^{-4}$, has been added to the bare Hamiltonians $H_{\rm
EDM}$.

The current analysis and further studies of the full model $H_{\rm cQED}$
below show that the formation of such subradiant ground states is not a
peculiarity of superconducting circuits, but rather a general property of non-
perturbative cavity QED. For an even number of dipoles, a possible way to
characterize these states is via the decoupling condition

$\frac{\partial}{\partial g}\langle a^{\dagger}a\rangle<0,$ (41)

which is used in Fig. 3 to mark the boundary between the normal and the
subradiant phase. It should be pointed out that such a light-matter decoupling
can already be predicted within the Holstein-Primakoff approximation, as
originally discussed in Ref. DeLiberato2014  for a multimode cavity QED
system. However, for the present single mode scenario, such linear decoupling
effects are not observable for the considered parameter range [see, for
example, Fig. 5(a)]. More specifically, for the current setting and within the
Holstein-Primakoff approximation the ground state photon number,

$\displaystyle\left.\lim_{g\rightarrow\infty}\langle
a^{\dagger}a\rangle\right|_{\rm
HP}=\dfrac{1+2\varepsilon-2\sqrt{\varepsilon(\varepsilon+1)}}{4\sqrt{\varepsilon(\varepsilon+1)}}>0,$
(42)

remains finite for very large couplings. Therefore, the suppression of the
photon number below this bound signifies the formation of highly entangled
anti-ferroelectric states, which decouple much more efficiently from the
cavity than the corresponding squeezed states of the linearized theory.

In Fig. 5(b) we plot the fluctuations of the observable voltage and flux
variables and find a very similar qualitative behavior as for the superradiant
transition. The voltage fluctuations show again a characteristic peak, which,
however, is much smoother and doesn’t sharpen when the number $N$ is
increased. Also the position of the maximum doesn’t vary significantly as a
function of $N$ or $\varepsilon$ and always occurs around $\alpha\approx 1$.
The absence of any significant even-odd effects make this peak in $\langle
U^{2}\rangle$ a robust signature for entering the non-perturbative coupling
regime. Interestingly, while the subradiant phase is characterized by a strong
decoupling of the dipoles from the cavity operator $a$, we find that the level
of voltage fluctuations is even higher than in the superradiant phase and also
the flux variance $\langle\Phi^{2}\rangle$ is substantially larger than for a
minimal uncertainty state.

Finally, a very interesting property which distinguishes the subradiant from
the normal and the superradiant phase, is the high degree of entanglement
between the dipoles, while being almost completely disentangled from the
cavity mode. This property can be visualized by introducing the two
entanglement entropies

$\mathcal{S}_{1}=-{\rm
Tr}\\{\rho_{1}\log_{2}(\rho_{1})\\},\qquad\mathcal{S}_{d}=-{\rm
Tr}\\{\rho_{d}\log_{2}(\rho_{d})\\}.$ (43)

Here $\rho_{d}={\rm Tr}_{c}\\{|G\rangle\langle G|\\}$ is the reduced density
operator of the dipoles, and $\rho_{1}={\rm Tr}_{N-1}\\{\rho_{d}\\}$ is the
reduced density operator of a single dipole. Therefore, $\mathcal{S}_{d}$
quantifies the entanglement between the dipoles and the cavity and
$\mathcal{S}_{1}$ the entanglement between a single dipole and the remaining
system. In Fig. 5(c) and (d) we plot the difference
$\Delta\mathcal{S}=\mathcal{S}_{1}-\mathcal{S}_{d}$ and the individual
entanglement entropies for different parameter regimes. The plots show that a
significant amount of ground-state entanglement occurs near the superradiant
phase transition, but also that this entanglement is established mainly
between the cavity and the dipoles. In contrast, when entering the subradiant
phase, $\mathcal{S}_{d}$ is strongly reduced, while the dipoles still remain
highly entangled among each other.

##  VI Non-perturbative cavity QED

The analysis in the previous section showed that for most parameter regimes
the ground state of $H_{\rm EDM}$ is either a normal vacuum state or a state
dictated by strong dipole-dipole interactions. From the perspective of cavity
QED, it is thus most interesting to consider the regime $\eta\approx 0$ and
$\alpha\gtrsim 1$, where dipole-dipole interactions play a minor role and the
influence of the cavity mode becomes important. As indicated in Fig. 3, in
this regime the superradiant and subradiant phases approach each other and a
new transition between these two very different phases emerges.

###  VI.1 Strong-coupling theory

For the following discussion we return to the full cavity QED Hamiltonian
$H_{\rm cQED}$ and focus on the regime $\omega_{0}\sim\omega_{c}$ and
$\alpha\geq 1$. In this case the coupling to the cavity $\sim g$ and the
dipole-dipole interactions $\sim g^{2}$ dominate over the bare energy
splitting of dipoles. It is thus useful to transform into a new basis, which
diagonalizes these two terms. This is achieved by a polaron transformation
$\tilde{H}_{\rm cQED}=\mathcal{U}H_{\rm cQED}\mathcal{U}^{{\dagger}}$, where
Jaako2016 ; Chen2008

$\mathcal{U}=e^{\frac{g}{\omega_{c}}S_{x}(a^{{\dagger}}-a)}.$ (44)

As a result we obtain

$\begin{split}\tilde{H}_{\rm
cQED}&=\omega_{c}a^{\dagger}a+\frac{g^{2}}{4\omega_{c}}\frac{N}{\nu}\sum_{i,j}\mathcal{D}_{ij}\sigma_{x}^{i}\sigma_{x}^{j}\\\
&+\frac{\omega_{0}}{2}\left(e^{\frac{g}{\omega_{c}}(a^{{\dagger}}-a)}\tilde{S}_{-}+e^{-\frac{g}{\omega_{c}}(a^{{\dagger}}-a)}\tilde{S}_{+}\right),\end{split}$
(45)

where $\tilde{S}_{\pm}=S_{z}\pm iS_{y}$ are collective ladder operators with
respect to $S_{x}$. Note that $\mathcal{U}$ is just the gauge transformation
(27) restricted to the two-level subspace. Therefore, Hamiltonian
$\tilde{H}_{\rm cQED}$ represents the appropriate USC cavity QED Hamiltonian
in the Coulomb gauge. By expanding the exponentials in the second line in Eq.
(45) up to first order in $\alpha$, we obtain

$\begin{split}\tilde{H}_{\rm
cQED}&\approx\omega_{c}a^{{\dagger}}a+\omega_{0}S_{z}-ig\frac{\omega_{0}}{\omega_{c}}(a^{{\dagger}}-a)S_{y}\\\
&+\frac{g^{2}}{2\omega_{c}^{2}}\omega_{0}(a^{\dagger}-a)^{2}S_{z}+\frac{g^{2}}{4\omega_{c}}\frac{N}{\nu}\sum_{i,j}\mathcal{D}_{ij}\sigma_{x}^{i}\sigma_{x}^{j},\end{split}$
(46)

which resembles very closely $H_{\rm cQED}^{\prime}$ given in Eq. (28) in the
limit of low excitation numbers. This correspondence is lost when highly
excited states or higher-order terms in the coupling parameter are taken into
account.

In the limit $\omega_{0}\rightarrow 0$, the first line of Eq. (45) is diagonal
in the photon number states $|n\rangle$ and the spin states $|s_{i}=\pm
1\rangle$, where $\sigma_{x}^{i}|s_{i}\rangle=s_{i}|s_{i}\rangle$. Therefore,
we obtain a set of eigenstates $|n,\\{s_{i}\\}\rangle$ with energies

$E^{0}_{n,\\{s_{i}\\}}=n\omega_{c}+\frac{g^{2}}{4\omega_{c}}\frac{N}{\nu}\sum_{i,j}\mathcal{D}_{ij}s_{i}s_{j}.$
(47)

Note that in the original basis the eigenstates represent displaced photon
number states,

$|\Psi\rangle_{n,\\{s_{i}\\}}=e^{-\frac{g}{\omega_{c}}S_{x}(a^{{\dagger}}-a)}|n,\\{s_{i}\\}\rangle,$
(48)

with a displacement amplitude $\beta=g/\omega_{c}\sum_{i}s_{x}^{i}$
proportional to the total spin projection along the $x$ axis.

For finite $\omega_{0}$ quantum fluctuations of the dipoles induce finite
couplings between different spin projections and different photon number
states. For $\alpha\gtrsim 1$ these couplings can be included in second order
perturbation theory following Ref. Jaako2016 . In the presence of dipole-
dipole interactions the result of such a calculation would still be very
involved, since the bare energy levels $E^{0}_{n,\\{s_{i}\\}}$ depend
explicitly on the spin configuration. For the purpose of this work we restrict
ourselves to $|\mathcal{D}_{ij}|N/\nu<g^{2}/\omega_{c}$, where this dependence
can be neglected. By projecting onto the $n=0$ sub-manifold we then obtain the
effective spin Hamiltonian

$H_{S}=\omega_{0}e^{-\frac{g^{2}}{2\omega_{c}^{2}}}S_{z}-\frac{\omega_{0}^{2}\omega_{c}}{2g^{2}}\left(\vec{S}^{2}-S_{x}^{2}\right)+\frac{g^{2}N}{4\omega_{c}\nu}\sum_{i,j}\mathcal{D}_{ij}\sigma_{x}^{i}\sigma_{x}^{j}.$
(49)

From this approximate model we see that the coupling to the cavity mode has
two main effects. First, due to the polaronic nature of the eigenstates
$|\Psi\rangle_{n,\\{s_{i}\\}}$, which contain both dipole and photonic
components, the transition frequency $\omega_{0}$ becomes exponentially
suppressed. Second, virtual excitations of higher photon-number states result
in collective dipole-dipole interactions, which favor states of maximal total
spin $S=N/2$, but minimal spin projection along $x$.

###  VI.2 Subradiant-to-superradiant phase transition

Given the effective spin Hamiltonian $H_{S}$ we can now investigate in more
detail the transition between the super- and the subradiant phase, which exist
for $\alpha\gtrsim 1$. In a first step we will consider again a collective
spin model where $\mathcal{D}_{ij}=\eta/N$. In this case the total spin is
conserved and we can restrict our analysis to states with $S=N/2$. In terms of
the effective finestructure constant we then obtain the LMG model,

$H_{S}=\omega_{0}e^{-\pi\alpha}S_{z}+\left(2\pi\alpha\varepsilon\omega_{c}+\frac{\omega_{0}^{2}}{4\pi\alpha\omega_{c}}\right)S_{x}^{2},$
(50)

with a renormalized frequency and a modified coupling term. By changing
$\varepsilon$ from positive to negative values, the expected transition from
the sub- into the superradiant phase can be determined from a Holstein-
Primakoff approximation for $H_{S}$ and we obtain the phase transition point

$\omega_{0}e^{-\pi\alpha}+N\left(2\pi\alpha\varepsilon\omega_{c}+\frac{\omega_{0}^{2}}{4\pi\alpha\omega_{c}}\right)=0.$
(51)

By omitting the second term in the brackets we obtain a condition for the
critical coupling parameter $\alpha_{c}$, which is analogous to Eq. (34), but
with a reduced dipole-frequency. This shows why for small $\varepsilon$ and
small $N$ the transition into the superradiant phase occurs at much smaller
couplings than predicted by the linearized theory (cf. Fig. 3). For $\alpha>1$
the dipole frequency is fully suppressed and the ground state phase is only
determined by the sign of the $S_{x}^{2}$-term. The resulting critical
interaction parameter $\varepsilon_{c}$ is independent of $N$ and given by

$\varepsilon_{c}\simeq-\frac{\omega_{0}^{2}}{8\pi^{2}\alpha^{2}\omega_{c}^{2}}.$
(52)

This result shows that cavity fluctuations stabilize the subradiant phase even
beyond $\varepsilon=0$ and that a small, but finite attraction between the
dipoles is required to push the system into the superradiant phase.

Figure 6: Subradiant-to-superradiant phase transition. (a) The adiabatic
potential $V_{\rm ad}(X)=X^{2}/2+E_{0}(X)$ for the cavity mode is plotted
together with the resulting ground state wavefunction for different values of
the interaction parameter $\varepsilon$. (b) Dependence of the mean cavity
field and the voltage and spin fluctuations as a function of $\varepsilon$ for
$\alpha=2$. (c) Zoom of the ground state phase diagram in the region
$|\varepsilon|\approx 0$. The color scale shows the ground state photon number
and the red solid lines indicate the same phase boundaries as in Fig. 3. The
dashed line is the critical value $\varepsilon_{c}$ given Eq. (51) and the
dotted line indicates the value of $\varepsilon_{c}$ obtained from the
classical transition point in Eq. (34). In (a) we have used $N=4$ qubits and
in (b) and (c) $N=8$. In all plots $\omega_{0}=\omega_{c}$ and a symmetry-
breaking bias field, $H_{\rm bias}=\lambda S_{x}$, where
$\lambda/\omega_{c}=10^{-3}$, has been assumed.

Figure 6(a) illustrates the subradiant-to-superradiant phase transition in
terms of the adiabatic Born-Oppenheimer potentials $V_{\rm
ad}(X)=X^{2}/2+E_{0}(X)$. Here $X=(a+a^{\dagger})/\sqrt{2}$ is the normalized
position quadrature of the cavity mode and $E_{0}(X)$ is the ground-state
energy of the spin part of the extended Dicke model,

$H_{\rm
EDM}(X)=\omega_{0}S_{z}+\sqrt{2}gXS_{x}+\frac{g^{2}}{\omega_{c}}\left(1+\varepsilon\right)S_{x}^{2},$
(53)

obtained for different values of $X$. The resulting potentials clearly show
the displaced quadratic lobes expected from the displaced oscillator states
given in Eq. (48) as well as the overall quadratic shift from the $S_{x}^{2}$
term in Eq. (50). This shift stabilized the subradiant state with $X\approx 0$
for $\varepsilon>\varepsilon_{c}$, while in the superradiant phase two minima
at $X\approx\pm Ng/\sqrt{2}\omega_{c}$ emerge. As shown in 6(b) the transition
point is indicated by a jump in $\langle a\rangle$ as well as in the voltage
fluctuations $\langle U^{2}\rangle$. This behavior is reminiscent of a first-
order phase transition, with the additional peculiarity that at the transition
point all the lobes in $V_{\rm ad}(X)$ become energetically degenerate.
Finally, Fig. 6(c) shows a zoom of the phase diagram in Fig. 3 with the
strongly modified phase boundaries in the non-perturbative regime.

###  VI.3 Beyond the collective spin approximation

In our analysis so far we have primarily focused on the collective spin model
where only the averaged dipole-dipole interaction strength appears. For
certain cavity QED implementation, in particular in the context of circuit
QED, this collective coupling arises naturally from the circuit design
Jaako2016 , in which case also $H_{\rm EDM}$ becomes exact. However, it is
clear that in general the approximation of an arbitrary coupling matrix
$\mathcal{D}_{ij}$ by a single parameter $\eta$ can lead to qualitatively very
different results. In the following we illustrate the relation between the
exact short-range and the collective spin model for two different settings.

In the first scenario shown in Fig. 7(a) the dipoles are arranged in a line
along the $x$-direction, but tilted by an angle $\theta$ with respect to the
$z$-axis. This slightly reduces the coupling to the cavity field, but changes
the dipole-dipole interactions from repulsive to attractive at a tilting angle
of about $\theta\approx 0.6$. In this case the short-range nature of the
interactions is taken into account, but the sign of all the nearest-neighbor
interactions is the same. For this configuration we evaluate the exact
coupling matrix $\mathcal{D}_{ij}(\theta)$ (including image charges) and
simulate the resulting full Hamiltonian $H_{\rm cQED}$ for different $\theta$.
We also construct the corresponding extended Dicke model by replacing
$g\longmapsto g\cos(\theta)$ and
$\eta\longmapsto\eta(\theta)/\cos^{2}(\theta)$, where $\eta(\theta)$ is
plotted in Fig. 7(b). In Fig. 7(c) we compare the results from the full model
and the corresponding $H_{\rm EDM}$ as we tune the system from repulsive to
attractive dipole-dipole interactions. We see that apart from small
quantitative differences, the qualitative features of the subradiant-to-
superradiant transition are in very good agreement.

Figure 7: Comparison between the EDM and full cavity QED Hamiltonian $H_{\rm
cQED}$. (a) Sketch of a setup with $N=4$ dipoles, where the sign of the
dipole-dipole interactions is varied by tilting the dipoles by an angle
$\theta$. The resulting averaged interaction parameter $\eta$ entering in
$H_{\rm EDM}$ is plotted in (b) as a function of $\theta$. (c) The dipole
fluctuations $\Delta S_{x}^{2}$ and the photon number $\langle
a^{{\dagger}}a\rangle$ are evaluated for the ground state of the full
Hamiltonian and for ground state of the corresponding EDM for different
$\theta$. In this plot $g/\omega_{c}=2$, $\omega_{c}=\omega_{0}$, $d=3r_{0}$
and $\nu=1/(4\pi)$. (d) Sketch of a setup, where two pairs of dipoles are
separated by $\Delta x$. For this configuration the parameter $\eta$ entering
the EDM is plotted in (e) as a function of $\Delta x$. The sign of $\eta$
changes at a value of around $\Delta x/r_{0}\approx 0.7$. (f) Plot of the
ground state correlations $\langle\sigma_{x}^{i}\sigma_{x}^{j}\rangle$
evaluated for the configuration shown in (d) using the EDM and the full model.
In all numerical simulations a symmetry-breaking bias field, $H_{\rm
bias}=\lambda S_{x}$, where $\lambda/\omega_{c}=10^{-4}$, has been added to
the bare Hamiltonians $H_{\rm EDM}$ and $H_{\rm cQED}$.

In a second scenario shown in Fig. 7(d) we consider two pairs of dipoles
placed on top of each other at a fixed distance $r_{0}$, but with a varying
separation $\Delta x$ along the $x$-direction. In this case there is a certain
distance, $\Delta x\approx 0.7r_{0}$, where the attractive interactions along
$z$ balances the repulsive interactions along $x$ and the parameter $\eta$
changes from a positive to a negative value [cf. Fig. 7(e)]. From a naive
application of the collective spin model $H_{\rm EDM}$ we would obtain around
this point a transition into a superradiant phase. However, as shown in Fig.
7(f), for the ground state of the full cavity QED Hamiltonian $H_{\rm cQED}$
this is not the case and it remains subradiant. As indicated by the values of
the dipole correlators $\langle\sigma_{x}^{i}\sigma_{x}^{j}\rangle$, this can
be understood from the fact that the two dipoles on top of each other align
and simply form a collective spin-1 particle. The two effective spin-1 dipoles
then anti-align in order to minimize the remaining attractive dipole-dipole
interactions as well as the collective coupling terms. From this basic example
we expect that in general the formation of subradiant rather than superradiant
ground states is more likely to occur.

##  VII Conclusions

In summary, we derived a minimal model for cavity QED, which is applicable in
the regime where the coupling between a single dipole and the field mode is
comparable to the bare photon energy. We discussed the physical parameters
which are required to achieve this condition in a generic setup of dipoles
coupled to the electric field of a lumped element resonator. This setting also
permitts a natural reinterpretation of the resulting dipole-field interactions
in terms an enhanced finestructure constant $\alpha$. For $\alpha\ll 1$ our
model differs from other commonly used cavity QED models mainly by the full
treatment of direct dipole-dipole interactions, which, however, is most
crucial for the correct prediction and interpretation of superradiant
instabilities. For $\alpha\gtrsim 1$ the hybridization of individual dipoles
and photons becomes relevant and leads to strong renormalization of the dipole
frequency and a cavity-induced anti-ferromagnetic ordering. This mechanism
favors highly entangled subradiant ground states, where the dipoles are almost
decoupled from the field.

While the analysis in this work was deliberately based on many idealizations
and approximations, the general findings are applicable for a large range of
different physical realizations of cavity QED systems. In traditional settings
with atoms in optical cavities the effects described in this work are not
directly accessible, since $\alpha\ll 1$ and also the required densities for
superradiant instabilities are so high that the cavity QED physics is masked
by solidification and other short-range interaction effects Griesser2016 . For
organic molecules or intersubband transitions in quantum wells the value of
$\alpha$ is still small, but ultrastrong collective couplings,
$G\sim\omega_{c},\omega_{0}$, and dipole-dipole interactions of similar
strength become possible. In these systems the interplay between dipole-field
and direct dipole-dipole interactions could be explored in more detail by
using differently structured samples, which either favor or suppress
ferroelectric order. This creates an interesting connection between
traditional studies of ferroelectric systems in confined geometries
Bratkovsky2008 , and the dynamical USC effects explored in cavity QED.

A value of $\alpha\sim 1$ can in principle be reached with superconducting
Cooper pair boxes or electrons in gate-defined quantum dots when coupled to an
$LC$ circuit with high impedance $Z\gg Z_{0}$ Devoret2007 ; Jaako2016 ;
Bosman2017 ; Manucharyan2017 . Such values are possible using superinductors
Manucharyan2009 , where $Z\sim R_{Q}\approx 26$ k$\Omega$ can become
comparable to the resistance quantum $R_{Q}$. Even higher values of $\alpha>1$
can be achieved with flux-coupled circuit QED systems, where a more favorable
scaling $g/\omega_{c}\sim\sqrt{R_{Q}/Z}$ is obtained Devoret2007 . While our
analysis has been restricted to electric systems, also the underlying
equations of motion for such flux-coupled circuits can be cast into the form
of Eqs. (2) and (6) and studied within the same theoretical framework. For
example, a serial coupling of flux qubits as considered in Ref. Jaako2016
corresponds to $\mathcal{D}_{ij}=0$ and a subradiant ground state is found.
For the parallel coupling considered in Ref. Bamba2016  we obtain
$\mathcal{D}_{ij}=\mathcal{D}<0$ and a superradiant ground state is expected.
Therefore, also quite abstract circuit geometries can be reinterpreted in
terms of interacting dipoles and described by $H_{\rm cQED}$.

Finally, let us emphasize that there are already many quantum simulation
platforms available, where $H_{\rm cQED}$ or $H_{\rm EDM}$ could potentially
be implemented as effective models Dimer2007 . For example, Rabi- and Dicke
models are currently studied with cold atoms Esslinger ; KlinderPNAS2015 ;
Schneeweiss2017  and trapped ions Lv2017 ; SafaviNaini2017 , or using digital
quantum simulation schemes with superconducting qubits Langford2016 . Similar
techniques can be used to engineer the additional $S_{x}^{2}$ terms required
for the simulation of $H_{\rm EDM}$. In Ref. Schneeweiss2017  it has been
discussed that such a collective spin term even appears for a single trapped
Rubidium atom, when its motion is coupled to Zeeman sublevels of the hyperfine
manifold via fictitious magnetic fields. For such effective models there are
in principle no constraints on the achievable parameter range and all the
different regimes of cavity QED can be explored.

###### Acknowledgements.

We thank Juraj Darmo, Karl Unterrainer and Philipp Schneeweiss for many
stimulating discussions and feedback on the manuscript. This work was
supported by the Austrian Science Fund (FWF) through the SFB FoQuS, Grant No.
F40, the DK CoQuS, Grant No. W 1210, and the START Grant No. Y 591-N16.

##  Appendix A Dipole-dipole interactions in the presence of metallic plates

To calculate dipole-dipole interactions in the presence of the capacitor
plates we follow a standard approach Takae2013  and solve the Poisson equation
$\nabla^{2}\phi(\vec{r})=-\rho(\vec{r})/\epsilon_{0}$, where $\phi(\vec{r})$
is the potential and $\rho(\vec{r})$ is the charge distribution of the dipole
ensemble, for metallic boundaries at $z=0$ and $z=d$ and with periodic
boundary conditions in the $(x,y)$-plane (with small differences, the same
calculation also holds in the case of a planar capacitor of infinite size).
This allows us to account for the overall dependence of system parameters on
the area $A=L^{2}$ and separation $d$, while avoiding a detailed numerical
simulation of the field distribution near the edges of the capacitor.

Here we consider the more general case in which all the dipoles are tilted by
an angle $\theta$ with respect to the $z$ axis. The dipole displacement is
thus given by $\vec{\xi}_{i}=\xi_{i}\vec{u}_{d}$, where
$\vec{u}_{d}=(\sin(\theta),0,\cos(\theta))$. For the evaluation of the
$\mathcal{D}_{ij}$ we need to calculate the field along the direction of each
dipole $E(\vec{r}_{i})=-\nabla\phi(\vec{r}_{i})\cdot\vec{u}_{d}$ produced by a
single dipole located at position $\vec{r}_{j}$ with charge distribution
$\rho(\vec{r})=q\xi_{j}\vec{u}_{d}\cdot\nabla\delta(\vec{r}-\vec{r}_{j})$. The
general result can be written as

$E(\vec{r})=-\frac{q}{\epsilon_{0}}\nabla\left[\vec{\xi}_{j}\cdot\nabla_{j}\left(G(\vec{r},\vec{r_{j}})\right)\right]\cdot\vec{u}_{d},$
(54)

where $G(\vec{r},\vec{r}^{\prime})$ is the Green’s function satisfying
$\nabla^{2}G(\vec{r},\vec{r}^{\prime})=-\delta^{(3)}(\vec{r}-\vec{r}^{\prime})$.
This equation can be solved by introducing the Fourier transform
$G_{\vec{k}}(z,z^{\prime})=1/(2\pi)^{2}\int
G(\vec{r},\vec{r}^{\prime})e^{i(k_{x}x+k_{y}y)}dxdy$, where
$\vec{k}=(k_{x},k_{y})$. For the boundary conditions specified above we obtain

$G_{\vec{k}}(z,z^{\prime})=g_{k}(z,z^{\prime})\Theta(z-z^{\prime})+g_{k}(z^{\prime},z)\Theta(z^{\prime}-z).$
(55)

Here $k=|\vec{k}|$, $\Theta(z)$ is the _Heaviside_ step-function and

$\begin{split}&g_{k}(z,z^{\prime})=\frac{e^{-k(z-z^{\prime})}}{2k}\\\
&-\frac{\sinh(kz)e^{-k(d-z^{\prime})}+\sinh(k(d-z))e^{-kz^{\prime}}}{2k\sinh(kd)}.\end{split}$
(56)

From this result we can immediately evaluate the total induced charge on the
capacitor plates $Q_{\rm in}=\int_{A}dxdy~{}\sigma_{\rm
in}(x,y)=\int_{A}dxdy~{}\epsilon_{0}\vec{e}_{z}\cdot\vec{E}(x,y,z=d).$ It
follows that $Q_{\rm
in}=-q\sum_{i=1}^{N}\xi_{i}\cos(\theta)\partial_{z}\partial_{z_{i}}G_{k=0}(z=d,z_{i}).$
To evaluate $G_{k=0}$ we use

$\lim_{k\rightarrow
0}G_{k}(z,z_{i})=\frac{(z+z_{i}-|z-z_{i}|)}{2}-\frac{zz_{i}}{d},$ (57)

and we obtain $Q_{\rm in}=-q\cos(\theta)\sum_{i}\xi_{i}/d$.

The full electric field in real space can be reconstructed by an inverse
Fourier transform of Eq. (56) (which can be performed in the cases of finite
size and periodic boundary conditions or infinite plane with the field
vanishing fast enough at infinity). The explicit expression of the local field
in the case of a finite system with periodic boundary conditions is

$\begin{split}\vec{E}(\vec{r})=-\frac{q}{4\pi\epsilon_{0}}\sum_{\vec{m}\in\mathbb{Z}^{3}}\nabla\cdot\nabla_{j}\left(\frac{\vec{\xi}_{j}}{\left|\vec{r}-\vec{r}_{j}-\vec{h}\right|}+\frac{\vec{\xi}_{*j}}{\left|\vec{r}-\vec{r}_{*j}-\vec{h}\right|}\right)\end{split}$
(58)

where $\vec{h}=(Lm_{x},Lm_{y},2dm_{z})$,
$\vec{\xi}_{*j}=(-\xi_{j}^{x},-\xi_{j}^{y},\xi_{j}^{z})$ and
$\vec{r}_{*j}=(x_{j},y_{j},-z_{j})$. This compact expression is nothing else
than the field generated by the $j$-th dipole, plus the field generated by the
infinite images of each dipole reflected by the metallic boundaries along $z$,
plus the field generated by the infinite copies of the system because of the
periodic boundaries along $(x,y)$. The same result holds for the infinite
capacitor, with the only difference that the summation over infinite copies of
the system disappears. Using the above definition of the tilted dipole moment,
and considering the case of an infinite planar size capacitor, we obtain

$\mathcal{D}_{ij}=\mathcal{D}_{ij}^{0}+\sum_{n\neq
0}(F^{n,S}_{ij}+F^{n,O}_{ij}).$ (59)

Here $\mathcal{D}_{ij}^{0}$ is the result given in Eq. (5) in the absence of
boundaries and

$\begin{split}F^{n,S}_{i,j}&=r_{0}^{3}\left(\frac{1}{|\vec{r}_{i}-\vec{r}_{*j}-\vec{h}_{n}|^{3}}\right.\\\
&\left.-\frac{3((z_{i}-z_{j}+2dn)\cos(\theta)+(x_{i}-x_{j})\sin(\theta))^{2}}{|\vec{r}_{i}-\vec{r}_{*j}-\vec{h}_{n}|^{5}}\right),\\\
F^{n,O}_{i,j}&=r_{0}^{3}\left(\frac{1}{|\vec{r}_{i}-\vec{r}_{*j}-\vec{h}_{n}|^{3}}\right.\\\
&\left.-\frac{3((z_{i}+z_{j}+2dn)^{2}\cos(\theta)^{2}-(x_{i}-x_{j})^{2}\sin(\theta)^{2})}{|\vec{r}_{i}-\vec{r}_{*j}-\vec{h}_{n}|^{5}}\right),\end{split}$
(60)

where $\vec{h}_{n}=(0,0,2dn)$.

##  Appendix B Polariton modes

For harmonically bound dipoles the general set of equations (2) and (6) can be
solved by decomposing the dipole variables as
$\xi_{i}=\sum_{n}c_{n}(i)\mathcal{Z}_{n}$, where the $c_{n}(i)$ are normalized
eigenmodes of the dipole-dipole interaction matrix, which obey

$\eta_{n}c_{n}(i)-\sum_{j}\mathcal{D}_{ij}c_{n}(j)=0.$ (61)

By introducing dimensionless variables $\tilde{\Phi}=\Phi/\bar{\Phi}$ and
$\tilde{\mathcal{Z}}_{n}=\mathcal{Z}/\bar{\mathcal{Z}}$, where
$\bar{\Phi}^{2}=\sqrt{L/C}m\omega_{0}\bar{\mathcal{Z}}^{2}$ and
$\bar{\mathcal{Z}}$ is an arbitrary length scale, we obtain the set of coupled
equations

$\begin{split}\ddot{\tilde{\mathcal{Z}}}_{n}+\left(\omega_{0}^{2}+\eta_{n}\omega_{\rm
p}^{2}\right)\tilde{\mathcal{Z}}_{n}=-\omega_{\rm
p}\sqrt{\frac{\omega_{0}}{\omega_{c}}\nu_{n}}\dot{\tilde{\Phi}},\\\
\qquad\ddot{\tilde{\Phi}}+\omega_{c}^{2}\tilde{\Phi}=\sum_{n}\omega_{\rm
p}\sqrt{\frac{\omega_{c}}{\omega_{0}}\nu_{n}}\dot{\tilde{\mathcal{Z}}}_{n}.\end{split}$
(62)

Here we defined the parameters
$\nu_{n}=r_{0}^{3}\left[\sum_{i}c_{n}(i)\right]^{2}/V$, which characterize the
relative coupling strength between the resonator and each dipole mode. In the
limit where the resonator is dominantly coupled to a single collective mode,
i.e., $\nu_{0}\simeq\nu$ and $\nu_{n\neq 0}\simeq 0$, we recover Eqs. (7) and
(8). The resulting eigenvalue equation is given by

$(\Omega^{2}-\omega_{c}^{2})\prod_{n}(\Omega^{2}-\omega_{n}^{2})\left[1-\sum_{n}\frac{\nu_{n}\omega_{\rm
p}^{2}\Omega^{2}}{(\Omega^{2}-\omega_{c}^{2})(\Omega^{2}-\omega_{n}^{2})}\right]=0,$
(63)

where $\omega_{n}^{2}=\omega_{0}^{2}+\eta_{n}\omega_{\rm p}^{2}$. Since we are
interested in the spectrum of coupled modes, we can assume
$\Omega\neq\omega_{c},\omega_{n}$ and look for the solutions of

$\sum_{n}\frac{\omega_{\rm
p}^{2}\nu_{n}}{\Omega^{2}-\omega_{n}^{2}}=1-\frac{\omega_{c}^{2}}{\Omega^{2}}.$
(64)

The appearance of an unstable mode is indicated by a solution of this equation
for which $\Omega\rightarrow 0$. This is only possible if one of the mode
frequencies of the dipole ensemble vanishes, i.e., for $\omega_{n}\rightarrow
0$.

##  Appendix C Single-mode approximation in the electric dipole gauge

Starting from Hamiltonian (17) in the electric dipole gauge, we can decompose
the operator for the displacement field into a set of orthogonal modes,

$\vec{D}(\vec{r})=\sum_{k}\sqrt{\frac{\hbar\omega_{k}\epsilon_{0}}{2}}\left[a_{k}\vec{f}_{k}(\vec{r})+\vec{f}^{*}_{k}(\vec{r})a_{k}^{\dagger}\right],$
(65)

where the $a_{k}$ ($a_{k}^{\dagger}$) are annihilation (creation) operators
for a mode of frequency $\omega_{k}$ and mode function $\vec{f}_{k}(\vec{r})$.
The $\vec{f}_{k}(\vec{r})$ are solutions of the Helmholtz equation for the
geometry under consideration and they are normalized to $\int
d^{3}r\vec{f}_{k}^{*}(\vec{r})\vec{f}_{k^{\prime}}(\vec{r})=\delta_{k,k^{\prime}}$.
Using this decomposition the Hamiltonian reads

$\begin{split}H_{\rm HM}&=\sum_{k}\hbar\omega_{k}a^{{\dagger}}_{k}a_{k}+H_{\rm
matter}+\frac{1}{2\epsilon_{0}}\int d^{3}rP^{2}(\vec{r})\\\ &+\sum_{k}\int
d^{3}r\sqrt{\frac{\hbar\omega_{k}}{2\epsilon_{0}}}\left[a_{k}\vec{f}_{k}(\vec{r})\cdot\vec{P}(\vec{r})+{\rm
H.c.}\right],\end{split}$ (66)

which at this stage is still exact.

We are now interested in the situation where only the lowest frequency mode is
resonant with the dipoles, i.e.
$\omega_{k_{0}}\equiv\omega_{c}\approx\omega_{0}$, while all other modes are
far separated in frequency, $\omega_{k}\gg\omega_{0}$. By looking, for
example, at the Heisenberg equations of motion for these modes,

$\partial_{t}a_{k}=-i\omega_{k}a_{k}-i\sqrt{\frac{\omega_{k}}{2\hbar\epsilon_{0}}}\int
d^{3}r\vec{f}_{k}(\vec{r})\cdot\vec{P}(\vec{r}),$ (67)

we can adiabatically eliminate their dynamics by approximating $a_{k\neq
k_{0}}(t)\simeq-\sqrt{\frac{1}{2\hbar\omega_{k}\epsilon_{0}}}\int
d^{3}r\vec{f}_{k}(\vec{r})\cdot\vec{P}(\vec{r},t)$. The resulting dynamics for
the dipoles and the remaining cavity mode can then be modelled by an effective
low-frequency Hamiltonian

$\begin{split}H_{\rm SMH}&=\hbar\omega_{c}a^{{\dagger}}_{c}a_{c}+H_{\rm
matter}+H_{\rm dd}\\\
&+\sqrt{\frac{\hbar\omega_{c}}{2\epsilon_{0}}}(a_{c}+a_{c}^{{\dagger}})\int
d^{3}r\vec{f}_{c}(\vec{r})\cdot\vec{P}(\vec{r})\\\
&+\frac{1}{2\epsilon_{0}}\left[\int
d^{3}r\vec{f}_{c}(\vec{r})\cdot\vec{P}(\vec{r})\right]^{2},\end{split}$ (68)

where we replaced the index $k_{0}$ by the index $c$ to be consistent with the
notation used in the main text. Here we have introduced the dipole-dipole
interaction term

$H_{\rm dd}=\frac{1}{2\epsilon_{0}}\sum_{i,j}\int
d^{3}rd^{3}r^{\prime}P_{i}(\vec{r})K_{ij}(\vec{r},\vec{r}^{\prime})P_{j}(\vec{r}^{\prime}),$
(69)

where
$K_{ij}(\vec{r},\vec{r}^{\prime})=\delta_{ij}\delta(\vec{r}-\vec{r}^{\prime})-\sum_{k}f_{k,i}(\vec{r})f_{k,j}(\vec{r}^{\prime})$.
Note that the sum in $K_{ij}$ runs over _all_ modes. To compensate for the
$k_{0}$ term, which has not been adiabatically eliminated, the nonlocal
$P^{2}$-term in Eq. (68) has been introduced.

Equation (68) represents a generic _single-mode_ version of the Hopfield model
in the dipole gauge (for a similar calculation in the Coulomb gauge see Ref.
Tufarelli2015 ). It shows that high frequency modes cannot be just omitted,
but they contribute in second-order perturbation theory to relevant
interactions terms between the dipoles. To illustrate this point let us
consider the limiting case $V\rightarrow\infty$, where the mode functions are
plane waves,

$\vec{f}_{\vec{k},\lambda}(\vec{r})=\vec{\varepsilon}_{\lambda}(k)\frac{1}{\sqrt{V}}e^{i\vec{k}\cdot\vec{r}},$
(70)

labeled by the wavevector $\vec{k}$ and the polarization vector
$\vec{\varepsilon}_{\lambda=1,2}\perp\vec{k}$. The kernel matrix for the
complex mode functions gives us the dipole-dipole interaction

$\begin{split}&K_{ij}(\vec{r},\vec{r}^{\prime})=\delta_{ij}\delta(\vec{r}-\vec{r}^{\prime})-\frac{1}{V}\sum_{\vec{k},\lambda}\varepsilon^{\lambda}_{i}(\vec{k})\varepsilon^{\lambda}_{j}(\vec{k})e^{i\vec{k}\cdot(\vec{r}-\vec{r}^{\prime})}\\\
&=\frac{1}{4\pi}\left[\frac{\delta_{i,j}}{\left|\vec{r}-\vec{r}^{\prime}\right|^{3}}-\frac{3(\vec{r}-\vec{r}^{\prime})_{i}(\vec{r}-\vec{r}^{\prime})_{j}}{\left|\vec{r}-\vec{r}^{\prime}\right|^{5}}\right]-\frac{1}{3}\delta(\vec{r}-\vec{r}^{\prime}),\end{split}$
(71)

where we made use of the transversality of the electro-magnetic field
PhotonsAndAtoms  and replaced the sums over the $k$-vectors by integrals. For
a finite volume and metallic boundaries, a similar calculation would reproduce
the modified dipole-dipole couplings $\sim\mathcal{D}_{ij}$, as used in our
model.

## References

  * (1) S. Haroche, and J.-M. Raimond, Exploring the Quantum: Atoms, Cavities and Photons (Oxford University Press, Oxford, 2006).
  * (2) M. H. Devoret, S. Girvin, and R. Schoelkopf, Circuit-QED: How strong can the coupling between a Josephson junction atom and a transmission line resonator be?, Ann. Phys. (NY) 16, 767 (2007).
  * (3) Note that the bound derived in Refs. HarocheCQED  and Devoret2007  assumes a resonant coupling between a microwave cavity and a Rydberg atom, which results in a different scaling with $\alpha_{\rm fs}$.
  * (4) R. H. Dicke, Coherence in Spontaneous Radiation Processes, Phys. Rev. 93, 99 (1954).
  * (5) T. Brandes, Coherent and collective quantum optical effects in mesoscopic systems, Physics Reports 408, 315 (2005).
  * (6) B. M. Garraway, The Dicke model in quantum optics: Dicke model revisited, Phil. Trans. R. Soc. 369, 1137 (2011).
  * (7) K. Hepp, and E. H. Lieb, On the superradiant phase transition for molecules in a quantized radiation field: the Dicke maser model, Ann. Phys. 76, 360 (1973).
  * (8) Y. K. Wang, and F. T. Hioe, Phase Transition in the Dicke Model of Superradiance, Phys. Rev. A 7, 831 (1973).
  * (9) T. Hioe, Phase Transitions in Some Generalized Dicke Models of Superradiance, Phys. Rev. A 8, 1440 (1973).
  * (10) K. Rzazewski, K. Wodkiewicz, and W. Zakowicz, Phase Transitions, Two-Level Atoms, and the $A^{2}$ Term, Phys. Rev. Lett. 35, 432 (1975).
  * (11) Y. A. Kudenko, A. P. Slivinsky, and G. M. Zaslavsky, Interatomic Coulomb interaction influence on the superradiance phase transition, Phys. Lett. A 50, 411 (1975).
  * (12) V. L Emelanov, and Yu. I.I Qimontovioh, Appearance of collective polarisation as a result of phase transition in an ensemble of two-level atoms, interacting through electromagnetic field, Phys. Lett. A 59, 366 (1976).
  * (13) J. M. Knight, Y. Aharonov, and G. T. C. Hsieh, Are super-radiant phase transitions possible?, Phys. Rev. A 17, 1454 (1978).
  * (14) I. Bialynicki-Birula, and K. Rzazewski, No-go theorem concerning the superradiant phase transition in atomic systems, Phys. Rev. A 19, 301 (1979).
  * (15) C. M. Bowden and C. C. Sung, Effect of spin-spin interaction on the super-radiant phase transition for the Dicke Hamiltonian, J. Phys. A: Math. Gen. 12, 2273 (1979).
  * (16) M. Yamanoi, On polariton instability and thermodynamic phase transition in a photon-matter system, J. Phys. A 12, 1591 (1979).
  * (17) J. Keeling, Coulomb interactions, gauge invariance, and phase transitions of the Dicke model, J. Phys: Cond. Mat. 19, 295213 (2007).
  * (18) D. Hagenmüller, S. De Liberato, and C. Ciuti, Ultrastrong coupling between a cavity resonator and the cyclotron transition of a two-dimensional electron gas in the case of an integer filling factor, Phys. Rev. B 81, 235303 (2010).
  * (19) Y. Todorov, and C. Sirtori, Intersubband polaritons in the electrical dipole gauge, Phys. Rev. B 85, 045304 (2012).
  * (20) A. Vukics, and P. Domokos, Adequacy of the Dicke model in cavity QED: A counter-no-go statement, Phys. Rev. A 86, 053807 (2012).
  * (21) S. De Liberato, Light-Matter Decoupling in the Deep Strong Coupling Regime: The Breakdown of the Purcell Effect, Phys. Rev. Lett. 112, 016401 (2014).
  * (22) A. Vukics, T. Griesser, and P. Domokos, Elimination of the $A$-Square Problem from Cavity QED, Phys. Rev. Lett. 112, 073601 (2014).
  * (23) M. Bamba, and T. Ogawa, Stability of polarizable materials against superradiant phase transition, Phys. Rev. A 90, 063825 (2014).
  * (24) A. Vukics, T. Grießer, and P. Domokos, Fundamental limitation of ultrastrong coupling between light and atoms, Phys. Rev. A 92, 043835 (2015).
  * (25) T. Tufarelli, K. R. McEnery, S. A. Maier, and M. S. Kim, Signatures of the $A^{2}$ term in ultrastrongly coupled oscillators, Phys. Rev. A 91, 063840 (2015).
  * (26) T. Grießer, A. Vukics, and P. Domokos, Depolarization shift of the superradiant phase transition, Phys. Rev. A 94, 033815 (2016).
  * (27) P. Nataf, and C. Ciuti, No-go theorem for superradiant quantum phase transitions in cavity QED and counter-example in circuit QED, Nature Commun. 1, 72 (2010).
  * (28) P. Nataf, and C. Ciuti, Vacuum Degeneracy of a Circuit QED System in the Ultrastrong Coupling Regime, Phys. Rev. Lett. 104, 023601 (2010).
  * (29) O. Viehmann, J. von Delft, and F. Marquardt, Superradiant Phase Transitions and the Standard Description of Circuit QED, Phys. Rev. Lett. 107, 113602 (2011).
  * (30) M. Leib, and M. J. Hartmann, Synchronized Switching in a Josephson Junction Crystal, Phys. Rev. Lett. 112, 223603 (2014).
  * (31) T. Jaako, Z.-L. Xiang, J.J. Garcia-Ripoll, and P. Rabl, Ultrastrong coupling phenomena beyond the Dicke model, Phys. Rev. A 94, 033850 (2016).
  * (32) M. Bamba, K. Inomata, and Y. Nakamura, Superradiant Phase Transition in a Superconducting Circuit in Thermal Equilibrium, Phys. Rev. Lett. 117, 173601 (2016).
  * (33) M. Bamba, and N. Imoto, Circuit configurations which can/cannot show super-radiant phase transitions, Phys. Rev. A 96, 053857 (2017).
  * (34) C. Ciuti, G. Bastard, and I. Carusotto, Quantum vacuum properties of the intersubband cavity polariton field, Phys. Rev. B 72, 115303 (2005).
  * (35) T. Schwartz, J. A. Hutchison, C. Genet, and T. W. Ebbesen, Reversible Switching of Ultrastrong Light-Molecule Coupling, Phys. Rev. Lett. 106, 196405 (2011).
  * (36) S. Kena-Cohen, S. A. Maier, and D. D. C. Bradley, Ultrastrongly Coupled Exciton-Polaritons in Metal-Clad Organic Semiconductor Microcavities, Adv. Opt. Mater. 1, 827 (2013).
  * (37) M. Mazzeo, A. Genco, S. Gambino, D. Ballarini, F. Mangione, O. Di Stefano, S. Patan , S. Savasta, D. Sanvitto, and G. Gigli, Ultrastrong light-matter coupling in electrically doped microcavity organic light emitting diodes, Appl. Phys. Lett. 104, 233303 (2014).
  * (38) Y. Todorov, A. M. Andrews, R. Colombelli, S. De Liberato, C. Ciuti, P. Klang, G. Strasser, and C. Sirtori, Ultrastrong Light-Matter Coupling Regime with Polariton Dots, Phys. Rev. Lett. 105, 196402 (2010).
  * (39) M. Geiser, F. Castellano, G. Scalari, M. Beck, L. Nevou, and J. Faist, Ultrastrong Coupling Regime and Plasmon Polaritons in Parabolic Semiconductor Quantum Wells, Phys. Rev. Lett. 108, 106402 (2012).
  * (40) A. Benz, S. Campione, S. Liu, I. Montano, J.F. Klem, A. Allerman, J. R. Wendt, M. B. Sinclair, F. Capolino, and I. Brener, Strong coupling in the sub-wavelength limit using metamaterial nanocavities, Nature Commun. 4, 2882 (2013).
  * (41) D. Dietze, A. M. Andrews, P. Klang, G. Strasser, K. Unterrainer, and J. Darmo, Ultrastrong coupling of intersubband plasmons and terahertz metamaterials, Appl. Phys. Lett. 103, 201106 (2013).
  * (42) A. Vasanelli, Y. Todorov, and C. Sirtori, Ultra-strong light-matter coupling and superradiance using dense electron gases, Comptes Rendus Physique 17, 8, 871-873 (2016).
  * (43) G. Scalari, C. Maissen, D. Turcinkova, D. Hagenmüller, S. De Liberato, C. Ciuti, C. Reichl, D. Schuh, W. Wegscheider, M. Beck, and J. Faist, Ultrastrong Coupling of the Cyclotron Transition of a 2D Electron Gas to a THz Metamaterial, Science 335, 1323 (2012).
  * (44) C. Maissen, G. Scalari, F. Valmorra, M. Beck, S. Cibella, R. Leoni, C. Reichl, C. Charpentier, W. Wegscheider, and J. Faist, Ultrastrong coupling in the near field of complementary split-ring resonators, Phys. Rev. B 90, 205309 (2014).
  * (45) Q. Zhang, M. Lou, X. Li, J. L. Reno, W. Pan, J. D. Watson, M. J. Manfra, and J. Kono, Collective non-perturbative coupling of 2D electrons with high-quality-factor terahertz cavity photons, Nature Phys. 12, 1005 (2016).
  * (46) J. Keller, G. Scalari, F. Appugliese, C. Maissen, J. Haase, M. Failla, M. Myronov, D. R. Leadley, J. Lloyd-Hughes, P. Nataf, and J. Faist, Critical softening of cavity cyclotron polariton modes in strained germanium 2D hole gas in the ultra-strong coupling regime, arXiv:1708.07773 (2017).
  * (47) A. Wallraff, D. I. Schuster, A. Blais, L. Frunzio, R. S. Huang, J. Majer, S. Kumar, S. M. Girvin, and R. J. Schoelkopf, Strong coupling of a single photon to a superconducting qubit using circuit quantum electrodynamics, Nature (London) 431, 162 (2004).
  * (48) A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin, and R. J. Schoelkopf, Cavity quantum electrodynamics for superconducting electrical circuits: An architecture for quantum computation, Phys. Rev. A 69, 062320 (2004).
  * (49) J. Q. You, and F. Nori, Atomic physics and quantum optics using superconducting circuits, Nature (London) 474, 589 (2011).
  * (50) X. Mi, J. V. Cady, D. M. Zajac, P. W. Deelman, and J. R. Petta, Strong Coupling of a Single Electron in Silicon to a Microwave Photon, Science 355, 156 (2017).
  * (51) A. Stockklauser, P. Scarlino, J. V. Koski, S. Gasparinetti, C. K. Andersen, C. Reichl, W. Wegscheider, T. Ihn, K. Ensslin, and A. Wallraff, Strong Coupling Cavity QED with Gate-Defined Double Quantum Dots Enabled by a High Impedance Resonator, Phys. Rev. X 7, 011030 (2017).
  * (52) L. E. Bruhat, T. Cubaynes, J. J. Viennot, M. C. Dartiailh, M. M. Desjardins, A. Cottet, and T. Kontos, Strong Coupling Between an Electron in a Quantum Dot Circuit and a Photon in a Cavity, arXiv:1612.05214.
  * (53) A. Cottet, M. C. Dartiailh, M. M. Desjardins, T. Cubaynes, L. C. Contamin, M. Delbecq, J. J. Viennot, L. E. Bruhat, B. Doucot, and T. Kontos, Cavity QED with hybrid nanocircuits: from atomic-like physics to condensed matter phenomena, J. Phys.: Condens. Matter 29, 433002 (2017).
  * (54) T. Niemczyk, F. Deppe, H. Huebl, E. P. Menzel, F. Hocke, M. J. Schwarz, J. J. Garcia-Ripoll, D. Zueco, T. Hü mmer, E. Solano, A. Marx, and R. Gross, Circuit quantum electrodynamics in the ultrastrong-coupling regime, Nature Phys. 6, 772 (2010).
  * (55) P. Forn-Diaz, J. Lisenfeld, D. Marcos, J. J. Garcia-Ripoll, E. Solano, C. J. P. M. Harmans, and J. E. Mooij, Observation of the Bloch-Siegert Shift in a Qubit-Oscillator System in the Ultrastrong Coupling Regime, Phys. Rev. Lett. 105, 237001 (2010).
  * (56) A. Baust, E. Hoffmann, M. Haeberlein, M. J. Schwarz, P. Eder, J. Goetz, F. Wulschner, E. Xie, L. Zhong, F. Quijandría, D. Zueco, J.-J. García-Ripoll, L. García-Álvarez, G. Romero, E. Solano, K. G. Fedorov, E. P. Menzel, F. Deppe, A. Marx, and R. Gross, Ultrastrong coupling in two-resonator circuit QED, Phys. Rev. B 93, 214501 (2016).
  * (57) P. Forn-Diaz, J. J. García-Ripoll, B. Peropadre, M. A. Yurtalan, J.-L. Orgiazzi, R. Belyansky, C. M. Wilson, and A. Lupascu, Ultrastrong coupling of a single artificial atom to an electromagnetic continuum, Nature Phys. 13, 39 (2017).
  * (58) F. Yoshihara, T. Fuse, S. Ashhab, K. Kakuyanagi, S. Saito, and K. Semba, Superconducting qubit-oscillator circuit beyond the ultrastrong-coupling regime, Nature Phys. 13, 44 (2017).
  * (59) Z. Chen, Y. Wang, T. Li, L. Tian, Y. Qiu, K. Inomata, F. Yoshihara, S. Han, F. Nori, J. S. Tsai, and J. Q. You, Multi-photon sideband transitions in an ultrastrongly-coupled circuit quantum electrodynamics system, Phys. Rev. A 96, 012325 (2017).
  * (60) S. J. Bosman, M. F. Gely, V. Singh, D. Bothner, A. Castellanos-Gomez, and G. A. Steele, Approaching ultrastrong coupling in transmon circuit QED using a high-impedance resonator, Phys. Rev. B 95, 224515 (2017).
  * (61) J. A. Hutchison, T. Schwartz, C. Genet, E. Devaux, and T. W. Ebbesen, Modifying Chemical Landscapes by Coupling to Vacuum Fields, Angew. Chem., Int. Ed. Engl. 51, 1592 (2012).
  * (62) S. Wang, A. Mika, J. A. Hutchison, C. Genet, A. Jouaiti, M.W. Hosseini, and T.W. Ebbesen, Phase Transition of a Perovskite Strongly Coupled to the Vacuum Field, Nanoscale 6, 7243 (2014).
  * (63) J. Galego, F. J. Garcia-Vidal, and J. Feist, Cavity-Induced Modifications of Molecular Structure in the Strong-Coupling Regime, Phys. Rev. X 5, 041022 (2015).
  * (64) F. Herrera, and F. C. Spano, Cavity-Controlled Chemistry in Molecular Ensembles, Phys. Rev. Lett. 116, 238301 (2016).
  * (65) J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, Atoms and Molecules in Cavities: From Weak to Strong Coupling in QED Chemistry, PNAS USA 114, 3026 (2017).
  * (66) G. Romero, D. Ballester, Y. M. Wang, V. Scarani, and E. Solano, Ultrafast Quantum Gates in Circuit QED, Phys. Rev. Lett. 108, 120501 (2012).
  * (67) F. Armata, G. Calajo, T. Jaako, M. S. Kim, and P. Rabl, Harvesting Multiqubit Entanglement from Ultrastrong Interactions in Circuit Quantum Electrodynamics, Phys. Rev. Lett. 119, 183602 (2017).
  * (68) This regime is also known as the deep strong coupling regime Casanova2010 .
  * (69) J. Casanova, G. Romero, I. Lizuain, J. J. García-Ripoll, and E. Solano, Deep Strong Coupling Regime of the Jaynes-Cummings Model, Phys. Rev. Lett. 105, 263603 (2010).
  * (70) C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, _Photons and Atoms_ (Wiley, New York, 1997).
  * (71) S. S. Scheel, L. Knöll, and D.-G. Welsch, _Coherence and Statistics of Photons and Atoms_ (Wiley, New York, 2001).
  * (72) J. J. Hopfield, Theory of the Contribution of Excitons to the Complex Dielectric Constant of Crystals, Phys. Rev. 112, 1555 (1958).
  * (73) J. W. Perram, Simulations at conducting interfaces: Boundary conditions for electrodes and electrolytes, J. Chem. Phys 104, 5174 (1996).
  * (74) L. D. Landau, E. M. Lifshitz, _Electrodynamics of Continuous Media_ , (Pergamon Press, Oxford, 1960).
  * (75) K. Takae, and A. Onuki, Applying electric field to charged and polar particles between metallic plates: Extension of the Ewald method, J. Chem. Phys. 139, 124108 (2013).
  * (76) A. M. Bratkovsky, and A. P. Levanyuk , Continuous Theory of Ferroelectric States in Ultrathin Films with Real Electrodes, Journal of Computational and Theoretical Nanoscience 6, 10.1166/jctn.2009.1058 (2008).
  * (77) A. P. Levanyuk, and I. B. Misirlioglu, Strong influence of non-ideality of electrodes on stability of single domain state in ferroelectric-paraelectric superlattices, J. App. Phys. 119, 024109 (2016).
  * (78) B. R. A. Nijboer, and F. W. De Wette, The internal field in dipole lattices, Physica 24, 442 (1958).
  * (79) L. Dell’Anna, and M. Merano, Clausius-Mossotti Lorentz-Lorenz relations and retardation effects for two-dimensional crystals, Phys. Rev. A 93, 053808 (2016).
  * (80) E. Lipparini, _Modern Many-Particle Physics: Atomic Gases, Nanostructures and Quantum Liquids_ (World Scientific Publishing Company, Singapore, 2008).
  * (81) D. De Bernardis, P. Pilar, T. Jaako, S. De Liberato, and P. Rabl _in preparation_.
  * (82) Y. Todorov, and C. Sirtori, Few-Electron Ultrastrong Light-Matter Coupling in a Quantum LC Circuit, Phys. Rev. X 4, 041031 (2014).
  * (83) T. Holstein and H. Primakoff, Field Dependence of the Intrinsic Domain Magnetization of a Ferromagnet, Phys. Rev. 58, 1098 (1940).
  * (84) C. Ciuti, and I. Carusotto, Input-output theory of cavities in the ultrastrong coupling regime: The case of time-independent cavity parameters, Phys. Rev. A 74, 033811 (2006).
  * (85) S. Fedortchenko, S. Huppert, A. Vasanelli, Y. Todorov, C. Sirtori, C. Ciuti, A. Keller, T. Coudreau, and P. Milman, Output squeezed radiation from dispersive ultrastrong light-matter coupling, Phys. Rev. A 94, 013821 (2016).
  * (86) F. Dimer, B. Estienne, A. S. Parkins, and H. J. Carmichael, Proposed realization of the Dicke-model quantum phase transition in an optical cavity QED system, Phys. Rev. A 75, 013804 (2007).
  * (87) H. J. Lipkin, N. Meshkov, and A. J. Glick, Validity of many-body approximation methods for a solvable model: (I). Exact solutions and perturbation theory, Nucl. Phys. 62, 188 (1965).
  * (88) Q.-H. Chen, Y.-Y. Zhang, T. Liu, and K.-L. Wang, Numerically exact solution to the finite-size Dicke model, Phys. Rev. A 78, 051801(R) (2008).
  * (89) V. E. Manucharyan, A. Baksic, and C. Ciuti, Resilience of the quantum Rabi model in circuit QED, J. Phys. A: Math. Theor. 50, 294001 (2017).
  * (90) V. E. Manucharyan, J. Koch, L. I. Glazman, and M. H. Devoret, Fluxonium: Single Cooper-Pair Circuit Free of Charge Offsets, Science 326, 113 (2009).
  * (91) K. Baumann, C. Guerlin, F. Brennecke, and T. Esslinger, Dicke quantum phase transition with a superfluid gas in an optical cavity, Nature 464, 1301 (2010).
  * (92) J. Klinder, H. Keßler, M. Wolke, L. Mathey, and A. Hemmerich, Dynamical phase transition in the open Dicke model, PNAS 112, 3290 (2015).
  * (93) P. Schneeweiss, A. Dareau, and C. Sayrin, Cold-atom based implementation of the quantum Rabi model, arXiv:1706.07781.
  * (94) D. Lv, S. An, Z. Liu, J.-N. Zhang, J. S. Pedernales, L. Lamata, E. Solano, and K. Kim, Quantum simulation of the quantum Rabi model in a trapped ion, arXiv:1711.00582.
  * (95) A. Safavi-Naini, R. J. Lewis-Swan, J. G. Bohnet, M. Garttner, K. A. Gilmore, E. Jordan, J. Cohn, J. K. Freericks, A. M. Rey, and J. J. Bollinger, Exploring adiabatic quantum dynamics of the Dicke model in a trapped ion quantum simulator, arXiv:1711.07392.
  * (96) N. K. Langford, R. Sagastizabal, M. Kounalakis, C. Dickel, A. Bruno, F. Luthi, D. J. Thoen, A. Endo, and L. DiCarlo, Experimentally simulating the dynamics of quantum light and matter at ultrastrong coupling, Nature Commun. 8, 1715 (2017).