Photon Condensation and Enhanced Magnetism in Cavity QED

A system of magnetic molecules coupled to microwave cavities ($LC$ resonators) undergoes the equilibrium superradiant phase transition. The transition is experimentally observable. The effect of the coupling is first illustrated by the vacuum-induced ferromagnetic order in a quantum Ising model and then by the modification of the magnetic phase diagram of ${\rm Fe_8}$ dipolar crystals, exemplifying the cooperation between intrinsic and photon-induced spin-spin interactions. Finally, a transmission experiment is shown to resolve the transition, measuring the quantum electrodynamical control of magnetism.

A system of magnetic molecules coupled to microwave cavities (LC resonators) undergoes the equilibrium superradiant phase transition. The transition is experimentally observable. The effect of the coupling is first illustrated by the vacuum-induced ferromagnetic order in a quantum Ising model and then by the modification of the magnetic phase diagram of Fe8 dipolar crystals, exemplifying the cooperation between intrinsic and photon-induced spin-spin interactions. Finally, a transmission experiment is shown to resolve the transition, measuring the quantum electrodynamical control of magnetism.
In 1973, Hepp and Lieb showed that N → ∞ polar molecules located inside a resonant electromagnetic cavity undergo a second order transition from a normal to a superradiant phase. The Z 2 (parity) symmetry is spontaneously broken leading to a ferroelectric-like state in the "matter" and to a nonzero population of photons in the cavity at equilibrium [1][2][3]. However, 47 years and a global pandemic later, this quantum phase transition is yet to be measured [4]. During this time, the community has enjoyed a winding succession of proposals on how to achieve the superradiant phase transition (SPT), each shortly matched by its corresponding no-go theorem [5][6][7][8][9][10]. Nowadays, the non-observation of the phase transition is well understood. The crux of the matter resides in the approximations used to derive the Hamiltonian solved by Hepp and Lieb, the so-called Dicke model. On one side, there is the treatment of the A 2 -term (the diamagnetic term) [11]. On the other, matter truncations in different gauges must be done consistently as pointed out by Keeling, showing that the phase transition, if any, is completely attributable to matter interactions [12]. Consequently, the matter phase diagram remains unaltered despite it being immersed in a cavity and no photonic population develops. The same conclusion has been recently revisited [13]. A final step for closing the debate is found in the work by Andolina and collaborators [14]. In the thermodynamic limit N → ∞, their no-go theorem illustrates how Gauge invariance inherently prohibits the SPT for electric dipoles in the long-wavelength limit. The latter implies that matter cannot respond to a static and uniform electromagnetic field. Therefore, to make the SPT observable, either the nature of the coupling or the spatial field distribution must be, in some way, modified.
Theoretical proposals consider systems in the ultrastrong light-matter coupling regime [15] or use electron gases that either possess a Rashba spin-orbit coupling [16] or are subjected to a spatially-varying electromagnetic field [16][17][18] or have the role of photons be played by magnons [19]. Here, we propose an alternative set-up, based on magnetic molecules that couple to superconducting microwave resonators via the Zeeman interaction [20]. Ar- Figure 1. Schematic picture of a CPW resonator coupled to a spin ensemble. Emw and Bmw are the cavity's microwave electric and magnetic fields, respectively. Bext is the external magnetic field that induces a Zeeman splitting between the spin energy levels.
tificial magnetic molecules [21,22], designed and synthesized by chemical methods, consist of a high-spin cluster core surrounded and stabilized by a cloud of organic molecular ligands. The ability to chemically tune their relevant properties, such as the ground state spin, magnetic anisotropy and mutual interactions, combined with their stability as isolated molecular units, confer them a potential interest as magnetic memories in spintronic devices [23] and as qubits for scalable quantum information schemes [24][25][26]. Besides, they tend to organize forming crystals, which makes them model systems to study pure magnetic dipolar order and quantum phase transitions [27][28][29][30][31]. This work explores the realization of the Dicke model (and generalizations of it), which undergoes the equilibrium SPT, in a crystal of molecular nanomagnets coupled to a on-chip microwave cavity. We compute the critical condition that triggers the SPT at zero and finite temperatures under rather general conditions, e.g. for a spatially-varying magnetic field or in the presence of direct molecule-molecule interactions, and discuss a feasible method to detect it. Besides, we build the phase diagram for the purely cavity-driven ferromagnetism and study how the spin-photon coupling enhances the intrinsic ferromagnetic order of a crystal of Fe 8 molecular clusters [31].
Magnetic cavity QED.-Hybrid platforms coupling electron [32] and, particularly, spin [33,34] ensembles to superconducting resonators or cavities complement circuit QED. Here, the "spins" are superconducting qubits [35,36]. Different magnetic species have been studied in this context, including impurity spins in semiconductors [37][38][39], lanthanide ions [40,41] and magnetic molecules [42][43][44][45]. To observe the SPT, set ups hosting molecular crystals offer the crucial advantage of coupling a macroscopic number of identical and perfectly organized spins to a single cavity mode (cf Fig. 1). A vast majority of molecular nanomagnets are both neutral and exhibit a close to zero electric dipole. Their response to external stimuli is then accurately described by a simple "giant"-spin effective Hamiltonian H S , which includes the effects of magnetic anisotropy and the couplings to magnetic field B [20,21]. The latter enter 2m is the Bohr magneton and g e the Landé factor (= 2 for an electron spin). The diamagnetic response, which arises mainly from the molecular ligands surrounding the magnetic core, is much smaller that the paramagnetic one and can be safely neglected, especially at sufficiently low temperatures. These results provide the experimental background for the following discussion. Let us summarize the main steps to model the moleculescavity system in the single mode case. For the multimode case, see [46]. The electromagnetic field is quantized, yielding the cavity Hamiltonian H c = Ωa † a, with a (a † ) the photonic annihilation (creation) operators, [a, a † ] = 1. The resonance frequency Ω ranges typically between 1 and 10 GHz. The local quantized magnetic field generated by the superconducting currents can be written as B mw (r) = B rms (r)(a † + a), with B 2 rms (r) = 0|B 2 mw (r)|0 its zero-point fluctuations. The resulting gauge-invariant cavity-spin Hamiltonian can be written as follows [20,47] with N the number of spins, the ladder operators S ± = S x ± iS y and the coupling constants Here, r j is the position vector of the j-th spin. The phases θ j in Eq. (1) are defined through B rms,x (r j ) + iB rms,y (r j ) = |B rms (r j )|e iθj . If the molecules are S = 1/2 non-interacting spins, like in a sufficiently diluted free radical sample, H S = ωz 2 j σ z j and the Hamiltonian (1) matches exactly the Dicke model. Notice that this model does not suffer from the "A 2 -issue" because the coupling is of Zeeman kind, rather than minimal (electric). Besides, the truncation of the electronic degrees of freedom to a finite dimensional Hilbert space is not an approximation but follows from the fact that we are dealing with "real" spins obeying the angular momentum commutation relations. Both properties combined permit to avoid the no-go theorems for the SPT [14].
Exact results at N → ∞.-In order to study specific spin models and how their properties are affected by the coupling to light and vice versa, it is convenient to obtain an effective spin Hamiltonian where the light degrees of freedom have been traced out. Following Hepp and Lieb's original derivation [2], this effective Hamiltonian is defined by the following expression, exact in the N → ∞ limit, where H(α) = α|H|α and H is the total Hamiltonian given by Eq. (1). The resulting (See [46]) effective Hamiltonian is It is apparent that the light-matter coupling translates into an effective Ising-type ferromagnetic interaction among all spins that drives the quantum phase transition. This formulation is particularly handy for studying spin models that are well captured by a mean field approach, as we show below. The critical point can be obtained by noticing that Then, in the thermodynamic limit, system and cavity factorize. Generalizing the prescription in [14] to the finite temperature case (fully developed in [46]), the critical condition can be written in terms of the static response function R(T ) of the bare spin model, where with |ψ m , m the eigenstates and eigenenergies of the bare spin Hamiltonian H S and ∆ mn = m − n . In the case of uniform coupling, λ i ≡ λ, the static response function is proportional to the magnetic susceptibility R = ( λ) 2 χ ⊥ in the xy-plane, perpendicular to the external dc magnetic field (see Fig. 1). A similar condition is found in three-dimensional electronic systems when the spin degrees of freedom are considered [18]. Both approaches show a relation between light and matter observables given by indicating that the phase transition is marked by the onset of both a macroscopic population of photons in the cavity and a spontaneous magnetization in the spin system.
In the case of independent spins coupled to a field pointing along x, we obtain the critical condition of the Dicke model for a spin S which depends on the coupling only through the collective parameterλ 2 ≡ N −1 j λ 2 j , i.e. the root mean square of the spin-dependent couplings where V sample is the cavity volume occupied by spins and ρ = N/V sample is the density of spins. In the second equality, we have used Eq. (2) and the Virial theorem and introduced the filling factor ν = I(V sample )/I(V total ) with I(V ) = V |B rms (r)| 2 dV . Vacuum-fluctuations-driven ferromagnetism.-To illustrate the effects of the coupling between the cavity and the magnetic molecules as well as the interplay between intrinsic and light-induced ferromagnetic interactions, we showcase two examples using realistic experimental parameters. From a theoretical perspective, Fig.  2 shows the alterations to the zero-temperature phases of the paradigmatic quantum Ising chain for S = 1/2 spins, defined by . Note that according to Eq. (7) ferromagnetic and superradiant phases are synonymous, i.e. ferromagnetic ordering is always accompanied by a finite photon population in the cavity. In particular, regions I, II and III are superradiant, as evidenced by the finite values of magnetization | σ x | and photon number a † a /N shown in the insets of Fig. 2. Remarkably, the light-induced effective ferromagnetic interaction overpowers the intrinsic antiferromagnetic interaction extending region I into the J < 0 sector. In Region II it is the synergy between in- trinsic and induced ferromagnetism that gives rise to superradiance. Finally, region III is intrinsically ferromagnetic even in the absence of the cavity and thus becomes superradiant when coupled to one. The use of the mean field approximation is validated in [46] by comparing the superradiant phase boundary obtained with mean field against the one obtained using exact diagonalization to compute the response function R (5).
Next, we consider a more realistic model. It corresponds to a specific molecular material, a crystal of Fe 8 clusters with S = 10, which shows a ferromagnetic phase transition purely induced by dipolar interactions below a critical temperature T c (B ⊥ = 0) ∼ = 0.6 K and a zerotemperature critical magnetic field B ⊥,c 2.65 T [31].
The magnetic phase diagram measured for a magnetic field perpendicular to the magnetic anisotropy axis is shown in Fig. 3. As expected for a system dominated by long-range dipolar interactions, it agrees very well with the predictions of a mean-field Hamiltonian with parameters given in Fig. 3. The combination of high spin, thus high susceptibility, and negligible exchange in-  (10) and is verified against experimental data (dots) obtained for a crystal of Fe8 molecular clusters [31], whose structure [48] is shown in the inset. The external field, B, is applied perpendicular to the easy magnetization axis x: The magnetic anisotropy parameters are D/kB = 0.294 K and E/kB = 0.046 K, and the coupling J/kB = 2.85×10 −3 K. The remaining lines are obtained with the same method taking into account the coupling to a microwave cavity for a range of filling factors ν. The parameters used are ρ = 5.1 · 10 20 cm −3 , taken from the crystal lattice of Fe8 and Ω = 1.4 × 10 9 s −1 .
teractions, which lead to a quite low T c even for a densely concentrated spin lattice, makes this material well suited to obtain and measure the SPT.
This expectation is borne out by calculations that include the coupling to a microwave cavity, whose results are also included in Fig. 3. They show that light-matter interaction has a remarkable effect on the equilibrium phase diagram, enhancing both B ⊥,c and T c , the latter by almost as much as a factor six, depending on the filling of the cavity. This can be understood by noting that the effective Hamiltonian is given by Eq. (10) but with the coupling J replaced by J eff = J + λ2 /Ω. This enhancement evidences that the cavity induces quite strong ferromagnetic correlations, a characteristic signature of the SPT, which in this case co-operate with the intrinsic interactions between the spins in the crystal. Achieving filling factors well above 0.1 seems quite reasonable provided that one achieves a sufficiently good interface between the chip and the magnetic material. It can also be seen from these results that, even after the introduction of direct spin-spin interactions, the transition depends only on the filling factor, as shown in [46]. We emphasize that the two previous examples show how, with magnetic coupling, the light-induced and intrinsic interactions add up, modifying the bare matter phase diagram. This is the signature of photon condensation and, thus, of the occurrence of the SPT [12,14].
Transmission experiment for resolving the transition.-Finally, we discuss how to measure a signature of the phase transition. A direct route would be to measure the order temperature of the magnetic material inside and outside the cavity. However, conventional methods to measure T c (or B c ), based on magnetic susceptibility or neutron diffraction [31], do not lend themselves easily to include a superconducting cavity. Besides, they require very large crystals, often much larger than the typical cavity volumes. Therefore, we envision here a more accessible way: a transmission experiment, where the cavity is coupled to a microwave transmission line. A signal, sent through it, interacts with the cavity-spins system and the transmitted signal is recorded. Since this signal is proportional to the dynamical response of the system, we expect to observe a signature near the transition. Technically, the calculation involves computing the dynamical susceptibility of the whole system (cavity plus spins) using a quantum master equation. Typically, the dissipation for both the spins and the cavity is added phenomenologically, inserting the terms as if the cavity and the spins were not coupled. However, this, so-called, local approach has been criticized. A rigorous derivation involves taking into account the coupling between the two subsystems, obtaining the global master equation. The differences between both approaches are relevant in strongly coupled systems [49][50][51]. Since we are interested in what happens near the transition, we explore the differences between local and global approaches to rule out that the signature is an artifact of the approximations taken. We notice, however, that in the global approach both eigenvectors and eigenvalues for the coupled system are needed. These are impossible to obtain in a full treatment. On the other hand, if we consider a spin 1/2 Dicke model with homogeneous coupling to the cavity a Holstein-Primakoff transformation allows us to write the total Hamiltonian as two coupled oscillators. Then, the system is exactly solvable and the global master equation can be obtained. By doing so, we can compare both local and global approaches for resolving the transition.
The explicit formulas are rather involved and thus sent to [46]. In figure 4 we summarize our findings. We show 2D transmission plots for temperatures larger or smaller than the zero-field critical temperature T c (0). In the former limit, the system remains disordered, regardless of ω z and only shows the well-known avoided crossing at ω = Ω = ω z . As temperature is decreased below T c (0) the magnetic and photon states depend on ω z . The phase transition, for ω z ∼ ω z,c (0), is then characterized by the appearance of a new resonance, i.e. a new transmission channel, at ω z ∼ ω z,c (0). On physical grounds, this is a consequence of the vanishing frequency of the lower mode at the transition, which increases the response at equilibrium (ω → 0). Besides, the resonance appears in both the local and global approaches, supporting our findings. Quantitative differences appear, obviously, but this is expected since one of the pitfalls of the local approach is not reproducing the correct equilibrium states. This behaviour constitutes a clear signature of the SPT and shows that this transition can be detected in a standard temperature-dependent transmission experiment, provided that the spin-photon coupling is large enough.
Conclusions.-We have shown that the coupling of a macroscopic number of spins, in particular crystals of molecular nanomagnets, to the quantum vacuum fluctuations of cavities or LC resonators generates ferromagnetic spin-spin interactions that lead to an equilibrium superradiant phase transition. Using realistic parameters, we find that it also gives rise to detectable signatures in the transmission of microwaves through such a hybrid set-up, thus providing a feasible solution to the long-standing problem of measuring the SPT. Our results also present these systems as ideal for exploring the quantum electrodynamical control of matter, baptized as cavity QED materials [52]. Recent studies have shown that quantum light fluctuations can modify properties such as excitonic transport [53][54][55], chemical reactivity [56,57], superconductivity [58][59][60][61] and the ferroelectric phase in quantum paramagnetic materials [62][63][64]. Here, we show that they can also generate, modify and control long-range ordered magnetic phases.
In this Supplemental Material we provide additional information on the modelling of the interaction between light and molecules (Sec. I), the effective spin Hamiltonian (Sec. II), the critical condition in terms of the spin static response function (Sec. III), the finite temperature stiffness theorem (Sec. IV), the calculation of the modified Ising phase diagram (Sec. V), the dependence of the effective spin-spin interaction onλ (Sec. VI) and the transmission experiment (Sec. VII).

