Thermodynamics of ultrastrongly coupled light-matter systems

We study the thermodynamic properties of a system of two-level dipoles that are coupled ultrastrongly to a single cavity mode. By using exact numerical and approximate analytical 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 analysis we identify the lowest-order cavity-induced corrections to those quantities in the collective ultrastrong coupling regime and show that for even stronger interactions the presence of a single cavity mode can strongly modify extensive 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 temperature and a characteristic broadening and collapse of the black-body spectrum of the cavity mode. Apart from a purely fundamental interest, these general insights will be important for identifying potential applications of ultrastrong-coupling effects, for example, in the field of quantum chemistry or for realizing quantum thermal machines.


I. INTRODUCTION
Undoubtedly, the interplay between statistical physics and the theory of electromagnetic (EM) radiation played a very important role in the history of modern 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 microscopically. 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 surrounding 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 molecular 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 potentials [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 systems 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 possible, one usually resorts to simplified descriptions, for example, based on the Dicke [25,26] or the Hopfield [27] model. However, such reduced models often do not represent the complete energy of the system [28][29][30][31][32][33][34][35] or contain gauge artefacts [33,[36][37][38][39] that prevent their applicability in the USC regime. More generally, while in weakly coupled cavity QED systems the role of static dipole-dipole inter-actions can often be neglected or modelled independently of the dynamical EM mode, this is no longer the case in the USC regime [33,[40][41][42]. An inconsistent treatment of static and dynamical fields can thus very easily lead to wrong predictions or a misinterpretation of results. A prominent example in this respect is the superradiant phase transition of the Dicke model [43][44][45], which is often described as cavity-induced, but which can be understood as a regular ferroelectric instability in a system of strongly attractive 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 systems so far.
In this paper we study the thermodynamics of cavity 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 remains consistent with basic electrodynamics at arbitrary interaction strengths and distinguishes explicitly between static and dynamical electric fields. It thus allows us to evaluate the FIG. 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. free energy of the most 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 . Therefore, when taking the limit N → ∞ for a fixed value of G, the changes in the free energy per particle, F g /N , vanish. This result contradicts common claims that a large collective coupling to a quantized field mode can induce substantial modifications of material properties. This prediction is, however, consistent with similar studies about molecular properties in the ground state of cavity QED systems [46,47], and can be intuitively explained by a simple polariton 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 material. 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 we increase g, the spectrum evolves from the usual polariton doublet into a broad and disordered set of lines and, finally, collapses 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 suppress the total radiated power. Therefore, the analysis 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. II and Sec. III 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. IV we then analyze in more detail the cavity-induced modifications for the cases of paraelectric and ferroelectric ensembles of dipoles. Finally, in Sec. V we evaluate the black-body spectrum of the cavity mode in different coupling regimes and we conclude our work in Sec. VI. The appendices A-D contain additional details about different approximation methods for the free energy and the derivation of the emission spectrum.

II. MODEL
We consider a generic cavity QED setup, where N two-level dipoles are coupled to a single electromagnetic mode. However, since we are interested in both thermal und USC effects, we can restrict our discussion to cavity and circuit QED setups in the GHz to THz regime, where these effects are experimentally most relevant. In this case the confined electromagnetic field can be represented by the fundamental mode of a lumped-element LC resonator [33,48] with capacitance C and inductance L (see Fig. 1). This configuration also ensures that all higher EM excitations are well separated in frequency and that the electric field is approximately homogeneous across the ensemble of dipoles. The dipoles themselves are modeled as effective two-level systems with states |0 and |1 . The two states are separated by an energy ω 0 and they are coupled via an electric transition dipole moment µ to the electric field.

A. Hamiltonian
The Hamiltonian of the whole cavity QED system can be written as 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 where the σ i k are the usual Pauli operators for the ith dipole. The couplings J ij account for the effect of static dipole-dipole interactions as well as possible other types of non-electromagnetic couplings between the twolevel 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 [49] where the S k = 1/2 i σ i k are collective spin operators. 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.
In the lumped-element limit, the energy of the EM mode is given by where V is the voltage difference across the capacitor 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 quantized theory. For a sufficiently homogeneous field, the charge is given by Q = CV − P/d, where 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 where Q 0 = /(2Z), Φ 0 = Z/2 and Z = L/C is the cavity impedance. Altogether, we obtain the canonical cavity QED Hamiltonian [33] where ω c = 1/ √ LC and g = µQ 0 /( Cd) is the coupling strength. Note that within this model the regular Dicke model is recovered as the special case of a ferroelectric dipolar system with J ij = −g 2 /ω c . A potential realization of this scenario in circuit QED systems is discussed in Ref. [50].

B. 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 observables. For the dipoles, quantities of interest are the population imbalance, S z , or the polarization along the cavity field, S x ∼ P , etc. Since the operator Q for the total charge includes the induced charges from the dipoles, its value is typically not directly measurable. Therefore, for the cavity mode the relevant observables are the magnetic flux Φ and the voltage drop V and it is convenient to introduce the displaced photon annihilation operator With the help of this definition we obtain [33,51] where V 0 = Q 0 /C, and the total Hamiltonian can be written in a compact form as By comparing with Eq. (1), we see that the expectation value of A † A can be interpreted as the energy of the dynamical cavity mode in units of ω c . This is in contrast to the conventional photon number a † a , which depends on the chosen gauge [33,51] and has no simple interpretation in a strongly coupled cavity QED system. Note, however, that despite this more intuitive physical interpretation, A and A † do not commute with all the dipole operators and we must still use a and a † to represent the independent cavity degree of freedom. While we focus here on a lumped-element realization 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 systems using the replacements [33,41,48,51,52] V → E, 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 justification of Hamiltonian (6) in dipolar cavity QED and circuit QED settings, see Refs. [32,33,36].

III. THE FREE ENERGY IN CAVITY QED
By assuming that the cavity QED system is weakly coupled to a large reservoir of temperature T , the resulting equilibrium state of the system is where β = 1/(k B T ) and Z = Tr{e −βHcQED } is the partition function. In this case the central quantity of interest is the free energy, which we divide into the free energies F 0 c and F 0 dip of the decoupled subsystems and a remaining contribution 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 properties of the bare subsystems.
For small and moderately large ensembles the partition 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 Here Z s is the partition function of H cQED constrained to a total spin s and accounts for the multiplicity of the respective multiplet due to permutation symmetry. For small and moderate temperatures, we use this approach to evaluate the exact free energy for systems with N ≈ 1 − 100 dipoles.

A. Mean-field theory
In the analysis of the regular Dicke model [43][44][45] with N 1, a frequently applied approximation for evaluating the free energy is based on the mean-field decoupling of the dipole-field interaction, where the expectation values α = a and Σ x = S x must be determined self-consistently. Under this meanfield approximation Hamiltonian (6) can be written as the sum of two independent parts, where x /ω c . The first two terms describe the energy 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 and (17), the first factor is the partition function of the bare cavity and Z MF dip is the partition function of an ensemble of dipoles with effective Hamiltonian (which includes the constant energy shift from the displaced oscillator) By applying a second mean-field decoupling for the dipoles (see Appendix A), which is expected to give accurate results for the non-interacting or ferroelectric LMG model [53][54][55], the effect of the cavity vanishes completely and F MF To be more precise we can evaluate the partition function of the spin system,Z MF dip , exactly. In this case we see from Eq. (18) that in the paraelectric phase, where Σ x = 0, the cavity induces a huge renormalization of the interaction term, J → J + g 2 N/ω c . This renormalization becomes a substantial modification of the dipolar system in the collective USC regime, g √ N ∼ √ ω 0 ω c , and would even prevent the ferroelectric instability for J < 0. However, as can be seen from a comparison with the exact free energy in Fig. 2 and other examples discussed below, such a strong modification is not observed in this regime. This discrepancy can be traced back to the fact that the mean-field decoupling neglects contributions which are second order in g(a + a † )S x and also scale approximately as Therefore, depending on the precise parameters, the mean-field decoupling either over-or underestimates the influence of the cavity mode and the approximation becomes uncontrolled. We conclude that a mean-field theory, as frequently employed to study ground states and thermal phases of the Dicke-or the LMG model, is not appropriate to describe the coupling of dipoles to a dynamical field mode in realistic cavity QED systems.

B. Collective USC regime
Many cavity QED experiments are operated in the regime G = g √ N ω c and N 1, where the collective 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, as a small perturbation and expand the free energy in powers of g. As a result of this calculation, which is detailed in Appendix B, we obtain the lowest-order correction to the bare free energy, which can be written in the form 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 resulting expression for f g still involves non-trivial correlation functions of spin operators, which for interacting dipoles must be evaluated numerically. For non-interacting dipoles this calculation can be carried out analytically and we obtain the explicit .
(22) In Fig. 2, this prediction from perturbation theory is compared to the exact free energy and we find that the cavityinduced 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. (22) provides us directly with the lowest-order correction to the ground state energy of a cavity QED system [1], which agrees with Hamiltonian perturbation theory. In the opposite high-temperature limit we obtain Therefore, the cavity-induced corrections to the free energy vanish quadratically with increasing temperature. 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 interacting dipoles (see Appendix B) For conventional dipolar systems with finite correlation length (in particular away from any critical points) one typically finds (∆S x ) 2 N/4. Therefore, by taking the limit N → ∞ with G kept fixed, we obtain This result confirms our basic intuition that the coupling of many dipoles to a single mode should not affect extensive thermodynamic properties of the matter system. While systems with stronger fluctuations, i.e.
(∆S x ) 2 ∼ N 2 , could in principle violate this conclusion, such systems would already exhibit very atypical thermodynamic properties in the absence of the cavity mode.

C. Non-perturbative regime
The physics of cavity QED changes drastically once g ∼ ω c and the light-matter interactions become nonperturbative at the level of individual dipoles. To analyze this regime, it is usually more convenient [56][57][58] to change to a polaron frame,H cQED = U H cQED U † , via the unitary transformation U = e g ωc Sx(a † −a) . In this frame the cavity QED Hamiltonian can be written as [32,33,38,59] where the interaction part now takes the form An immediate benefit of the polaron representation is that the interaction term is proportional to ω 0 . This shows that for ω 0 → 0 the coupling to the dynamical mode vanishes and we recover the electrostatic limit, 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] This effective model captures two important signatures 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-toall anti-ferroelectric coupling ∼ ω 2 0 . In principle, we can again apply perturbation theory to evaluate F g up to second order in H int and extend the results from Sec. III B 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 B. 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, Here ρ * is the thermal state and F * the corresponding free energy for the trial Hamiltonian H * . Based on the discussion above we choose which describes a non-interacting cavity QED system, but with a variable frequencyω 0 . By minimizing F V with respect toω 0 for each g we obtain (see Appendix C) where N th = 1/(e β ωc − 1). While from the comparison in Fig. 2 we see that overall 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.

IV. PARA-AND FERROELECTRICITY IN THE USC REGIME
While the free energy contains all the relevant information about the cavity QED system, we are usually interested in derivative quantities such the susceptibility, the specific heat, etc., or the existence of different 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.
A. USC modifications of the Curie law As a first example it is instructive to consider the simplest case of non-interacting dipoles, where 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 S z on the level splitting ω 0 . Specifically, we consider the limit ω 0 → 0, where we obtain the zero-field susceptibility from the second derivative of the free energy. In the limit of a vanishing bias field the susceptibility follows the usual Curie law 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 decouple from the cavity mode for ω 0 = 0, the zero-field susceptibility 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 become 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 behavior we approximate the susceptibility by the two dominant contributions, where p n and E n are the thermal occupation probabilities and the energy of the n-th eigenstate, respectively. The first term emerges from the change of S z due to small changes of the thermal populations 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 variational free energy F V discussed in Sec. III C, from which we obtain 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. (36) 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 systems. However, as evident from the effective spin model in Eq. (29), for couplings g > ω c the energy levels show indeed a quadratic scaling, E n ∼ ω 2 0 . From this effective model and by assuming that all spin levels are equally populated, p n ≈ 1/2 N , we obtain the approximate result 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. (29) breaks down. Beyond this point, which for large g is indicated by a small pump, 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 B). Although the resulting 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.

B. USC modifications of the specific heat
A second quantity of general interest in statistical physics is the heat capacity, For a decoupled system, the heat capacities of the cavity 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 negligible contribution to the specific heat C/N of the combined system. This can be shown explicitly in the collective USC regime, where for non-interacting dipoles 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 plotted 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 system 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. (40). In the non-perturbative regime we observe substantial modifications. On a qualitative level, these corrections can be understood from a cavityinduced suppression of ω 0 , but overall we find that the dependence of C g is not very accurately reproduced by the variational ansatz in Eq. (31) 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 and shifts, but does not vanish with increasing N . Combined with the behavior of the susceptibility discussed above, this finding demonstrates that in the nonperturbative regime the coupling to a single dynamical field degree of freedom can have a substantial influence on extensive thermodynamic quantities of a large ensemble of dipoles.

C. 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 equilibrium states of the standard Dicke model [43][44][45]. While in more accurate models for light-matter interactions 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 → 2ΣS x − Σ 2 (see Appendix A), from which one can derive the general relation between the critical coupling strength J c and the critical temperature T c [53][54][55], For T → 0 there is a critical coupling strength J c = −ω 0 , beyond which the dipoles enter a ferroelectric phase with S x = 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 classical 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 S x = 0 cannot be used to characterize the ferroelectric phase in exact numerical calculations. In Fig. 5(a) we show instead the quantitym = S 2 x , which provides a good indicator for the ferroelectric phase of the LMG model [55]. However, for the rather small numbers of dipoles assumed in the simulations of the full cavity QED model below, the variation ofm around the phase transition line is still rather smooth. Therefore, for the following analysis we consider instead the probability distribution p(m x ) = Tr{P mx ρ th }, where P mx = s P s,mx and P s,mx is the projector on all states with S x |ψ = m x |ψ 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, 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 S x = 0 is naturally associated with a non-vanishing expectation value of a −g/ω c S x , similar to what is expected for the superradiant phase in the Dicke model. Note, however, that this expectation value is real and corresponds to a finite charge (or displacement field) Q ∼ a + a † . The relevant cavity observables, V ∼ A + A † and Φ ∼ i A † − A , are not affected by this transition [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 dynamical 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 different values of g/ω c . In this plot, the exact analytic results are compared with a modified mean-field theory, where in Eq. (42) the bare dipole frequency ω 0 is replaced by the renormalized frequencỹ ω 0 given in Eq. (32). From this comparison we find that the variational free energy F V captures the overall trend, although 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 couplings 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.

V. BLACK-BODY RADIATION
The emission spectrum of a hot body was one of the first examples that could not be explained by combining the otherwise very successful theories of statistical mechanics and electromagnetism. In the correct quantum statistical derivation of the black-body spectrum it is assumed that the EM field thermalizes through weak interactions with the material, but that it can be treated otherwise as a set of independent harmonic modes. Therefore, it is particularly interesting to see how the thermal emission spectrum of a cavity mode changes under strong and USC conditions [62][63][64].

A. 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. [51]). By assuming that the transmission 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 D) where ω nm = (E n − E m )/ are the transition frequencies between the eigenstates |E n of H cQED with energies E n . In Eq. (43), κ denotes the decay rate of the bare cavity into the transmission line. In addition, we have introduced a phenomenological rate γ to account for a small but finite thermalization rate with the surrounding bath. For consistency we require κ γ 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 coupling strength g and for different temperatures. For small couplings, g ω c , we see the expected splitting of the unperturbed cavity resonance into two polaritonic resonances at frequencies ω ± ≈ ω c ±G/2. Although the lower polariton mode has a higher thermal occupation, the upper branch is slightly brighter. This observation can be partially explained by the 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. [61]) 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 frequency. 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.
Note that in previous studies of the absorption and emission spectra of the EDM [48,64] or the thermal radiation spectrum of the Rabi model (N = 1) [63] 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 statistics of the emitted photons can be sub-Poissonian for a certain range of couplings and temperatures [62]. We don't find such a behavior for larger numbers of two-level dipoles.  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. (46) and (47)] 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 system (P rad ) and from the bare cavity (P 0 rad ) as a function of temperature. The solid lines are obtained from exact numerical calculations for N = 6 and the dashed lines show the corresponding results predicted by Eq. (48) based on a Holstein-Primakoff approximation.

