Effective many-body Hamiltonians of qubit-photon bound states

Quantum emitters (QEs) coupled to structured baths can localize multiple photons around them and form qubit-photon bound states. In the Markovian or weak coupling regime, the interaction of QEs through these single-photon bound states is known to lead to effective many-body QE Hamiltonians with tuneable but yet perturbative interactions. In this work we study the emergence of such models in the non-Markovian or strong coupling regime in different excitation subspaces. The effective models for the non-Markovian regime with up to three excitations are characterized using analytical methods, uncovering the existence of doublons or triplon states. Furthermore, we provide numerical results for systems with multiple excitations and demonstrate the emergence of polariton models with optically tuneable interactions, whose many-body ground state exhibits a superfluid-Mott insulator transition.


Introduction
Quantum systems coupled to a common environment experience interactions mediated by the bath excitations [1]. In quantum electrodynamics, this is the basic mechanism behind the forces between electrons, atoms, or molecules. Those interactions can be tuned if one controls the coupling to the bath, which opens up exciting possibilities in quantum information. One of the most prominent examples in this context is that of two-level quantum emitters (QEs) coupled to the free-space electromagnetic field, resulting in the well known dipole-dipole interactions between QEs [2,3]. Unfortunately, these interactions are generally accompanied by spontaneous emission, which limits their applications. The latter can be avoided if the density of modes at the transition frequency vanishes, since the spontaneous decay rate is proportional to that quantity. This occurs, for instance, if one embeds QEs in a cavity such that QEs are far-off resonance from the cavity modes [4,5].
Another way of cancelling spontaneous emission while obtaining exchange interactions among QEs is by endowing the bath with a periodic structure, which strongly influences the density of states [6]. In fact, band gaps where the density of states vanishes can emerge, so that spontaneous emission in that bath can be prevented, but still interactions between the emitters can be mediated by virtual processes via the common bath [7,8]. Experiments with atoms, quantum dots, or superconductors interacting with structured bath has renewed the interest in investigating these phenomena [9,10,11,12,13,14,15]. Other experimental scenarios involving cold atoms in optical lattices with state-dependent potentials are amenable to the same description, and thus the appearance of analogous phenomena have been predicted [16,17] and have recently been observed [18].
The theoretical description of the interactions mediated by a bath is relatively simple in the so-called Markovian limit [2,3,19]. There, it is possible to derive a master equation for the quantum emitters only, where the bath degrees of freedom are traced out. This effective description contains both Hamiltonian and dissipative Lindblad terms. The latter vanishes if the QEs transition frequency lies in the band gap, whereas the first one describes the interactions between the QEs mediated by the bath, where one of them is excited when another one is de-excited. The emergence of dipole-dipole interactions in this scenario, which opens the door to investigate spins models with long range interactions [20,21,22], can be attributed to the existence of a single-photon bound states [7]. An intuitive picture [21] is that the single-photon bound state acts as an off-resonant cavity mode that mediates interaction between the QEs. The strength and functional form of these interactions depend on the QE-bath coupling strength, detuning as well as the band-dispersion relation. Although, these interactions can be made relatively strong to overcome other dissipative mechanism, their predicted strength is ultimately limited by the Markovian conditions under which these effective description has been derived.
In this work, we study the effective many-body Hamiltonians emerging when QEs couple to structured baths beyond the Markovian limit, and investigate to which point the dipole-dipole description survives in this regime. Furthermore, we study the consequences on the effective QE interactions of the emergence of multi-photon bound states, which were recently predicted in the single QE regime [23,24,25,26,15], but which impact in the multi QE situation has not yet been fully considered. Our analysis allows us to uncover qualitatively different interaction Hamiltonians, in which multiphoton bound states (doublons/triplons) hop from QE to the other, and analytically characterize them up to three excitations. With more excitations, we numerically characterize the emergent polariton models using density matrix renormalization group (DMRG), and show that their interactions can be controlled optically through the QEbath interactions, allowing us to probe the quantum phase transition between superfluid and Mott insulator. This manuscript is organized as follows. In Sec. 2, we explain the model that will be used throughout the paper and review the results in the Markovian limit to have them as reference for the next Sections. In Sec. 3, we study the single excitation subspace, deriving an effective Hamiltonian to describe the dipole-dipole interaction for two and many QEs. In Secs. 4-6, we study the two-excitation subspace for both the two and many QE regimes. Since the phenomenology in this regime differs significantly from what is expected in a Markovian description, we use Sec. 4 to explain the analytical tools we use to characterize it, and describe qualitatively the main features that emerge in this subspace, namely, the scattering of two polariton modes and the hopping of doublon states. This emergent dynamics will be discussed in detail in Secs. 5 and 6, respectively. In Sec. 7, we analytically characterize the three excitation subspace, and then go to the many-excitations regime to numerically explore a superfluid to Mottinsulator phase transition appearing in the ground state of the systems. Finally, we summarize our results and conclude in Sec. 8.

