Superradiance, charge density waves and lattice gauge theory in a generalized Rabi-Hubbard chain

We investigate a one-dimensional Rabi-Hubbard type of model, arranged such that a qdot is sandwiched between every cavity. The role of the qdot is to transmit photons between neighboring cavities, while simultaneously acting as a photon non-linearity. We consider three-level qdots in the $\Lambda$ configuration, where the left and right leg couples exclusively to the left or right cavity. This non-commuting interaction leads to two highly entangled incompressible phases, separated by a second order quantum phase transition: the qdot degrees-of-freedom act as a dynamical lattice for the photons and a Peierls instability breaks a second $\mathbb{Z}_2$ symmetry which leads to a dimerization in entanglement and photon number. We also find a normal insulating phase and a superfluid phase that acts as a quantum many-body superradiant phase. In the superradiant phase, a $\mathbb{Z}_2$ symmetry is broken and the phase transition falls within the transverse field Ising model universality class. Finally, we show that a limit of the model can be interpreted as a $\mathbb{Z}_2$ lattice gauge theory.

We investigate a one-dimensional Rabi-Hubbard type of model, arranged such that a qdot is sandwiched between every cavity. The role of the qdot is to transmit photons between neighboring cavities, while simultaneously acting as a photon non-linearity. We consider three-level qdots in the Λ configuration, where the left and right leg couples exclusively to the left or right cavity. This non-commuting interaction leads to two highly entangled incompressible phases, separated by a second order quantum phase transition: the qdot degrees-of-freedom act as a dynamical lattice for the photons and a Peierls instability breaks a second Z2 symmetry which leads to a dimerization in entanglement and photon number. We also find a normal insulating phase and a superfluid phase that acts as a quantum many-body superradiant phase. In the superradiant phase, a Z2 symmetry is broken and the phase transition falls within the transverse field Ising model universality class. Finally, we show that a limit of the model can be interpreted as a Z2 lattice gauge theory.