B. 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, H em . Here, the total power is obtained from Eq. (43) by integrating over all frequencies, For an empty cavity, this expression reduces to P 0 rad = ω c κN th , which we use to normalize the power values. Interestingly, for moderate temperatures we find a rather counterintuitive behavior: While the average energy that is stored in the mode increases for intermediate couplings, the emitted power decreases at the same time.
To explain this behavior we consider moderate values of G ω c and low temperatures. In this case we can use a Holstein-Primakoff approximation [60] and replace the spin operators σ i − by bosonic annihilation operators c i . The resulting linearized Hamiltonian can then be diagonlized and written as H cQED H HP − ω c /2 − N ω 0 , where 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 decoupled from the cavity due to symmetry. In terms of these polariton op-erators we can write 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 We see that there are various competing effects. With increasing coupling G, the frequency ω + and the matrix 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. (48) indeed predicts the observed decrease in P rad for increasing values of G. However, for higher and lower temperatures the dependence can also be reversed. In particular, for k B T /( ω c ) 1 the occupation number of the bare cavity mode is exponentially suppressed. However, for G ∼ ω c , also the lower polariton frequency is strongly reduced and N th (ω − ) ≈ k B T /( ω − ). Under 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 We see that the prefactors for the thermal contributions in the first line of this equatoin 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. (49). This part does not contribute to the radiated power such that overall P rad and H em display a very different dependence on G.