Model
As shown schematically in Fig. 1(a), we consider N b > 1 QEs with two energy levels {|g j , |e j } (j = 1, . . . , N b ) coupled to a structured bath. The total Hamiltonian of the system reads (with = 1): where H A = ω e N b j=1 σ j ee describes an array of N b two-level QEs with transition frequency ω e . Periodic boundary conditions (PBCs) for the bath are used, so we can label its modes in terms of the quasi-momentum k. Here, a k (a † k ) is the annihilation (creation) operator of the bath mode with quasi-momentum k, ε k is the energy dispersion relation, N is the total number of bath modes, and Ω is the QE-bath coupling strength. We will assume that the bath has a single band of width W , although the results can be easily extended to other situations. To obtain analytical results, we will go into the continuum limit, where N → ∞ so that k becomes a continuous vector, and we can replace sums by integrals. Finally, σ j αβ = |α j β| are the spin operators for the j-th QE and n j its corresponding vector position.
In optical and microwave implementations, we require that Ω ω e , minε k , so the counter-rotating terms for the QE-bath coupling can be neglected. For convenience, we will work in the rotating frame at the frequency minε k , which amounts to taking (2) where ∆ = ω e − min ε k . A crucial feature of the Hamiltonian in Eq. (1) is that the number of excitations, defined as is conserved. This allows us to derive effective models separately in the subspaces with different numbers of excitations. In this work, we will concentrate on the few-body scattering and bound-state behaviors in the subspaces with N exc = 1, 2, 3, and the quantum phase transitions in the ground states of the subspaces with many excitations. Though most of the expressions we derive are valid for an arbitrary energy dispersion ε k (see Appendices), in the main text we focus on the results for a 1D tightbinding model, where the energy dispersion reads: with J being the hopping strength and k = 0, 2π/N, . . . , 2π(N − 1)/N .

QEs as hard-core bosons
For the calculations performed in this paper, it is convenient to describe QEs using hardcore bosons. In this representation, we replace σ j ge → b j , where b j is an annihilation operator fulfilling bosonic commutation relations, and we restrict the Hilbert space to the states with (b † j ) 2 = 0. In practical terms, this can be done by writing: and taking the U → ∞ limit, which forbids double occupation of the b j modes.
In the case of many QEs, we assume that they are equally spaced, which allows us to work in Fourier space by defining b j = p b p e ip·n j . The quadratic part of the that commute with each other (i.e., [H p , H p ] = 0), where the photon momentum k = p + K is given by the QE quasi-momentum p and the reciprocal momentum K of the sublattice, and z = N/N b is the number of bath modes per unit cell. The hardcore interaction part H hc becomes

Markovian limit
Let us here remind the results obtained in the Markovian limit when the QEs transition frequency lies within the band-gap, which means that ∆ < 0, and |∆|, W Ω. In that limit, one can eliminate the bath degrees of freedom using standard quantum optics techniques and obtain an effective Hamiltonian for the QEs [2,3]. For two QEs with relative position d = n 2 − n 1 , one obtains: where V dd (d) is the dipole-dipole interaction strength, which in the limit N → ∞, reads: for a general D-dimensional bath and σ + j = (σ − j ) † ≡ σ j eg = (σ j ge ) † . For many QEs, one obtains Note that the dipole-dipole interaction for many QEs is the same as for two QEs under the Markov approximation. For the 1D tight-binding model, one obtains [20] with the decay length ξ = ln −1 (|∆|/J). We will use these effective Hamiltonians as a baseline to compare with the results of the next sections. In particular, we will see what is the regime of validity, and how it has to be modified outside that regime. To do that, we will first study analytically up the three-excitation manifold, and finally perform DMRG calculations for the case with many excitations.

Single Excitation
In this section we study the physics of the single excitation subspace. This regime has been extensively studied in the literature (see, for instance, Ref. [26] and references therein). Here, we will review results that are relevant for the other sections, and also derive simple formulas for the emergent effective models. We divide this section in two parts: in Section 3.1 we study the situation when only two QEs are coupled to the bath, deriving an effective exchange interaction Hamiltonian valid in a particular region of the (∆/J, Ω/J) parameter regime that we will define. In Section 3.2, we study the situation with many QEs and derive an effective hopping model in the lowest band, which can be defined in all parameter regimes. In both regimes, when |∆| Ω we recover the effective spin models predicted for the Markovian limit in Eqs. 8-10. However, in the strong-coupling limit, Ω J, ∆, an effective spin model can still be derived to characterize the hopping of the strong hybridized QE and photon, i.e., a polariton, where the effective hopping strengths are dramatically enhanced compared to those in the Markovian regime.

