Thermodynamics of ultrastrongly coupled light-matter
systems
Philipp Pilar
1
, Daniele De Bernardis
1
, and Peter Rabl
1
1
Vienna Center for Quantum Science and Technology, Atominstitut, TU Wien, 1040 Vienna, Austria
We study the thermodynamic properties of
a system of two-level dipoles that are coupled
ultrastrongly to a single cavity mode. By us-
ing exact numerical and approximate analyt-
ical methods, we evaluate the free energy of
this system at arbitrary interaction strengths
and discuss strong-coupling modifications of
derivative quantities such as the specific heat
or the electric susceptibility. From this analy-
sis we identify the lowest-order cavity-induced
corrections to those quantities in the collec-
tive ultrastrong coupling regime and show that
for even stronger interactions the presence of
a single cavity mode can strongly modify ex-
tensive thermodynamic quantities of a large
ensemble of dipoles. In this non-perturbative
coupling regime we also observe a significant
shift of the ferroelectric phase transition tem-
perature and a characteristic broadening and
collapse of the black-body spectrum of the cav-
ity mode. Apart from a purely fundamen-
tal interest, these general insights will be im-
portant for identifying potential applications
of ultrastrong-coupling effects, for example, in
the field of quantum chemistry or for realizing
quantum thermal machines.
1 Introduction
Undoubtedly, the interplay between statistical physics
and the theory of electromagnetic (EM) radiation
played a very important role in the history of mod-
ern physics. Discrepancies between the predicted and
the measured spectrum of black-body radiation led
to the birth of quantum mechanics. Based on purely
thermodynamic arguments, Einstein introduced his
A-coefficient and postulated the effect of spontaneous
emission, long before it was understood microscop-
ically. Investigations of photon-photon correlations
from thermal and coherent sources of light stood at
the beginning of the field of quantum optics, and so
on. In most of these and related examples the EM
field can be treated as an independent subsystem,
which thermalizes via weak interactions with the sur-
rounding matter. This assumption breaks down in the
so-called ultrastrong coupling (USC) regime [1, 2, 3],
where the interaction energy can be comparable to the
bare energy of the photons. Such conditions can be
reached in solid-state [4, 5, 6, 7, 8, 9, 10] and molecu-
lar cavity QED experiments [11, 12, 13, 14, 15], where
modifications of chemical reactions [16, 17] or phase
transitions [18] have been observed and interpreted
as vacuum-induced changes of thermodynamic poten-
tials [19]. Together with the ability to realize even
stronger couplings between artificial superconducting
atoms and microwave photons [20, 21, 22, 23, 24],
these observations have led to a growing interest [2, 3]
in the ground and thermal states of light-matter sys-
tems under conditions where the coupling between the
individual parts can no longer be neglected.
Since an exact theoretical treatment of light-matter
systems in the USC regime is in general not pos-
sible, one usually resorts to simplified descriptions,
for example, based on the Dicke [25, 26] or the
Hopfield [27] model. However, such reduced mod-
els often do not represent the complete energy of
the system [28, 29, 30, 31, 32, 33, 34, 35] or con-
tain gauge artefacts [33, 36, 37, 38, 39] that prevent
their applicability in the USC regime. More gener-
ally, while in weakly coupled cavity QED systems
the role of static dipole-dipole interactions can of-
ten be neglected or modelled independently of the
dynamical EM mode, this is no longer the case in
the USC regime [33, 40, 41, 42, 43, 44]. An inconsis-
tent treatment of static and dynamical fields can thus
very easily lead to wrong predictions or a misinter-
pretation of results. A prominent example in this re-
spect is the superradiant phase transition of the Dicke
model [45, 46, 47], which is often described as cavity-
induced, but which can be understood as a regular
ferroelectric instability in a system of strongly attrac-
tive dipoles [33, 41]. In the past, these and other
subtle issues have led to many controversies in this
field and prevented a detailed understanding of the
ground- and thermal states of USC light-matter sys-
tems so far.
In this paper we study the thermodynamics of cav-
ity and circuit QED systems within the framework of
the extended Dicke model (EDM) [32, 33]. Although
based on several simplifications, such as the two-level
and the single-mode approximation, this model re-
mains consistent with basic electrodynamics at arbi-
trary interaction strengths and distinguishes explic-
itly between static and dynamical electric fields. It
thus allows us to evaluate the free energy of the most
Accepted in Quantum 2020-09-15, click title to verify. Published under CC-BY 4.0. 1
arXiv:2003.11556v5 [quant-ph] 22 Sep 2020
relevant dynamical degrees of freedom in cavity QED
and to study the thermal equilibrium states of the
combined system for a macroscopic number of N 1
dipoles and in various coupling regimes.
Our analysis shows, first of all, that in the collective
USC regime, where G = g
N is comparable to the
cavity frequency ω
c
, but where the coupling g between
the cavity and each individual dipole is still small, the
coupling-induced corrections to the free energy scales
as F
g
~g
2
N
c
> 0. This very generic result, which
also holds for arbitrary dipolar systems, shows that
the coupling to a cavity mode leads to a positive shift
of the free energy and when taking the limit N
for a fixed value of G, the changes in the free energy
per particle, F
g
/N, vanish. Both findings contradict
the common intuition built upon the analysis of the
Dicke model, which predicts negative corrections to
the free energy [19] and that a large collective cou-
pling to a quantized field mode can induce substantial
modifications of material properties [45, 46, 47]. Our
results are, however, consistent with similar conclu-
sions obtained in studies about molecular properties
in the ground state of cavity QED systems [48, 49, 50],
and can be intuitively explained by a simple polari-
ton picture: In the collective USC regime the cavity
field only couples to a single collective dipole mode
while the other N 1 orthogonal excitation modes
remain unaffected. Therefore, the presence of a single
cavity mode should not have a considerable impact
on the thermodynamics of a macroscopic ensemble of
dipoles.
Surprisingly, in the regime g ω
c
this intuition
is no longer true and we find that the coupling to
the cavity can indeed influence quantities such as the
electric susceptibility or the specific heat, or even the
phase transition temperature of a ferroelectric mate-
rial. This creates a highly unusual situation where
the addition of a single degree of freedom changes the
thermodynamics of a macroscopic system. Further,
we show that the different coupling regimes of cavity
QED result in very distinct features in the black-body
spectrum of the cavity. As g in increased, the spec-
Figure 1: Sketch of a cavity QED setup where an ensemble
of dipoles is coupled to the electric field of a lumped-element
LC resonator. The system is in thermal contact with a bath
of temperature T . The black-body spectrum of the cavity
mode, S
bb
(ω), can be measured through a weak capacitive
link to a cold transmission line. See text for more details.
trum evolves from the usual polariton doublet into
a broad and disordered set of lines and, finally, col-
lapses again to a single resonance. At the same time
we find that already at moderate coupling strengths,
the light-matter interaction can either enhance or sup-
press the total radiated power. Therefore, the anal-
ysis of the EDM already provides many conceptually
important predictions, which can serve as a basis for
more detailed investigations of thermal effects in real
and artificial light-matter systems.
The remainder of the paper is structured as follows:
In Sec. 2 and Sec. 3 we first introduce the EDM and
discuss some general properties of the free energy of
a cavity QED system in different coupling regimes.
In Sec. 4 we then analyze in more detail the cavity-
induced modifications for the cases of paraelectric and
ferroelectric ensembles of dipoles. Finally, in Sec. 5
we evaluate the black-body spectrum of the cavity
mode in different coupling regimes and we conclude
our work in Sec. 6. The appendices A-D contain addi-
tional details about different approximation methods
for the free energy and the derivation of the emission
spectrum.
2 Model
We consider a generic cavity QED setup, where N
two-level dipoles are coupled to a single electromag-
netic mode. However, since we are interested in both
thermal und USC effects, we can restrict our discus-
sion to cavity and circuit QED setups in the GHz
to THz regime, where these effects are experimen-
tally most relevant. In this case the confined electro-
magnetic field can be represented by the fundamental
mode of a lumped-element LC resonator [33, 51] with
capacitance C and inductance L (see Fig. 1). This
configuration also ensures that all higher EM exci-
tations are well separated in frequency and that the
electric field is approximately homogeneous across the
ensemble of dipoles. The dipoles themselves are mod-
eled as effective two-level systems with states |0i and
|1i. The two states are separated by an energy ~ω
0
and they are coupled via an electric transition dipole
moment µ to the electric field.
2.1 Hamiltonian
The Hamiltonian of the whole cavity QED system can
be written as
H
cQED
= H
em
+ H
dip
, (1)
where the two terms represent the energies of the EM
mode and the dipoles, respectively. We model the
bare dynamics of the dipoles by a spin Hamiltonian
of the form
H
dip
=
~ω
0
2
N
X
i=1
σ
i
z
+ ~
N
X
i,j=1
J
ij
4
σ
i
x
σ
j
x
, (2)
Accepted in Quantum 2020-09-15, click title to verify. Published under CC-BY 4.0. 2
where the σ
i
k
are the usual Pauli operators for the
i-th dipole. The couplings J
ij
account for the effect
of static dipole-dipole interactions as well as possi-
ble other types of non-electromagnetic couplings be-
tween the two-level systems. For all of the explicit
calculations below we will consider the special case
of all-to-all interactions, J
ij
= J/N. In this limit
the spin Hamiltonian reduces to the Lipkin-Meshkov-
Glick (LMG) model [52]
H
dip
= ~ω
0
S
z
+
~J
N
S
2
x
H
LMG
, (3)
where the S
k
= 1/2
P
i
σ
i
k
are collective spin opera-
tors. For the current purpose, this model is sufficient
to capture the qualitative features of non-interacting
(J = 0), ferroelectric (J < 0) and anti-ferroelectric
(J > 0) dipolar systems, while still being simple
enough to allow for exact numerical simulations for
moderate system sizes. However, we emphasize that
none of the general conclusions and theoretical ap-
proaches in this work depend on the assumption of
all-to-all interactions and can be extended to arbitrary
dipolar systems using more sophisticated numerical
techniques [53].
In the lumped-element limit, the energy of the EM
mode is given by
H
em
=
CV
2
2
+
Φ
2
2L
=
(Q + P/d)
2
2C
+
Φ
2
2L
, (4)
where V is the voltage difference across the capaci-
tor plates and Φ the magnetic flux. After the second
equality sign we have expressed the capacitive energy
in terms of the total charge Q, which is the variable
conjugate to Φ and obeys , Q] = i~ in the quan-
tized theory. For a sufficiently homogeneous field, the
charge is given by Q = CV P/d, where P =
P
i
µσ
i
x
is the total polarization and d is the distance between
the capacitor plates. As usual we express Φ and Q in
terms of annihilation and creation operators a and a
as
Q = Q
0
(a + a
), Φ = iΦ
0
(a
a), (5)
where Q
0
=
p
~/(2Z), Φ
0
=
p
~Z/2 and Z =
p
L/C
is the cavity impedance. Altogether, we obtain the
canonical cavity QED Hamiltonian [33]
H
cQED
= ~ω
c
a
a + ~g(a + a
)S
x
+
~g
2
ω
c
S
2
x
+ H
dip
,
(6)
where ω
c
= 1/
LC and g = µQ
0
/(~Cd) is the cou-
pling strength.
The form of H
cQED
given in Eq. (6) allows a clear
distinction between electrostatic and dynamical ef-
fects. Here the terms J
ij
σ
i
x
σ
j
x
represent the electro-
static energy of the ensemble with a fixed orientation
of the dipoles. This energy might be modified in the
presence of metallic cavity mirrors [33, 42, 53], but it
is independent of the frequency or the vacuum field
amplitude of the dynamical mode. The coupling to
the dynamical field is then proportional to g and in-
cludes the collective dipole-field coupling as well as
the so-called P
2
-term S
2
x
[32, 33, 51]. This dis-
tinction shows that the regular Dicke model, which
is recovered for J
ij
= g
2
c
, describes a very spe-
cial case of a dipolar system with attractive all-to-all
dipole-dipole interactions. Although such a scenario
can be realized in circuit QED [54], the analysis of
this specific model does not provide much insights on
the behavior of more general cavity QED systems.
2.2 Observables
Apart from H
cQED
, which determines the dynamics
and the equilibrium states of the system, it is also
important to identify the relevant measurable observ-
ables. For the dipoles, quantities of interest are the
population imbalance, hS
z
i, or the polarization along
the cavity field, hS
x
i hPi, etc. Since the operator Q
for the total charge includes the induced charges from
the dipoles, its value is typically not directly measur-
able. Therefore, the relevant observables for the cav-
ity mode are the magnetic flux Φ and the voltage drop
V (see Fig. 1) and it is convenient to introduce the
displaced photon annihilation operator
A = a +
g
ω
c
S
x
, [A, A
] = 1. (7)
With the help of this definition we obtain [33, 55]
V = V
0
(A + A
), Φ = iΦ
0
(A
A), (8)
where V
0
= Q
0
/C, and the total Hamiltonian can be
written in a compact form as
H
cQED
= ~ω
c
A
A + H
dip
. (9)
By comparing with Eq. (1), we see that the expecta-
tion value of hA
Ai can be interpreted as the energy
of the dynamical cavity mode in units of ~ω
c
. This is
in contrast to the conventional photon number ha
ai,
which depends on the chosen gauge [33, 55] and has
no simple interpretation in a strongly coupled cavity
QED system. Note, however, while A + A
, A
A, etc.
represent physical properties of the cavity mode only,
the operators A and A
do not commute with all the
dipole operators and on a formal level we must still
use a and a
to represent the independent cavity de-
gree of freedom.
While we focus here on a lumped-element realiza-
tion of the EM mode as an explicit example, the model
in Eq. (4) and all the results discussed in this work
can be readily applied to arbitrary cavity QED sys-
tems using the replacements [33, 41, 42, 51, 55, 56]
V E, Q D, Φ B. (10)
Here E, D and B are the operators for the electric
field, the displacement field and the magnetic field,
respectively. For a detailed derivation and justifica-
tion of Hamiltonian (6) in dipolar cavity QED and
circuit QED settings, see Refs. [32, 33, 36].
Accepted in Quantum 2020-09-15, click title to verify. Published under CC-BY 4.0. 3
3 The free energy in cavity QED
By assuming that the cavity QED system is weakly
coupled to a large reservoir of temperature T , the re-
sulting equilibrium state of the system is
ρ
th
=
1
Z
e
βH
cQED
, (11)
where β = 1/(k
B
T ) and Z = Tr{e
βH
cQED
} is the
partition function. In this case the central quantity
of interest is the free energy,
F = k
B
T ln(Z) = F
0
c
+ F
0
dip
+ F
g
, (12)
which we divide into the free energies F
0
c
and F
0
dip
of
the decoupled subsystems and a remaining contribu-
tion F
g
. In the following discussion we will mainly
focus on the coupling-induced part of the free energy,
F
g
, which allows us to separate the influence of light-
matter interactions from the thermodynamic proper-
ties of the bare subsystems.
For small and moderately large ensembles the par-
tition function Z and the resulting free energy can be
evaluated by diagonalizing H
cQED
numerically. For
the LMG model, which conserves the total spin s,
this can be done for each spin sector separately and
we obtain
Z =
X
s
ζ
s,N
Z
s
. (13)
Here Z
s
is the partition function of H
cQED
con-
strained to a total spin s and [45]
ζ
s,N
=
N!(2s + 1)
(
N
2
s)!(
N
2
+ s + 1)!
(14)
accounts for the multiplicity of the respective mul-
tiplet due to permutation symmetry. For small and
moderate temperatures, we use this approach to eval-
uate the exact free energy for systems with N
1 100 dipoles. More details about the numerical
calculations are given in Appendix A.
3.1 Mean-field theory
In the analysis of the regular Dicke model [45, 46, 47]
with N 1, a frequently applied approximation for
evaluating the free energy is based on the mean-field
decoupling of the dipole-field interaction,
(a + a
)S
x
(α + α
)S
x
+ (a + a
x
(α + α
x
,
(15)
where the expectation values α = hai and Σ
x
=
hS
x
i must be determined self-consistently. Under
this mean-field approximation Hamiltonian (6) can be
written as the sum of two independent parts,
H
MF
= ~ω
c
a
a + ~g(a + a
x
+ H
MF
dip
(α, Σ
x
), (16)
where H
MF
dip
(α, Σ
x
) = H
dip
+ ~g(α + α
)(S
x
Σ
x
) +
~g
2
S
2
x
c
. The first two terms describe the energy
0.01 0.1
1
0.01 0.1
1
10
0.01 0.1
1
10
0
0.1
0.2
0.3
0.4
0.01 0.1
1 10
0.2
0.1
0
0.2
0.1
0
0.2
0.1
0
exact
1
2
3
10
-3
10.0 1.0 10.0 1.0
10
-3
0
2
4
-3
10
10.0 1.0
0
1
3
2
-3
10
10.0 1.0
0
1
3
2
10
Figure 2: Dependence of the coupling-induced part of the
free energy, F
g
, on the cavity-dipole coupling strength, g.
This dependence is shown in the individual plots for different
temperatures and dipole-dipole coupling strengths, J, and for
ω
c
= ω
0
. In each plot the exact numerical results for N = 20
dipoles are compared with approximate results obtained from
mean-field theory (F
MF
g
), second-order perturbation theory
(F
(2)
g
) and from a variational calculation (F
V
g
).
of a displaced oscillator, which is minimized for α =
(g
c
x
. With the help of this relation between
α and Σ
x
, the total partition function in mean-field
approximation is given by
Z
MF
x
) = Z
0
c
×
¯
Z
MF
dip
x
). (17)
Here, the first factor is the partition function of the
bare cavity and
¯
Z
MF
dip
x
) = Tr{exp[β
¯
H
MF
dip
x
)]} is
the partition function of an ensemble of dipoles with
effective Hamiltonian (which includes the constant en-
ergy shift from the displaced oscillator)
¯
H
MF
dip
x
) = H
dip
+
~g
2
ω
c
(S
x
Σ
x
)
2
. (18)
The free energy for the whole system in mean-field
approximation is then given by
F
MF
= min
Σ
x
{−k
B
T ln [Z
MF
x
)]}, (19)
and F
MF
g
= F
MF
F
0
c
F
0
dip
are the correspond-
ing coupling-induced corrections. Note that the min-
imization of the free energy also ensures that the self-
consistency condition hS
x
i = Σ
x
is satisfied.
Equation (18) shows that cavity-induced correc-
tions to the thermodynamic properties of a dipolar
system are only affected by fluctuations, but not by
the mean orientation of the dipoles. Therefore, by ap-
plying a second mean-field decoupling for the dipoles
(see Appendix B), the effect of the cavity vanishes
completely and F
MF
g
= 0. We conclude that a full
Accepted in Quantum 2020-09-15, click title to verify. Published under CC-BY 4.0. 4
mean-field treatment, as frequently employed to study
the ground states and thermal phases of the Dicke
model or of collective spin models [57, 58, 59], can-
not be used to analyze the thermodynamics of actual
cavity QED systems.
To take fluctuations of the dipoles into account we
can evaluate the partition function of the spin sys-
tem,
¯
Z
MF
dip
x
), exactly. In this case we find that
in the paraelectric phase, where Σ
x
= 0, the cav-
ity induces a renormalization of the interaction term,
J J + g
2
N
c
. This renormalization becomes a
substantial modification of the dipolar system already
in the collective USC regime, g
N
ω
0
ω
c
, and
could, at first sight, even prevent the ferroelectric in-
stability for J < 0. While such a shift of the phase
transition point is not observed when a proper min-
imization over Σ
x
is carried out, a comparison with
the exact free energy in Fig. 2 shows that mean-field
theory systematically overestimates the influence of
the cavity mode, in particular at higher temperatures.
This somewhat counterintuitive trend can be traced
back to the fact that the mean-field decoupling in
Eq. (15) neglects contributions which are second order
in ~g(a + a
)S
x
and scale approximately as
H
(2)
g
~g
2
ω
c
S
2
x
. (20)
Therefore, the mean-field decoupling neglects an es-
sential contribution from the light-matter interaction
and the approximation becomes uncontrolled.
As shown in the examples in Fig. 2, at larger cou-
pling parameters g
c
and low temperatures, the
mean-field predictions agree reasonably well with the
exact results. However, this agreement seems to
be accidental since at higher temperatures there are
again substantial deviations and the limit g
c
1
is not reproduced correctly. Although not shown ex-
plicitly, a very similar trend is also found in the anti-
ferroelectric case, J > 0. In summary, these results
for small and large couplings indicate that also on a
qualitative level Hamiltonian
¯
H
MF
dip
x
) does not cor-
rectly capture the influence of the cavity mode.
3.2 Collective USC regime
Many cavity QED experiments are operated in the
regime G = g
N . ω
c
and N 1, where the collec-
tive coupling G can become comparable to the photon
frequency, but the coupling of each individual dipole
to the cavity mode is still very small, g ω
c
. In this
regime, we can treat the dipole-field interaction,
H
g
= ~g(a + a
)S
x
+
~g
2
ω
c
S
2
x
, (21)
as a small perturbation and expand the free energy
in powers of g. As a result of this calculation, which
is detailed in Appendix C, we obtain the lowest-order
correction to the bare free energy. It can be written
in the form
F
(2)
g
= N
~g
2
4ω
c
f
g
. (22)
The dimensionless function f
g
O(1) contains two
contributions, one arising from the average value of
the S
2
x
term and a second-order contribution from
the linear coupling term, ~g(a
+ a)S
x
. The result-
ing expression for f
g
still involves non-trivial correla-
tion functions of spin operators, which for interact-
ing dipoles must be evaluated numerically. For non-
interacting dipoles this calculation can be carried out
analytically and we obtain the explicit result
f
g
(J = 0) =
ω
2
0
ω
0
ω
c
tanh
~ω
0
2k
B
T
coth
~ω
c
2k
B
T
(ω
2
0
ω
2
c
)
.
(23)
In Fig. 2, this prediction from perturbation theory is
compared to the exact free energy and we find that
the cavity-induced corrections to the free energy are
very accurately reproduced by F
(2)
g
at low and high
temperatures, even for collective coupling strengths
as large as G ω
c
.
By taking the limit T 0, Eq. (23) provides us
directly with the lowest-order correction to the ground
state energy of a cavity QED system [1],
E
(2)
0
= F
(2)
g
(T 0, J = 0) = N
~g
2
4ω
c
ω
0
ω
0
+ ω
c
, (24)
which agrees with Hamiltonian perturbation theory.
In the opposite high-temperature limit we obtain
F
(2)
g
(T ) ' N
~
3
g
2
ω
2
0
48ω
c
k
2
B
T
2
. (25)
Therefore, the cavity-induced corrections to the free
energy vanish quadratically with increasing tempera-
ture. Importantly, we find that for all temperatures
F
(2)
g
0, which is in stark contrast to the negative
correction terms obtained within the framework of the
regular Dicke model [19]. More generally, one can
show that also for an arbitrary system of interacting
two-level dipoles (see Appendix C)
0 f
g
<
4(∆S
x
)
2
N
, (26)
where (∆S
x
)
2
= hS
2
x
i
0
hS
x
i
2
0
. For a regular dipolar
system away from a critical point, the spatial extent of
the individual correlations, hσ
i
x
σ
j
x
ihσ
i
x
ihσ
j
x
i, is finite
and the collective fluctuations scale as (∆S
x
)
2
N.
Therefore, under very generic conditions, by taking
the limit N with G kept fixed, we obtain
lim
N→∞
F
(2)
g
N
= 0. (27)
This result confirms our basic intuition that the cou-
pling of many dipoles to a single mode should not af-
fect extensive thermodynamic properties. Note that
Accepted in Quantum 2020-09-15, click title to verify. Published under CC-BY 4.0. 5
this conclusion does not necessarily hold for a stronger
scaling of fluctuations, i. e. (∆S
x
)
2
N
2
. This scal-
ing can be found, for example, in the LMG model for
J < ω
0
when the symmetry of the ferroelectric phase
is not explicitly broken. However, even in this special
case our exact numerical calculations confirm that the
bound f
g
< 1 still holds.
3.3 Non-perturbative regime
The physics of cavity QED changes drastically once
g ω
c
and the light-matter interactions become non-
perturbative at the level of individual dipoles. To
analyze this regime, it is usually more convenient to
change to a polaron frame,
˜
H
cQED
= UH
cQED
U
, via
the unitary transformation U = e
g
ω
c
S
x
(a
a)
[60, 61,
62]. In this frame the cavity QED Hamiltonian can
be written as [32, 33, 38, 63]
˜
H
cQED
= ~ω
c
a
a + H
dip
+ H
int
, (28)
where the interaction part now takes the form
H
int
= ~ω
0
US
z
U
S
z
. (29)
An immediate benefit of the polaron representation
is that the interaction is proportional to ω
0
. This
shows that for ω
0
0 the coupling to the dynamical
mode vanishes and we recover the electrostatic limit,
˜
H
cQED
(ω
0
0) = ~ω
c
a
a +
P
i,j
~J
ij
4
σ
i
x
σ
j
x
.
For finite ω
0
the effects of H
int
are more involved.
For T = J = 0 it can be shown that up to second order
in H
int
and for g
c
& 2, the low-energy behavior of
the dipolar system is well-described by the effective
spin Hamiltonian [32]
H
eff
= ~ω
0
e
g
2
2ω
2
c
1
S
z
+
~ω
2
0
ω
c
2g
2
S
2
x
~
S
2
.
(30)
This effective model captures two important signa-
tures of non-perturbative light-matter interactions,
which will be relevant for the discussions below.
Firstly, there is a strong suppression of the dipole
oscillation frequency when g
c
& 1. Secondly, the
cavity mediates an all-to-all anti-ferroelectric coupling
ω
2
0
.
In principle, we can again apply perturbation the-
ory to evaluate F
g
up to second order in H
int
and ex-
tend the results from Sec. 3.2 into the strong-coupling
regime. However, such an approach is only reliable
when ω
0
ω
c
and the resulting expressions are much
more involved. Therefore, this method is only briefly
summarized in Appendix C. As a less accurate, but
more intuitive approach we can use the variational
principle of Bogoliubov to derive an upper bound F
V
for the free energy,
F F
+ h
˜
H
cQED
H
i
ρ
F
V
. (31)
Here ρ
is the thermal state and F
the corresponding
free energy for the trial Hamiltonian H
. Based on
(a) (b)
10
0
10
-2
10
-4
0.1
1
10
3
0
0.1
0.2
exact
variational
0 1 2
Figure 3: (a) Plot of the zero-field susceptibility χ
z
(solid
lines) for different coupling parameters g
c
. The dashed
lines indicate the predictions from the approximate formula
given in Eqs. (37)-(39). The x-markers show the results ob-
tained from the perturbation theory discussed in Appendix
C.3. (b) Dependence of the Curie constant α
C
(g) on the
dipole-field coupling strength. The exact numerical results
are in perfect agreement with the analytic scaling derived in
Eq. (38) from the variational free energy, F
V
. For all plots
N = 20, ω
c
= ω
0
and J = 0 have been assumed.
the discussion above we choose
H
= ~ω
c
a
a + ~˜ω
0
S
z
+
~J
N
S
2
x
, (32)
which describes a non-interacting cavity QED sys-
tem, but with a variable frequency ˜ω
0
. By minimizing
F
V
with respect to ˜ω
0
for each g we obtain (see Ap-
pendix D)
˜ω
0
(g) = ω
0
e
g
2
2ω
2
c
(1+2N
th
)
, (33)
where N
th
= 1/(e
β~ω
c
1). While from the compari-
son in Fig. 2 we see that overall F
V
g
= F
V
F
0
c
F
0
dip
does not reproduce the quantitative behavior of F
g
very accurately, we will see in the following that there
are still many cavity-induced modifications that can
be directly explained by this simple renormalization
of the dipole frequency.
4 Para- and ferroelectricity in the USC
regime
While the free energy contains all the relevant infor-
mation about the cavity QED system, we are usually
interested in derivative quantities such the suscepti-
bility, the specific heat, etc., or the existence of dif-
ferent phases and the transitions between them. To
understand in which way the coupling to a quantized
cavity mode can influence such quantities, we discuss
in this section three elementary examples.
4.1 USC modifications of the Curie law
As a first example it is instructive to consider the
simplest case of non-interacting dipoles, where
H
dip
= ~ω
0
S
z
. (34)
Accepted in Quantum 2020-09-15, click title to verify. Published under CC-BY 4.0. 6
This means, that in the absence of the cavity, the
dipoles form an ideal paraelectric material. For this
system, we are interested in the dependence of the
population imbalance hS
z
i on the level splitting ω
0
.
Specifically, we consider the limit ω
0
0, where we
obtain the zero-field susceptibility
χ
z
=
1
N
hS
z
i
ω
0
ω
0
=0
=
1
N~
2
F
ω
2
0
ω
0
=0
(35)
from the second derivative of the free energy. In the
limit of a vanishing bias field the susceptibility follows
the usual Curie law
χ
z
|
g0
(T ) =
α
C
T
, (36)
with a Curie constant α
C
= ~/(4k
B
). In the context
of cavity QED, the behavior of this quantity for finite
g is of particular importance. Since the dipoles decou-
ple from the cavity mode for ω
0
= 0, the zero-field sus-
ceptibility captures the lowest-order deviations from
the electrostatic limit.
In Fig. 3(a) we plot χ
z
(T ) for a cavity QED system
with J = 0, N = 20 and different coupling strengths
g. For small g we still recover the 1/T behavior with
a small reduction of the Curie constant. In the non-
perturbative regime, g
c
& 1, the modifications be-
come more significant. Although in this regime the
susceptibility still diverges for T 0, there appears
an additional plateau for an intermediate range of
temperatures, T < ~ω
c
/k
B
. To understand this be-
havior we approximate the susceptibility by the two
dominant contributions,
χ
z
α
C
(g)
T
1
N~
X
n
p
n
2
E
n
ω
2
0
ω
0
=0
, (37)
where p
n
and E
n
are the thermal occupation prob-
abilities and the energy of the n-th eigenstate, re-
spectively. The first term emerges from the change
of hS
z
i due to small changes of the thermal popula-
tions when ω
0
is varied. Since we are interested in
the limit ω
0
0, this change results in the same
high-temperature scaling 1/T as in the case of free
dipoles. This effect is already captured by the vari-
ational free energy F
V
discussed in Sec. 3.3, from
which we obtain
α
C
(g) '
~
4k
B
e
g
2
ω
2
c
(1+2N
th
)
. (38)
In Fig. 3(b) we show that this analytic result is in
perfect agreement with the low-temperature limit of
χ
z
obtained from exact numerical simulations.
The second term in Eq. (37) is the contribution to
the susceptibility, which arises from quadratic changes
of the energy eigenstates with varying ω
0
. For free
dipoles, H
dip
ω
0
and therefore this contribution is
absent in regular paramagnetic and paraelectric sys-
tems. However, as evident from the effective spin
(b)
0.2
-0.2
-0.6
-0.4
0
(a)
1
2
0 2 4 6 8
0
0.2
0.4
-0.2
Figure 4: (a) Plot of the dimensionless quantity c
g
defined
in Eq. (41), which determines the corrections to the specific
heat, C
g
/N, in the collective USC regime. (b) Dependence of
C
g
on the coupling strength g for a fixed temperature and an
increasing number of dipoles, N = 10, 20, 30, 40, 50. For this
plot ω
0
= 0.5ω
c
and a photon cutoff number of N
ph
= 100
has been used.
model in Eq. (30), for couplings g > ω
c
the energy lev-
els show indeed a quadratic scaling, E
n
ω
2
0
. From
this effective model and by assuming that all spin lev-
els are equally populated, p
n
1/2
N
, we obtain the
approximate result
1
N~
X
n
p
n
2
E
n
ω
2
0
ω
0
=0
ω
c
2g
2
. (39)
As shown in Fig. 3(a), this estimate is in very good
agreement with the value of the plateau of χ
z
(T )
found in exact numerical simulations. For T &
~ω
c
/k
B
, the thermal population of the photon states
is no longer negligible and the approximate model in
Eq. (30) breaks down. Beyond this point, which for
large g is indicated by a small bump, the regular Curie
law is approximately recovered.
Note that since the susceptibility χ
z
is evaluated at
ω
0
= 0, it can be calculated exactly from a second-
order perturbation theory in H
int
ω
0
in the polaron
representation (see Appendix C). Although the result-
ing expressions must still be evaluated numerically,
this method allows us to obtain χ
z
from correlation
functions of the dipolar system only. Therefore, this
approach can be very useful for performing similar
calculations for more complicated dipolar systems.
4.2 USC modifications of the specific heat
A second quantity of general interest in statistical
physics is the heat capacity,
C = T
2
F
T
2
= C
0
c
+ C
0
dip
+ C
g
. (40)
For a decoupled system, the heat capacities of the cav-
ity and the dipoles, C
0
c
and C
0
dip
, are both bounded and
scale as C
0
c
k
B
and C
0
dip
N k
B
, respectively. This
scaling suggests that for a large ensemble of dipoles
and similar energy scales, ω
c
ω
0
, the presence of
a single cavity degree of freedom should have a neg-
ligible contribution to the specific heat C/N of the
Accepted in Quantum 2020-09-15, click title to verify. Published under CC-BY 4.0. 7
combined system. This can be shown explicitly in
the collective USC regime, where for non-interacting
dipoles
C
g
k
B
N
'
~g
2
4ω
c
(k
B
T )
× c
g
(T, ω
c
, ω
0
). (41)
Here, c
g
= (k
B
T )
2
2
f
g
/∂(k
B
T )
2
is a dimensionless
function, which is independent of N and which is plot-
ted in Fig. 4(a) for different frequency ratios, ω
0
c
.
Therefore, taking the limit N for a fixed G, the
corrections to the specific heat vanish.
In Fig. 4(b) we plot C
g
/N for a cavity QED sys-
tem with an increasing number of N non-interacting
dipoles. For small couplings, G/ω
c
. 0.5, the
correction is accurately reproduced by the analytic
weak-coupling result given in Eq. (41). In the non-
perturbative regime we observe substantial modifica-
tions. On a qualitative level, these corrections can be
understood from a cavity-induced suppression of ω
0
,
but overall we find that the dependence of C
g
is not
very accurately reproduced by the variational ansatz
in Eq. (32) or any of the other approximation schemes.
However, from the exact numerical results plotted in
Fig. 4(b) we see that the maximal correction to the
specific heat is
|C
g
|
k
B
N
O(1), (42)
and shifts, but does not vanish with increasing N.
Combined with the behavior of the susceptibility dis-
cussed above, this finding demonstrates that in the
non-perturbative regime the coupling to a single dy-
namical field degree of freedom can have a substantial
influence on extensive thermodynamic quantities of a
large ensemble of dipoles.
4.3 USC modifications of the ferroelectric
phase transition
A central topic of interest in the field of USC cavity
QED is the so-called superradiant phase transition,
which is predicted for the ground and thermal equilib-
rium states of the standard Dicke model [45, 46, 47].
While in more accurate models for light-matter in-
teractions this transition does not occur for non-
interacting dipoles [28, 29, 30, 31, 32, 33, 34, 35],
the system can still undergo a regular ferroelectric
phase transition in the case of attractive electrostatic
interactions, J
ij
< 0. Within the framework of the
LMG model, such a transition is well-described by
a mean-field decoupling of the collective interaction
term, S
2
x
S
x
Σ
2
(see Appendix B), from which
one can derive the general relation between the criti-
cal coupling strength J
c
and the critical temperature
T
c
[57, 58, 59],
tanh
~ω
0
2k
B
T
c
=
ω
0
J
c
. (43)
(a)
(b)
-3
-2
-1
0
3
4
5
6
7
8
9
10
0
0.1
(I)
(II)
0-1 -0.5 0.5
0.6
0.7
0.8
(c)
10.1 10
-10 -5 50
0
0.1
10-10 -5 50
(I)
(II)
-2
-1
0
0-1 -0.5 0.5
Figure 5: (a) Phase diagram of the LMG model without
cavity, where the color scale shows the value of the parameter
¯m =
p
hS
2
x
i [59] for N = 20 dipoles. For each point, we also
evaluate the probability distribution p(m
x
) for the projection
quantum number m
x
, which exhibits a single maximum in
the paraelectric phase (I) and two maxima in the ferroelectric
phase (II). The transition between the single- and bi-modal
distribution is indicated by the solid line, while the dotted
line depicts the phase boundary obtained from mean-field
theory, see Eq. (43). The same boundaries are shown in (b)
for different coupling strengths g, where for the mean-field
results ω
0
has been replaced by ˜ω
0
. (c) Dependence of the
critical temperature T
c
on the coupling parameter g
c
for
a fixed inter-dipole coupling strength of J
c
= 1.5. In all
plots ω
0
= ω
c
and N = 20.
For T 0 there is a critical coupling strength J
c
=
ω
0
, beyond which the dipoles enter a ferroelectric
phase with hS
x
i 6= 0. For ω
0
0 this phase only
exists below a critical temperature T
c
= ~J/(2k
B
),
which is just the transition temperature of the classi-
cal Ising model. For arbitrary ω
0
, the phase boundary
of the LMG model in the limit N is indicated
by the dashed line in Fig. 5(a).
Since symmetry breaking does not occur for finite
N, the criterion hS
x
i 6= 0 cannot be used to charac-
terize the ferroelectric phase in exact numerical cal-
culations. In Fig. 5(a) we show instead the quantity
¯m =
p
hS
2
x
i, which provides a good indicator for the
ferroelectric phase of the LMG model [59]. However,
for the rather small numbers of dipoles assumed in
the simulations of the full cavity QED model below,
the variation of ¯m around the phase transition line
is still rather smooth. Therefore, for the following
analysis we consider instead the probability distribu-
tion p(m
x
) = Tr{
m
x
ρ
th
}, where
m
x
=
P
s
P
s,m
x
and P
s,m
x
is the projector on all states with S
x
|ψi =
m
x
|ψi and total spin s. In this case, we can define
the phase boundary as the line where this function
changes from a single to a bi-modal distribution, as
illustrated in Fig. 5(a). For the bare LMG model,
Accepted in Quantum 2020-09-15, click title to verify. Published under CC-BY 4.0. 8
this approach reproduces very accurately the phase
boundary derived from mean-field theory, even for a
small number of N = 20 dipoles.
When the dipoles are coupled to the cavity mode,
a finite polarization hS
x
i 6= 0 is naturally associ-
ated with a non-vanishing expectation value of hai '
g
c
hS
x
i, similar to what is expected for the su-
perradiant phase in the Dicke model. Note, however,
that this expectation value is real and corresponds to
a finite charge (or displacement field) hQi ha + a
i.
The relevant cavity observables, hV i hA + A
i and
hΦi ihA
Ai, are not affected by this transi-
tion [33, 41]. For the ground state, it has further
been shown that in the collective USC regime, i.e.
when G
ω
0
ω
c
but g ω
c
, also the transition
point is not influenced by the coupling to the dynam-
ical cavity mode [33]. This is no longer the case when
g ω
c
.
In the current study we are primarily interested
in USC effects beyond the ground state and show
in Fig. 5(b) the coupling-induced modification of the
phase boundary in the whole T J plane for differ-
ent values of g
c
. In this plot, the exact analytic
results are compared with a modified mean-field the-
ory, where in Eq. (43) the bare dipole frequency ω
0
is replaced by the renormalized frequency ˜ω
0
given in
Eq. (33). From this comparison we find that the vari-
ational free energy F
V
captures the overall trend, al-
though the actual phase transition line deviates from
the exact results, in particular for g
c
> 1 and for
low temperatures. In Fig. 5(c) we fix the value of J
and plot the dependence of the critical temperature
on the coupling strength g. Consistent with the other
examples above, we observe only minor corrections
for G . ω
c
, but a substantial increase of T
c
for cou-
plings g
c
& 1. This means that in this coupling
regime the presence of the cavity mode stabilizes the
ferroelectric phase against thermal fluctuations. This
behavior is qualitatively reproduced by the modified
mean-field ansatz.
5 Black-body radiation
The emission spectrum of a hot body was one of the
first examples that could not be explained by com-
bining the otherwise very successful theories of sta-
tistical mechanics and electromagnetism. In the cor-
rect quantum statistical derivation of the black-body
spectrum it is assumed that the EM field thermal-
izes through weak interactions with the material, but
that it can be treated otherwise as a set of indepen-
dent harmonic modes. Therefore, it is particularly
interesting to see how the thermal emission spectrum
of a cavity mode changes under strong and USC con-
ditions [64, 65, 66].
(a)
(b)
1
0
2
3
0 2 8
1
0
2
3
1
0
2
3
4 6 0 2 84 6 0 2 84 6
0 2 84 6 0 2 84 6 0 2 84 6
1
0
2
1
0
2
1
0
2
-1
0
1
2
Figure 6: (a) The black-body spectrum S
bb
(ω)/(~κ) is plot-
ted on a logarithmic scale as a function of g for three dif-
ferent temperatures. The green dashed lines indicate the
frequencies ω
±
of the two polariton modes obtained from
a Holstein-Primakoff approximation. (b) Plot of the total
emitted power, P
rad
, and the average value of the EM en-
ergy, hH
em
i/(~ω
c
) = hA
Ai+ 1/2, for the same parameters.
Note that for better visibility we have included for this plot
the offset ~ω
c
/2 (indicated by the dashed line) into the def-
inition of H
em
. For all plots N = 6, J = 0 and γ
c
= 0.04
have been assumed.
5.1 Power spectral density
In the setup shown in Fig. 1, the black-body spectrum
can be measured, for example, by coupling the cavity
via a weak capacitive link to a cold transmission line.
The emitted power will then be proportional to the
fluctuations of the voltage operator V = V
0
(A + A
)
(see also Ref. [55]). By assuming that the transmis-
sion line can be modeled as an Ohmic bath and that
the capacitive link is sufficiently weak, we can write
the spectrum of the emitted black-body radiation as
(see Appendix E)
S
bb
(ω) =
~κγ
2πω
c
X
n>m
e
βE
n
Z
ω
2
nm
|hE
n
|A + A
|E
m
i|
2
(ω ω
nm
)
2
+ γ
2
/4
,
(44)
where ω
nm
= (E
n
E
m
)/~ are the transition fre-
quencies between the eigenstates |E
n
i of H
cQED
with
energies E
n
. In Eq. (44), κ denotes the decay rate of
the bare cavity into the transmission line. In addi-
tion, we have introduced a phenomenological rate γ
to account for a small but finite thermalization rate
with the surrounding bath. For consistency we re-
quire κ γ and γ |E
n
E
m
|/~, but otherwise the
precise values of κ and γ are not important.
In Fig. 6(a) we plot S
bb
(ω) as a function of the cou-
pling strength g and for different temperatures. For
small couplings, g ω
c
, we see the expected split-
ting of the unperturbed cavity resonance into two po-
laritonic resonances at frequencies ω
±
ω
c
± G/2.
Although the lower polariton mode has a higher ther-
mal occupation, the upper branch is slightly brighter.
This observation can be partially explained by the
Accepted in Quantum 2020-09-15, click title to verify. Published under CC-BY 4.0. 9
scaling of the emission rate ω
2
±
, but a more detailed
analysis is presented below. At intermediate coupling
strengths and temperatures the spectrum becomes
rather complex. This is related to the large spread
of the eigenenergies E
n
for these coupling values (see,
for example, Fig. 1 in Ref. [67]) and the fact that the
dipoles and photons are still strongly hybridized. At
very large interactions, the spectrum collapses again
to a single line centered around the bare cavity fre-
quency. This collapse is a striking consequence of the
approximate factorization of the eigenstates at very
large interactions [32] and provides a clear signature
of the non-perturbative coupling regime, which can
be detected in the emitted radiation field. While
in Fig. 6(a) we plot the spectrum only for the case
of non-interacting dipoles, the same overall trend is
found independently of the dipole-dipole interaction
strength J.
Note that in previous studies of the absorption and
emission spectra of the EDM [51, 66] or the thermal
radiation spectrum of the Rabi model (N = 1) [65]
only moderate values of g have been considered, where
this spectral collapse does not yet occur. In the case of
the Rabi model, it has also been shown that the statis-
tics of the emitted photons can be sub-Poissonian for
a certain range of couplings and temperatures [64].
We don’t find such a behavior for larger numbers of
two-level dipoles.
5.2 Radiated power and EM energy
In Fig. 6(b) we also plot the total radiated power,
P
rad
, and compare it with the equilibrium value of
the EM field energy, hH
em
i. Here, the total power
is obtained from Eq. (44) by integrating over all fre-
quencies,
P
rad
=
~κ
ω
c
X
n>m
e
βE
n
Z
ω
2
nm
|hE
n
|A + A
|E
m
i|
2
. (45)
For an empty cavity, this expression reduces to P
0
rad
=
~ω
c
κN
th
, which we use to normalize the power val-
ues. Interestingly, for moderate temperatures we find
a rather counterintuitive behavior: While the average
energy that is stored in the mode increases for inter-
mediate couplings, the emitted power decreases at the
same time.
To explain this behavior we consider moderate val-
ues of G . ω
c
and low temperatures. In this case
we can use a Holstein-Primakoff approximation [68]
and replace the spin operators σ
i
by bosonic annihi-
lation operators c
i
. The resulting linearized Hamilto-
nian can then be diagonlized and written as H
cQED
'
H
HP
~ω
c
/2 N~ω
0
, where
H
HP
=
X
η=±
~ω
±
c
±
c
±
+
1
2
+
N1
X
k=1
~ω
0
c
k
c
k
+
1
2
.
(46)
(b)(a)
1 20 0 0.20.5 1.5
1
2
0
0.5
1.5
10.4 0.6 0.8
0.5
1
1.5
Figure 7: (a) The dimensionless matrix elements V
±
and Φ
±
,
which determine the decomposition of the voltage and the
magnetic flux operators in terms of the polariton operators
c
±
[see Eqs. (47) and (48)] are plotted as a function of the
collective coupling strength, G. (b) Plot of the ratio between
the total power emitted from the coupled cavity QED sys-
tem (P
rad
) and from the bare cavity (P
0
rad
) as a function
of temperature. The solid lines are obtained from exact nu-
merical calculations for N = 6 and the dashed lines show
the corresponding results predicted by Eq. (49) based on a
Holstein-Primakoff approximation.
Here the c
±
are bosonic operators for the two bright
polariton modes with frequencies ω
±
. The other
bosonic operators c
k
represent dark polaritons, i.e.,
collective excitations of the dipoles, which are decou-
pled from the cavity due to symmetry. In terms of
these polariton operators we can write
A
+ A = V
+
c
+
+ c
+
+ V
c
+ c
, (47)
(A
A) = Φ
+
c
+
c
+
+ Φ
c
c
, (48)
where the dimensionless matrix elements V
±
and Φ
±
are plotted in Fig. 7(a) as a function of the collective
coupling strength G.
Within the Holstein-Primakoff approximation only
the bright polariton modes contribute to the emission
spectrum and we obtain
P
rad
P
0
rad
' V
2
+
ω
2
+
ω
2
c
N
th
(ω
+
)
N
th
(ω
c
)
+ V
2
ω
2
ω
2
c
N
th
(ω
)
N
th
(ω
c
)
.
(49)
We see that there are various competing effects. With
increasing coupling G, the frequency ω
+
and the ma-
trix element V
2
+
for the upper polariton mode goes
up, while at the same time the corresponding mode
occupation, N
th
(ω
+
), gets exponentially suppressed.
The opposite is true for the lower polariton mode. As
shown in Fig. 7(b), this competition leads to a non-
monotonic influence of the light-matter coupling on
the radiated power. For temperatures k
B
T/(~ω
c
)
0.2 0.5, as considered in Fig. 6(b), Eq. (49) indeed
predicts the observed decrease in P
rad
for increasing
values of G. However, for higher and lower tempera-
tures the dependence can also be reversed. In particu-
lar, for k
B
T/(~ω
c
) 1 the occupation number of the
bare cavity mode is exponentially suppressed. How-
ever, for G ω
c
, also the lower polariton frequency is
strongly reduced and N
th
(ω
) k
B
T/(~ω
). Under
Accepted in Quantum 2020-09-15, click title to verify. Published under CC-BY 4.0. 10
such conditions we observe a huge coupling-induced
enhancement of the radiated power, P
rad
/P
0
rad
1.
By expressing also the EM energy, H
em
= ~ω
c
A
A,
in terms of the mode operators for the bright polariton
modes we obtain
hH
em
i
~ω
c
=
(V
2
+
+ Φ
2
+
)
2
N
th
(ω
+
) +
(V
2
+ Φ
2
)
2
N
th
(ω
)
+
V
2
+
+ Φ
2
+
+ V
2
+ Φ
2
4
1
2
.
(50)
We see that the prefactors for the thermal contribu-
tions in the first line of this equation have a much
weaker dependence on G. Further, we find that for
G & ω
c
and for temperatures k
B
T/(~ω
c
) . 0.5,
the main contribution to the EM energy comes from
a positive vacuum term given in the second line of
Eq. (50). This part does not contribute to the radi-
ated power such that overall P
rad
and H
em
display a
very different dependence on G.
6 Conclusions
In summary, we have analyzed the basic thermal prop-
erties of cavity QED systems in the USC regime. By
using various analytic approximations and exact nu-
merical results for a moderate number of two-level
dipoles, we have derived the coupling-induced correc-
tions to the free energy and some of its derivative
quantities. In the collective USC regime our analytic
results confirm the basic intuition that the coupling to
a single cavity mode cannot significantly change the
properties of a larger ensemble of dipoles. While a
large collective coupling strength G has a substantial
influence on the emission spectrum and also on the to-
tal radiated power, the corrections to material prop-
erties are small and scale only with the single-dipole
coupling strength as g
2
c
. This means that ma-
jor modifications of ground-state chemical reactions
or cavity-induced shifts of ferroelectric phase transi-
tions cannot be simply explained by a strong collec-
tive coupling to a single quantized mode. To identify
the detailed origin of such effects further experimen-
tal and theoretical investigations are still required, in
particular also on the influence of the metallic bound-
aries on electrostatic interactions [33, 42, 53, 69] and
other modifications of the background EM environ-
ment [15].
The behaviour of cavity QED systems changes
completely once the single-dipole coupling parame-
ter, g
c
, becomes of order unity. For this regime
we have shown in terms of several explicit exam-
ples that the coupling to a single cavity mode can
strongly modify material properties such as the spe-
cific heat, the susceptibility or the ferroelectric phase
transition temperature. While this regime is cur-
rently not accessible with atomic, molecular or ex-
citonic cavity QED systems, such coupling conditions
can be reached in superconducting quantum circuits,
where also the effect of temperature is more rele-
vant than in the optical regime. Apart from real-
izations with non-interacting [32] or collective ferro-
electric [54] qubits, such setups would also allow to
engineer H
cQED
with arbitrary local qubit-qubit in-
teractions by adding to those circuits additional ca-
pacitances or SQUID loops [70]. Therefore, such arti-
ficial light-matter systems constitute an ideal testbed
for studying strongly-coupled quantum systems with
unconventional thermodynamical properties. Poten-
tially, this can also lead to more accurate descrip-
tions of fundamental thermodynamical processes or
the optimization of quantum thermal machines, for
which cavity QED and collective spin models have al-
ready been considered as simple toy systems in the
past [71, 72, 73, 74, 75, 76].
Acknowledgements
This work was supported by the Austrian Academy
of Sciences (
¨
OAW) through a DOC Fellowship (D.D.)
and a Discovery Grant (P.P., P.R.) and by the Aus-
trian Science Fund (FWF) through the DK CoQuS
(Grant No. W 1210) and Grant No. P31701 (UL-
MAC).
References
[1] C. Ciuti, G. Bastard, and I. Carusotto, Quantum
vacuum properties of the intersubband cavity po-
lariton field, Phys. Rev. B 72, 115303 (2005).
[2] P. Forn-D´ıaz, L. Lamata, E. Rico, J. Kono, and
E. Solano, Ultrastrong coupling regimes of light-
matter interaction, Rev. Mod. Phys. 91, 025005
(2019).
[3] A. F. Kockum, A. Miranowicz, S. De Liberato,
S. Savasta, and F. Nori, Ultrastrong coupling be-
tween light and matter, Nat. Rev. Phys. 1, 19
(2019).
[4] Y. Todorov, A. M. Andrews, R. Colombelli, S.
De Liberato, C. Ciuti, P. Klang, G. Strasser,
and C. Sirtori, Ultrastrong Light-Matter Cou-
pling Regime with Polariton Dots, Phys. Rev.
Lett. 105, 196402 (2010).
[5] G. Scalari, C. Maissen, D. Turcinkova, D. Ha-
genm¨uller, 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, Sci-
ence 335, 1323 (2012).
[6] D. Dietze, A. M. Andrews, P. Klang, G. Strasser,
K. Unterrainer, and J. Darmo, Ultrastrong cou-
pling of intersubband plasmons and terahertz
metamaterials, Appl. Phys. Lett. 103, 201106
(2013).
Accepted in Quantum 2020-09-15, click title to verify. Published under CC-BY 4.0. 11
[7] C. R. Gubbin, S. A. Maier, and S. K´ena-Cohen,
Low-voltage polariton electroluminescence from
an ultrastrongly coupled organic light-emitting
diode, App. Phys. Lett. 104, 233302 (2014).
[8] Q. Zhang, M. Lou, X. Li, J. L. Reno, W. Pan,
J. D. Watson, M. J. Manfra, and J. Kono, Col-
lective non-perturbative coupling of 2D electrons
with high-quality-factor terahertz cavity pho-
tons, Nature Phys. 12, 1005 (2016).
[9] A. Bayer, M. Pozimski, S. Schambeck, D. Schuh,
R. Huber, D. Bougeard, and C. Lange, Terahertz
Light-Matter Interaction beyond Unity Coupling
Strength, Nano Lett. 17, 6340 (2017).
[10] B. Askenazi, A. Vasanelli, Y. Todorov, E. Sakat,
J.-J. Greffet, G. Beaudoin, I. Sagnes, and C. Sir-
tori, Midinfrared Ultrastrong Light-Matter Cou-
pling for THz Thermal Emission, ACS Photonics
4, 2550 (2017).
[11] T. Schwartz, J. A. Hutchison, C. Genet and T.
W. Ebbesen, Reversible Switching of Ultrastrong
Light-Molecule Coupling, Phys. Rev. Lett. 106,
196405 (2011).
[12] J. George, T. Chervy, A. Shalabney, E. Devaux,
H. Hiura, C. Genet, and T. W. Ebbesen, Mul-
tiple Rabi Splittings under Ultrastrong Vibra-
tional Coupling, Phys. Rev. Lett. 117, 153601
(2016).
[13] J. Flick, M. Ruggenthaler, H. Appel, and A.
Rubio, Atoms and Molecules in Cavities: From
Weak to Strong Coupling in QED Chemistry,
Proc. Natl. Acad. Sci. 114, 3026 (2017).
[14] R. F. Ribeiro, L. A. Martinez-Martinez, M. Du,
J. Campos-Gonzalez-Anguloand, and J. Yuen-
Zhou, Polariton chemistry: controlling molecu-
lar dynamics with optical cavities, Chem. Sci. 9,
6325 (2018).
[15] V. N. Peters, S. Prayakarao, S. R. Koutsares,
C. E. Bonner, and M. A. Noginov, Control of
Physical and Chemical Processes with Nonlocal
MetalDielectric Environments, ACS Photonics 6,
3039 (2019).
[16] J. A. Hutchison, T. Schwartz, C. Genet, E. De-
vaux, and T. W. Ebbesen, Modifying Chemi-
cal Landscapes by Coupling to Vacuum Fields,
Angew. Chem., Int. Ed. 51, 1592 (2012).
[17] A. Thomas, J. George, A. Shalabney, M.
Dryzhakov, S. J. Varma, J. Moran, T. Chervy,
X. Zhong, E. Devaux, C. Genet, J. A. Hutchi-
son, and T. W. Ebbesen, Ground-State Chemi-
cal Reactivity under Vibrational Coupling to the
Vacuum Electromagnetic Field, Angew. Chem.,
Int. Ed. 55, 11462 (2016).
[18] S. Wang, A. Mika, J. A. Hutchison, C. Genet,
A. Jouaiti, M. W. Hosseini, and T. W. Ebbe-
sen, Phase Transition of a Perovskite Strongly
Coupled to the Vacuum Field, Nanoscale 6, 7243
(2014).
[19] A. Canaguier-Durand, E. Devaux, J. George, Y.
Pang, J. A. Hutchison, T. Schwartz, C. Genet,
N. Wilhelms, J.-M. Lehn, and T. W. Ebbesen,
Thermodynamics of Molecules Strongly Coupled
to the Vacuum Field, Angew. Chem., Int. Ed.
52, 10533 (2013).
[20] M. H. Devoret, S. Girvin, and R. Schoelkopf,
Circuit-QED: How strong can the coupling be-
tween a Josephson junction atom and a trans-
mission line resonator be?, Ann. Phys. (NY) 16,
767 (2007).
[21] T. Niemczyk, F. Deppe, H. Huebl, E. P. Menzel,
F. Hocke, M. J. Schwarz, J. J. Garcia-Ripoll, D.
Zueco, T. H¨ummer, E. Solano, A. Marx, and R.
Gross, Circuit quantum electrodynamics in the
ultrastrong-coupling regime, Nature Phys. 6, 772
(2010).
[22] 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 Ultra-
strong Coupling Regime, Phys. Rev. Lett. 105,
237001 (2010).
[23] P. Forn-D´ıaz, J. J. Garc´ıa-Ripoll, B. Peropadre,
J.-L. Orgiazzi, M. A. Yurtalan, R. Belyansky, C.
M. Wilson, and A. Lupascu, Ultrastrong cou-
pling of a single artificial atom to an electromag-
netic continuum, Nature Phys. 13, 39 (2017).
[24] F. Yoshihara, T. Fuse, S. Ashhab, K.
Kakuyanagi, S. Saito, and K. Semba, Super-
conducting qubit-oscillator circuit beyond the
ultrastrong-coupling regime, Nature Phys. 13, 44
(2017).
[25] R. H. Dicke, Coherence in Spontaneous Radia-
tion Processes, Phys. Rev. 93, 99 (1954).
[26] T. Brandes, Coherent and collective quantum op-
tical effects in mesoscopic systems, Physics Re-
ports 408, 315 (2005).
[27] J. J. Hopfield, Theory of the Contribution of
Excitons to the Complex Dielectric Constant of
Crystals, Phys. Rev. 112, 1555 (1958).
[28] K. Rzazewski, K. Wodkiewicz, and W. Zakowicz,
Phase Transitions, Two-Level Atoms, and the A
2
Term, Phys. Rev. Lett. 35, 432 (1975).
[29] O. Viehmann, J. von Delft, and F. Marquardt,
Superradiant Phase Transitions and the Stan-
dard Description of Circuit QED, Phys. Rev.
Lett. 107, 113602 (2011).
[30] Y. Todorov and C. Sirtori, Intersubband polari-
tons in the electrical dipole gauge, Phys. Rev. B
85, 045304 (2012).
Accepted in Quantum 2020-09-15, click title to verify. Published under CC-BY 4.0. 12
[31] M. Bamba, and T. Ogawa, Stability of polariz-
able materials against superradiant phase transi-
tion, Phys. Rev. A 90, 063825 (2014).
[32] T. Jaako, Z.-L. Xiang, J. J. Garcia-Ripoll, and
P. Rabl, Ultrastrong coupling phenomena be-
yond the Dicke model, Phys. Rev. A 94, 033850
(2016).
[33] D. De Bernardis, T. Jaako, and P. Rabl,
Cavity quantum electrodynamics in the non-
perturbative regime, Phys. Rev. A 97, 043820
(2018).
[34] V. Rokaj, D. M. Welakuh, M. Ruggenthaler, and
A. Rubio, Light-matter interaction in the long-
wavelength limit: no ground-state without dipole
self-energy, J. Phys. B: At. Mol. Opt. Phys. 51,
034005 (2018).
[35] G. M. Andolina, F. M. D. Pellegrino, V. Gio-
vannetti, A. H. MacDonald, and M. Polini, Cav-
ity quantum electrodynamics of strongly corre-
lated electron systems: A no-go theorem for pho-
ton condensation, Phys. Rev. B 100, 121109(R)
(2019).
[36] D. De Bernardis, P. Pilar, T. Jaako, S. De Liber-
ato, and P. Rabl, Breakdown of gauge invariance
in ultrastrong-coupling cavity QED, Phys. Rev.
A 98, 053819 (2018).
[37] A. Stokes and A. Nazir, Gauge ambiguities im-
ply Jaynes-Cummings physics remains valid in
ultrastrong coupling QED, Nat. Commun. 10,
499 (2019).
[38] O. Di Stefano, A. Settineri, V. Macri, L.
Garziano, R. Stassi, S. Savasta, and F. Nori,
Resolution of gauge ambiguities in ultrastrong-
coupling cavity QED, Nature Phys. 15, 803
(2019).
[39] M. Roth, F. Hassler, and D. P. DiVincenzo, Op-
timal gauge for the multimode Rabi model in cir-
cuit QED, Phys. Rev. Research 1, 033128 (2019).
[40] Y. A. Kudenko, A. P. Slivinsky, and G. M. Za-
slavsky, Interatomic Coulomb interaction influ-
ence on the superradiance phase transition, Phys.
Lett. A 50, 411 (1975).
[41] J. Keeling, Coulomb interactions, gauge invari-
ance, and phase transitions of the Dicke model,J.
Phys: Cond. Mat. 19, 295213 (2007).
[42] A. Vukics and P. Domokos, Adequacy of the
Dicke model in cavity QED: A counter-no-go
statement, Phys. Rev. A 86, 053807 (2012).
[43] T. Grießer, A. Vukics, and P. Domokos, Depolar-
ization shift of the superradiant phase transition,
Phys. Rev. A 94, 033815 (2016).
[44] A. Stokes and A. Nazir, Uniqueness of the phase
transition in many-dipole cavity QED systems,
arXiv:1905.10697 (2019).
[45] K. Hepp and E. H. Lieb, On the superradiant
phase transition for molecules in a quantized ra-
diation field: the Dicke maser model, Ann. Phys.
76, 360 (1973).
[46] Y. K. Wang, and F. T. Hioe, Phase Transition in
the Dicke Model of Superradiance, Phys. Rev. A
7, 831 (1973).
[47] H. J. Carmichael, C. W. Gardiner, and D. F.
Walls, Higher order corrections to the Dicke su-
perradiant phase transition, Phys. Lett. A 46, 47
(1973).
[48] 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).
[49] J. A. Cwik, P. Kirton, S. De Liberato, and J.
Keeling, Excitonic spectral features in strongly
coupled organic polaritons, Phys. Rev. A 93,
033840 (2016).
[50] L. A. Martinez-Martinez, R. F. Ribeiro, J.
Campos-Gonzalez-Angulo, and J. Yuen-Zhou,
Can Ultrastrong Coupling Change Ground-State
Chemical Reactions?, ACS Photonics 5, 167
(2018).
[51] Y. Todorov and C. Sirtori, Few-Electron Ultra-
strong Light-Matter Coupling in a Quantum LC
Circuit, Phys. Rev. X 4, 041031 (2014).
[52] H. J. Lipkin, N. Meshkov, and A. J. Glick, Va-
lidity of many-body approximation methods for
a solvable model: (I). Exact solutions and per-
turbation theory, Nucl. Phys. 62, 188 (1965).
[53] M. Schuler, D. De Bernardis, A. M. auchli, and
P. Rabl, The Vacua of Dipolar Cavity Quantum
Electrodynamics, arXiv:2004.13738 (2020).
[54] M. Bamba, K. Inomata, and Y. Nakamura, Su-
perradiant Phase Transition in a Superconduct-
ing Circuit in Thermal Equilibrium, Phys. Rev.
Lett. 117, 173601 (2016).
[55] A. Settineri, O. Di Stefano, D. Zueco, S.
Hughes, S. Savasta, and F. Nori, Gauge
freedom, quantum measurements, and time-
dependent interactions in cavity and circuit
QED, arXiv:1912.08548 (2019).
[56] C. Cohen-Tannoudji, J. Dupont-Roc, and G.
Grynberg, Photons and Atoms (Wiley, New
York, 1997).
[57] A. Das, K. Sengupta, D. Sen, and B. K.
Chakrabarti, Infinite-range Ising ferromagnet
in a time-dependent transverse magnetic field:
Quench and ac dynamics near the quantum crit-
ical point, Phys. Rev. B 74, 144423 (2006).
[58] H. T. Quan and F. M. Cucchietti, Quantum fi-
delity and thermal phase transitions, Phys. Rev.
E 79, 031101 (2009).
Accepted in Quantum 2020-09-15, click title to verify. Published under CC-BY 4.0. 13
[59] J. Wilms, J. Vidal, F. Verstraete, and S. Dusuel,
Finite-temperature mutual information in a sim-
ple phase transition, J. Stat. Mech. P01023
(2012).
[60] E. Irish, Generalized Rotating-Wave Approxima-
tion for Arbitrarily Large Coupling, Phys. Rev.
Lett. 99, 173601 (2007).
[61] 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).
[62] A. Le Boite, Theoretical methods for ultrastrong
light-matter interactions, Adv. Quantum Tech-
nol. 3, 1900140 (2020).
[63] M. Aparicio Alcalde, M. Bucher, C. Emary, and
T. Brandes, Thermal phase transitions for Dicke-
type models in the ultrastrong-coupling limit,
Phys. Rev. E 86, 012101 (2012).
[64] A. Ridolfo, S. Savasta, and M. J. Hartmann,
Nonclassical Radiation from Thermal Cavities
in the Ultrastrong Coupling Regime, Phys. Rev.
Lett. 110, 163601 (2013).
[65] A. Ridolfo, M. Leib, S. Savasta, and M. J.
Hartmann, Thermal emission in the ultrastrong-
coupling regime, Phys. Scr. 2013, 014053 (2013).
[66] T. Chervy, A. Thomas, E. Akiki, R. M. A. Ver-
gauwe, A. Shalabney, J. George, E. Devaux, J. A.
Hutchison, C. Genet, and T. W. Ebbesen, Vibro-
Polaritonic IR Emission in the Strong Coupling
Regime, ACS Photonics 5, 217 (2018).
[67] F. Armata, G. Calajo, T. Jaako, M. S. Kim, and
P. Rabl, Harvesting Multiqubit Entanglement
from Ultrastrong Interactions in Circuit Quan-
tum Electrodynamics, Phys. Rev. Lett. 119,
183602 (2017).
[68] T. Holstein and H. Primakoff, Field Dependence
of the Intrinsic Domain Magnetization of a Fer-
romagnet, Phys. Rev. 58, 1098 (1940).
[69] A. M. Bratkovsky and A. P. Levanyuk, Con-
tinuous Theory of Ferroelectric States in Ul-
trathin Films with Real Electrodes, Journal of
Computational and Theoretical Nanoscience 6,
10.1166/jctn.2009.1058 (2008).
[70] T. Jaako, J. J. Garcia-Ripoll, and P. Rabl,
Ultrastrong-coupling circuit QED in the radio-
frequency regime, Phys. Rev. A 100, 043815
(2019).
[71] L. Fusco, M. Paternostro, and G. De Chiara,
Work extraction and energy storage in the Dicke
model, Phys. Rev. E 94, 052122 (2016).
[72] Y. Ma, S. Su, and C. Sun, Quantum thermo-
dynamic cycle with quantum phase transition,
Phys. Rev. E 96, 022143 (2017).
[73] N. Cottet, S. Jezouin, L. Bretheau, P.
Campagne-Ibarcq, Q. Ficheux, J. Anders, A.
Auffeves, R. Azouit, P. Rouchon, and B. Huard,
Observing a quantum Maxwell demon at work,
Proc. Natl. Acad. Sci. 114, 7561 (2017).
[74] M. Naghiloo, J. J. Alonso, A. Romito, E. Lutz,
and K. W. Murch, Information Gain and Loss for
a Quantum Maxwell’s Demon, Phys. Rev. Lett.
121, 030604 (2018).
[75] Y. Masuyama, K. Funo, Y. Murashita, A.
Noguchi, S. Kono, Y. Tabuchi, R. Yamazaki,
M. Ueda, and Y. Nakamura, Information-to-work
conversion by Maxwells demon in a supercon-
ducting circuit quantum electrodynamical sys-
tem, Nat. Commun. 9, 1291 (2018).
[76] M. A. Alcalde, E. Arias, Quantum Heat En-
gine and Quantum Phase Transition: through
Anisotropic LMG and Full Dicke models,
arXiv:1906.00292 (2019).
A Numerics
To perform an exact numerical evaluation of the par-
tition function and of other quantities derived from
it we make use of the fact that [H
cQED
,
~
S
2
] = 0.
This means that the Hamiltonian is block-diagonal
in the eigenbasis of the collective spin operator,
~
S =
(S
x
, S
y
, S
z
), and we can evaluate the partition func-
tions Z
s
for each subspace separately, taking the mul-
tiplicity of each angular-momentum sector into ac-
count [see Eqs. (13) and (14)]. This reduces the di-
mension of the spin part from 2
N
to N + 1 states.
For the actual calculation of Z
s
we first change into
the polaron frame, as explained in Sec. 3.3. The
Hamiltonian
˜
H
cQED
is then projected onto the col-
lective spin eigenstates |s, mi, where m = s, . . . , s,
and the Hilbertspace of the photon mode is trun-
cated at a maximal photon number N
ph
. For this
truncated Hamiltonian we compute the full set of
eigenvalues
i
s
to obtain Z
s
=
P
i
e
β
i
s
. The cutoff
number N
ph
is increased until the results converge.
Since the polaron transformation absorbes any mean
displacements of the photon field, the actual cutoff
numbers remain moderately large, even in the fer-
roelectric phase. For the parameter regimes consid-
ered in this paper we find that the heuristic choice
N
ph
= max{40, 20 × dk
B
T
max
/(~ω
c
)e} already pro-
duces very accurate results. If not mentioned other-
wise, this value of N
ph
is used for the plots presented
in this paper.
B Mean-field theory for the LMG
model
For J < 0 the LMG model constitutes a simple model
for ferroelectricity with all-to-all interactions. It is
Accepted in Quantum 2020-09-15, click title to verify. Published under CC-BY 4.0. 14
thus expected that for N 1 the phase transition
point is accurately predicted by the mean-field Hamil-
tonian
H
MF
LMG
= ~ω
0
S
z
+ ~
2J
N
Σ
x
S
x
~
J
N
Σ
2
x
, (51)
where Σ
x
= hS
x
i. Under this approximation, the re-
sulting free energy is given by
F
MF
LMG
x
) = k
B
T N ln
2 cosh
~
2k
B
T