I. MODELING THE INTERACTION OF LIGHT AND MAGNETIC MOLECULES
If we consider the source of the electromagnetic field to be an optical or superconducting cavity, we quantize the electromagnetic field in the Coulomb gauge, where A(r i ) = k A k (r i )(a k + a † k ), B(r i ) = ∇ × A(r i ) and a k , a † k are the bosonic annihilation and creation operators. The interaction of the spins with the magnetic field through Zeeman coupling H Z = −g e µ B S · B can be split into two contributions, one from the external field inducing level splitting on the spin energy levels which is included in H S , and another from the microwave field generated by the cavity which constitutes H I . The local multimode quantized magnetic field generated by the superconducting currents can be written as B mw (r) = k B k,rms (r)(a † k + a k ), where we have used the orthogonality of the electromagnetic modes to write B 2 rms (r) = 0|B 2 mw (r)|0 = k B 2 k,rms (r). The gauge invariant full cavity-spin Hamiltonian is in general where

II. COMPUTING THE EFFECTIVE SPIN HAMILTONIAN
Despite the fact that the bounds for the partition function were essentially proven by Hepp and Lieb [1,2], let us prove it again here for self-consistency. We want to show that where H(α) = α|H|α . Consider a Hamiltonian of the form In essence, the only relevant details are that the light-matter coupling is linear in the bosonic operators and that the only non linear contribution is trough the photonic term Ωa † a. We will start with the lower bound. For that, we make use of Jensen's inequality for convex functions, which states that, for a convex function f Because the exponential is a convex function, it follows that The upper bound requires some more algebra. We will make use of the Golden-Thompson inequality: let A and B be Hermitian operators, then We are concerned with obtaining an upper bound for Z, which we can express as for m → ∞. Before advancing further, it is convenient to express a † a = aa † − 1 in terms of a sum of projectors of coherent states, using the closure relation π −1 d 2 α |α α|, a † a = 1 π d 2 α a |α α| a † − |α α| = 1 π d 2 α |α| 2 − 1 |α α| .
We can now write Defining the multiplication operator A we can write the previous expression as At this point we use the Golden-Thompson inequality, noticing that B n = I which finally gives us the upper bound It can be readily seen that the generalization to a multimode field simply yields Z ≤ exp[β k Ω k ]Z. Provided the number of modes is finite or alternatively that k Ω k is, the bounds for Z become exact equalities in the thermodynamic limit.
At this point, we define the effective spin Hamiltonian through the following relation Where, for a general spin Hamiltonian, we define H(α) as Which yields We have used the fact that [H S /N, j λj √ N e iθj S + j + h.c. /N ] → 0 in the thermodynamic limit, provided that H S is extensive, in order to extract the exponential of H S for integration and then to put together H S,eff .