Two QEs
The single-excitation eigenstates for the system with two QEs can be generally written The coefficients u 1,λ , u 2,λ , and f λ (k) for all the eigenstates (including bound states and scattering states) can be obtained by solving the Schrödinger equation H|Ψ 1λ = E 1λ |Ψ 1λ (see Appendix A). The probability weights of symmetric and anti-symmetric In general, the two lower eigenstates, |Ψ 1λ=± = β † ± |0 , represent the symmetric and anti-symmetric superpositions of two local single-excitation bound states around the QEs. However, as we will show below (and already derived in Ref. [26]), only the symmetric bound state survives in certain parameter regimes. In the regime where both bound states exist, the Hamiltonian can be projected into the subspace spanned by these two bound states, which gives rise to the following effective model for the low energy part of the spectrum: A basis transformation converts this Hamiltonian into a hopping model for the two localized Wannier modes β 1,2 = (β + ± β − )/ √ 2, where E 0 = (E 1+ + E 1− )/2 is the effective chemical potential and t eff = (E 1+ − E 1− )/2 is the effective hopping strength.
If the distance d between two neighboring QEs is sufficiently large, the two nearly degenerate ground states with E 1+ ∼E 1− describe the single-excitation bound states localized around two individual QEs. These two bound states begin to hybridize with each other as d decreases. When the distance d is much smaller than the localization length of the bound states, the strong mixing of the bound states induces the energy level splitting 2t eff around E 0 . If the splitting is large enough, the anti-symmetric bound state energy merges in the continuum and the state is not bound anymore. For example, for the 1D bath with the tight-binding dispersion relation, the symmetric bound state always exists, but the anti-symmetric bound state can only be found in the regime ∆ < Ω 2 d/2 [26]. Obviously, the effective hopping model of Eq. 14 can not be defined when only the symmetric bound state exists. In this paper, we mainly focus on the regime where both bound states exist, such that the low-energy physics of the single excitation is described by the effective Hamiltonian (14), where the relevant parameter is the hopping strength t eff .
As shown in Fig. 2, in the Markovian |∆| Ω and strong coupling Ω J, ∆ regimes, a parameter regime that we denote as the arc region, the photon is strongly localized around QEs. Thus, the effective chemical potential E 0 ∼ E 1B tends to the energy E 1B of the bound state around a single QE. The effective hopping strength when t eff E 0 can be approximated by t eff ∼ − Z 1B Ω 2 e −d/ξ / E 1B (E 1B − 4J) (see Appendix A), where the single-particle weight and the decay length are In the Markovian limit, the effective hopping strength t eff ∼ − J eff |∆/J| 1−d decays exponentially with the distance d, where J eff = JΩ 2 /∆ 2 , such that it reproduces the result of Eq. (11). In the strong coupling limit, the effective hopping strength t eff = −J eff (Ω/J) 1−d also decays exponentially, however, J eff = J/2, which means that the coupling strengths are significantly enhanced due to the strong hybridization between QE excitation and photon.
In Fig. 2, the effective hopping strength t eff is shown in the ∆ − Ω plane for two QEs with d = 1, 2 coupled to the 1D tight-binding bath. For ∆ < 0, |t eff | increases monotonically and saturates to J/2 in the strong coupling limit. For ∆ > 0, |t eff | increases to the maximal value at Ω max slightly above the boundary Ω = 2∆/d, and  decreases to J/2 in the strong coupling limit Ω J, ∆. The Rabi frequency Ω max maximizing |t eff | is highlighted by the dashed red lines in Fig. 2.

QE array
Let us now consider the situation where we have now many QEs coupled to every z bath lattice sites. By imposing PBCs for the QE array, the excitations will have z + 1 bands for the QE propagation, although here we focus on the lowest energy band. Using the intuition from the previous Section (see also [26]), we expect that the bath mediates QE interactions giving rise to a polariton propagating in the QE lattice. Let us now explain how to characterize this emergent behaviour.
Since we are restricting in this Section to the single excitation subspace, the hardcore interaction Hamiltonian plays no role, and therefore the single-particle modes with different quasi-momenta p in the first Brillouin zone are decoupled. Thus, the QE-bath Hamiltonian is quadratic and can be diagonalized, H 0 = pλ E 1λ (p)β † pλ β pλ , by the annihilation (creation) operator β pλ (β † pλ ) of the single polariton with momentum p and dispersion relation E 1λ (p) in the λ-th band, where λ ∈ [0, 1, 2, · · · , z] labels the different energy bands. In the polariton mode β † pλ |0 , the probability weight of QE mode of quasimomenta p being in the excited state is given by: The dynamics in the lowest band (i.e., λ = 0) is given by: where P is the projector into the states β † p |0 ≡β † pλ=0 |0 in the lowest band λ = 0, and E 1 (p)≡E 1λ=0 (p) its energy dispersion relation. The Wannier modes can be obtained via Nearest neighbours Nearest neighbours Next-Nearest neighbours Next-Nearest neighbours -0.050 0 0.050 Fourier transform as β j = p β p e ip·n j / √ N b , which describes the local single-excitation bound state around the j-th QE. In terms of these local Wannier modes, the effective Hamiltonian in the coordinate space can be written as with hopping strengths In the arc region, the dominant hopping constant t j ∼ −t 1 δ j1 can be deduced from the first order degenerate perturbation theory, where t 1 = J(1 − Z 1B ). In the Markovian and large Rabi coupling limit, t 1 = JΩ 2 /∆ 2 and J/2 agrees with J eff in the two-QE case.
In the large d limit, the vanishing t j →0 indicates that the local Wannier modes reduce to the single-excitation bound states localized around different individual QEs.
As the lattice spacing decreases, t j becomes finite and the local Wannier modes hybridize with each other to form the dispersive polariton band. Here, we note that since PBC for the QE array is applied, the translational symmetry ensures that the number of states in the lowest band is the same as the QE number due to the Bloch theorem. However, if the QE array has finite size and we have open boundary condition, like it occurs for two QEs coupled to the bath, some polariton modes may vanish in certain parameter regimes due to the boundary effect, as we showed in previous section.
In Fig. 3, the nearest neighbor (NN) and next-nearest neighbor (NNN) hopping strengths t l=1,2 in the ∆−Ω plane are shown for the 1D tight-binding dispersion relation. In the arc region of the ∆ − Ω plane, the Wannier mode in the lowest band is the single excitation bound state strongly localized around each individual QE, and the overlap of two Wannier modes exponentially decays with the localization length ξ≤1. As a result, the NN hopping strength t 1 ∼t eff can be determined by the hopping strength in the two-QE case, and the long range hopping strength t l>1 ∼0. In the Markovian limit |∆| Ω, the single particle weight Z 1λ=0 (p)≡Z 1 (p)∼1 shows that the polariton state in the lowest band is mostly composed of QE excitations, and the band is only slightly deformed from a completely flat one. In the strong coupling limit Ω J, ∆, the reduced Z 1 (p)∼1/2 exhibits the strong hybridization of QE excitation and bath photon modes, which gives rise to the finite width∼ 2J of the polariton band.
Non-Markovian effects emerge in the intermediate regime, where the long range hopping strengths t l>1 in general do not vanish, as shown in Figs. 3c and 3d. In the small Ω/J limit, the hopping strengths saturate to t 1 /J = −0.424 and t 2 /J = 0.085 for d = 2 when the detuning ∆/J > 2. This saturation can be understood using the single-excitation band structure in the limit of Ω→0. If 0 < ∆ < ∆ c ≡ε k=π/d and Ω/J is small, the lowest band states are mostly composed of photons in the bath and the lowest band has a width ∼ ∆. If ∆ > ∆ c , the width of the lowest band saturates to ∆ c at small Ω, and t l →2 Summing up, we have derived and characterized the effective hopping model emerging in the single-excitation subspace when many QE are coupled periodically to the bath modes. In particular, we have shown it reduces to an effective tight-binding model in the arc region of the ∆ − Ω plane. In the Markovian regime |∆| Ω [strong coupling limit Ω J, ∆], the effective model describes the propagation of the bare QE excitation [the hybridized QE-photon polariton excitation with the single particle weight 1/2]. In the latter, the effective NN hopping is significantly enhanced compared to that in the Markovian regime.