I. INTRODUCTION
Light has since always served as a tool for detecting states of matter; whether it is our own eyes registering our surrounding, or photon detectors measuring the wavelengths of light from distant stars. Light can also be used to control and change the states of matter. In modern times, for example, we can cool and trap individual particles with designed laser light. At the quantum level of single photons, strongly coupled light and matter can form novel states with no counterparts in other branches in physics [1]; new interacting quantum manybody models derive which may host exotic phases. The study of phase transitions (PTs) in light-matter systems dates back to the early days of the laser, when the onset of lasing with increasing pump power was identified as a non-equilibrium continuous PT [2].
Predating the laser, in 1954 Dicke showed how the rate of spontaneous emission for a set of N two-level atoms could be enhanced by a factor √ N [3]. This phenomenon, arising due to collective multi-partite interference, has been termed superradiance and can be derived from the Dicke model describing the coupling of N identical twolevel systems with a single photon mode [4]. In 1973, first by Hepp and Lieb [5] and shortly afterwards by Wang and Hioe [6], it was demonstrated that the Dicke model supports a second-order PT from a 'normal' to a 'superradiant' phase as the light-matter coupling is raised above a critical value. The corresponding PT is accompanied by a spontaneous breaking of a Z 2 symmetry. In more recent times, there has been a strong interest in realizing PTs, like the Dicke one, in quantum optical lattice systems [7]. This development was spurred by the increased experimental control over many-body quantum systems and in the wake of quantum simulators [8].
One manifestation of a quantum PT is a non-analytic behavior of the ground-state [9]. For the Dicke PT one finds such a non-analyticity, and as such it has often been referred to as a quantum PT (QPT). Nevertheless, the transition is of the mean-field type, and the role of quantum fluctuations becomes irrelevant in the thermodynamic limit [10]. In this respect, the normal-superradiant transition cannot be driven by quantum fluctuations, and depending on ones personal taste one may question referring to the Dicke PT as a proper QPT.
In a 'true' many-body quantum normal-superradiant PT, as the light-matter coupling g is increased the system would enter a superradiant phase while quantum fluctuations remain extensive in the thermodynamic limit. This is more in the vein of paradigm quantum critical models like the transverse field Ising or Bose-Hubbard models [11]. As a way to enhance the role of quantum fluctuations one may instead consider multimode Dicke models [12] or cavity arrays [13]. In the latter, arrays of cavities/resonators are manufactured on microchips such that photons can tunnel between neighboring cavities. Each transmission line resonator is equipped with a quantum dot that acts as an artificial two-level atom [14]. This produces a Jaynes-Cummings nonlinearity [15], which acts as an effective photon-photon interaction. The resulting physics (neglecting photon losses) is essentially the same as for the Bose-Hubbard model with insulating and superfluid phases of polaritons [16].
In a more recent work, the light-matter terms of the model were used to produce an effective interaction as well as the kinematics [17]. The authors studied a onedimensional array of resonators, with a two-level system (qubit) placed between neighboring resonators, such that photon tunneling is mediated by a Jaynes-Cummings interaction. It was demonstrated that the low energy physics is described by the transverse field Ising model. While not discussed in [17], we may think of the emerging Ising transition as a normal-superradiant PT that is driven by quantum fluctuations; in the normal phase, all qubits are to a good approximation in their lower state and the cavity modes in vacuum, while if the lightmatter coupling is increased beyond a critical value g c , the photon modes as well as the excited qubit states are populated.
In this work we also consider a cavity array with mediated photon tunneling, but instead of qubits we take three-level systems, i.e. qutrits with three internal states |1 , |2 , and |3 . We assume strong light-matter couplings g ∼ 1. The novel features of our model are due to the non-commuting structure of the interaction terms; every other resonator couples to the |1 ↔ |3 transitions of its left and right qutrit, while every other instead addresses the |2 ↔ |3 transitions [18].
It may be helpful to picture the photons as 'site variables' and the qutrits as 'bond variables', forming a 'dynamical lattice' for the photons as in quantum link models [19]. A salient result that may arise from a dynamical lattice is an altered periodicity, leading to novel phases such as charge density waves [20] first predicted by Peierls in one dimension, and supersolids [21]. Indeed, for not too large coupling strengths we find a phase (CDW) with a charge density wave order. Contrary to a traditional Peierls instability for free fermions, the filling in our model is not set by an external chemical potential but is determined mainly by the light-matter coupling g, and there are no strict constraints on the filling for the appearance of the CDW. This is akin the 'photonic Peierls transition' recently discussed in an extended Bose-Hubbard model with a similar dynamical lattice as the one considered here [22]. Beyond the CDW, as the coupling is increased further, there is a phase transition to another phase (NE), characterized by unbroken translational symmetry and high pairwise entanglement over the bonds.
Neither of these phases are superradiant, i.e. polaronic superfluids. In order to open up for the aforementioned quantum normal-superradiant PT we assume a static dipolar field driving the |1 ↔ |2 transition. Such a field has a 'resetting' effect on the qutrits; a qutrit making the transition of say |1 → |3 → |2 by transferring one photon between two resonators can be reset by the field to its original state |1 without involving further photons from the resonators. Above some critical drive amplitude s c and coupling strength g c we find a superradiant phase which spontaneously breaks the Z 2 symmetry. For weaker couplings, on the other hand, a symmetric normal phase emerges. In addition, our numerical results indicate a multi-critical point where all four phases meet. All transitions are found to fall within the universality class of the transverse field Ising model, i.e. not of the mean-field type but truly quantum by nature.
The outline of the paper is as follows. In the next section we introduce our model and find its symmetries. In Sec. III we present the zero temperature phase diagram as obtained from DMRG calculations, and characterize the different phases. We also discuss the possible transitions between the phases. Finally, in Sec. IV we conclude FIG. 1. (Color online) Schematic picture of the system set-up (a) and the qutrit-cavity coupling (b). An array of identical cavities supporting a single mode, with frequency ω, are coupled via three-level Λ-systems on the bonds between the cavities (a). These qutrit Λ-systems host two degenerate bare ground states |1 and |2 , and an excited state |3 separated from the lower states by Ω, see (b). One of the legs of the Λ-system is coupled to the cavity to the left and the other leg is coupled to the cavity to the right. Selecting the transitions in this particular way can be achieved by alternating the polarization, σ ± , between the cavities. Photon transport between consecutive cavities are only accompanied by transitions within the corresponding Λ-system. In the rotating wave approximation we would have a qutrit transition |1 → |2 for a photon moving to the right from cavity i to cavity i + 1, and vice versa for the other direction. The counter rotating terms, neglected within the rotating wave approximation, break this property. In addition, we add a direct coupling between the |1 and |2 states with an amplitude s, as depicted in (b). As explained in the main text, a polaronic superfluid (superradiant phase) is only possible for a non-zero coupling s.
with a summary.

