The Vacua of Dipolar Cavity Quantum Electrodynamics

The structure of solids and their phases is mainly determined by static Coulomb forces while the coupling of charges to the dynamical, i.e., quantized degrees of freedom of the electromagnetic field plays only a secondary role. Recently, it has been speculated that this general rule can be overcome in the context of cavity quantum electrodynamics (QED), where the coupling of dipoles to a single field mode can be dramatically enhanced. Here we present a first exact analysis of the ground states of a dipolar cavity QED system in the non-perturbative coupling regime, where electrostatic and dynamical interactions play an equally important role. Specifically, we show how strong and long-range vacuum fluctuations modify the states of dipolar matter and induce novel phases with unusual properties. Beyond a purely fundamental interest, these general mechanisms can be important for potential applications, ranging from cavity-assisted chemistry to quantum technologies based on ultrastrongly coupled circuit QED systems.


Introduction
In QED, the relevant dimensionless coupling parameter is the finestructure constant α fs = E C /E ph . It can be expressed as the ratio between the Coulomb energy, E C = e 2 /(4πε 0 d), of two electrons at a distance d and the energy E ph = cħ h/d, which is needed to create a photon confined approximately to the same region in space. The small value of α fs 1/137 already suggests that the quantized modes of the electromagnetic field play a minor role in the physics of atoms, molecules and solids, as confirmed by more rigorous calculations. However, this argument does not necessarily hold in structured electromagnetic environments, such as nanoplasmonic systems or LC circuits, where the energy of a photon can be tuned independently of its wavelength. In this case the coupling between an electric dipole and a quantized field mode is characterized by an effective parameter α = α fs (Z/Z 0 ) [1,2], which can be considerably enhanced by increasing the impedance of the mode, Z, compared to the value in free space, Z 0 . This raises an important fundamental question: Can the properties of matter be influenced by such an artificially boosted coupling to the quantized field, and, if so, how would the properties change?
In the following analysis we address this open theoretical problem by considering the conceptually simplest scenario of a lattice of interacting two-level dipoles coupled to a single cavity mode. While this prototypical cavity QED system has been the primary workhorse for modelling light-matter interactions in quantum optics and solid-state physics for many decades, its ground state properties are, surprisingly, still unknown. Here we use exact numerical calcula- tions to evaluate, without any approximations, the non-perturbative effect of vacuum fluctuations on the ground states of strongly-correlated dipolar systems. This analysis confirms, first of all, the existence of so-called superradiant [25][26][27] and subradiant [2,23] phases, which already appear in the analysis of simple collective spin models. However, we also observe the formation of completely new phases of dipolar matter and cavity-induced ordering mechanisms, which have not been discussed in the literature before. By that, we are able to provide a first complete phase diagram for the 'vacua' of dipolar cavity QED in elementary lattice geometries. This study also reveals that there is still a wealth of unexplored phenomena in cavity and circuit QED, which may soon become accessible with further experimental advances in these fields.