General features of two excitation subspace
The dynamics with more than one QE excitation can be characterized by the effective spin model of Eq. 10 only in the Markovian limit. To the best of our knowledge, the use of that effective model beyond the Markovian regime is not justified, and how to characterize the dynamics in the whole ∆-Ω plane is not clear yet. In this section, we pedagogically review the tools on how to treat both the two and many QE situation in the two-excitation sector, and explain qualitatively the emergent features. In particular, we will see that two different types of excitations appear in the two-excitation subspace, namely, the scattering states of two single polaritons and the doublon states formed by two bound polaritons, as we show schematically in Fig. 1(b). Then, we discuss them in detail in Sections 5 and 6, respectively.

Two QEs
When only two QEs are coupled to the bath, the general two-excitation eigenstate has the form The eigenvalue E 2λ of |Ψ 2λ , the probability weight Z 2λ , and the wavefunctions ϕ 1λ (k), ϕ 2λ (k), ϕ 2λ (k, k ) can be obtained from the analytical structure of the Green function G 2 (ω) = dtG 2 (t)e iωt in the frequency domain [25], where The ground state energy corresponds to the smallest isolated pole of G 2 (ω) in the first Riemann surface. As we shall discuss in Sec. 5A, the ground state describes two repulsive polaritons localized around two different QEs, where the localization behavior is analyzed by the ground state wavefunctions in Sec. 5A. We also construct a variational wavefunction and an effective Hamiltonian to describe the low energy physics in the twoexcitation subspace.
Apart form the pole corresponding to the ground state, one can also find two additional isolated poles corresponding to higher excited states in the certain parameter regime (see Sec. 6.1 and Appendix B). These two states can be illustrated in the following way: for a single QE with two photons, it has been demonstrated that the two photons can be both localized around the QE and form a two-excitation bound state [23,25,26], which is referred to as the doublon state, schematically depicted in Fig. 1(b). Here, the two higher excited states represent the symmetric and anti-symmetric superpositions of doublon states around different QEs. The properties of the doublon state is studied in Sec. 6.1, where an effective hopping model for the doublon is derived.

QE array
Here, two new types of states appear compared to the single-excitation regime, i.e., the scattering state of two polaritons and the propagating doublon state. For a system with two excitations in the QE array, the eigenstate at quasi-momentum q has the general form where the momenta p and q are restricted to the first Brillouin zone of the QE's sublattice. One can introduce the four-point Green function in the time domain, where the operators α 1 ∈ (b q/2+p , b q/2+p , a q/2+p+K ) and α 2 ∈ (b q/2−p , a q/2−p+K , a q/2−p+K ). Its Fourier transform G 2 (q, ω) = dtG 2 (q, t)e iωt gives the dispersion relation E 2λ (q) of the state |Ψ 2λ (q) in the band λ, the wavefunctions f b (p), f ba (p, K), and f a (p, K, K ). The energy band structure in the two-excitation subspace can be determined from the analytical structure of G 2 (q, ω) (see Appendix C).
In Fig. 4, we show the three lowest bands of the two excitation subspace for different sets of system parameters with 1D tight-binding dispersion relation. The bands corresponding to the scattering of two polaritons are identified by the continuum in Fig. 4 (in gray shading), where the lowest band belongs to them. The scattering properties will be studied in detail in Sec. 5, where an effective two-body Hamiltonian is established to describe the low energy dynamics of two polaritons. Between the scattering continuum, the single isolated band appears in the midgap, which corresponds to the propagating doublon, which can be described with an effective hopping model as we show in Section 6.