III. COMPUTING THE CRITICAL CONDITION IN TERMS OF THE SPIN STATIC RESPONSE FUNCTION
Let us decompose the Hamiltonian as H = H ph + H S + H int . Where H S is model dependent and We found that a generalization of Andolina and colleagues' approach [3] is convenient for dealing with general Hamiltonians of this sort. Let us first consider the zero temperature case, and then advance to finite temperature.
A. Zero temperature T = 0 In the thermodynamic limit, N → ∞, it can be shown that [H ph /N, H int /N ] = [H S /N, H int /N ] = 0, thus light and spin degrees of freedom are effectively decoupled and the states of the system are separable. In such scenario a mean-field treatment becomes exact, and we can propose for the photons in the cavity a coherent state a |α = α |α . Taking α ∈ R and a spin state |ψ , we find Minimizing with respect to α yields which gives a constraint in the search for the ground state of H S . Introducing Eq. (22) in Eq. (21) gives the minimum energy for a state with non-zero α For the state |Ψ = |ψ ⊗ |α with α = 0 to be the true minimum requires |ψ to be the least energetic state obeying the constraint (22) (termed constrained ground state) and, of course, Where |ψ 0 is the unconstrained ground state of H S . The stiffness theorem [4] allows us to expand ψ|H S |ψ − ψ 0 |H S |ψ 0 in powers of α.
withR being the static response functioñ Note that the static response function used in the main text R is defined as R = ( Ω) 2R , whereR is the response function defined here above.
Firstly, let us show that minimizing the energy is either sufficiently or completely equivalent to minimizing the free energy at low temperatures. The free energy is defined as F = E − T S. Consider now a variation in F at constant temperature Clearly, at T = 0, ∆F = ∆E, so minimizing the free energy is equivalent to minimizing the energy at T = 0. For higher temperatures, we can write Taking the low temperature limit we have Applying L'Hopital's rule and Nernst postulate lim T →0 ∆S = 0, we find Therefore, not only are ∆E and ∆F equal at T = 0 but so are their slopes, ensuring that ∆E and ∆F remain equal in a range of low temperatures, and not only at strictly T = 0. This implies that the minimization of F and of E become equivalent at low temperatures. Now, if we consider that our thermal spin ensemble state is captured by the density matrix ρ, the energy of the system is which we can minimize with respect to α to yield a constrained minimum condition for the spin system Reintroducing it in our expression for the energy gives With that, the critical condition becomes E ρ (α) ≤ E ρ0 (0), which amounts to We can now make use of the extended stiffness theorem (See Sec. IV) to write the left hand side of Eq. (35) in terms of powers of α and then obtain a critical condition analogous to the zero temperature case, but with a static response function that depends on temperature Considering the temperature dependent static response function we obtain the critical condition Where |ψ m and m are respectively the eigenstates and eigenenergies of H S and ∆ mn = m − n .

