(Almost) Everything is a Dicke model – Mapping correlated light-matter systems to the exactly solvable Dicke model

We investigate classes of interacting quantum spin systems in a single-mode cavity with a Dicke coupling, as a paradigmatic example of strongly correlated light-matter systems. Coming from the limit of weak light-matter couplings and large number of matter entities, we map the relevant low-energy sector of a broad class of models onto the exactly solvable Dicke model. We apply the outcomes to the Dicke-Ising model as a paradigmatic example [ 1,2 ] , in agreement with results obtained by mean-field theory [ 2 ] . We further accompany and verify our findings with finite-size calculations, using exact diagonalization and the series expansion method pcst ++ [ 3 ]


Introduction
Over the last decades, advancements in the field of cavity quantum electrodynamics as well as circuit quantum electrodynamics paved the way to study systems of matter strongly and collectively coupled to a light mode.These experimental breakthroughs made it possible to realize and study paradigmatic theoretical models like the Rabi, Tavis-Cummings, and Dicke model with strong light-matter interactions in the lab [4][5][6][7][8][9][10][11].With these tools at hand, a fundamental question is how the interactions between light and matter do influence each other, altering properties of the separated (potentially complicated) individual parts, like observables, local interactions, or the location of phase transitions [12][13][14][15][16][17][18][19][20][21][22].
One of the paradigmatic light-matter systems is the Dicke model, with a minimal setup both on the light and the matter part [23,24].The model consists out of N individual spin-1/2 particles that are individually coupled to a single cavity mode.Hepp and Lieb showed for the thermodynamic limit N → ∞ that this system can be solved analytically by a Bogoliubov transformation and features a second-order phase transition from a normal to a superradiant phase with a non-vanishing photon density in the ground state [24].While the matter part of the Dicke model consists out of an arbitrary number of spins, it breaks down to a noninteracting problem in the case of no light-matter interaction, as the local degrees of freedom are only coupled through the cavity, making it trivially solvable.
In order to make the composite system more interesting, various generalizations for the Dicke model were proposed and discussed, like more complex local spin structures [25], multimode cavities [24,26,27], non-Hermitian generalizations [28], open systems [29,30], altered light-matter interactions [31,32], non-equilibrium systems [33], and added matter-matter interactions between the spins [2,34,35].One prime example is the Dicke-Ising model, where an additional Ising interaction between nearest-neighbor spins is present.Using mean-field theory and a classical approximation of the spin degrees of freedom, Zhang et al. found a rich phase diagram for antiferromagnetic interactions including superradiant phases, with both an antiferromagnetic and a paramagnetic phase in the matter part [2].However, using quantitative numerical techniques, deviations were found for the phase transitions both in the location as well as in the order in 1D [1,36].
In this work we elaborate on this line of thinking by considering a more generalized setting for the matter part, consisting of long-range hopping and correlated processes, and coupling it to a single light mode.This enables us to investigate the interplay between the correlations and effects induced by the light-matter coupling and the matter-matter coupling.Restricting ourselves to the non-superradiant phases that are adiabatically connected to the case of vanishing light-matter interaction, we establish an analytical solution of the low-energy part of this model by mapping it to an effective Dicke model.This enables us to study the lowlying excitations of this generalized Dicke model in the non-superradiant phases analytically, including the closing of the gap, potentially inducing second-order phase transitions.
The structure of this paper is as following.First, in Sec. 2 we introduce the general framework, including the generalized model and the prerequisites in order to derive the effective Dicke model.The latter is done in the Subsecs.2.2 and 2.3, first giving some physical intuition how the system can be solved, followed by a general derivation on operator level.In Sec. 3 we apply our general findings onto the Dicke-Ising model as an exemplary case.We compare our outcomes, obtained in the thermodynamic limit, with results from exact diagonalization (ED) and the series expansion method pcst ++ [3] on finite systems to strengthen the validity of the effective model.In Sec. 4 we draw conclusions and give an outlook for potential research directions.

Derivation of the effective Dicke model
In this section we will derive the main finding of this work on a general level.First, we will define the framework including the light-matter model and further conditions we assume.Thereafter, we present our two approaches to show the decoupling of the Dicke part and the rest of the Hamiltonian.While the second one is more rigorous and general, the first one helps understanding the decoupling in a more intuitive way being less technical.