Two-QE ground state
In this subsection, we study the ground state of the two excitation subspace when two QEs coupled to the photonic bath. In particular, i) we compute the exact ground state wavefunction; ii) we construct a variational ansätz to reveal the properties of the ground state; and iii) we derive an effective Hamiltonian to describe the low-energy dynamics in the so-called arc region.
The ground energy E G is determined by the position of the isolated pole in the first Riemann surface of G 2 (ω), more precisely, E G is the smallest solution of the equation where E 1λ and Z σ 1λ are defined in Sec. 3.1. The ground state configuration is visualized by the wavefunctions ϕ jG (k) and ϕ 2G (k, k ) (their analytical expressions are given in Appendix B). In the left and right panels of This repulsive interaction between the two excitations can also be identified by the energy difference δE G ≡ E G − 2E 1B the ground state energy E G and the energy of the excitations localized around two individual QEs far apart from each other. In Figs. 6a  and 6c, we show δE G in the ∆ − Ω plane for d = 1 and 2, respectively. In the arc region, the energy difference δE G ∼ 0. This is because the two single-excitation bound states are strongly localized around different QEs, which suppresses their interaction. In the regime ∆ > 0 and Ω/J 1, the non-interacting photon excitations dominate in the bound state, and δE G ∼0. In the vicinity of the boundary Ω = 2∆/d, where the anti-symmetric bound state β † − |0 vanishes, there is a still a considerable probability for the QEs to be excited, while having a large overlap between the two single-excitation bound states. This induces the interaction between the two Wannier modes such that δE G < 0, as shown by the dark (blue) regions in Figs. 6a and 6c.
The low energy dynamics of the single excitation subspace is governed by the effective Hamiltonian H (1) eff projected into the subspace spanned by the polariton modes β † ± |0 . Therefore, we expect that the two-excitation ground state describing the two interacting polaritons is mostly composed of the low energy excitations β †2 ± |0 / √ 2. In Figs. 6c and 6d, we show the probability p = σ=± | 0|β 2 σ |Ψ 2G | 2 /2 to find two (anti-) symmetric excitations β †2 ± |0 / √ 2 in the exact ground state |Ψ 2G . Here, we note that in the regime ∆ > Ω 2 d/2 the anti-symmetric mode β † − |0 vanishes, and the probability is defined as p = 0|β 2 + |Ψ 2G 2 /2, namely, only the symmetric mode is taken into account. The large probability p > 0.85 even in the non-Markovian regime shows that the excitations β †2 ± |0 / √ 2 dominate the ground state, which indicates that the low energy dynamics in the single and two-excitation subspaces can be described by some interacting Hamiltonian for the polariton modes β † ± |0 . A remarkable feature of the two excitation ground state is that the first order correlation matrix M ij = α † i α j , for operators α j = b 1,2 , a k , and the two-photon wavefunction ϕ 2G (k, k ) at most have two dominating singular values in the whole ∆ − Ω plane. This fact inspires us to construct a variational ansätz |Ψ 2v = N −1/2 2 γ †2 + − γ †2 − |0 /2 for the ground state by two deformed single-excitation modes defined by: where v ± are real numbers, ϕ ± (k) are orthonormal wavefunctions, and the normalization factor N 2 = σ (1 + v 2 σ ) 2 /2. Under the normalization constraint k |ϕ σ (k)| 2 = 1, the minimization of the variational energy E = Ψ 2v |H|Ψ 2v with respect to ϕ σ (k) indicates that the variational function should have the form where e σ is a variational parameter, η kσ = (1 + σe −ik·d )/ √ 2, and the normalization factor In terms of the variational parameters v σ and e σ , the variational ground state energy can be written as where One can minimize the ground state energy with respect to v σ and e σ to find the optimal variational state |Ψ v2 . Figures 7a and 7b show the overlap p v between |Ψ v2 and the exact ground state for 1D tight-binding dispersion relation at d = 1, 2. The fact that p v > 0.98 in all parameter regimes underscores the accuracy of our variational state |Ψ v2 .
If we define the normalized symmetric and anti-symmetric modes γ N,± = γ ± / 1 + v 2 ± of the variational state, we can calculate the overlap between these modes and the single-particle ones, β † ± |0 , given by p ± = 0|β ± γ † N,± |0 . The symmetric mode γ † N,+ |0 only deviates slightly from the bare single-particle mode β † + |0 , i.e., p + > 0.99, in the whole parameter space, including the non-Markovian regimes (not shown). On the contrary, the mode γ † N,− |0 is very different from the anti-symmetric single-particle mode β † − |0 in the vicinity of the boundary ∆ = Ω 2 d/2, as shown in Figs. 7c and 7d for d = 1, 2. Furthermore, in the regime ∆ > Ω 2 d/2, the anti-symmetric bound state vanishes and the γ †2 N,+ |0 / √ 2 dominates the ground state. This is why the probability p, which measures the overlap between the exact and the variational wavefunction, is still very large in this non-Markovian regime.
In the so-called arc region, p ± ∼1 and v + = v − indicates that the ground state , where the small component of the double occupation states β †2 j=1,2 |0 in |Ψ 2v shows the hard-core nature of β † j |0 . As a result, we can construct the effective spin model for the single and two excitations subspaces, where J z = E G −2E 0 and the spin operators σ − j=1,2 denote the local Wannier modesβ 1,2 . In the Markovian regime the spin operators σ − j become the bare QE transition operator σ − 1,2 , while in the strong coupling regimeσ − j denotes the annihilation operator of the strongly hybridized polariton mode.