IV. EXTENDING THE STIFFNESS THEOREM TO FINITE TEMPERATURES.
Let us extend the Stiffness theorem [4] to deal with finite temperatures. We consider a constrained minimization problem, with the energy given by where ρ A indicates that A = Â = T r(Âρ A ). We now construct the HamiltonianĤ A =Ĥ + F AÂ such that By the definition of the static response function, we impose . (40) With that, we can express which we now wish to compute. For that purpose, consider a parametrized Hamiltonian with parameter 0 ≤ λ ≤ 1, Extending the equations presented thus far, we find with Applying the Hellman-Feynman theorem again we find Substituting back into Eq. (41) we finally obtain So the result is analogous to the zero temperature case with the only difference that the static response function R AA (T ) is now a function of T . We have assumed that Â 0 = 0 throughout

V. COMPARISON BETWEEN EXACT DIAGONALIZATION AND THE MEAN FIELD APPROXIMATION FOR COMPUTING THE SUPERRADIANT PHASE BOUNDARY OF AN ISING SPIN MODEL
In Figure 2 of the main text, we study the different phases that a hybrid spins-cavity system can host when a direct Ising nearest-neighbours interaction between spins 1/2 in a 1D model is considered. We also assume for simplicity that the coupling of the spins to the cavity is uniform, λ i ≡ λ.
Two approaches are then possible to compute the phase boundary between the superradiant and normal phases. From a practical viewpoint, it is convenient to perform the mean field approximation on the direct Ising interaction part of the Hamiltonian. Since the cavity-mediated Ising interaction is all-to-all, it can be exactly treated with mean field. Proceeding in this fashion from Eq. (4) in the main text (with = 1), we obtain Which is the mean field effective Hamiltonian of a Dicke model with a modified light-matter coupling From Eq. 47 one obtains an analytical critical condition via free energy minimization with respect to the variational parameter m x An alternative calculation simply resorts to the critical condition given by Eq. (5) in the main text, which requires the computation of the response function R. In the case of homogeneous coupling, the latter is proportional to the magnetic susceptibility in the xy-plane, χ ⊥ . This can be computed numerically by introducing a perturbative field along the x-direction in the Ising Hamiltonian The susceptibility is then defined as σ x = σ x 0 + χ ⊥ . With a Lanczos routine, we perform exact diagonalization of the Hamiltonian to compute σ x and σ x 0 . In Fig. 1 we compare the phase boundaries obtained with the two methods here described. We find that the mean field approximation provides an excelent qualitative description of the sinergy between intrinsic and cavity-mediated Ising interactions.