II. MODEL SYSTEM AND ANALYSIS
We investigate a one-dimensional array of cavities, where photon tunneling between each pair of neighbor cavities is mediated by quantum dots. The system is schematically pictured in Fig. 1 (a). The idea of interconnecting cavities via superconducting qubits dates back to the works of [23], where the qubits were considered as controllable switches for the tunneling of photons between the cavities. Based on these ideas, it was possible to experimentally connect three resonators with the help of qubits [24]. The extension to arrays of multiple cavities was also later considered [17]. The system is then truly in the quantum many-body regime and the low energy physics can be seen as bosons living on a 'dynamical lattice' [22]. In this work we extend the model of [17] to three-level quantum dots, where the transitions between the qdot levels can be addressed separately.

A. Full Hamiltonian
We assume ideal qutrits (three-level systems) in the Λ configuration: the two states |1 and |2 are degenerate, while the state |3 is of a higher energy. Angular momentum must be preserved, so if, e.g., the lower states are angular momentum states with m = ±1 and the excited state has m = 0, then the qutrit transition 3 → 1 (3 → 2) emits a σ + -polarized (σ − -polarized) photon. One possible method of realizing our model would be by using polarization-selective cavities. Due to these selection rules, the qutrit transition 3 → 1 (3 → 2) can only emit a photon into the left (right) cavity. We return to the question about experimental realization in the conclusion. As mentioned in the introduction, we also introduce a static dipolar field between the two lower qutrit states |1 and |2 . Note that such a field breaks the angular momentum conservation, which is possible to realize by external two photon driving [25]. The level configuration and system parameters are shown in Fig. 1 (b). The full model system, as shown in Fig. 1 (a), consists of an infinite chain of cavities interspersed with qutrits. We label each qutrit and the resonator to its right by the same site index i. Cavities with even (odd) index select σ + -polarized (σ − -polarized) modes and are coupled to the 3 ↔ 1 (3 ↔ 2) transition with a Rabi interaction of strength g. We do not restrict the analysis to moderate couplings g ω, in which the rotating wave approximation is applicable, but allow for g ∼ ω, i.e. the deep strong coupling regime [26].
The Hamiltonian can be split in a bare and an interaction part The bare Hamiltonian is comprised of the harmonic oscillators representing the single light modes of the resonators, the bare energies of the qutrits, as well as the dipolar field on the two lower qutrit states. Expressed in terms of Gell-Mann matrices [27] (see Appendix A for the definition of the Gell-Mann matrices λ (α) i ), the bare Hamiltonian reads ( = 1 throughout) where ω is the resonant angular frequency of the cavities, Ω is the frequency (energy) of the upper qutrit states, s is the strength of the dipolar field on the lower qutrit states and a † i (a i ) are the photon creation (annihilation) operators for the i'th cavity, i.e. for a photon Fock state of the i'th cavity, a i |n i = √ n i |n−1 i , a † i |n i = √ n i + 1|n+1 i , and for the number operator n i |n i = a † i a i |n i = n i |n i .
The interaction Hamiltonian can be written as where g is the coupling strength.
We note that there is a duality which maps the Hamiltonian H(s, g, ω, Ω) → H(−s, g, ω, Ω). The unitary transformation is Hence, it is sufficient to consider s > 0, and we remark that for any superradiant phase for s < 0 there exists a corresponding superradiant phase for s > 0, but with a staggered ordered parameter a i ∝ (−1) i .