General framework
For this work we concentrate on a general Hamiltonian H of the form where H D denotes the Dicke Hamiltonian [24] and H matter represents an Hamiltonian with additional matter-matter interactions.For the Dicke Hamiltonian we use the following notation where the first term describes the coupling of the spins to a magnetic field ω 0 /2 > 0 in z direction.The second term defines the energy of the single light mode with frequency ω > 0.
The last term induces a coupling between light and matter degrees of freedom with a coupling constant g/ N ≥ 0, with N being the number of spins in the system.We restrict ourselves to the non-superradiant phase with ω 0 , ω ≫ g inducing two quasiparticle types, namely photons and magnons.The photons are given as excitations of the light field, the magnons are defined as spin-flips out of the fully polarized state induced by the magnetic field.We can rewrite the Dicke Hamiltonian with hardcore bosonic creation and annihilation operators b † j , b j via the Matsubara-Matsuda transformation [37] as with E 0 = −ω 0 N /2 being the ground-state energy of H D for g = 0 and n j = b † j b j the number operator.
In addition, we define a pure matter Hamiltonian of the form similarly to generalizations done before [38][39][40][41].The first term describes (long-range) hopping processes of magnonic excitations, while the second one describes (long-range) interactions including density-density terms and correlated hopping terms.The sums over δ i describe all possible distances on a chosen lattice.To keep the sums over all distances finite, we restrict the coefficients to stay finite in the limit of N → ∞ as For our findings we further restrict ourselves onto the low-energy subspace demanding a finite expectation value for the number of magnonic and photonic excitations in the thermodynamic limit as Figure 1: Visualization of the general light-matter Hamiltonian H from Eq. ( 1) on a square geometry.The model consists out of a large number of spin-1/2 particles (white spheres), with a local energy splitting caused by a uniform magnetic field, which interact with each other via (long-range) hopping and interaction terms, as introduced in Eq. ( 4), given in purple.The interacting matter system is put into a cavity with a single light mode of frequency ω.The spins and the light mode interact via a Dicke coupling proportional to g.
In Fig. 1 we sketch the constituents of the model on a square lattice geometry, while we will not restrict ourselves on any particular lattice structure in the following.For parts of the discussion it is convenient to express the Hamiltonian in momentum basis via Fourier transformation.We write the operators in momentum space as bk = 1 with k being the momentum of the annihilated or created magnon.The Hamiltonian can thus be written as with δ k 1 +k 2 ,p 1 +p 2 being the Kronecker delta preserving momentum conservation.See App.A for an in-depth derivation of Eq. ( 8).
For the case of vanishing c δ 1 ,δ 2 ,δ 3 , we obtain again the Dicke model with a rescaled magnetic field, depending on the momentum for c δ ̸ = 0.In the thermodynamic limit we can show that the momentum operators fulfill the bosonic commutation relation for the low-energy spectrum (see App. A).Thus, the Dicke model can be solved analytically by performing a Bogoliubov transformation in the k = 0 subspace, as described in [42].All other k ̸ = 0 sectors are already in diagonal form.Therefore, we can read off the energies directly.
The general case c δ 1 ,δ 2 ,δ 3 ̸ = 0 is not block diagonal in momentum space, as the mattermatter interactions induce scattering of magnon pairs, coupling different single-magnon momenta.For the further discussion, we will define where HD is the rescaled Dicke Hamiltonian from the first line of Eq. ( 8) and H MM all matter processes including k ̸ = 0 magnons.In the following, we will motivate and show that HD and H MM decouple in the thermodynamic limit, given the prerequisites in Eqs. ( 5) and (6) to stay in the non-superradiant phase, so that we can solve the two parts of the Hamiltonian independently from each other.

Physical intuition
This subsection is not meant as a rigorous proof of our findings, but shall rather motivate on a physical level why the proof in Subsec.2.3 works.Apart from the correlated processes proportional to c δ 1 ,δ 2 ,δ 3 , the general Hamiltonian H in Eq. ( 1) is solely made up of uncorrelated one-magnon terms in the matter part, including number operators and hopping processes.This means there is no interaction between magnonic excitations, apart from the hardcore constraint coming from the spin-1/2 sites.Therefore it is beneficial to transform the Hamiltonian into momentum space, as done in Eq. ( 8), obtaining effectively free magnon modes with different momenta.In this basis, the light-matter interaction only couples to the k = 0 sector.As the b( †) k operators obey bosonic commutation relations in the thermodynamic limit for the low-energy spectrum (see App. A), it is intuitive that only the k = 0 part of the effective hopping in H matter contributes to the Dicke dynamics, as the other magnon modes with momenta k ̸ = 0 do not couple to the light mode.Thus, the eigenstates of HD in Eq. ( 9) are hybridized excitations of the k = 0 magnon mode and the photonic cavity mode.
Moving on to the interaction processes H MM in Eq. (10), the terms differ fundamentally to all other terms in H as they couple two magnons with each other.This results in effective scattering processes of the form coupling magnons with different momenta with each other.Only considering these correlated processes, we obtain eigenstates with (multiple) magnon pairs with a fixed distance to each other depending on the coefficients c δ 1 ,δ 2 ,δ 3 , in contrast to the free magnon modes from the rest of the Hamiltonian.This leads to a localization of the states in real space, directly corresponding to a broadening in momentum space.Thus, the weight of the states on a single frequency is reduced, leading to a vanishing weight in the thermodynamic limit, as can be shown.As only the k = 0 mode couples to the light field, it is reasonable that Eq. ( 11) will decouple from the light part in H D for N → ∞ in the low-energy subspace.So, the low-energy sector of H features three kinds of dynamics in the matter part in the thermodynamic limit, which we will prove in subsection 2.3.The corresponding low-energy spectrum is sketched in Fig. 2 1) for small light-matter couplings.The parameters of the exemplary model are chosen appropriately to show the typical features of this class of systems.The energy of the perturbed photon states are omitted for clarity.The one-magnon subspace consists out of a conventional free mode (green line) with an energy jump at k = 0, where the magnon mode couples to the light mode (blue dot).The two-magnon subspace consists out of a conventional continuum (green area) complemented by composite two-magnon states with one hybridized magnon at k = 0 and a magnon with k ̸ = 0 that does not couple to the light mode (blue area).For suitable H MM one or more bound states can form (orange line) that are not affected by the light-matter coupling.
1. a free magnon state with momentum k = 0 forming a hybridized state with photons from the cavity due to the light-matter coupling term in HD (blue dot), 2. free magnon states with non-vanishing momentum in H MM that do not couple to the light mode (green line), 3. bound magnon-pair states with a fixed distance in H MM that do not couple to the light mode (orange line).
While the first one commutes with the latter two, the second and third part may only be described together, as it is common for correlated matter Hamiltonians.The first part can be solved analogously to the Dicke model by applying a Bogoliubov transformation because this part does not involve any matter-matter interactions (see App. C).As H MM conserves the number of magnons, this part of the model is block diagonal with respect to the number of particles and can be solved, e.g., with exact diagonalization.