Many QEs
For two QEs, the two-excitation ground state can be described by two interacting single-excitations β †2 ± |0 / √ 2. In this subsection, we investigate the scattering of two excitations in the periodic QE array coupled to the photonic bath. As discussed in Sec. 4, the two-excitation spectrum is obtained by the four-point connected Green function (see Appendix C), which displays a band structure due to the PBC imposed. To study the low energy scattering of two individual polaritons, we derive the energy-dependent two-body interaction in the whole parameter plane using the Feshbach treatment [27]. We demonstrate that the effective interaction strength between two polaritons can be tuned by the detuning ∆ and the Rabi frequency Ω. In particular, the hardcore nature of two excitations justifies the validity of the spin model in the arc region.
Let us first derive the effective Hamiltonian describing the lowest band of twoexcitation spectrum using the Feshbach treatment. For the incident state β † p β † q−p |0 of two polaritons with total momentum q and energy E, the scattering state |Ψ 2sc (q) = |Ψ 2P (q) + |Ψ 2Q (q) can be written as the superposition of |Ψ 2P (q) = P |Ψ 2sc (q) and |Ψ 2Q (q) = Q|Ψ 2sc (q) , where P and Q = 1 − P are the projectors into the lowest and higher bands, i.e., the "open" and "closed" channels in the Feshbach resonance.
By eliminating the higher energy band in the Schrödinger equation, the component |Ψ 2P (q) in the lowest channel obeys the secular equation The effective Hamiltonian H eff (E) = H eff + P H int (E)P describes the scattering of two excitations in the lowest band by where the effective two-body interaction depends on the incident energy E, and the summation does not include the contribution from the lowest band. The component in the closed channel follows from the Schrödinger equation, which is a bound state of two excitations in the higher energy bands. Thus, the closed channel component will not contribute to the scattering wavefunction in the asymptotic limit. We characterize the scattering process of two low-energy polaritons at the band bottom by the effective interaction U eff = Z 2 1 (0)U eff (0, E 0 ), which we plot in Fig. 8 for d = 1, 2 in the ∆ − Ω plane. In the ∆ > 0 regime, the polariton mode β † p |0 is mostly composed of bath photons for small Ω, so the interaction strength U eff is extremely softened. As Ω increases, the weight of QE excitation in the mode β † p |0 becomes larger, and the effective interaction U eff increases accordingly. This means a wide range of U eff can be achieved by tuning the Rabi frequency Ω and ∆.
In the Markovian regime, the single-particle weight Z 1 ∼1 and a large gap (compared with the bandwidth) separates the lowest band and the higher energy bands, which result in the divergent interaction U eff and the small component ||Ψ 2Q (q) | 2 in the closed channel. Thus, the low-energy dynamics can be described by the effective spin model H eff = jj t j−j σ + jσ − j projected in the lowest band, where the spin operator σ − j ∼ σ − j can be approximated by the bared QE transition operator. In the strong coupling regime, Z 1 ∼ 1/2 and the two-body interaction U eff between strongly hybridized polaritons is divergent. Therefore, the spin model is also valid in this region, andσ − j ∼ β j denotes the local Wannier mode with the strong on-site interaction.

Two-QE doublon state
As it occurred for the single-photon bound states, only for a given region of the ∆ − Ω parameter space, two bound states |Ψ 2± with higher energies E 2± than E G appear in the spectrum. The exact condition for the existence of two bound states is determined via the Green function approach in Appendix B. The localization behavior of two bound states can be understood by inspecting their wavefunctions ϕ j± (n) and ϕ 2± (n, m), which we plot in Fig.. 9, for d = 2, Ω/J = 2, and ∆/J = −1. There, we observe that the two bound states are the symmetric and anti-symmetric superpositions of the two-photon bound states [25] localized around different QEs.
When both doublons exist an effective hopping model for them can be constructed, which reads: In Fig. 10 we plot the hopping strength t D in the ∆ − Ω plane for d = 1, 2. For fixed distance d and small Rabi coupling, the large localization length of the doublon Wannier modes gives rise to a strong hybridization between them and a large energy level splitting 2t D . As a result, the two doublon states vanish by merging into the continuum of scattering band. As the Rabi frequency or the distance d increases, the localization length is shorter than d, and the overlap between two doublon Wannier modes becomes smaller. Therefore, the energy level splitting 2t D is reduced, so one can find two doublon states and define the effective hopping strength t D . This intuitive picture also explains why the regime where these bound states vanish shrinks as the distance d increases, as shown in Fig. 10(b).