B. Effective qutrit Hamiltonian
It is possible to eliminate the photon degrees-offreedom to derive an effective model of interacting SU (3) spins in one dimension. To obtain such an effective model, we employ the polaron transformation (also called Lang-Firsov transformation) [28,29], which is an ansatz for the ground state relying on a time-scale separation between fast and slow variables. The method has proven efficient in describing phonon physics, where the states of the ions in a crystal are 'displaced' according to the electronic state; phonons (fast variables) dress the electrons (slow variables) to form a polaronic excitation [30]. If the phonons are displaced, the idea is to find the new minima, i.e. the displaced phononic vacuum [29]. The variational ansatz of the polaron transformation for our model is given by where we have defined | α = i |α i . Here, {c i } are the 3 N coefficients of the entire many-body state of the qutrits, α i ∈ R is the amplitude of a photonic coherent state |α i , and γ is a real parameter which is later chosen such that the photonic and qutrit variables are decoupled. While this ansatz is a mean-field ansatz for the photonic variables, it allows for general quantum correlations in the qutrit variables. We define the unitary polaron transformation as where we have defined Taking the expectation value with respect to the photonic variables, α|H| α , yields the effective qutrit Hamiltonian As motivated in the introduction, one reason for studying the present model is to find a truly quantum normalsuperradiant PT. Similar to the Dicke model but with extensive quantum fluctuations in the thermodynamic limit. Our Hamiltonian is invariant under a π-rotation with respect to the total excitation number i.e. the corresponding unitary is Since the eigenvalues of N ex are integers, we clearly have Π 2 = I as required for a Z 2 symmetry. That the Hamiltonian is symmetric under the action of Π follows from noticing Apart from a few special exceptions, continuous phase transitions occur when the ground state of a system spontaneously breaks the symmetries of its Hamiltonian [31]. We can then find a (local) order parameter, which is zero in the symmetric phase and non-zero in the symmetry broken phase. For the normal-superradiant PT, which is connected to the breaking of a similar Z 2 symmetry, a proper choice is φ = a [32]. Given a real light-matter coupling g, the order parameter will be real; the sign of φ determines the parity [10]. Similarly in our model, breaking of the parity symmetry results in a real non-zero value of where sign(s) i ensures that the order parameter works whether the parity breaks ferro-magnetically (for s > 0) or anti-ferromagnetically (for s < 0). The magnitude of φ is the same for every lattice site, but the sign can alternate between neighboring sites; constant signs represent ferromagnetic order, while alternating signs antiferromagnetic order. Since the symmetry Π involves both qutrits and resonators, there are two other possible order parameters which should show the same (anti)ferromagnetic order: As a one-dimensional spin chain supporting a Z 2 symmetry, the critical point should belong to the universality class of the transverse field Ising model. To elucidate this, we investigate the limit of s → ∞ and extrapolate to finite s. By a Hadamard rotation of the qutrit states the Hamiltonian can be rewritten as The bare Hamiltonian has eigenvalues E 1n = s + ωn, E 2n = −s + ωn, E 3n = Ω + n. (17) The population in the state with bare energy E 2n will be small if s is large enough and the mean photon number n of each resonator is small. In this limit, we project the qutrit degrees of freedom of the Hamiltonian by P = where σ z i = |1 i i 1| − |3 i i 3| and σ x i = |1 i i 3| + |3 i i 1|. In this limit, the model is equivalent to that considered by [17], who already concluded that the phase transition falls within the Ising universality class. By the duality (4), we conclude that the Hamiltonian is 'antiferromagnetic' for s < 0.