VI. CONCLUSIONS
In summary, we have analyzed the basic thermal properties of cavity QED systems in the USC regime. By using various analytic approximations and exact numerical results for a moderate number of two-level dipoles, we have derived the coupling-induced corrections 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 total radiated power, the corrections to material properties are small and scale only with the singledipole coupling strength as ∼ g 2 /ω c . This means that major modifications of ground-state chemical reactions or cavity-induced shifts of ferroelectric phase transitions cannot be simply explained by a strong collective coupling to a single quantized mode. To identify the detailed origin of such effects further experimental and theoretical investigations are still required, in particular also on the influence of the metallic boundaries on electrostatic interactions [33,65] and other modifications of the background EM environment [15].
The behaviour of cavity QED systems changes completely once the single-dipole coupling parameter, g/ω c , becomes of order unity. For this regime we have shown in terms of several explicit examples that the coupling to a single cavity mode can strongly modify material properties such as the specific heat, the susceptibility or the ferroelectric phase transition temperature. While this regime is currently not accessible with atomic, molecular or excitonic cavity QED systems, such coupling conditions can be reached in superconducting quantum circuits, where also the effect of temperature is more relevant than in the optical regime. Therefore, such artificial light-matter systems constitute an ideal testbed for studying strongly-coupled quantum systems with unconventional thermodynamical properties. Potentially, this can also lead to more accurate descriptions of fundamental thermodynamical processes or the optimization of quantum thermal machines, for which cavity QED and collective spin models have already been considered as simple toy systems in the past [66][67][68][69][70][71].
where Σ x = S x . Under this approximation, the resulting free energy is given by In the paraelectric phase the free energy has only a single minimum at Σ x = 0, while the ferroelectric phase is characterized by the appearance of two degenerate minima at a finite value of Σ x . The transition between the two phases is given by Eq. (42).
Equation (A2) allows us to derive some useful relations for the spin expectation values of the LMG model. From the condition ∂F MF LMG /∂Σ x = 0 we find For J < J c a solution Σ x = 0 exists. Subsequently, an expression for Σ z = Tr{S z ρ MF LMG } can be derived via Eq. (A3) and Eq. (A4) are transcendental and therefore, in general, no explicit expressions for Σ x and Σ z can be found.
Appendix B: 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, Here, · 0 denotes the average with respect to the thermal state of H 0 , H g (τ ) = e τ H0 H g e −τ H0 and T is the time-ordering operator along the imaginary time axis. In systems where H g ∼ g is only a small perturbation to H 0 we can use a cumulant expansion of the expectation value in terms of the F (n) g ∼ g n .