General case
To show on a general level that HD and H MM from Eqs. ( 9) and ( 10) decouple, we prove in this concluding derivation section that the commutator [ HD , H MM ] vanishes in the thermodynamic limit for the given prerequisites.We therefore perform calculations on finite systems with N sites and afterwards take the limit N → ∞.Contrary to intuition, the terms HD and H MM with different distinct momenta p ̸ = k do not commute for finite systems as [ bp , b † k ] ̸ = 0 (see App.A for details).Thus, we have to include these terms in the calculation to rigorously show that the commutator vanishes in the thermodynamic limit.In the following we split the calculation of the commutator into the two terms of H MM , starting with the free-particle term.First, we find We can express both commutators in Eq. ( 12) in terms of the commutators using the product identity scale with N −1/2 and thus are zero in the thermodynamic limit.As the sum over all |c δ | has a finite value and 〈a † + a〉 , 〈 bk 〉 , 〈 b † k 〉 < ∞ for the low-energy subspace, we conclude that the commutator of Eq. ( 12) vanishes in the thermodynamic limit.
For the correlated hopping term in H MM we build the commutator analogously -using a mixed representation of momentum and real-space operators to keep notation short -as Again, we reexpress the first commutator of Eq. ( 15) by the second one by using the product identity.In App.B we show that the norm of the commutators scale with N −1/2 .With Eq. ( 5) and 〈a † + a〉 , 〈 bk 〉 , 〈 b † k 〉 < ∞ we conclude analogous to above that the commutator of Eq. ( 15) approaches zero in the thermodynamic limit for the lowenergy spectrum.
Thus, the commutator [ HD , H MM ] vanishes in the thermodynamic limit.Therefore, we can discuss the rescaled Dicke model HD and the correlated processes from H MM , separately.This induces that the potential bound states evolving out of H MM do not couple to the light mode in the thermodynamic limit, agreeing on the physical intuition gained in Subsec.2.2.

Application onto the Dicke-Ising model
After motivating and proving the main finding of this work that HD and H MM commute in the thermodynamic limit, we will now apply it onto a concrete model, namely the Dicke-Ising model.First, we introduce the model and give an overview of the research done on it in the last years.Afterwards, we bring it into the form of Eq. ( 1), followed by a discussion of the two non-superradiant phases of the model, which obey the prerequisites in Eq. ( 6), namely the paramagnetic normal phase and the antiferromagnetic normal phase.In the end, we will compare our results, which were derived in the thermodynamic limit, with calculations on finite systems to check for the convergence towards the effective model for larger system sizes.).The four phases PNP, ANP, ASP, and PSP are given in different colors.The phase transitions obtained by mean-field theory in [2] are plotted in solid lines.For the PNP and ANP we discuss effective models for the low-energy spectrum becoming exact in the thermodynamic limit.The phase transitions out of the two normal phases coincide with the closing of lowest gap in the presented effective theory, denoted with black solid lines.For J = 0, we recover the standard Dicke model, indicated with a faint dashed line.

Introduction to the model
The Dicke-Ising model is defined as a combination of the Dicke model with an additional Ising interaction in the direction of the external magnetic field.It can be written as [2] with ferromagnetic (antiferromagnetic) Ising interactions for J < 0 (J > 0) and N being the number of sites.Due to the competition of the Ising interactions with the magnon-photon coupling, the model features a variety of phases.In [2] the phase diagram for antiferromagnetic Ising interactions was investigated with mean-field theory and a classical approximation of the matter part.This yields four distinct phases, as plotted in Fig. 3.According to this mean-field treatment, all quantum phase transitions between the four phases are of second order.In contrast, recent works using perturbation theory and exact diagonalization [1] showed that the phase transition at the limiting case ω 0 = 0 is of first order in 1D.A study on the full Dicke-Ising model using quantum Monte Carlo [36] has found that the first order transition survives around the limit ω 0 = 0 in 1D and 2D.For larger ω 0 the gap shrinks at the phase-transition point until changing to a continuous phase transition, as predicted by the mean-field study.We refer to [36] for a quantitative discussion of the nature of the phase transitions.
The different phases can be characterized by the photon density 〈a † a〉 /N and the xmagnetization of the spins.For small g, the 'paramagnetic normal phase' (PNP) and the 'antiferromagnetic normal phase' (ANP) have zero photon density.For larger g, the remaining two phases, namely the 'paramagnetic superradiant phase' (PSP) and the intermediate 'antiferromagnetic superradiant phase' (ASP), have a non-vanishing photon density.The photon density is therefore the relevant order parameter.Because of our prerequisites in Eq. ( 6), we concentrate on the former two phases, staying adiabatically connected to the case of no lightmatter interaction g = 0. We discuss our findings for the two phases in separate Subsecs.3.2 and 3.3, as the derivation of Eq. (1) from Eq. ( 17) depends on the phase.It turns out that the phase transitions out of the PNP and ANP obtained by mean-field calculations coincide with the closing of the lowest excitation level in the effective theory that we derived (as plotted in Fig. 3).
For g = 0, the Dicke-Ising Hamiltonian boils down to a longitudinal-field Ising model ignoring the light part of the Hamiltonian, as it behaves trivially in this limit.The phase transition between the paramagnetic and antiferromagnetic phase in this limit is of first order.The critical point J crit depends on the connectivity of the chosen lattice, yielding on a bipartite lattice with c being the number of bonds per site, e.g., c = 2 for a chain geometry.
For the limiting case ω 0 = 0, one recovers the case studied in [1] for 1D with the Hamiltonian In the work a first-order phase transition between the normal and superradiant phase was found at by comparing the ground-state energies of the two phases [1,36].In the weak light-matter coupling phases, the ground-state energy per site stays unaffected by perturbations in g [1], thus having the constant energy for PNP and ANP on bipartite lattices, respectively.In contrast, the excitations are affected by the light-matter coupling in these low-coupling phases, as we discuss in the following.