D. Quasi-translation symmetry
Assuming periodic boundary conditions, there is a second global symmetry of the model which we call quasitranslation symmetry: a permutation of the qutrit states |1 , |2 , followed by a translation with one site. Calling itT , the action on the operators is for α = 1, 8. Since this symmetry involves no change in parity of the photonic operators, no superradiance is involved in a corresponding symmetry breaking phase transition. It can still be called a parity symmetry, sincẽ T 2 = 1, but the parity change involves the qutrit operator λ E. Local gauge symmetries for s = 0 In the limit that s = 0, we find an infinite set of local conserved quantities where B k i = e iπ|k i i k| (k = 1, 2). We note that the global parity symmetry Π = i Π i and that Π i have the form of generators of a Z 2 gauge theory, where the photons play the role of the 'matter field' and the qutrits are bond variables [19]. Gauge theories can lead to interesting phenomena such as confinement, but in one-dimensional quantum systems the gauge theory is 'trivial' [33]. We discuss further aspects of gauge theory in the conclusion.
It is unclear whether the presence of these gauge symmetries implies that the model is integrable along the s = 0 line. For the quantum Rabi model, comprising a single spin-1/2 particle interacting with one boson mode, the Z 2 symmetry implies that the model is integrable [34]. The present model cannot, however, be thought of as interconnected Rabi models since we rather have an SU (3) algebra instead of an SU (2) one.