Weak coupling limit
We first apply this perturbation theory to the coupling Hamiltonian H g as defined in Eq. (20). Since in this case the thermal average of the linear term vanishes, a + a † 0 = 0, the lowest-order correction to the free energy is given by where For interacting dipoles the spin expectation value must be evaluated numerically. In terms of the energies E n and eigenstates |E n of H dip we obtain For noninteracting dipoles this calculation can be carried out analytically using and S ± S ∓ 0 = N [1 ∓ tanh ( ω 0 /2k B T )] /2. Altogether and writing F (2) g = N g 2 /(4ω c )f g we obtain where After some further simplifications the expression for f g reduces to the result given in Eq. (22).

Bound on the free energy correction
From the general expression for the correlation function given in Eq. (B4) one can show that S x (τ 1 )S x (τ 2 ) 0 ≤ S 2 x 0 for τ 1 ≥ τ 2 . This inequality can be used together with (N th + 1)I(−ω c ) + N th I(ω c ) = 1 in the second line of Eq. (B2) in order to derive the upper and lower bounds 0 ≤ f g ≤ 4 S 2 x 0 /N . To improve the upper bound we can repeat the whole perturbation calculation in a displaced frame,H cQED = D † (α)H cQED D(α) = H 0 +H g (α), where D(α) = e αa † −α * a is the displacement operator. Specifically, by choosing α = −g S x 0 /ω c we obtaiñ where ∆S x = S x − S x 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. (25).

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. III C and decompose the total Hamiltonian asH cQED = H 0 + H ω0 . Here . According to this partitioning, the bare Hamiltonian H 0 is diagonal in the σ x basis and H ω0 0 = 0. Thus, up to second order in ω 0 the free energy is given by The correlation function for the photons can be calculated exactly using the properties of the displacement operator. Specifically, given a pair of complex number z 1 and z 2 the following formula holds The two-point correlation functions can then be expressed in terms of an infinite series, K rq e (q−r)ωc(τ1−τ2) (B12) and Q rq e (q−r)ωc(τ1−τ2) , with coefficients g ω c 2(r+q) (1 + N th ) r N q th r! q! , g ω c 2(r+q) (1 + N th ) r N q th r! q! . (B14) Altogether we obtain × K rq S z (τ 1 )S z (τ 2 ) 0 + Q rq S y (τ 1 )S y (τ 2 ) 0 , whereω 0 = ω 0 exp[−g 2 /(2ω 2 c )(1 + 2N th )]. In Fig. 3(a) this result is used to evaluate χ z for noninteracting dipoles and, as expected, we find that it is in perfect agreement with the full numerical simulations 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 processes with multiphoton scattering, emission and absorption. These processes become more and more important as the coupling strength increases.