Paramagnetic normal phase
The paramagnetic normal phase is adiabatically connected to the limit J = g = 0 so that photons and spin flips are elementary excitations in this phase.Thus, it is suitable to introduce hardcore bosons in the same way as done in Eq. ( 3), resulting in the Hamiltonian n j n l (23) with E 0 = −ω 0 N /2 + J cN /2 and c being the number of bonds per site.We can now map H DI to Eq. ( 1) up to a constant offset by defining with δ •,• being the Kronecker delta, effectively setting all distances in Eq. ( 4) to zero or nearest neighbors (by defining a distance δ i with |δ i | = 1 as nearest neighbor) and fulfilling the prerequisites of Eq. ( 5).
As the low-energy subspace fulfills Eq. ( 6), we can apply the finding of the last section onto this model.The rescaled Dicke Hamiltonian in momentum space from Eq. ( 9) has the form and can be diagonalized in the thermodynamic limit by using a Bogoliubov transformation (see App. C for the derivation).The diagonalized Hamiltonian can be written as with b -being a linear combination of the a ( †) , b( †) 0 operators and the corresponding eigenenergies defining h ≡ ω 0 −2cJ as the coefficient of the rescaled magnetic-field term.The remaining part of the Hamiltonian does not couple to the light field and thus stays unaffected by varying g.We can write H MM in the original form with a correction term for k = 0 as Note that H MM is no longer diagonal in the real-space basis due to the correction term.Nonetheless, the interaction term and the correction term commute in the thermodynamic limit for the low-energy subspace, as can be shown analogously to Subsec.2.3.So, we obtain a flat dispersion for k ̸ = 0 in the one-magnon sector with energy ϵ k̸ =0 = ω 0 − 2cJ.The higher particle channels involve (anti)-bound states for neighboring magnons due to the interaction term in H MM .The energy of a free magnon pair is given directly as twice the one-magnon energy, as can be found using the free particle approximation [43].The (anti)-bound magnon-pair state is given by two magnons next to each other to (maximize) minimize the Ising term in Eq. (28).In the thermodynamic limit the energy of the (anti)-bound state is given exactly as ϵ bound = 2 • (ω 0 − 2cJ) + 4J, taking twice the free magnon energy plus the energy offset caused by the nearest-neighbor Ising interaction.
In Fig. 4 we plot the excitation energies for k = 0 varying the light-matter interaction g.For the plot, we have fixed the parameters of the model to be ω/h = 9/8, J/h = −1/20 with h ≡ ω 0 − 2cJ and a connectivity of c = 2, as present for a chain geometry.Note that other values for the parameters give qualitatively similar results.To keep the visualization light, we only plot the energies of states with up to two quasiparticles.For the unperturbed case (g = 0), we have the one-magnon (in blue), one-photon (in red), bound two-magnon (in orange), free two-magnon (in green and blue), and two-photon state (in red) from lower to higher excitation energies.As argued, the bound magnon state and the free magnon states, where all individual magnons have non-vanishing momenta, stay constant in energy as they do not couple with the light mode (in orange and green).The energies of the two-magnon and two-photon states are given by the doubled value of the respective one-particle energies, as they form bosonic modes in the thermodynamic limit.For the critical value the gap of the magnon mode closes with a critical exponent of 1/2 indicating a second-order phase transition.This critical value coincides with the one calculated by [2] with mean-field methods.Nonetheless, it does not rule out a potential first-order phase transition for smaller 0.0 0.2 0.4 g values, as we do not have any information about the energetics of the superradiant phase.Therefore, the first-order nature of phase transition lines around ω 0 = 0 observed by perturbation theory, exact diagonalization [1] and quantum Monte Carlo calculations [36] is in no contradiction to our findings.As can be seen in Eq. ( 27) whether the renormalized one-magnon or one-photon gap closes depends on which of the two have a higher unperturbed energy.For deeper insights of the excited states in low-energy space, like corresponding observables, we can use the eigenvectors obtained by the Bogoliubov transformation, as presented in App. C.