Cavity QED of interacting dipoles
We consider a prototypical cavity-QED system as depicted in Fig. 1(a), where N anharmonic dipoles are coupled to a single quantized mode of an LC resonator [29] with frequency ω c . We approximate the dipoles by two-level systems with transition frequency ω d , located at fixed lattice positions r i (in units of the lattice spacing a). The dipoles couple to the electric field E of the cavity with a transition dipole moment µ E and among each other via static dipole-dipole interactions. Under these assumptions the quantized dynamics of the system is described by the Hamiltonian [2] (ħ h = 1) where σ i α denote the Pauli operators at site i, S α = i σ i α /2 are collective spin operators, and a is the annihilation operator for the field mode. Note that H cQED represents light-matter interactions in the dipole gauge, where gauge-dependent artefacts in the two-level truncation can be avoided [48]. For other recent contributions on this topic, see [49][50][51][52].
The cavity affects the dynamics of the dipoles in two different ways. First, in Eq. (1) the dipole frequency, ω d , and the dipole-dipole couplings, J i j , already include screening effects from the metallic boundaries and can differ considerably from their bare values ω 0 d and J 0 i j in free space. This behaviour is illustrated in Fig. 1(b-d), which shows that the usual dipoledipole interactions, J 0 i j = J 0 /| r i − r j | 3 , become short-ranged and suppressed as the distance D between the plates is decreased. This boundary effect can strongly modify the properties of para-and ferroelectric systems [2,[53][54][55][56], but it is of electrostatic origin and not the main focus of this study. Therefore, we simply treat ω d and J i j as arbitrary model parameters and investigate the additional modifications induced by the coupling to the dynamical field mode with frequency ω c ∼ ω d . For a sufficiently homogeneous mode, these effects are described by the collective dipole-field coupling g a † + a S x , with a single-dipole coupling strength g, and the associated depolarization shift ∼ S 2 x . We emphasize that the Hamiltonian H cQED represents a minimal consistent extension of the more commonly used Dicke model in quantum optics. In particular, the inclusion of the so-called P 2 -term ∼ S 2 x [2,17,23,29,57] ensures that for ω d → 0 we recover the correct electrostatic limit, H cQED ω c a † a+ i< j , which is not the case in the aforementioned Dicke model or variants thereof [32,33,41]. Therefore, although being based on several simplifying assumptions, H cQED still respects all Maxwell's equations and allows us to treat electrostatic and electrodynamical interactions in a consistent and transparent manner. For a more detailed derivation and justification of this model in both cavity and circuit QED systems see Refs. [2,23].

The ground states of cavity QED
The physics of H cQED and variants thereof has been studied extensively in quantum optics and solid-state physics, but primarily in the regime g/ω c 1. In this limit, the system may still feature huge collective Rabi-splittings of Ω R = g N ∼ ω c in the excitation spectrum [3,5,6], but qualitative changes in the ground and equilibrium states are still only perturbative [2,16,[45][46][47]. In turn, preceding studies of H cQED in the non-perturbative coupling regime, g/ω c 1, have been restricted to very few spins or the special case of all-to-all interactions J i j = J [2,23,47]. This strongly reduces the computational complexity of the problem, but also ignores all non-collective correlations, the influence of the lattice geometry and other essential effects. Here we perform exact, large-scale numerical simulations of finite-sized dipole systems to obtain the ground states of H cQED without any further approximations [see Appendix A for details].
Phase diagram. In Fig. 2(a) we first show the ground state phase diagram for N = 26 dipoles on a square lattice with nearest-neighbor only couplings J i j = Jδ 〈i, j〉 and ω d = ω c . For g = 0 the cavity is completely decoupled and Eq. (1) reduces to the familiar transverse field Ising (TFI) model. In this limit we observe the expected transition from a disordered paraelectric to the ordered ferroelectric or antiferroelectric Néel phase when |J| exceeds a critical value of |J * | ≈ 0.7ω d , which agrees within a few percent with the transition point of the infinite system [58].
Although on finite-size systems symmetries cannot be broken spontaneously and the order parameters are strictly zero, the ordered phases can be uniquely identified through the correlations 〈σ i x σ j x 〉 between the spins at positions i and j. For that, we define the (normalized) structure factor for a momentum k in the Brillouin zone of the lattice, and r i0 = r i − r 0 . In an ordered phase Σ x ( k) shows single peaks at specific momenta k * , which identify the ordering pattern, and the value of Σ x ( k * ) can be used to define the ordering strength. To identify the ferroelectric and staggered Néel ordered phases, we introduce which are nonzero in the corresponding phases, but vanishingly small elsewhere, as shown in Fig. 2(b). These correlation measures are related to the second moments of the order parameters of the respective phases, in particular where we have introduced the standard order parameters for the ferroelectric and Néel phases, respectively, Here, we have decomposed the lattice into two sublattices A and B to capture the staggered structure of the Néel phase and introduced the sublattice polarizations p I = i∈I σ x i /2. For finite g/ω c 1 the phases do not change considerably, except that in the ferroelectric state now also the photon number acquires a large expectation value, 〈a † a〉 (g N /(2ω c )) 2 . Eventually, in the thermodynamic limit N → ∞ the 2 symmetry of the model [see Appendix A] is spontaneously broken in the ferroelectric regime leading to a finite spin order parameter 〈S x 〉 = 0 and a coherent photonic state with amplitude α ≡ 〈a〉 = g/ω c 〈S x 〉 [see also Appendix B]. In the quantum optics literature one commonly refers to such a phase as 'superradiant' [25][26][27], a notation that we adopt in the following. The anti-aligned Néel phase, on the other hand, shows a staggered arrangement of dipoles with 〈S x 〉 = 0 in the thermodynamic limit. While such a state breaks the 2 symmetry (in a non-trivial way together with lattice symmetries), the photon sector obeys 〈a〉 = 0, such that we refer to it as 'subradiant'. On finite-size systems, this leads to a strongly reduced photon number 〈a † a〉 0 compared to the fluctuating paraelectric phase [see Fig. 2 Collective subradiant phase. The paraelectric phase gradually evolves into a new 'collective subradiant' phase with unusual properties for g/ω c 3. This phase exhibits no order and . At the same time, also the photon number 〈a † a〉 vanishes, indicating that all the dipoles are still anti-aligned.
These seemingly contradicting properties can be understood by looking at the limit J ≈ 0 and g/ω c 1, where we can eliminate the photons using strong-coupling perturbation theory. The remaining low-energy physics of H cQED is then approximately described by the effective spin model [2,23] where S = (S x , S y , S z ). This Hamiltonian describes Ising interactions with J i j = Jδ 〈i, j〉 subject to a renormalized 'transverse field' h z = ω d exp −g 2 /(2ω 2 c ) and with an additional cavitymediated collective coupling of strength In the considered regime where the collective subradiant state is observed the term ∝ J c dominates over the exponentially suppressed h z . Thus, for J = 0, the Hamiltonian is minimized by a perfectly anti-aligned state |ψ cs 〉, which obeys S x |ψ cs 〉 = 0 and has maximal total spin S = N /2. Compared to a Néel-ordered configuration, this highly entangled state represents an equal superposition of all possible combinations with half of the dipoles pointing along x and the other half into the opposite direction, without any spatial order. Surprisingly, this peculiar collective phase survives to a high degree in the presence of competing short-range interactions (J = 0) and it is separated from both ordered phases by a sharp transition.
Although both the collective and the Néel ordered phases are subradiant, a crucial difference between them is visualized in Fig. 2(c-d). These two plots compare the photon number 〈a † a〉 and 〈σ x σ x 〉 stag for fixed J/ω d , but increasing g/ω c . For a value of J = 0.5 ω d the system then directly transitions from the paraelectric into the Néel phase. This is signified by an increase of 〈σ x σ x 〉 stag and simultaneously, or better to say as a consequence of that, also the photon number vanishes. In the evolution from the paraelectric to the collective subradiant states, as shown exemplarily for J = 0, the staggered polarization is always vanishingly small, but the dipoles still completely decouple from the photons for large g/ω c . The formation of such a collective subradiant phase is thus much more subtle than an interaction induced spatially ordered anti-alignment of dipoles.
Estimating phase boundaries. After establishing the different phases which appear as vacua of H cQED on the square lattice, let us in the following estimate their boundaries in the phase diagram. The ordered phases (ferroelectric and Néel) are separated from the disordered phases as shown in Fig. 3. The finite-size results for g = 0 agree within a few percent with the known phase boundaries in the thermodynamic limit. With increasing g we observe a significant cavity-induced reduction of the critical coupling strength |J * (g)| for both ordered phases. Within the limited range of available system sizes we still see a weak dependence of the phase boundaries on the particle number N [see Fig. 2(a)], but also for larger N no significant qualitative changes of the phase transitions are expected. Compared to the g = 0 results, where the transitions are known to be continuous and in the (2+1)D Ising universality class, with increasing g/ω c the width of the peaks in the fluctuations shrinks, which indicates a narrowing of the critical region. Moreover, for the transition between the ferroelectric and disordered phases, the shape of the fluctuations changes when g/ω c 4, while it remains the same (up to rescaling) for the transition between the Néel and the disordered phases, as shown in Fig. 3(c). Although numerical evidence is still limited, this behavior indicates a cavity-induced change from a continuous to a first order phase transition in the ferroelectric case, while the transition into the Néel phase remains continuous for all g/ω c .
The evolution from the paraelectric to the collective subradiant regimes, on the other hand, does not show any trend towards a non-analytic behaviour in our analysis. Since also both of these regimes do not break any symmetries, we expect this evolution to be better described by a smooth crossover, instead of a sharp phase transition. As such there is no well-defined transition line and in our plots we use the characteristic peak in the polaron photon number 〈a † a〉 polaron to define a crossover boundary [see Appendix F]. For this crossover we observe a non-negligible dependence on N , which in detail depends on the observable that is used to define the boundary. In the non-interacting case J = 0 the evolution of the crossover bound-ary can be numerically estimated also for much larger values of N [see Appendix F]. These simulations further support a smooth crossover between the paraelectric and the subradiant phase, but do not yet provide a conclusive picture about the behaviour of this crossover in the limit N → ∞.
Finally, let us emphasize that H cQED describes the coupling of all dipoles to a single cavity mode with fixed properties, in particular, a fixed interaction region. Simply increasing N does not represent a well-defined thermodynamic limit, while a rescaling of the coupling constant, g → g/ N , would render the dipole-field coupling non-perturbative. Therefore, the finitesize phase diagrams discussed so far and in the following are already representative for the practically relevant scenarios, where small or mesoscopic ensembles of dipoles are coupled to a field mode localized within a tiny mode volume.

Order and fluctuations
From the analysis above we can extract two basic cavity-induced many-body effects, namely the stabilization of phases with pre-existing order and the suppression of fluctuations in the disordered phase through the formation of highly entangled collective states. One thus expects that also in general cavity-induced modifications will be most significant in situations where strong fluctuations occur already in the bare system. As a prototypical example we now consider the ground state phases of H cQED for repulsive dipoles on a triangular lattice [see also Fig. 6], where we again assume nearest-neighbor interactions J i j = Jδ 〈i, j〉 . In this configuration the dipole-dipole interactions are strongly frustrated and lead, for g = 0, to large fluctuations (in S x ) even in the ordered ground states at J > 0. As shown in the corresponding full phase diagram in Fig. 4(a), under this condition completely new cavity-induced phases appear at sufficiently large g/ω c .
To understand these observations, let us first summarize the established results of the frustrated TFI model at vanishing coupling g = 0. In the classical limit J/ω d → ∞, the strong frustration prevents the spins from ordering even at zero temperature and the model exhibits an exponentially large (in N ) ground-state manifold with a finite T = 0 entropy density [59]. However, quantum fluctuations from a transverse field, i.e. the term ω d S z in Eq. (1), select an ordered subset of states in an "order-by-disorder" (OBD) process [60]. As shown in the inset in Fig. 4(a), the selected three-sublattice (3SL) antiferroelectric state [61,62] can be depicted as an arrangement of anti-aligned dipoles on two of the sublattices (in the S x direction), while on the third sublattice the dipoles align (paraelectrically) with the transverse field and do not point in any particular direction according to S x . This phase is thus characterized by a long-range 3SL order, while it still exhibits strong fluctuations in S x [see Fig. 4(c), left panel]. When the interaction strength J is decreased below a critical value, the 3SL phase eventually becomes unstable towards a disordered, paraelectric phase. The phase transition is continuous and features an emergent O(2) symmetry, such that the universality class is not of Ising, but of XY type [61,62].
To investigate the properties of this system for g/ω c > 0, we compute the correlations corresponding to 3SL order and find an extended region above a critical line J * (g) where they become large. This indicates the stability of the 3SL phase also in the presence of the cavity mode [see Fig. 4(a-b)]. Similar to the square lattice case, increasing the coupling to the cavity reduces the critical value J * (g), which we estimate by maxima in the fluctuations 〈∆|p 3SL |〉 for the complex 3SL order parameter where the lattice was decomposed into the three sublattices A, B, C.
For J < J * (g) we observe the crossover between the paraelectric and the collective subradiant phase also for the triangular lattice, since the geometry becomes irrelevant in this regime. While the formation of such a homogeneous state is hindered by the 3SL order above the transition line J > J * (g), we discover a new type of '3SL subradiant' regime for g/ω c 3. This regime is characterized by a finite order, 〈σ x σ x 〉 3SL > 0, and is thus separated from the collective phase by a sharp transition line [see Fig. 4(a-b)]. At the same time it differs from the normal 3SL phase in terms of its vanishing photon number, 〈a † a〉 0, which indicates strongly reduced polarization fluctuations. This difference can be clearly seen in the ground state distribution of S x -values in the two regimes, as shown in Fig. 4(c): while the polarization distribution is broad in the normal 3SL phase, it is pinned to a single value of S x = 0 deep in the subradiant regime. This behavior can be intuitively explained, by adopting again the simplified picture of a 3SL state, where the fluctuating dipoles on one sublattice participate in the formation of a collective subradiant configuration, similar to the state |ψ cs 〉, while the two S x -polarized sublattices remain unaffected. Note, however, that this is only an oversim-plified picture of the actual state, where correlations among different sublattices do not vanish completely.

Order by cavity-induced disorder
A very surprising finding in the case of a triangular lattice is the appearance of an additional superradiant phase for antiferroelectric dipole interactions [blue region in Fig. 4(a)]. As shown by the histogram in Fig. 4(c), also in this phase the polarization is well-defined, but assumes non-zero values S x = ±1, 2, . . . , and consequently, 〈a † a〉 = (g/ω c S x ) 2 1. Although this value is much smaller than in the regular superradiant phase, this property is completely unexpected, since, at first sight, in this regime both the direct dipole-dipole as well as the cavityinduced interactions would favor a fully anti-aligned configuration. As shown in Fig. 4(d) the transition into this phase is associated with a sharp peak in the photon-number fluctuations, 〈∆a † a〉, which also remain finite within this phase.
To further investigate the properties of this new type of superradiant states we focus on the relevant regime g/ω c 1, where, as shown above, we can eliminate the photons using strong-coupling perturbation theory to obtain the effective spin model H S [Eq. (6)]. Although it is derived under the assumption J i j → 0, a comparison with full numerical simulations up to N = 24 shows that H S accurately reproduces the qualitative features of H cQED for large J i j , as long as g/ω c 3 [see Fig. 4(a)].
As already discussed above, the regular OBD process on the frustrated, antiferroelectric Ising interactions (for g/ω c 3) is driven by the perturbation with a transverse field ∝ S z [i.e. H S with artificially set J c = 0], which stabilizes the normal 3SL phase. On top of that, the cavity-mediated collective term ∝ J c can induce a crossover into the 3SL subradiant regime. The appearance of a superradiant phase suggests that this hierarchy no longer holds for large J/ω d and g/ω c 3. In this regime, the transverse field h z is already strongly (exponentially) suppressed and the cavity-mediated collective term ∝ J c is the dominant perturbation on the Ising interactions. Its quantum fluctuations select a distinctly ordered subset of states compared to the regular OBD process and yield a superradiant phase with a different 3SL ordering pattern.
A common way to analyze different ordering patterns on the triangular lattice is to look at the distribution of the 3SL order parameter p 3SL = p 3SL e iθ (10) in the complex plane [see also Appendix E]. In Fig. 5(a) we use this method to represent the different ground-state structures in the normal, the subradiant and the superradiant 3SL phases. The large values of 〈|p 3SL |〉 show that all three phases exhibit 3SL order, as also seen from 〈σ x σ x 〉 3SL in Fig. 4(b), but with different levels of fluctuations. Even more, the pattern for the superradiant phase differs qualitatively from the other two plots, in particular the positions of the largest peaks are shifted, and indicate a configuration where dipoles in one sublattice are (almost) fully polarized in S x , while dipoles on the other two sublattices are equally, but only partially polarized along the opposite direction. For N = 24 this ordering leaves a residual net polarization S x = ±1.
Within the effective spin model we can investigate the superradiant phase also for larger lattices and find that this residual polarization increases with the system size and leads, with increasing ratio J/J c , to a whole series of superradiant phases, characterized by |S x | = 1, 2, 3, . . . , S max From this analysis we can extract a linear scaling for the maximal ground state polarization S max x = δ · N and the photon number 〈a † a〉 = (gδ · N /ω c ) 2 with a polarization density of δ 0.07. Interestingly, very similar distributions of p 3SL and a finite net polarization (although at a much smaller value of δ) have been previously discussed in connection with supersolidity in frustrated spin systems [63][64][65][66][67], where magnetic and superfluid order parameters coexist. While outside the scope of the current study, this connection between superradiance and supersolidity in cavity QED is a particularly exciting direction to explore further. Finally, we want to point out that the interplay between short-range and cavity-mediated collective interaction, such as in H S , has been a source to observe unexpected superradiant and supersolid phases in experiments with cold atoms in cavities [68,69]. While the interactions and the nature of the phases in these systems are different, interesting connections between the different models might be explored in future works.

Conclusions
In summary, we have addressed the many-body problem in cavity QED, which arises from the interplay between short-range electrostatic interactions and the non-perturbative coupling to a common cavity mode. Based on exact numerical calculations, we have obtained a first complete phase diagram for the 'vacua' of cavity QED covering the full range of dipole-field interaction strengths, from vanishingly small to ultrastrong. By taking the influence of short-range dipole-dipole interactions in different lattice geometries fully into account, we have shown how the competition between conventional and cavity-induced correlations can lead to the formation of several novel phases, which have no direct counterparts in the collective Dicketype models [25][26][27]30] usually studied in quantum optics, nor in regular solid-state spin systems. Although we focused here on the conceptually simplest case of two-level dipoles, the basic mechanisms identified in this work, i.e. the cavity-induced reduction of fluctuations, extended subradiant states without order, or a cavity-induced OBD process, will also be rele-vant in various other strongly interacting cavity QED systems [36][37][38][40][41][42], where so far the analysis of ground states has been constrained to mean-field methods or perturbative coupling regimes.
While the most interesting regime of light-matter interactions, g/ω c > 1, is not accessible in cavity QED experiments with atoms and molecules today, recent advances in the fabrication of electromagnetic resonances with Z/Z 0 > 50 [13,70] show that this regime is by no means out of reach. Of immediate relevance are our findings for the field of circuit QED, where equivalent models can also be obtained with alternative galvanic coupling schemes, and values of g/ω c > 1 [11,12] have already been demonstrated for single superconducting two-level systems. The extension of these experiments to larger arrays of superconducting qubits coupled directly and via microwave modes will provide a natural platform to explore the new phases and physical mechanisms identified in this work. Such systems are currently developed for quantum simulation and quantum annealing schemes [71], where ultrastrong coupling effects, similar to what we have analyzed here, can find direct practical applications [72].

Acknowledgements
We thank Alessandro Toschi for stimulating discussions and feedback on the manuscript.

A Numerical simulations
The numerical results in this manuscript have been achieved by Exact Diagonalization using a Lanczos algorithm [73,74] on regular clusters with a finite number of N two-level systems [see Fig. 6]. To reduce finite-size effects we use periodic boundary conditions along both directions of the square and triangular lattices and, to study genuine two-dimensional properties we only use clusters with an aspect ratio ε = 1, i.e. the loops around both periodic directions have equal length. To fit the antiferroelectric phases with a two (three) sublattice structure, we only consider square (triangular) clusters with N mod 2 = 0 (N mod 3 = 0). Here, it is worth mentioning that the subradiant states discussed in this work cannot exist on clusters with odd N , where the possible total polarizations S x are half-integers, so that a subradiant state with fixed S x = 0 cannot be obtained.
The Hilbertspace H is kept finite by additionally introducing a photon-number cutoff n max ph for the cavity mode in H cQED , such that a † |n max ph 〉 ≡ 0 and dim[H] = 2 N n max ph + 1 . n max ph has to be chosen large enough to achieve accurate results throughout the different regimes of the external parameters [see Appendix C for a thorough discussion]. It can also be favorable to transform H cQED into the polaron frame, which describes a distinct basis, where, in some regimes, a much smaller cutoff than in the original frame can be sufficient [see Appendix B].
To further reduce the Hilbertspace dimension, we use the 2 symmetry of H cQED , given by the operator S = e −iπ(a † a+S z ) , together with the lattice translational and point-group symmetries to block-diagonalize H.
since a ( †) → a ( †) − g/ω c S x is displaced proportional to S x . Within this formulation, it becomes obvious that the correct electrostatic limit H cQED ω c a † a + i< j J i j 4 σ x i σ x j is achieved for ω d → 0, since the depolarization shift in Eq. (1) exactly cancels the additional terms ∝ S 2 x from the transformation of a † a.
The Hamiltonian Eq. (11) can also be advantageous to study superradiant phases, because the photon number 〈a † a〉 polaron in this polaron frame ignores the part from the direct coupling to the polarization S x and remains much smaller than in the standard frame, in particular for superradiant phases. Therefore, a substantially lower photon number cutoff n max ph can be sufficient for precise numerical simulations, with the disadvantage of having to deal with a dense photonic Hamiltonian, when ω d = 0.
Also, the polaron photon number, which can be computed by 〈a † a〉 polaron = (a † − α)(a − α) with α = g/ω c S x in the standard frame, can be a useful observable. In particular, we use characteristic peaks in the polaron photon number to identify the crossover regime between the paraelectric and collective subradiant phases [see Appendix F].
Using strong-coupling perturbation theory for g/ω c 1 and projecting onto the lowestenergy sector without polaronic excitations |0〉 polaron ph , the last term in H cQED can be approximated as [2,23] where we have introduced the total spin operator S = (S x , S y , S z ). Within this approximation, we thus obtain the effective spin model H S given in Eq. (6). It is important to note that the eigenstates in the original basis |Ψ〉 = e −g/ω c S x (a † −a) |Ψ〉 spin ⊗|0〉 polaron ph , and the photon number 〈a † a〉 = g 2 /ω 2 c 〈|S x |〉 2 is generally non-zero in the standard basis. Note that the approximate expression in Eq. (12) has been derived for non-interacting dipoles J i j = 0.

C Photon number cutoff
In this appendix we discuss the implications of introducing a photon number cutoff n max ph to obtain a finite Hilbert space. This cutoff has to be chosen large enough such that the true ground state in the full Hilbert space only shows negligible deviations (up to some defined precision) when it is projected into the restricted Hilbert space. Appropriate values for n max ph strongly depend on the chosen external parameters, i.e., for parameters belonging to a superradiant phase much larger cutoffs have to be chosen than for parameters belonging to a subradiant phase.
For the simulations in this work we choose n max ph large enough, such that doubling this cutoff does not change the measured observables. In Fig. 7 we show an analysis of the dependence of observables on n max ph for a square lattice configuration with N = 16 dipoles and a constant coupling g/ω c = 2. For antiferroelectric interactions, J/ω c > 0, a small cutoff n max ph = 64 is already sufficient to obtain converged results in all observables, since the average photon number 〈a † a〉 remains small. Contrarily, for ferroelectric interactions, J/ω c < 0, the average photon number and its fluctuations become large [see Fig. 7(a)] when the superradiant regime is entered. Then, a too small cutoff yields false results not only for 〈a † a〉, but also for pure dipole observables [see Fig. 7(c)-(f)]. The distribution of the photon number a † a in the ground state [see inset in Fig. 7(a)] reveals, that a cutoff n max ph 400 would be sufficient for a simulation with these particular external parameters.
In Fig. 7 we also show results obtained from the polaron frame with Hamiltonian H cQED . In this formulation, a much smaller cutoff is already enough to obtain converged results in the ferroelectric regime, since the polaron photon number 〈a † a〉 polaron remains small even in the superradiant phase [see Fig. 7(b)].

D OBD simulations
The classical Ising model, which is obtained from H cQED or H S for J → ∞, is strongly frustrated on the triangular lattice with nearest-neighbor interactions and does not order even at zero temperature T = 0. It features an exponentially large (in N ) ground-state manifold, with an extensive T = 0 entropy S ≈ 0.323k B N [75]. This ground-state manifold can be destabilized by quantum fluctuations from other interaction terms, such as a transverse field or the cavity-mediated collective coupling ∝ J c in H S , when a particular subset of states, with the softest response to the fluctuations, is selected. If the set of the selected states is ordered, this process is termed "order by disorder" [60], and for a perturbation with a transverse field this mechanism is known to induce a 3SL ordered phase [61,62].
To study the OBD process from the collective coupling ∝ J c in Eq. (6) in the J/J c → ∞ regime, we restrict the Hilbertspace of the system to the degenerate, classical ground-state manifold, and define the effective Hamiltonian Here, P J is the projector onto the classical ground-state manifold. The low-energy eigenvectors of H OBD yield the states stabilized by the OBD mechanism (in the sense of a first-order degenerate perturbation theory), from which we can compute the observables, as shown in Fig. 5. The advantage of this approach is that, compared to a full simulation of H S , larger system sizes N can be simulated.

E Histograms of the three-sublattice order parameter
In this appendix we illustrate the properties of the 3SL order parameter histograms in the complex plain. A (classical) state with sublattice polarizations p = (p A , p B , p C ) gives a single peak in the histogram according to the axes defined in Fig. 8(a). While this identification is not unique for any p, the strength of the 3SL ordering p 3SL is given by the distance of the peak from the center, and the angle θ (in the complex plain) indicates different types of 3SL ordering patterns. In particular, states of type p = (0, 1, −1) with two fully, but oppositely polarized, and a non-polarized sublattice (zero net polarization), give peaks at the centers of the hexagonal boundaries of the histogram. The six different peaks correspond to the possible permutations of the three sublattices. States of type p = (1, 1, −1), where all sublattices are fully polarized, one of them oppositely to the others (non-zero net polarization) have peaks at the vertices of the hexagonal boundary. The six peaks correspond to the possible permutations of the sublattices and the inversion of the polarization p → − p.
More generally, patterns of type p = (0, m, −m) yield peaks at angles θ l = π/6 + lπ/3, l ∈ {0, . . . , 5}, with a radius proportional to m, as shown in Fig. 8(b). A comparison with the full histograms for the normal and 3SL subradiant regimes [c.f. Fig. 5(a)] shows that the maxima for those phases correspond to such a pattern with maximal m = 1, as indicated by the red dots. Furthermore, we want to note that fluctuations of the non-polarized sublattice lead to a broadening of the six single peaks parallel to the edges of the hexagonal boundary, as seen for the normal 3SL phase. The strength of these fluctuations can be further used to distinguish the normal and the fluctuation-free subradiant 3SL regimes [c.f. Fig. 5(a)].

F Estimating crossover boundaries
In contrast to phase transitions with an abrupt change in the behavior of the order parameter (in the thermodynamic limit), crossovers between two regimes of states with different physical properties show a smooth change (if any) in all of the observables, since the ground state evolves smoothly. Therefore, the crossover region, or a 'boundary' between the two regimes, can typically not be determined uniquely, but depends on the chosen observable and the feature used to estimate the boundary.
To estimate a boundary between the paraelectric regime and the collective subradiant regime, we use maxima in the polaron photon number 〈a † a〉 polaron , as shown in Fig. 9(a). In comparison, the standard photon number 〈a † a〉 is a bad estimator for the crossover when J/ω c < 0, where the proximity to the superradiant phase spoils its characteristic features in the narrow collective subradiant regime.
Based on this definition, we observe a shift of the boundary to larger g/ω c with increasing  Fig. 2(a) and Fig. 4(a)]. (b-d) System size dependence of the crossover boundary for J = 0. The polaron photon number (b) and normal photon number (c) do not show any trend towards a non-analytic behaviour with system size N , which would be a sign for a phase transition. (d) We utilize the maxima in (b, c) to estimate the dependence of the crossover boundary g * /ω c with N . Characteristic for a crossover, the boundary depends on the choice of the observable used to define it.
the system size N [see also Fig. 2(a)]. The achievable system sizes in exact diagonalization are too small to make reliable predictions about the fate of the collective subradiant phase in the limit N → ∞ in general. However, for the non-interacting case J = 0 much larger system sizes can be analyzed, as shown in Fig. 9(b-d). Both the polaron and normal photon number increase with system size, but their characteristic shape remains unchanged, so that we do not observe any trend towards a non-analytic behaviour (also for other observables) when increasing N , making a phase transition scenario very unlikely. Again, we can use the maxima in these observables to estimate the size dependence of the crossover boundary as shown in Fig. 9(d). While this boundary depends on the choice of observable used to define it, we find a stable collective subradiant phase for all considered system sizes at large enough g/ω c , consistent with the analysis in [23].
Nevertheless, these numerical results do not provide a conclusive picture about the fate of the collective subradiant phase in the limit N → ∞. In this respect it is important to emphasize, that this thermodynamic limit is also not properly defined for the single-mode model used in this work, where H cQED is super-extensive. For finite systems, the single-mode approximation is, however, expected to capture the main results and our analysis can be directly applied to most near-term experiments, where intermediate-scale systems, far away from the thermodynamic limit, will be realized.
Because of the similarity with the evolution from the paraelectric to the collective subradiant regimes, we also expect the evolution from the normal to the subradiant 3SL regime to be described by a crossover instead of a sharp phase boundary. We characterize the regimes by a distinct strength of the polarization fluctuations 〈∆|p|〉, since the polarization fluctuates strongly in the 3SL normal regime and becomes pinned to S x = 0 in the 3SL subradiant regime [see Fig. 4(b)]. We, therefore, define the boundary by a rather sharp drop in 〈∆|p|〉 and estimate its location for constant J/ω d by a peak in the negative gradient of 〈∆|p|〉 with respect to the coupling g/ω c [see Fig. 10].