III. PHASE DIAGRAM
In this section we present the numerical results obtained from DMRG simulations using the TenPy library [35]. In principle, for gapped systems, DMRG will  (13) and (20), calculated with DMRG for the full boson-qutrit model (a) and the effective qutrit model (b). In the simulations, the maximal bond dimension χ = 100 was used and the bosonic Hilbert spaces were truncated to dimension 4. On the two axes we have s and g. For larger values of the light-matter coupling g, the photon number ni increases and we expect a breakdown of the adiabatic elimination. This is seen by the quantitative differences between the phase diagrams. Nevertheless, the structure of the phases are the same. This supports the claim that the two models share identical universal properties. We distinguish between four different phases; normal (N), charge density wave (CDW), superradiant (SR), and normal-entangled (NE). The remaining (dimensionless) parameters used for the calculation are ω = 1 and Ω = 1. In figure (b), the sweeps made to calculate figures 4 and 5 are indicated with white lines. be able to reproduce the zero temperature phase diagram to arbitrary accuracy provided one keeps large enough bond dimension [36]. In practice, to extract the critical exponents can be computationally costly. For bosons, DMRG relies on truncating the infinite local Hilbert space dimension on every site. In our model, to reach numerical convergence close to the critical point may imply including up to 10 Fock states per site, making the simulations rather time-consuming. To circumvent this slow convergence, we perform some of the DMRG simulations for the effective qutrit model (8) which shares the symmetries and qualitative features of the full model. This is justified when we focus on universal critical features, and not quantitative results.
In Fig. 2 we show the sum of the two order parameters, (13) and (20), for the full model and the effective qutrit model in the (s, g)-plane. The characteristics of the four phases, normal (N), charge density wave (CDW), superradiant (SR), and normal-entangled (NE), will be specified below. We were not able to determine conclusively whether the four phases merge in a single multi-critical point involving all four phases, or in two tri-critical points involving only three phases each. Nevertheless, within the precision of the numerics, we lean towards the former scenario. In obtaining the figures, The bosonic Hilbert spaces were truncated to dimension 4. The entanglement entropy is identical between resonators of odd index and the left and right adjacent qutrits, but differs from the entanglement between resonators of even index and its adjacent resonators. (b) Expectation value of the photon number operator a † i ai , for odd (dashed line) and even (dotted line) indices, calculated from the same DMRG simulations. The spontaneous breaking of theT symmetry (the CDW in photon number) is evident in the "charge" dimerization in (b) and in the staggered entanglement entropy around odd and even resonators in (a). Note that entanglement entropy is high inside both phases, such that a non-trivial bond dimension χ is needed.
we used a maximal bond dimension of χ = 100. It was checked that setting a higher bond dimension did not improve the quality of the figure. We remark that, while the qualitative structure of the phase diagrams are the same, the polaron ansatz seems to give a bad quantitate approximation in most of the phase diagram. This is in contrast to ( [17]) who find quantitative agreement with a polaron ansatz.
There are two symmetry-breaking phases in the model: φ = 0 in the SR phase, and ϕ = 0 in the CDW phase. We mentioned above that the SR phase can alternatively be viewed as a superfluid, in contrast to the CDW phase which is an insulator in the sense that the superfluid order parameter φ = 0. We find no phase with both symmetries broken simultaneously, i.e. a supersolid. The N and NE phases are both symmetric with respect to both symmetries.
The normal (N) and the normal-entangled (NE) phase are both fully symmetric, and thus considered identical in the Landau classification of phases. One qualitative difference between them is the entanglement entropy [37]. We consider the von Neumann entropy where ρ 1 is the reduced density operator of subsystem 1, obtained by tracing the full density operator over the degrees-of-freedom of subsystem 2. It is understood that the trace in the above expression is over the degreesof-freedom for subsystem 1. The subsystems 1 and 2 are defined from splitting the full system into two equal halves. We expect an 'area law' behavior of the entropy in gapped phases, i.e. up to logarithmic corrections the entropy is proportional to the length of the boundary between the two subsystems [37], which in one dimension is a single bond. For gapless phases, on the other hand, we anticipate a 'volume law' where the entropy instead grows linearly with the subsystem size. The same linear growth in entanglement is expected at the critical points.
In all phases, we find that the entanglement entropy saturates In Fig. 3, we plot the expectation value of the photon number as well as the entanglement entropy across the CDW-NE phase transition for the full model of qutrits and resonators. In obtaining the figures, the matrix product states were truncated to a bond dimension of χ = 100 and the bosonic Hilbert spaces were truncated to dimension 4. The entanglement entropy in the two phases is shown in Fig. 3 (a); in the CDW phase, the entanglement is higher between resonators of odd/even index and its adjacent qutrits. The state chooses a higher entanglement around odd or even sites. This staggered structure of the entanglement entropy is another indication of the breaking of the Z 2 quasi-translation symmetry. Indeed, in Fig. 3 we see that also the photon number breaks quasitranslation symmetry. The CDW is an entangled state which is best described by a MPS (matrix product state) of moderate bond dimension. Observe that the NE phase is fully symmetric. This is reflected in that entanglement is constant in all bonds, and no superradiance is observed. It is a gapped phase since the DMRG simulations converge for a moderate but quite high bond dimension. In the entanglement structure and MPS description, the NE phase is similar to the AKLT state. This paradigm many-body spin state has been found to be a symmetryprotected topological state for spin S of odd integer values. One telltale sign of a symmetry-protected state is a two-fold degeneracy in the entanglement spectrum of the MPS. We have found no such degeneracy in the entanglement spectrum, which indicates that the NE phase is not a symmetry-protected topological state [39,40]. This raises one issue: if the Landau paradigm is not enough to distinguish the N phase from the NE phase, and no topological [41] feature distinguishes them, there is the possibility that the N phase could be smoothly deformed into the NE phase.
Let us briefly return to the discussion about the noncommuting interactions. Assume no dipolar field on the qutrits, i.e. s = 0, and assume that all qutrits are initially prepared in states |1 i (|2 i ). Then, photons are only able to tunnel between sites i, i + 1 for i odd (even), precisely because of the non-commuting structure of the interactions. This 'blockade' causes the ground state to be an insulator, either of CDW or NE type (the s = 0 line in the phase diagram). These two phases survive for small non-zero couplings s, but for some critical s c the two phases terminate. If the light-matter coupling g is small, the photon number remains small and the system stays an insulator even for large |s| values -the normal phase. In order to find a superfluid state, g has to be increased and the system enters into the SR phase. This is further evidence of a charge density wave and not the aforementioned supersolid. One could possibly imagine that the vanishing order parameters, φ = 0 and ϕ = 0, in the NE phase is some numerical artifact and the two symmetries are indeed broken. However, this is physically unlikely since for at least s = 0 we should find an insulator with φ = 0.
Having identified the phases, we also wish to determine the types of transitions. With the two Z 2 symmetries broken in the SR and CDW phases, one expects that the PTs fall within the universality class of the transverse field Ising model in one dimension [11]. Indeed, we find that, close to the phase transitions, our numerically calculated order parameters can be fitted with excellent agreement to a power law with the exponent β = 1/8, i.e.
where Θ(x) is the Heaviside step function. This is demonstrated in figs. 4 and 5 displaying the order parameters (13) and (20) respectively. The numerical data are complemented with the expected fits (23). We cross the critical points by varying g, while keeping s fixed. The symmetry broken phases SR and CDW can be reached from either the N or the NE phases, and in all cases we have numerically verified that the exponent β = 1/8 is obtained. In fig. 4, a maximal bond dimension of χ = 100 was used, while in 5, the maximal bond dimension was χ = 200.