Antiferromagnetic normal phase
For the antiferromagnetic normal phase we will restrict ourselves to bipartite lattices (denoting A, B as the two subsets of sites) to avoid any kind of geometric frustration caused by the Ising interaction for J > 0. For models with geometric frustration, a straightforward mapping to the generalized Dicke model in Eq. ( 1) is not possible, as the model features an extensive ground-state degeneracy in the Ising limit, which is not contained in our model.The phase is adiabatically connected to the case of ω 0 = g = 0, with an unperturbed ground state with alternating spin orientations and an empty cavity mode.To bring the Hamiltonian closer to the form of Eq. ( 1), we first have to apply a rotation on one of the sublattices (here sublattice A).Therefore, we apply a sublattice rotation U about the x-axis giving for j ∈ A, altering the eigenfunctions but not the eigenvalues of the Hamiltonian.The rotated Dicke-Ising Hamiltonian possesses a fully polarized ground state for ω 0 = g = 0, as the sign in front of the Ising interaction is flipped by the transformation.The expression (−1) j is the shorthand notation for inspired by the usual notation in the case of a chain geometry.We perform a mapping to hardcore bosons, as done in Eq. ( 3), resulting in with E 0 = −J cN /2.While we can map parts of H RDI to H matter in Eq. ( 4) up to a constant offset with the Dicke Hamiltonian H D can not directly be recovered due to the staggered magnetic field in H RDI .Instead, we define two separate magnon modes in momentum space for the A and B sublattice of the form bk,s = 1 The form is analogous to the general case in Eq. ( 1), with the difference that there are two bosonic modes which are coupled via the light mode and the Ising interaction.We can prove analogously to before that the Ising term and the k ̸ = 0 modes decouple from the k = 0 magnon modes which are coupled to the light field for the prerequisites in Eq. ( 6) in the thermodynamic limit.Thus, we have again a simplified Hamiltonian at k = 0 describing the complete dynamics induced by the light field, analogous to the renormalized Dicke model HD , reading We can solve this system again with a Bogoliubov transformation, as shown in App.C, yielding the Hamiltonian  5: Low-energy spectrum of the Dicke-Ising model in the antiferromagnetic normal phase.The parameters are chosen as ω/h = 9/8, (2cJ + ω 0 )/h = 13/8, and c = 2 (chain geometry) with h ≡ 2cJ − ω 0 at vanishing momentum k = 0.For the spectrum only dressed states with one or two quasiparticles are considered, omitting the mixed 2QP states with one dressed magnon and photon for clarity of the plot.The three unperturbed magnon (photon) states hybridize for finite g leading to a decreased (increased) energy.At g crit , one of the elementary magnon gaps closes, indicating the end of validity of the found eigenvalues.The free magnons and bound magnons in the 2QP sector are not coupled to the light mode and therefore stay constant in energy when varying g. with three elementary excitation types, which are built out of the photon mode and the two magnon modes for the even and odd sites.While the eigenenergies ε + , ε -, ε ⋆ can be calculated analytically, we omit them here due to complexity of the analytic expressions (see the supplementary data for the full expressions) [44].
For the rest of H RDI we can argue, analogous to the paramagnetic phase with H MM , that all non-diagonal corrections for k = 0 will commute with the rest of the matter-matter interactions.So, we again obtain flat modes for k ̸ = 0 in the one-magnon sector and bound states in the higher particle channels.
We plot the low-energy spectrum for vanishing momentum varying the light-matter interaction in Fig. 5.We define h ≡ 2cJ − ω 0 as the lowest unperturbed excitation energy and fix the parameters to representative values ω/h = 9/8, (2cJ + ω 0 )/h = 13/8, and c = 2.As before, we only plot energies of states with up to two quasiparticles.The spectrum offers a richer structure than for the paradigmatic phase in Fig. 4 as we have two independent dressed magnon modes (both plotted in blue).This also results in three flat free-magnon bands in the 2QP sector (in green), coming from the combinations of the two elementary magnon modes.Using the free particle approximation, the energy of the flat bands is given as the sum of the energies of the 1QP modes.Note that the 1QP and 2QP sectors overlap for all values of g, as the bound magnon state (in orange), with two excitations next to each other, has a lower energy than the upper 1QP magnon band.In the thermodynamic limit the energy of the bound state is given as Figure 6: Comparison of the results for the 1QP sector obtained by the effective theory in the thermodynamic limit and exact diagonalization of finite systems.The parameters of the model are the same as in Fig. 4. In (a), the finite-system results are plotted as crosses for various system sizes 4 ≤ N ≤ 14, the results from the effective model are given as lines.In (b), the absolute energy difference between the effective model ϵ ∞ and the finite systems ϵ N is plotted for varying system sizes N .The different lines correspond to various g values and the magnon and photon mode.
critical value again coinciding with the one found within mean field [2].We are again not capable of detecting the first-order phase transition found with numerical methods [1,36] around the limit ω 0 = 0, as expected.