Many QEs
For many QEs, we already saw in Fig. 4, that a doublon band (E 2B (q), in red triangles) appears in the midgap of the scattering bands, whose dispersion depends strongly on the parameter regime. For example, for ∆/J = 1 and Ω/J = 1 in Fig. 4f, the doublon band has a visible curvature which implies that the doublons have a large hopping strength, whereas for ∆/J decreases to −1 in Fig. 4d, the doublon band is very flat which tell us that the doublons are mostly localized. To gain more intuition of the features of this doublon band, we plot the coordinate space structure of the doublon in Fig. 11 using the Fourier transforms f b (r), f ba (r, m), and f a (r, n, m) of the functions f b (p), f ba (p, K), and f a (p, K, K ) defined in Eq. 21.
Here, f b (r) is the amplitude of having two excited QEs separated by a distance r, f ba (r, m) is the amplitude of finding one photon at the position m and one excited QE separated from the photon by a distance r, and f a (r, n, m) is the amplitude of detecting two photons at positions n and m. In Fig. 11, we observe that for d = 2, ∆/J = 0 and Ω/J = 2, the square norms |f b (r)| 2 , |f ab (r, m)| 2 , and |f a (r, n, m)| 2 of the state with q = 0 in the first two isolated bands display that the two polaritons attract each other to form a propagating doublon. One can also see that the two polaritons are tightly bound in the lowest doublon band but loosely bound in the higher doublon band.
From the doublon state, |Ψ 2B (q) , and its dispersion relation E 2B (q), we can construct the effective doublon Hamiltonian as follows: which describes the hopping of the doublon mode d † n |0 . In the single doublon space, d n = |0 Ψ 2B (n)| is defined by the Wannier state |Ψ 2B (n) = q e −iq·n |Ψ 2B (q) / √ N b localized around QE at the position n. The effective hopping strength of the doublon is where the NN hopping strength t D n−m=1 is t D (the effective hopping strength in the two-QE case). We note that in the dilute gas limit, where the density of doublons 1, the effective hopping model H doublon is still valid.

Many Excitations
In this Section, we study the spectrum for the N -excitation subspaces (with N > 2), focusing on the QE-array situation. In the first subsection, we apply the Greenfunction approach to characterize analytically the emergent dynamics in the N = 3 subspace, finding a continuum band that describe two types of scattering processes, i.e., the scattering of three individual polaritons and that between one polariton and one doublon. There exists also an isolated band in which the three excitations form a bound state (referred to as the triplon state) and co-propagate along the QE array. Even though the Green-Function method can be in principle extended for larger excitation number, it becomes very challenging due to the emergence of many topologically inequivalent Feynman diagrams. Thus, in the second subsection we use the intuition developed by in the previous sections to numerically characterize the ground state properties of the system using DMRG [28,29]. In particular, we will be able to show that one can go from a regime where the system behaves as a Mott-insulator to a superfluid behaviour, just by tuning the system parameters.

Three excitations
In this subsection, we study the properties of three excitations using the Green function method (see Appendix D). In Fig. 12, we show the band structure for three excitations for a situation with ∆ = 0 and lattice spacing d = 1. The lowest continuum band describes the scattering between three individual polaritons, while the second continuum band is composed of the scattering states between one polariton and one doublon. In the three-polariton scattering band, the two-body interaction between polaritons can be tuned by Ω and ∆, as shown in Sec. 5. With relatively large Rabi frequencies Ω/J = 5, 6, the midgap opens between the two lowest scattering bands, and a triplon band (denoted by the red triangles) with dispersion relation E 3B (q) appears. In the triplon band, the explicit analytic form of the triplon state |Ψ 3B (q) reads: In Figs. 13(a-d), we plot the triplon wavefunctions in real coordinate space, by Fourier transforming the functions , respectively. As we observe in the figure, the wavefunction has a bound-state behaviour, where the QE excitations and photons are localized around each other. The effective hopping model for the triplon has the same form as Eq. (34). In the single triplon subspace, T n = |0 Ψ 3B (n)| is defined by the three-excitation bound state |Ψ 3B (n) = q e −iq·n |Ψ 3B (q) / √ N b localized around QE at the position n, and the effective hopping strength is

