Fermat's Library | q-2020-09-28-335 annotated/explained version.

Thermodynamics of ultrastrongly coupled light-matter

systems

Philipp Pilar

1

, Daniele De Bernardis

1

, and Peter Rabl

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 modiﬁcations of

derivative quantities such as the speciﬁc 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 signiﬁcant

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 eﬀects, for example, in

the ﬁeld 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-coeﬃcient and postulated the eﬀect 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 ﬁeld of quantum optics, and so

on. In most of these and related examples the EM

ﬁeld 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

modiﬁcations 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 artiﬁcial 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 simpliﬁed descriptions,

for example, based on the Dicke [25, 26] or the

Hopﬁeld [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 ﬁelds 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

ﬁeld 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 simpliﬁcations, 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 ﬁelds. 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, ﬁrst 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 ﬁxed value of G, the changes in the free energy

per particle, F

g

/N, vanish. Both ﬁndings 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 ﬁeld mode can induce substantial

modiﬁcations 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

ﬁeld only couples to a single collective dipole mode

while the other N − 1 orthogonal excitation modes

remain unaﬀected. 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 ﬁnd that the coupling to

the cavity can indeed inﬂuence quantities such as the

electric susceptibility or the speciﬁc 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 diﬀerent 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 ﬁeld 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, ﬁnally, col-

lapses again to a single resonance. At the same time

we ﬁnd 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 eﬀects in real

and artiﬁcial light-matter systems.

The remainder of the paper is structured as follows:

In Sec. 2 and Sec. 3 we ﬁrst introduce the EDM and

discuss some general properties of the free energy of

a cavity QED system in diﬀerent coupling regimes.

In Sec. 4 we then analyze in more detail the cavity-

induced modiﬁcations 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 diﬀerent coupling regimes and we conclude

our work in Sec. 6. The appendices A-D contain addi-

tional details about diﬀerent 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 eﬀects, we can restrict our discus-

sion to cavity and circuit QED setups in the GHz

to THz regime, where these eﬀects are experimen-

tally most relevant. In this case the conﬁned electro-

magnetic ﬁeld 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

conﬁguration also ensures that all higher EM exci-

tations are well separated in frequency and that the

electric ﬁeld is approximately homogeneous across the

ensemble of dipoles. The dipoles themselves are mod-

eled as eﬀective 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 ﬁeld.

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 eﬀect

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 suﬃcient

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

2L

=

(Q + P/d)

2

2C

+

Φ

2

2L

, (4)

where V is the voltage diﬀerence across the capaci-

tor plates and Φ the magnetic ﬂux. 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 suﬃciently homogeneous ﬁeld, 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 ﬁxed orientation

of the dipoles. This energy might be modiﬁed in the

presence of metallic cavity mirrors [33, 42, 53], but it

is independent of the frequency or the vacuum ﬁeld

amplitude of the dynamical mode. The coupling to

the dynamical ﬁeld is then proportional to g and in-

cludes the collective dipole-ﬁeld 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 speciﬁc 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 ﬁeld, 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 ﬂux Φ 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 deﬁnition 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

ﬁeld, the displacement ﬁeld and the magnetic ﬁeld,

respectively. For a detailed derivation and justiﬁca-

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 inﬂuence 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-ﬁeld 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-ﬁeld

decoupling of the dipole-ﬁeld 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-ﬁeld 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 ﬁrst 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 diﬀerent

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-ﬁeld 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-ﬁeld

approximation is given by

Z

MF

(Σ

x

) = Z

0

c

×

¯

Z

MF

dip

(Σ

x

). (17)

Here, the ﬁrst 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

eﬀective 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-ﬁeld

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 satisﬁed.

Equation (18) shows that cavity-induced correc-

tions to the thermodynamic properties of a dipolar

system are only aﬀected by ﬂuctuations, but not by

the mean orientation of the dipoles. Therefore, by ap-

plying a second mean-ﬁeld decoupling for the dipoles

(see Appendix B), the eﬀect 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-ﬁeld 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 ﬂuctuations 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 ﬁnd 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 modiﬁcation of the dipolar system already

in the collective USC regime, g

√

N ∼

√

ω

0

ω

c

, and

could, at ﬁrst 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-ﬁeld

theory systematically overestimates the inﬂuence of

the cavity mode, in particular at higher temperatures.

This somewhat counterintuitive trend can be traced

back to the fact that the mean-ﬁeld 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-ﬁeld 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-ﬁeld 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 inﬂuence 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-ﬁeld 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 ﬁnd 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 ﬁnd 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

i−hσ

i

x

ihσ

j

x

i, is ﬁnite

and the collective ﬂuctuations scale as (∆S

x

)

2

∼ N.

Therefore, under very generic conditions, by taking

the limit N → ∞ with G kept ﬁxed, we obtain

lim

N→∞

F

(2)

g

N

= 0. (27)

This result conﬁrms 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 ﬂuctuations, 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 conﬁrm 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 beneﬁt 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 ﬁnite ω

0

the eﬀects 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 eﬀective

spin Hamiltonian [32]

H

eﬀ

= ~ω

0



e

−

g

2

2ω

2

c

− 1



S

z

+

~ω

2

0

ω

c

2g

2



S

2

x

−

~

S

2



.

(30)

This eﬀective 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 brieﬂy

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-ﬁeld susceptibility χ

z

(solid

lines) for diﬀerent 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-ﬁeld 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 modiﬁcations 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 speciﬁc 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 inﬂuence such quantities, we discuss

in this section three elementary examples.

4.1 USC modiﬁcations of the Curie law

As a ﬁrst 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

.

Speciﬁcally, we consider the limit ω

0

→ 0, where we

obtain the zero-ﬁeld 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 ﬁeld the susceptibility follows

the usual Curie law

χ

z

|

g→0

(T ) =

α

C

T

, (36)

with a Curie constant α

C

= ~/(4k

B

). In the context

of cavity QED, the behavior of this quantity for ﬁnite

g is of particular importance. Since the dipoles decou-

ple from the cavity mode for ω

0

= 0, the zero-ﬁeld 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 diﬀerent 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 modiﬁcations be-

come more signiﬁcant. 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 ﬁrst 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 eﬀect 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 eﬀective 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

deﬁned

in Eq. (41), which determines the corrections to the speciﬁc

heat, C

g

/N, in the collective USC regime. (b) Dependence of

C

g

on the coupling strength g for a ﬁxed temperature and an

increasing number of dipoles, N = 10, 20, 30, 40, 50. For this

plot ω

0

= 0.5ω

c

and a photon cutoﬀ 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 eﬀective 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 modiﬁcations of the speciﬁc 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 speciﬁc 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 diﬀerent frequency ratios, ω

0

/ω

c

.

Therefore, taking the limit N → ∞ for a ﬁxed G, the

corrections to the speciﬁc 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 modiﬁca-

tions. On a qualitative level, these corrections can be

understood from a cavity-induced suppression of ω

0

,

but overall we ﬁnd 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

speciﬁc 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 ﬁnding demonstrates that in the

non-perturbative regime the coupling to a single dy-

namical ﬁeld degree of freedom can have a substantial

inﬂuence on extensive thermodynamic quantities of a

large ensemble of dipoles.

4.3 USC modiﬁcations of the ferroelectric

phase transition

A central topic of interest in the ﬁeld 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-ﬁeld decoupling of the collective interaction

term, S

2

x

→ 2Σ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-ﬁeld

theory, see Eq. (43). The same boundaries are shown in (b)

for diﬀerent coupling strengths g, where for the mean-ﬁeld

results ω

0

has been replaced by ˜ω

0

. (c) Dependence of the

critical temperature T

c

on the coupling parameter g/ω

c

for

a ﬁxed 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 ﬁnite

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 deﬁne

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-ﬁeld theory, even for a

small number of N = 20 dipoles.

When the dipoles are coupled to the cavity mode,

a ﬁnite 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 ﬁnite charge (or displacement ﬁeld) hQi ∼ ha + a

†

i.

The relevant cavity observables, hV i ∼ hA + A

†

i and

hΦi ∼ ihA

†

− Ai, are not aﬀected 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 inﬂuenced 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 eﬀects beyond the ground state and show

in Fig. 5(b) the coupling-induced modiﬁcation of the

phase boundary in the whole T − J plane for diﬀer-

ent values of g/ω

c

. In this plot, the exact analytic

results are compared with a modiﬁed mean-ﬁeld 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 ﬁnd 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 ﬁx 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 ﬂuctuations. This

behavior is qualitatively reproduced by the modiﬁed

mean-ﬁeld ansatz.

5 Black-body radiation

The emission spectrum of a hot body was one of the

ﬁrst 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 ﬁeld 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-Primakoﬀ 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 oﬀset ~ω

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

ﬂuctuations 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 suﬃciently 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 ﬁnite 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 diﬀerent 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 ﬁeld. 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 ﬁnd 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 ﬁeld 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 ﬁnd

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-Primakoﬀ 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



+

N−1

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 ﬂux 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-Primakoﬀ 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-Primakoﬀ 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 eﬀects. 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 inﬂuence 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 ﬁrst line of this equation have a much

weaker dependence on G. Further, we ﬁnd 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 diﬀerent 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 conﬁrm the basic intuition that the coupling to

a single cavity mode cannot signiﬁcantly change the

properties of a larger ensemble of dipoles. While a

large collective coupling strength G has a substantial

inﬂuence 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 modiﬁcations 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 eﬀects further experimen-

tal and theoretical investigations are still required, in

particular also on the inﬂuence of the metallic bound-

aries on electrostatic interactions [33, 42, 53, 69] and

other modiﬁcations 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-

ciﬁc 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 eﬀect 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-

ﬁcial 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 ﬁeld, 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. Greﬀet, 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 artiﬁcial 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 eﬀects in mesoscopic systems, Physics Re-

ports 408, 315 (2005).

[27] J. J. Hopﬁeld, 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 inﬂu-

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 ﬁeld: 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 Modiﬁcations 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. L¨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, Inﬁnite-range Ising ferromagnet

in a time-dependent transverse magnetic ﬁeld:

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 ﬁ-

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 ﬁnite-

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. Primakoﬀ, 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.

Auﬀeves, 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 ﬁrst 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 cutoﬀ

number N

ph

is increased until the results converge.

Since the polaron transformation absorbes any mean

displacements of the photon ﬁeld, the actual cutoﬀ

numbers remain moderately large, even in the fer-

roelectric phase. For the parameter regimes consid-

ered in this paper we ﬁnd 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-ﬁeld 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-ﬁeld 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 ﬁnite 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 ﬁnd

∂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

dτ 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 ﬁnite

number of F

(n)

g

∼ g

n

.

C.1 Weak coupling limit

We ﬁrst apply this perturbation theory to the cou-

pling Hamiltonian H

g

as deﬁned 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

dτ

1

Z

τ

1

0

dτ

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

dτ

1

Z

τ

1

0

dτ

2

e

~∆(τ

1

−τ

2

)

. (61)

After some further simpliﬁcations 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.

Speciﬁcally, 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 ﬁrst 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

dτ

1

Z

τ

1

0

dτ

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. Speciﬁcally, 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 inﬁnite series,

C

c

(τ

1

, τ

2

) = e

−

g

2

ω

2

c

(1+2N

th

)

∞

X

r,q=0

K

rq

e

(q−r)ω

c

(τ

1

−τ

2

)

(66)

and

C

s

(τ

1

, τ

2

) = e

−

g

2

ω

2

c

(1+2N

th

)

∞

X

r,q=0

Q

rq

e

(q−r)ω

c

(τ

1

−τ

2

)

,

(67)

with coeﬃcients

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

dτ

1

Z

τ

1

0

dτ

2

e

(q−r)ω

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 ﬁnd 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

c−t

=

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 deﬁne

λ

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

−iω

k

t

− iλ

k

Z

t

0

dt

0

e

−iω

k

(t−t

0

)

X(t

0

).

(75)

Therefore, by assuming that initially the transmission

line is in the vacuum state, we ﬁnd that

hb

†

k

b

k

i(t) = λ

2

k

Z

t

0

Z

t

0

dt

0

dt

00

e

iω

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

n≥m

X

nm

e

iω

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

dτ hX

+

(τ)X

−

(0)i

0

e

−iω

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

dω 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

dτ hX

+

(τ)X

−

(0)i

0

e

−iωτ

. (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 ﬁnite broaden-

ing of all the spectral lines. A detailed investiga-

tion of such line-broadening eﬀects is beyond the

scope of the current analysis. Instead, we simply re-

place all resonances by a Lorentzian proﬁle 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

Comments

Products

Project