~
J
N
Σ
2
x
,
(52)
where =
p
ω
2
0
+ 4J
2
Σ
2
x
/N
2
. In the paraelectric
phase the free energy has only a single minimum at
Σ
x
= 0, while the ferroelectric phase is character-
ized by the appearance of two degenerate minima at
a finite value of Σ
x
. The transition between the two
phases is given by Eq. (43).
Equation (52) allows us to derive some useful re-
lations for the spin expectation values of the LMG
model. From the condition F
MF
LMG
/∂Σ
x
= 0 we find
F
MF
LMG
Σ
x
=
2~J
N
J
tanh
~
2k
B
T
+ 1
Σ
x
= 0.
(53)
For J < J
c
a solution Σ
x
6= 0 exists. Subsequently,
an expression for Σ
z
= Tr{S
z
ρ
MF
LMG
} can be derived
via
Σ
z
=
F
MF
LMG
(~ω
0
)
=
Nω
0
2Ω
tanh
~
2k
B
T
. (54)
Eq. (53) and Eq. (54) are transcendental and there-
fore, in general, no explicit expressions for Σ
x
and Σ
z
can be found.
C Perturbation theory
Consider a generic system with Hamiltonian H =
H
0
+ H
g
and free energy F = F
0
+ F
g
, such that F
0
is
the free energy of the unperturbed Hamiltonian H
0
.
In this case we obtain the following general expression
for the coupling-induced part of the free energy,
e
βF
g
= hT e
R
β
0
H
g
(τ)
i
0
= e
β
P
n
F
(n)
g
. (55)
Here, h·i
0
denotes the average with respect to the
thermal state of H
0
, H
g
(τ) = e
τH
0
H
g
e
τH
0
and T is
the time-ordering operator along the imaginary time
axis. In systems where H
g
g is only a small per-
turbation to H
0
we can use a cumulant expansion to
approximate the expectation value in terms of a finite
number of F
(n)
g
g
n
.
C.1 Weak coupling limit
We first apply this perturbation theory to the cou-
pling Hamiltonian H
g
as defined in Eq. (21). Since in
this case the thermal average of the linear term van-
ishes, ha + a
i
0
= 0, the lowest-order correction to the
free energy is given by
F
(2)
g
=
~g
2
ω
c
hS
2
x
i
0
~
2
g
2
β
Z
β
0
1
Z
τ
1
0
2
C(τ
1
, τ
2
)hS
x
(τ
1
)S
x
(τ
2
)i
0
,
(56)
where
C(τ
1
, τ
2
) = h[a(τ
1
) + a
(τ
1
)][a(τ
2
) + a
(τ
2
)]i
0
= (N
th
+ 1)e
~ω
c
(τ
1
τ
2
)
+ N
th
e
~ω
c
(τ
1
τ
2
)
.
(57)
For interacting dipoles the spin expectation value
must be evaluated numerically. In terms of the en-
ergies E
n
and eigenstates |E
n
i of H
dip
we obtain
hS
x
(τ
1
)S
x
(τ
2
)i
0
=
1
Z
0
dip
X
n,m
e
βE
n
+(τ
1
τ
2
)(E
n
E
m
)
× |hE
n
|S
x
|E
m
i|
2
.
(58)
For non-interacting dipoles this calculation can be
carried out analytically using
hS
x
(τ
1
)S
x
(τ
2
)i
0
=
1
4
e
(τ
1
τ
2
)~ω
0
hS
+
S
i
0
+
1
4
e
(τ
1
τ
2
)~ω
0
hS
S
+
i
0
(59)
and hS
±
S
i
0
= N [1 tanh (~ω
0
/2k
B
T )] /2. Alto-
gether and writing F
(2)
g
= N~g
2
/(4ω
c
)f
g
we obtain
Nf
g
= N (N
th
+ 1)hS
+
S
i
0
I(ω
0
ω
c
)
(N
th
+ 1)hS
S
+
i
0
I(ω
0
ω
c
)
N
th
hS
+
S
i
0
I(ω
0
+ ω
c
)
N
th
hS
S
+
i
0
I(ω
c
ω
0
),
(60)
where
I(∆) =
~ω
c
β
Z
β
0
1
Z
τ
1
0
2
e
~∆(τ
1
τ
2
)
. (61)
After some further simplifications the expression for
f
g
reduces to the result given in Eq. (23).
C.2 Bound on the free energy correction
From the general expression for the correlation
function given in Eq. (58) one can show that
hS
x
(τ
1
)S
x
(τ
2
)i
0
hS
2
x
i
0
for τ
1
τ
2
. This inequal-
ity can be used together with (N
th
+ 1)I(ω
c
) +
N
th
I(ω
c
) = 1 in the second line of Eq. (56) in or-
der to derive the upper and lower bounds 0 f
g
4hS
2
x
i
0
/N. To improve the upper bound we can re-
peat the whole perturbation calculation in a displaced
frame,
˜
H
cQED
= D
(α)H
cQED
D(α) = H
0
+
˜
H
g
(α),
Accepted in Quantum 2020-09-15, click title to verify. Published under CC-BY 4.0. 15
where D(α) = e
αa
α
a
is the displacement operator.
Specifically, by choosing α = ghS
x
i
0
c
we obtain
˜
H
g
α =
ghS
x
i
0
ω
c
= g(a + a
)∆S
x
+
g
2
ω
c
(∆S
x
)
2
,
(62)
where S
x
= S
x
hS
x
i
0
. This means that we can
repeat the analysis from above with S
x
being replaced
by S
x
. This leads to the stricter bound for f
g
given
in Eq. (26).
C.3 Low-frequency limit
We can use the same perturbation scheme to evaluate
the lowest-order corrections to the free energy when
ω
0
0. To do so we first change to the polaron
frame as described in Sec. 3.3 and decompose the total
Hamiltonian as
˜
H
cQED
= H
0
+ H
ω
0
. Here
H
0
= ~ω
c
a
a +
~
4
X
i,j
J
ij
σ
i
x
σ
j
x
,
H
ω
0
= ~ω
0
cos(
ˆ
θ)S
z
sin(
ˆ
θ)S
y
,
(63)
where
ˆ
θ = i(g
c
)(a
a). According to this parti-
tioning, the bare Hamiltonian H
0
is diagonal in the
σ
x
basis and hH
ω
0
i
0
= 0. Thus, up to second order in
ω
0
the free energy is given by F ' F
0
+ F
(2)
ω
0
, where
F
(2)
ω
0
=
~
2
ω
2
0
β
Z
β
0
1
Z
τ
1
0
2
× [C
c
(τ
1
, τ
2
)hS
z
(τ
1
)S
z
(τ
2
)i
0
+ C
s
(τ
1
, τ
2
)hS
y
(τ
1
)S
y
(τ
2
)i
0
],
(64)
where C
c
(τ
1
, τ
2
) and C
s
(τ
1
, τ
2
) are the correlation
functions of the operators cos(
ˆ
θ) and sin(
ˆ
θ), respec-
tively. The correlation function for the photons can be
calculated exactly using the properties of the displace-
ment operator. Specifically, given a pair of complex
number z
1
and z
2
the following formula holds
he
z
1
a
z
2
a
i
0
= e
z
1
z
2
(1+2N
th
)/2
. (65)
The two-point correlation functions can then be ex-
pressed in terms of an infinite series,
C
c
(τ
1
, τ
2
) = e
g
2
ω
2
c
(1+2N
th
)
X
r,q=0
K
rq
e
(qr)ω
c
(τ
1
τ
2
)
(66)
and
C
s
(τ
1
, τ
2
) = e
g
2
ω
2
c
(1+2N
th
)
X
r,q=0
Q
rq
e
(qr)ω
c
(τ
1
τ
2
)
,
(67)
with coefficients
K
rq
=
[1 + (1)
r+q
]
2
g
ω
c
2(r+q)
(1 + N
th
)
r
N
q
th
r! q!
,
Q
rq
=
[1 (1)
r+q
]
2
g
ω
c
2(r+q)
(1 + N
th
)
r
N
q
th
r! q!
.
(68)
Altogether we obtain
F
(2)
ω
0
=
~
2
˜ω
2
0
β
X
r,q=0
Z
β
0
1
Z
τ
1
0
2
e
(qr)ω
c
(τ
1
τ
2
)
×
K
rq
hS
z
(τ
1
)S
z
(τ
2
)i
0
+ Q
rq
hS
y
(τ
1
)S
y
(τ
2
)i
0
,
(69)
where ˜ω
0
= ω
0
exp[g
2
/(2ω
2
c
)(1 + 2N
th
)].
In Fig. 3(a) this result is used to evaluate χ
z
for
non-interacting dipoles and, as expected, we find that
it is in perfect agreement with the full numerical sim-
ulations for arbitrary g and arbitrary temperatures.
Note that in the series above one can interpret each
term as a process in which q-photons are absorbed
and r-photons are emitted. In this way it is clear
that in the ultrastrong coupling regime there are pro-
cesses with multiphoton scattering, emission and ab-
sorption. These processes become more and more im-
portant as the coupling strength increases.
D Variational LMG Hamiltonian
In the variational ansatz described in Sec. 3.3, the
trial Hamiltonian H
given in Eq. (32) is simply the
decoupled cavity QED system with ω
0
being replaced
by a renormalized frequency ˜ω
0
. Therefore, we obtain
F
= F
0
c
+ F
0
dip
(˜ω
0
) and
hH H
i
ρ
= ~ω
0
hcos(
ˆ
θ)S
z
i
0
~˜ω
0
hS
z
i
0
, (70)
where we have already assumed hS
y
i
0
= 0. The ther-
mal expectation value of the cosine operator can be
evaluated with the help of Eq. (65).
Altogether, the variational free energy F
V
can be
written as
F
V
= F
0
c
+F
0
dip
(˜ω
0
)+~
ω
0
e
g
2
ω
2
c
(N
th
+1/2)
˜ω
0
hS
z
i
0
.
(71)
Taking the derivative of F
V
with respect to ˜ω
0
yields
the extremal condition
F
V
˜ω
0
=
F
0
dip
(˜ω
0
)
˜ω
0
~hS
z
i
0
(72)
+ ~
ω
0
e
g
2
ω
2
c
(N
th
+1/2)
˜ω
0
hS
z
i
0
˜ω
0
!
= 0.
Using the general relation F
0
dip
(˜ω
0
)/∂(~˜ω
0
) = hS
z
i
0
we obtain
˜ω
0
= ω
0
e
g
2
2ω
2
c
(1+2N
th
)
. (73)
Note that this result is independent of the dipole-
dipole couplings J
ij
.
E Emission spectrum
For C
t
C, the capacitive coupling between the
cavity and the transmission line shown in Fig. 1 can
Accepted in Quantum 2020-09-15, click title to verify. Published under CC-BY 4.0. 16
be modeled by a Hamiltonian of the form H
ct
=
C
t
V V
t
. Here C
t
is the coupling capacitance and
V
t
=
P
k
V
k
(b
k
+ b
k
) is the voltage operator of the
transmission line, which we expressed in terms of
a set of free EM modes with bosonic operators b
k
and frequency ω
k
. We write V = V
0
X and define
λ
k
= C
t
V
0
V
k
/~ such that the total Hamiltonian of the
cavity QED system and the transmission line reads
H = H
cQED
+
X
k
~ω
k
b
k
b
k
+~
X
k
λ
k
(b
k
+b
k
)X. (74)
For a conventional transmission line we have λ
k
ω
k
.
We can formally integrate the equations of motion
for the mode operators b
k
,
b
k
(t) = b
k
(0)e
k
t
k
Z
t
0
dt
0
e
k
(tt
0
)
X(t
0
).
(75)
Therefore, by assuming that initially the transmission
line is in the vacuum state, we find that
hb
k
b
k
i(t) = λ
2
k
Z
t
0
Z
t
0
dt
0
dt
00
e
k
(t
00
t
0
)
hX(t
0
)X(t
00
)i.
(76)
To proceed we write X(t) = X
(t) + X
+
(t), such
that X
+
(t) =
P
nm
X
nm
e
nm
t
contains only con-
tributions that oscillate with a positive frequency,
ω
nm
= (E
n
E
m
)/~ 0, and X
(t) = X
+
(t) [64].
Then, for long enough times t γ
1
, where γ is the
characteristic thermalization rate of the cavity QED
system, we obtain
hb
k
b
k
i(t) ' λ
2
k
t × 2Re
Z
0
hX
+
(τ)X
(0)i
0
e
k
τ
.
(77)
The total power radiated into the transmission line
is
P
rad
=
X
k
~ω
k
t
hb
k
b
k
i(t) (78)
and by writing P
rad
=
R
0
S
bb
(ω) we obtain the
general expression for the black-body spectrum
S
bb
(ω) = ~ωJ(ω)C
X
(ω). (79)
Here J(ω) =
P
k
λ
2
k
δ(ω ω
k
) is the spectral density
of the transmission line and
C
X
(ω) = 2Re
Z
0
hX
+
(τ)X
(0)i
0
e
τ
. (80)
By assuming an Ohmic spectral density we can write
J(ω) = κω/(2πω
c
), where κ is the rate at which the
bare cavity decays into the transmission line. For a
completely isolated system the correlation function of
the system operator X is given by
C
X
(ω) = 2π
X
n>m
e
βE
n
Z
|hn|X|mi|
2
δ(ω ω
nm
), (81)
and we obtain the total emitted power given in
Eq. (45).
For the evaluation of the black-body spectrum
S
bb
(ω) one must keep in mind that the cavity QED
system is in constant interaction with the surround-
ing thermal bath, which induces a finite broaden-
ing of all the spectral lines. A detailed investiga-
tion of such line-broadening effects is beyond the
scope of the current analysis. Instead, we simply re-
place all resonances by a Lorentzian profile with a
phenomenological width γ, while still approximating
ωJ(ω) ' ω
nm
J(ω
nm
) for each transition. As a result,
we obtain the spectrum given in Eq. (44).
Accepted in Quantum 2020-09-15, click title to verify. Published under CC-BY 4.0. 17