Superfluid to Mott insulator transition
Let us finally consider the situation of many excitations in the QE array situation.
With the intuition we developed with the results of the previous Sections, we know that when the coupling is very strong, Ω J, or when we are in the deep Markov regime, |∆| J, Ω, we expect to have very localized bound states around the QEs. However, as one deviates from that conditions, the localization length of the bound states grows leading to strong hybridizing effects between the localized excitations. In this Section, we explore whether this localization length change leads to a superfluid-Mott insulating phase transition in the ground state of the system.
To obtain a detailed quantitative understanding, we study the system numerically using the DMRG method [28,29]. The DMRG algorithm is a variational method within the class of matrix product states and its nature imposes two constraints on the numerical studies. First, one should adopt open boundary conditions as this is more suitable for DMRG. The Bloch bands and momentum introduced for periodic systems cannot be defined for open systems but the physical properties should be the same if the system is sufficiently large. Secondly, the number of excitations on the bath sites should have an upper bound, that we denote as C, because the computational cost of DMRG is related to the Hilbert space dimensions of the lattice sites. C should be large enough such that the numerical results reflect the true physics. The ground states have been computed using DMRG in various cases and we find that C = 5 is sufficient because increasing it to 6 does not change the results significantly. The total number of excitations N exc is a conserved quantity so different N exc sectors can be studied separately. It is possible to access both the superfluid and the Mott insulating phases only if N exc is equal to the number of impurity sites N imp . We label the unit cells using Roman letters i, j etc. and the lattice sites (both the impurity and bath) within a unit cell using Greek letters α, β etc. The creation/annihilation operators for both all lattice sites can be expressed on a equal footing as a † iα /a iα . The first thing we would like to confirm is the existence of two phases in the system. The diagnostic tool we use is the correlation function matrix F iα,jβ = a † iα a jβ with iα (jβ) interpreted as the row (column) index. Fig. 14 shows the eigenvalues of F iα,jβ in the system with N imp = N exc = 80 at ∆ = 0.0, 0.2 and various different Ω. The single-excitation bound states in the large Ω regime are very localized so we expect to see Mott insulator behaviour, where F iα,jβ has multiple eigenvalues of similar magnitudes corresponding to the multiple modes occupied by the excitations. When the spatial extent of the singleexcitation bound states increases, the system transits to the superfluid phase where most excitations occupy the same mode so F iα,jβ has only one dominant eigenvalue. The low-energy effective theory for the superfluid phase is a Luttinger liquid theory, which predicts that F iα,jα ∼ |i − j| f to the first order [30]. Fig. 14 shows two examples of least where c is the central charge of the Luttinger liquid, g is a constant, and F is a nonuniversal oscillating term [31]. Fig. 15 shows two examples of least square fitting of S in the system with N imp = N exc = 80, where we choose 10≤L A ≤70 and discard those close to 0 or L to avoid edge effect. It turns out that the oscillating term is negligible and c = 1.0157 (c = 1.0595) for ∆ = 0.0J, Ω = 0.15J (∆ = 0.2J, Ω = 0.40J). This suggests that the superfluid state is a one-component Luttinger liquid with c = 1. In the Mott insulating phase at larger Ω, S is almost constant in the bulk of the system as one expects for a 1D gapped phase.

Conclusions
To sum up, we have studied the emergent dynamics of many QEs interacting with structured photonic reservoirs in the non-Markovian and many excitation regimes. In the two-and three-excitation subspaces, we provide analytical formulas for both the energies and wavefunctions of the relevant states governing the dynamics for arbitrary bath dimension and energy dispersion. We apply these formulas to study the case of a nearest-neighbour tight-binding one-dimensional bath uncovering several phenomena which are oblivious in perturbative descriptions. First, we show the emergence of effective hopping models in parameter regimes far from the Markovian ones, with the advantage of having stronger dipole-dipole couplings as compared to the perturbative regimes. Second, we also predict the emergence of new hopping models in the excited part of the spectrum, in which doublon/triplon states hop between the different QEs coupled to the bath. Finally, we numerically characterize the ground state in the many excitation sector of these quantum optical models, and show how the ground state can undergo an optically driven Mott-superfluid phase transition controlled by the localization length of the bound states. An interesting research direction is to apply the theoretical toolbox developed in the manuscript to study higher dimensional structured baths [32,33].

Appendix A. Single excitation
In this Appendix, we solve the Schrödinger equation to study the single excitation bound states of two-QEs with distance d = |d| and the QE array with the lattice spacing z = N/N b , where d = n 2 − n 1 is the vector connecting two quantum emitters (QEs). Without loss of generality, we choose n 1 = 0 and n 2 = d for two QEs. The parameters u 1,λ , u 2,λ , and the wavefunction f λ (k) are determined by the Schrödinger equation as By solving Eq. (A.1), we obtain the following equation to determine the bound state energies: and the vectors u ± = (u 1,± , u 2,± ) T , where the Green function is defined by the Pauli matrix σ x and the self-energies For the (anti-) symmetric bound state with energy E 1± , the parameters u 1,± = ±u 2,± ≡u ± / √ 2, and the wavefunction is determined by the normalization condition. The bound state energies E 1± and the parameters u 2 ± are the poles and the corresponding residues of the Green function As we show in the main text, there are certain parameter regimes in which the antisymmetric state merges into the continuum and only a single bound state exists. However, when both bound states exist a low-energy effective Hamiltonian can be written: obtained by projecting onto the subspace with the symmetric and anti-symmetric bound states.
In the Markovian limit |∆| Ω (∆ < 0), the single excitation bound state is extremely localized, such that the effective hopping |t eff | |∆|. Thus, the single particle bound state energies E 1± ∼E 1B ±t eff can be expanded around the single-excitation bound state energy in the presence of one QE. The secular equation determines the hopping strength By expanding the unitary evolution operator in Eq. (B.1), the Fourier transform G 2 (ω) = [G 2 (ω)] 0 + [G 2 (ω)] c can be written as the free propagation part and the connected part where · · · c denotes the connected Green function on the vacuum state.
In the interaction picture, the Green function G 3 (t) becomes G 3 (t) = −i 0|T α 3,I (t)α 2,I (t)α 1, Using the Dyson expansion, the connected part [G 3 (ω)] c reads × T 3 (q, p 1 , p 1 ; ω, ω 1 , ω 1 ) where P denotes the permutation of α 1,2,3 , and the three-excitation T -matrix satisfies the equation × G b (k, ω k )T (q − k, ω − ω k )T 3 (q, k, p 1 ; ω, ω k , ω 1 ). (D.5) The branch cuts of T 3 correspond to the scattering bands describing both the scattering of three individual polaritons and that between a single doublon and one polariton. The pole of T 3 determines the triplon band, where three polaritons form the bound state and co-propagate on the lattice.
The residues in the vicinity of the pole E 3B (q) leads to the wavefunctions