Comparison to results on finite systems
After presenting results for the Dicke-Ising model using the effective model in the thermodynamic limit, in this subsection we will compare the effective model with calculations on finite systems checking for the convergence towards the effective model when increasing the system size.We perform calculations using two methods: the first being exact diagonalization (ED) and the second being the perturbative series expansion method pcst ++ [3], treating the lightmatter interaction g as the perturbation parameter.To keep the discussion concise, we will limit ourselves to the chain geometry and the paramagnetic normal phase using the same parameters as in 3.2, while the same is possible for other geometries and the antiferromagnetic normal phase.
Exact diagonalization In Fig. 6 we compare our derived effective model with results on finite systems using ED.For that we look at the 1QP sector, namely the dressed one-magnon state (in blue) and the dressed one-photon state (in red).To determine the corresponding states in the ED calculations, we check the weight of the unperturbed 1QP states in the obtained eigenstates to choose the eigenstate which is closest to the desired unperturbed state.In Fig. 6(a) we see quite good agreement of the finite and the infinite systems for small light-matter couplings, while for couplings closer to the phase transition we obtain larger deviations.This is to be expected as the light-matter coupling, acting globally on the system, becomes non-negligible and therefore adds a stronger N dependence.Additionally, the ED results can be affected by the nearby phase transition, potentially happening before the closing of the gap by a first-order transition.
In Fig. 6(b) the absolute energy difference |ϵ N − ϵ ∞ | between the effective theory ϵ ∞ and the different finite systems ϵ N is plotted for various g values.All differences show the same overall exponent with respect to 1/N .Applying a fit to the data, we obtain an extrapolation of the form within the PNP for g/h < 0.35.In contrast, taking a point close to g crit , we obtain no clear exponent anymore (not shown).These deviations hint for a first-order phase transition, affecting the ED results in a way that is not taken into account in the effective model thus leading to a non-vanishing difference when taking the thermodynamic limit.We also calculate the operator norm of the commutator [ HD , H MM ], as defined in App.B, for the low-energy subspace of the finite systems (not shown).We find that the commutator norm clearly scales with N −1/2 , which fits perfectly to the derived scaling relation in Sec. 2. While we have no direct handle to the deviation of the energies of the finite systems in relation to the derived effective theory, we can motivate the scaling of N −1 in Eq. ( 41) with perturbation theory.As the two terms HD , H MM have a non-vanishing commutator, treating one as a perturbation to the other, yields a first potentially non-vanishing contribution in second order.As the commutator scales with N −1/2 -both shown in the analytic derivation and ED -we obtain the squared factor of N −1 in second order.This explains the found scaling of the energy in Eq. (41).It would be interesting to gain a deeper understanding between the interconnection of the energy and commutator scaling with N , which could be a future task.
pcst ++ In Fig. 7 we compare our effective model with results on finite systems using the perturbative series expansion method pcst ++ that was introduced in [3].The pcst ++ is a generalization of the perturbative continuous unitary transformation (pCUT) method [45,46], often used to treat quantum many-body systems with high-order linked-cluster expansions.In contrast to pCUT, pcst ++ is capable of having an arbitrary number of quasiparticle-types in the unperturbed part of the Hamiltonian, enabling us to apply it to the Dicke-Ising model, with g as the perturbation parameter.For more information about the method we refer to [3].
In Fig. 7(a) we plot the obtained series expansions in dashed lines, cutting the series at different orders, going from O2 to O6 for fixed N = 100.For both 1QP energies, the results match with the effective model up to around g/h ≈ 0.1 before diverging.This divergence, as explained in more detail in [3], comes from the energetic vicinity of the magnon and photon mode.Due to the perturbative treatment of the model, the series expansion's convergence radius depends on the energy differences between the individual unperturbed eigenenergies in the bare Hamiltonian giving rise to the obtained divergence.
While this issue can be overcome by modifying the underlying transformation, it is not needed for the following discussion.Instead, we calculate the series expansion of the effective model's energy around g = 0 and compare the coefficients of the individual orders with those of the series obtained by pcst ++ .The result is shown in Fig. 7(b).We plot the normalized total difference of the individual coefficients Figure 7: Comparison of the for the 1QP sector obtained by the effective theory in the thermodynamic limit and pcst ++ on finite systems.The parameters of the model are the same as in Fig. 4. In (a), the finite-system results are plotted as dashed lines for different orders, the results from the effective model are given as solid lines.In (b), the series expansion from pcst ++ is compared to the series expansion of the effective model in the thermodynamic limit.The relative difference of the coefficients of the orders is plotted vor varying system sizes, c N O denoting the coefficient of order O of the series for the system with N sites.obtained by the effective theory.The absolute error obeys again an extrapolation of the form for the given coefficients obtained by pcst ++ .It should be mentioned that the c N 2 for the onephoton mode already matches up to computer precision with c ∞ 2 and is therefore not shown in Fig. 7(b).This can be explained by the fact that there are no finite-size corrections in second order for the energy of the one-photon state, as there are no magnons in the finite system that would alter the possible processes in perturbation theory.As can be seen in the plot, this changes for higher orders.As for the ED case, we motivate the found scaling in Eq. ( 42) with the dominant second order perturbation theory between the two parts of the Hamiltonian.
All in all, the finite-size calculations are in perfect agreement with the findings of Sec. 2. While the results on finite systems differ from the effective model, the discrepancies shrink for larger N , allowing for an extrapolation in N with a vanishing difference in the thermodynamic limit.