IV. SUMMARY AND CONCLUDING REMARKS
In this work we have studied a generalized Rabi-Hubbard chain. A few aspects are crucially different between our model and that of the standard Rabi-Hubbard chain, which motivated our study. Let us summarize these below and explain how our work differs from earlier findings: 1. The qdots are not contained within the resonators but rather positioned between consecutive resonators. Photon tunneling between any two resonators occurs via the qdot, and as such the qdot acts as a 'bond variable' like in the quantum link models that have been analyzed for dynamical gauge theories [19]. Thus, we can picture the system as bosons in a dynamical lattice -the lattice itself constitutes quantum degrees-of-freedom. We know that dynamical lattices can give rise to novel phases and phenomena -the most well-known example being the Peierls instability [42]. It is not restricted to fermions (electrons) coupled to lattice bosons (phonons), but may occur also for bosons in a dynamical lattice [22] and in spin systems [43]. A similar transition was also found in our system, manifested as the appearance of a CDW phase of period 2.
2. While the qdots are necessary for any dynamics (they couple the bare photon oscillators), they also serve as a nonlinear medium, i.e. the finite local Hilbert space dimension implies a nonlinear spectrum and the occurrence of insulating phases.
3. Instead of qubits we considered qutrits. For qubits the physics is qualitatively described by a transverse field Ising model [17], with the insulating (normal) phase being the symmetric one, and the superfluid (superradiant) phase the symmetry broken one. In such a model, the more exotic CDW and NE phases do not appear. The specific coupling of our Λ qutrits breaks down the local gauge symmetries Π i to a global symmetry Π, as explained in the main text. Without the direct coupling of the lower qutrit states, i.e. s = 0, a qutrit acts as a single photon 'diode'; if one photon has tunneled to the left, a second one cannot tunnel. Thus, the phases for s = 0 are insulating and the order parameter φ of Eq. (12) vanishes.
4. For s = 0, the Hamiltonian is invariant under the infinite set of local symmetries (21). It has the structure of a Z 2 gauge theory. The Mermin-Wagner theorem tells us that in one dimension a continuous symmetry cannot be spontaneously broken to generate long range order (superfluidity) [44]. Elitzur's theorem [45] tells us that gauge symmetries can never be broken, and only gauge invariant observables can have non-zero expectation values. One may speculate if this is related to the fact that we find no superfluidity in the NE phase.
It is interesting to note that this phase extends to finite s, where the gauge theory description does not hold strictly. Note also that the superradiant phase with φ = 0 is not contradicting the Mermin-Wagner theorem since a discrete symmetry can be broken in one dimension at zero temperature.
While our generalized Rabi-Hubbard model is simply extended from the traditional one, a most relevant question is whether it is also accessible experimentally. The main experimental novelty is the selected coupling between a qutrit and its resonators to the left and right. Such selection is easily achieved by adjusting the photon frequencies to become resonant with their respective qutrit transition, but has the disadvantage of leading to a staggered self-energy ω of the resonators. While the symmetry Π connected to the N-SR phase transition survives, such a modification destroys the quasi-translation sym-metryT . By numerical simulation, we have confirmed that the SR phase is robust to such changes [46], while unfortunately the CDW and NE phases do not survive. Consequently, to access the full phase diagram we need to implement the selection rules by other means, for example employing polarizations as discussed in the main text. This is rather straightforward for Fabry-Prot cavities, but technically more involved for transmission line resonators. Another experimental aspect is how to reach the deep strong coupling regime. This is typically done with the help of external Raman driving [47], such that the effective light-matter coupling g is controlled by a classical field amplitude. Another option how to fullfil the specific selection of the transitions is to externally address different Raman transitions in the qutrit and so controlling how the frequencies of the drive fields hit the desired resonances.
The largest experimental hindrance is, however, the loss of photons. The N-SR PT in the Dicke model survives photon losses, even if the universality class is altered [10] (note that we consider the deep strong coupling regime such that the photon number in the ground state can be non-zero). It might well be that also our many-body N-SR PT would survive photon losses. How the properties of the CDW and NE phases are affected by photon losses is a most interesting issue. At the one hand, the CDW and NE phases are highly entangled and decoherence could demolish it, but on the other hand one could also argue that topology could stabilize in particular the excitations. Resolving this is certainly interesting but it lies outside the scope of the present work. Nevertheless, let us here assume that the experiment can be performed on time scales shorter than the typical photon life-time κ −1 , which can be of the order of ms in the microwave regime. This is of course experimentally challenging, but not completely unrealistic.
We end by mentioning the idea of studying similar qdot generated photon tunneling in two dimensions, which should not be experimentally much more demanding. In higher dimensions the lattice geometry may play an important role -one can imagine a plethora of different settings like lattices supporting flat bands, confinement, Dirac cones [48], or other more exotic lattices [49]. matrices is Note that the Gell-Mann matrices can be written as products of matrices of the Spin 1 representation of SU (2). However, only by using the Gell-Mann matrices can (1) be written as a quadratic Hamiltonian.