VI. EXPRESSING THE EFFECTIVE SPIN-SPIN INTERACTION IN TERMS OF THE ROOT MEAN SQUARE COUPLINGλ
Consider the cavity-mediated effective coupling that appears in the effective spin Hamiltonian H S , Without loss of generality, let us focus on the case that θ j = 0 and apply a mean field approximation to the spins The Cauchy-Schwarz inequality teaches us that the root mean square can be used as an upper bound for the arithmetic Considering that the λ j are definite positive, this upper bound is reasonably close to the actual value, allowing us to approximate H MF eff,int by Benefiting again from the bound, let us perform a mean field approximation on the λ j where we neglect fluctuations around the root mean square λ j =λ + δ(λ j ) λ . With this approximation, we can finally write Having shown thatλ 2 is a function of only the density of spins ρ and the filling factor ν, this approximation allows us to study the modification of magnetic properties depending exclusively on ρ and ν, as is exemplified in the main text.

VII. COMPUTING THE TRANSMISSION IN A DRIVEN DICKE MODEL
A. Phenomenological damping -Local master equation For simplicity, consider a spin 1/2 Dicke model coupled to a transmission line ( = 1), in which the coupling is modulated by a periodic driving of the cavity where f (t) = γ 1 e −iωt . The quantity to be measured in the experiment is the transmission, which is proportional to the susceptibility of the system to the periodic driving that is the excitation. The expression for the transmission is t = −iγ 1 χ a . So we are interested in computing χ a . For that, it is convenient to express the Hamiltonian in terms of total spin operators so we can write the equations of motion of the system in terms of these new observables. Presumably, the experiment would be carried out in a regime where the temperature can be approximated by zero, which simplifies the calculation. The equations of motion are the Heisenberg equations of motion (Ȯ = i[Ĥ,Ô]) plus a phenomenological damping originating from the coupling to the transmission line Here γ 1 and γ 2 are, respectively, the damping constants for the cavity and the spins. These can be formally obtained by considering each subsystem (spins and cavity) in isolation, and then postulating for each a more general model where the subsystem is coupled to a bosonic environment that acts as a heat reservoir. The dynamics of each each subsystem are then given by their so-called local master equations, in which the bath is traced-out, leaving the dissipation terms that we introduce here. The denomination local stems from the fact that the dissipation terms are computed for each subsystem independently but then introduced into the equations of motion of the joint system. Below we contrast this approach with the more consistent global master equation, in which a single master equation is written for the joint system. The dissipation terms are computed taking into account the complete system, instead of the isolated subsystems. In many instances, the phenomenological damping given by the local approach provides an accurate description of the dissipative dynamics. However, in cases where the intersubsystem coupling is strong, the significant modification to the energy levels warrants the use of the global master equation [5]. At the very least, to validate the results obtained with the simpler local master equation. If we consider f (t) as a perturbation, i.e. γ 1 Ω, we can write the average values of the operators in terms of their unperturbed values and the corresponding susceptibility: S µ = S µ 0 + γ 1 e −iωt χ µ , a = α 0 + γ 1 e −iωt χ a . Imposing the equilibrium condition:Ṡ µ =ȧ = 0, and assume a = a † to be real, we arrive to the system of equations that must be solved in order to obtain χ a       Solving for χ a yields Noting that S µ 0 = N σ µ 0 /2 allows us to write with The dependece on the normal and superradiant phases is through the equilibrium values of σ z 0 and σ x 0 .