Appendix C: Variational LMG Hamiltonian
In the variational ansatz described in Sec. III C, the trial Hamiltonian H * given in Eq. (31) 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 where we have already assumed S y 0 = 0. The thermal expectation value of the cosine operator can be evaluated with the help of Eq. (B11). Altogether, the variational free energy F V can be written as (C2) Taking the derivative of F V with respect toω 0 yields the extremal condition Note that this result is independent of the dipole-dipole couplings J ij .

Appendix D: Emission spectrum
For C t C, the capacitive coupling between the cavity and the transmission line shown in Fig. 1 can 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 = 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 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 t 0 dt e −iω k (t−t ) X(t ). (D2) Therefore, by assuming that initially the transmission line is in the vacuum state, we find that (D3) To proceed we write X(t) = X − (t) + X + (t), such that X + (t) = n≥m X nm e iωnmt contains only contributions that oscillate with a positive frequency, ω nm = (E n − E m )/ ≥ 0, and X − (t) = X † + (t) [62]. Then, for long enough times t γ −1 , where γ is the characteristic thermalization rate of the cavity QED system, we obtain The total power radiated into the transmission line is and by writing P rad = ∞ 0 dω S bb (ω) we obtain the general expression for the black-body spectrum S bb (ω) = ωJ(ω)C X (ω). (D6) Here J(ω) = k λ 2 k δ(ω − ω k ) is the spectral density of the transmission line and 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 and we obtain the total emitted power given in Eq. (44).
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 surrounding thermal bath, which induces a finite broadening of all the spectral lines. A detailed investigation of such line-broadening effects is beyond the scope of the current analysis. Instead, we simply replace 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. (43).