Conclusions
In this work we considered a class of Dicke models with a generalized matter part consisting of long-range hopping and correlated processes while keeping the simple form of the single light mode and the Dicke light-matter coupling.In the limit of weak light-matter interactions and large system sizes, we were able to establish an analytical solution of the low-energy physics by mapping the system to an effective conventional Dicke model without matter-matter interactions.To prove the presented mapping, we demonstrated that the commutator between the two parts of the Hamiltonian vanishes in the thermodynamic limit for the given prerequisites.With this mapping in place, we are able to determine low-lying excitations on general footing by performing a Bogoliubov transformation, including the closing of the gap indicating a potential second-order phase transition.This reasoning is exemplified for the paradigmatic Dicke-Ising chain, which is compared to mean-field [2] and quantum Monte Carlo [36] results as well as finite-size calculations using ED and pcst ++ [3], yielding a satisfying agreement with our found mapping.
Next, it would be interesting to generalize our findings further, both with a more complex matter or light part.This includes adding more light modes, be it in a discrete or continuous way [38,47,48], and matter-matter interactions that do not conserve the number of spin-flip excitations.The latter terms are naturally present in a variety of interacting quantum manybody systems relevant for quantum-optical and solid state platforms including paradigmatic examples like the transverse-field Ising or Heisenberg model.This generalization would give the mapping a wider range of potential applications to foresee the impact of light-matter interactions on well-discussed strongly correlated matter systems.
While these potential generalizations offer a broader range of systems to apply the mapping to, it would be also an interesting question to find the fundamental limits where this type of mapping is applicable.Introducing terms that prevent the mapping to hold would therefore allow to find emerging physical phenomena that can not be described by a non-interacting matter part in an effective model.These insights could be a building block toward tailored effective models by specific light-matter systems to induce new types of effective interactions.
Restricting ourselves to finite particle numbers according to Eq. ( 6), the second term vanishes in the thermodynamic limit.We are left with the bosonic commutation relation [ bp , b † k ] = δ p,k .The remaining commutation relations for bosonic operators

B Calculating the commutator norms
In this section we derive the result of [ HD , H MM ] vanishing in the thermodynamic limit in more detail.As mentioned in Subsec.2.3 we do so by calculating an upper bound of the operator norm ∥ • ∥ for the given commutators.The operator norm ∥O∥ of an operator O is defined as with |s〉 being any normalized state with 〈s|s〉 = 1.For our discussion, we will restrict the choice of |s〉 onto those states that fulfill the prerequisites given in Eq. ( 6), namely 〈s| j n j |s〉 < ∞.
We will start the derivation by showing that the norm of the hopping operator does not scale with the system size and thus stays finite in the thermodynamic limit for any chosen δ.First, we use the definition of the norm in Eq. (B.1) to write it as an expectation value for the given operator with m 1 +m 2 finite and arbitrary δ m i values.Even though the new operator offers a much more complicated structure, the above argumentation can be followed analogously, as the sums are again cut to a finite subset and ∥b ( †) j ∥ ≤ 1, leading to an upper bound of for the generalized operator.As 〈s| j n j |s〉 ≡ n and |m 1 − m 2 | are bound to finite values, we can again conclude that the expectation value stays finite.We now apply these findings to the commutators in Eqs. ( 14) and ( 16).To do so, we bring the expressions into the form of Eq. (B.9).We will show the derivation only for one of the terms each, as the other can be derived analogously.We use Eq. ( 7), to rewrite the commutator Eq. ( 14) in in terms of real-space operators: For above calculation we have used [b j , b † l ] = δ j,l (1 − 2n l ) and k e ik(l−m+δ) = N δ l+δ,m .When taking the norm of above commutator, we can use the identities ∥A B∥ = ∥A∥∥B∥ and ∥A + B∥ ≤ ∥A∥ + ∥B∥ to obtain As the operators l n l−δ b l and l n l fulfill the form of Eq. (B.9) their norms stay finite.For the remaining b0 the same can be shown.So, the overall commutator scales with N −1/2 and therefore vanishes in the thermodynamic limit.We move on with the commutator in Eq. ( 16).We write down the commutator in real-space operators and use the respective commutation relations: Eq. (B.16) can be expressed in terms of four operators of the form of Eq. (B.9), leading to the conclusion that the overall commutator scales with N −1/2 .To sum up this section, we have introduced the operator norm as a suitable measure to show the decoupling of HD and H MM in the thermodynamic limit by finding an upper bound for a general operator O in Eq. (B.9) and tracing back the commutator [ HD , H MM ] to terms of this form.