B. Global master equation
Introducing the total spin operators S z = j σ z j 2 , the Hamiltonian of the spins-cavity system is We use a Holstein-Primakoff transformation on the spins to switch to a bosonic representation, where the system consists of two coupled oscillators. When applying the transformation, we must distinguish two cases, when the system is in the normal phase and when it is superradiant. In the normal phase at zero temperature, where the ground state is an empty cavity with relaxed spins, we use the relations In the thermodynamic limit, N → ∞, these relations are exact. The new operators s † , s obey the bosonic commutation relations and hence constitute bosonic creation and annihilation operators, respectively. The resulting Hamiltonian describes two coupled harmonic oscillators. In the theory of open quantum systems, dissipation is introduced by expanding the system to include bosonic baths that couple to each of the system's degrees of freedom. In this case, to each of the oscillators. Henceforth, we follow the prescription of Breuer and Petruccione [6]. We assume the coupling operators to be a + a † and s + s † for each subsystem. The master equation formalism requires that we find the spectral decomposition of the coupling operators in terms of the eigenoperators of the Hamiltonian. We diagonalize the Hamiltonian with the use of the Bogoliubov transformation. It can be written in matrix form as H = 1 2 φ † M φ, with φ t = (a, s, a † , s † ) and To compute the susceptibility χ α , we recall the Bogoliubov decomposition of operator a (67), from which it is clear that χ α = g 1 χ 1 + g 2 χ 2 + g 3 χ † 1 + g 4 χ † 2 . Note that the susceptibilities are scalars and thus the † sign is not used indicate hermitian conjugacy, but simply to mark that χ † 1 (χ † 2 ) is the susceptibility of operator d † 1 (d † 2 ). Note as well that χ † 1 = χ * 1 . We are thus interested in the dynamics of the transformed operators. Focusing on d 1 , we find that the dynamics are completely uncoupled from d 2 , with master equatioṅ As for the terms that survive the secular approximation at the transition, we have They do not contribute to the dynamics. Assuming a flat bath γ j (ω) = 2γ j N (ω), with N (ω) the bosonic ocupation number at temperature T of the bosonic bath, we obtain a simplification for Eq. (73) valid at all temperatures, and particularly at T = 0,ḋ Analogusly, we obtain for d 2ḋ In order to compute χ α , we introduced a periodic perturbation in H of the form γ 1 e −iωt (a + a † ) in the local master equation formalism. Having switched to new bosonic operators, this translates to We define the susceptibilities as the response of each operator to this perturbation, for instance d 1 = d 1 0 + γ 1 e −iωt χ 1 , with d 1 0 the equilibrium value. Introducing this and the perturbed Hamiltonian (77) in Eq. 75 we obtain Proceeding analogously, we obtain the remaining susceptibilities Finally, the transmission is t = −iγ 1 χ α , or, written in terms of the new susceptibilities, t = −iγ 1 g 1 χ 1 + g 2 χ 2 + g 3 χ † 1 + g 4 χ † 2 .
It is convenient to recall the range of validity of the expression for the transmission just presented. The use of the Holstein-Primakoff transformation implies that we stay on the normal phase at zero temperature, where the equilibrium state consists entirely of relaxed spins. We can, however, extend it to finite temperatures with the use of a clever observation. It can be seen from the critical condition 4λ 2 = Ωω z coth (βω z /2) that we can define an effective temperature-dependent coupling λ 2 (T ) = λ 2 tanh (βω z /2). Realizing that | σ z 0 | = tanh (βω z /2), we see that as the In this case the result is not as clear as with the local approach as it depends on the coefficients of the Bogoliubov transformation. To gain further insight, in Figure 2 we plot the transmission in the low ω regime as a function of ω z . Note that this is at T = 0, where a phase transition is expected. What we find is that (g 1 − g 3 ), and subsequently the global transmission, goes to zero right at the transition, despite resonant behaviour right before and after the critical point. This contrasts with the result obtained for the local transmission, but in any case, both approaches show a signature either at (local) or immediately around (global) the phase transition.