Appendix B: Polaron ansatz
The goal of this appendix is to calculate the effective qutrit Hamiltonian (8), which is found by taking the expectation value with respect to the photonic variables only: where, as in the main text, we have defined | α = i |α i . It is useful to write the polaron unitary U γ in two ways: where we have defined the two-site operators and i p i+1 i even, Note that both P i−1,i and S i,i+1 are unitary, both involves operators from two sites and [P i−1,i , P i,i+1 ] = 0, [S i−1,i , S i,i+1 ] = 0. S i,i+1 only involves qutrit operators for one site, while a i − a † i P i−1,i only involves photonic operators from one site.
U γ acts like a displacement operator on the photonic variables, i.e. U γ a i U † γ = a i +γ a i − a † i , a i P i−1,i = a i +γP i−1,i , (B5) which allows us to calculate the contribution to H eff from the photonic part of H B immediately as i α|U γ a † i a i U † γ | α = i α 2 i + 2γα i P i−1,i + γ 2 P 2 i−1,i .
(B6) The transformation of the qutrit part of H B is more involved. Using the Hadamard lemma [50], the transformed qutrit operators can be expanded as power series of nested commutators as To compute coherent state expectation values, we use the formula (assuming α ∈ R): α|p n |α = 0 n odd Γ( n+1 2 )/ √ π n even (B8) Using (B8), we find where the functions of γ can be expanded as rapidly converging power series It remains to calculate the contribution from the interaction Hamiltonian, In the calculation, we will need to use the formula α|(a + a † )p n |α = in/ √ 2 α|p n−1 |α n odd 2α α|p n |α n even. (B12) Using (B12), we find Note that the denominators are the primorial numbers [51], so f g (γ) is the 'primorial version' of e −γ 2 . Furthermore, from (B8), we know that only coherent state expectation values of even powers of p n are nonzero. By expanding in a power series using the Hadamard lemma and setting all terms which are odd in p i−1 or p i to zero, we find 2gγ i Q i−1,i α|U γ P i−1,i U † γ | α = 2γf g (γ) i P 2 i−1,i .
(B15) The expression (B13) and the cross-term in (B6) are the only terms which couple the photonic mean-field and qutrit degrees of freedom. Both terms can be eliminated by choosing γ to be a root of the equation which has a solution γ ∼ −g/ω close to g/ω = 0. By this choice, the photonic and qutrit degrees of freedom are completely decoupled. Once decoupled, the photonic degrees of freedom can be minimized independently to α i = 0. In this sense, the polaron ansatz is a "dynamic displacement" of the photonic vacuum. The ansatz neglects any contribution to quantum fluctuations from the resonators around this vacuum state. We found that the photonic part of the effective qutrit Hamiltonian was minimized to zero, so the remaining part is a qutrit Hamiltonian: Assuming periodic boundary conditions, we find (the constant term is dropped): i−1 ) 2 + (λ (6) i−1 ) 2 + 2λ (6) i−1 λ (6) i