C Solving the Dicke model
In this section we derive the solution of the effective Dicke models occurring in the main text, namely Eqs. ( 9), (25), and (37).To keep this derivation concise, we will not cover all technical details needed for the transformations to work.For a more technical and detailed work on the Bogoliubov transformation itself, we recommend the paper by Xiao [55].For this section we will stick closely to the formulation in this paper.
The general starting point for all three Hamiltonians referenced above is their quadratic nature, namely that all terms consist out of exactly two operators apart from a constant term.All three Hamiltonians can therefore be rewritten in a general form as where b ( †) i are annihilation (creation) operators with pairwise bosonic commutation relations and α i, j , γ i j ∈ .Note that these commutation relations only hold in the thermodynamic limit for our models (see App. A) as well as for the bare Dicke model [42].Thus, the presented solution only applies in the limit N → ∞, which is nonetheless our point of interest.
For this class of models Bogoliubov introduced a linear transformation to diagonalize these systems [56].While the transformation was also generalized to particles with fermionic commutation relations [57][58][59], we can limit ourselves to the bosonic case here.It is convenient to write the Hamiltonian of Eq. (C.1) in matrix notation as where α, γ are matrices built out of the scalars α i, j , γ i, j and b, b † are vectors with the annihilation and creation operators b i , b † i , respectively, as their entries.The Hamiltonian is solved by introducing a linear transformation of the form b = Ad + Bd † (C.3) with d ( †) being vectors, analog to b ( †) , consisting out of new bosonic operators d i , d † i and coefficient matrices A, B. These matrices have to be chosen such that they diagonalize H and at the same time preserve the desired bosonic commutation relations for the d ( †) i operators, which is not guaranteed to be possible.In [55] an equivalence between the existence of such a transformation and the diagonalization of a so-called dynamical matrix D ≡ I 0 0 −I • M , with I being the identity matrix, (C.4) is proven, which is advantageous for actual calculations.If and only if D is diagonalizable with all its eigenvalues being real, the Bogoliubov transformation can be applied.To keep this section concise, we refer to Xiao for the derivation of this finding [55].If the transformation exists, the resulting Hamiltonian can be written as with ω i being the eigenenergy of the i-th elementary excitation and a constant term C. The eigenvalues ω i are the positive eigenvalues of D if the transformation exists.
For our models discussed in the main text, we find that the transformation is applicable for the case of small light-matter couplings.The transformation starts failing when the respective gap of the normal phase closes.The Hamiltonian in Eq. ( 25) can be written in matrix notation in terms of Eq. (C.2) as (C.6) Diagonalizing the dynamical matrix yields the positive eigenvalues of Eq. ( 27).The general Hamiltonian of Eq. ( 9) can be solved analogously.For the Hamiltonian of the antiferromagnetic phase in Eq. ( 37), we have to take into account three types of quasiparticles, resulting in the 6 × 6 matrix with a/ 2 = 2cJ − ω 0 , b/ 2 = 2cJ + ω 0 , which can be solved with a computer algebra program.As the eigenvalues can not be written down in a concise form, we omit to write them down here [44].

1 Figure 2 :
for a paradigmatic system.It contains Sketch of the low-energy spectrum of a paradigmatic Hamiltonian H in Eq. (

1 Figure 3 :
Figure 3: Phase diagram of the Dicke-Ising model for the fixed representative ratio ω = 4ω 0 /3 in 1D (c = 2).The four phases PNP, ANP, ASP, and PSP are given in different colors.The phase transitions obtained by mean-field theory in[2] are plotted in solid lines.For the PNP and ANP we discuss effective models for the low-energy spectrum becoming exact in the thermodynamic limit.The phase transitions out of the two normal phases coincide with the closing of lowest gap in the presented effective theory, denoted with black solid lines.For J = 0, we recover the standard Dicke model, indicated with a faint dashed line.

1 Figure 4 :
Figure4: Low-energy spectrum of the Dicke-Ising model in the paramagnetic normal phase.The parameters are chosen as ω/h = 9/8, J/h = −1/20, and c = 2 (chain geometry) with h ≡ ω 0 − 2cJ at vanishing momentum k = 0.For the spectrum only dressed states with one or two quasiparticles are considered, omitting the mixed 2QP states with one dressed magnon and photon for clarity of the plot.The unperturbed magnon (photon) states hybridize for finite g leading to a decreased (increased) energy.At g crit the magnon gap closes, indicating the end of validity of the found eigenvalues.The free magnons and bound magnons in the 2QP sector are not coupled to the light mode and therefore stay constant in energy when varying g.
with s, s ′ ∈ {A, B}.The new operators fulfill pairwise bosonic commutation relations in the low-energy sector as [ bp,s , b † k,s ′ ] N →∞ − −−− → δ p,k δ s,s ′ .With this we rewrite H RDI as ) being again the sum of the individual magnons (one from sublattice A, one from sublattice B) plus the energy correction from the Ising interaction.The lower magnon mode closes for the 0 where c N O is the coefficient of order O for a system of N sites obtained by pcst ++ and c ∞ O is the respective coefficient of the series 0 with |c s→s ′ | ≤ 1 and S s being a finite set of states s ′ depending on the initial state s.We can find an upper bound of |S s | ≤ n because of the annihilation operator b l+δ acting once on all j i in Eq. (B.3).As Eq. (B.5) holds for any state s in real-space basis, we can conclude thatO † O |s〉 = s ′ ∈S s s ′′ ∈S s ′ c s→s ′ c s ′ →s ′′ |s ′′ 〉 (B.6) consists out of a maximum of |S s | • |S s ′ | ≤ n 2 terms independent of the system size with |c s→s ′ c s ′ →sAs 〈 j n j 〉 has to stay finite according to our setup (see Eq. (6)), we can conclude that 〈O † O〉 stays finite, too.Last, we generalize this statement to arbitrary states, being a (potentially infinite) superposition of the states of Eq. (B.4).Assuming a normalized eigenstate of O as 3) While s may be any valid state, let us for the moment consider an n-magnon state in real space of the form |s〉 = | j 1 , j 2 , . . ., j n 〉 (B.4) where the j i denote a magnon at position j i .Applying O on this state leaves us with a state of the form O |s〉 = ′′ | ≤ 1.We can therefore conclude that 〈s|O † O|s〉 ≤ n = 〈s| j n j |s〉 .(B.7)