Bound states in ultrastrong waveguide QED

We discuss the properties of bound states in finite-bandwidth waveguide QED beyond the Rotating Wave Approximation or excitation number conserving light-matter coupling models. Therefore, we extend the \emph{standard} calculations to a broader range of light-matter strengths, in particular, in the so-called ultrastrong coupling regime. We do this using the Polaron technique. Our main results are as follows. We compute the spontaneous emission rate, which is renormalized as compared to the Fermi Golden Rule formula. We generalise the existence criteria for bound states, their properties and their role in the qubits thermalization. We discuss effective spin-spin interactions through both vacuum fluctuations and bound states. Finally, we sketch a perfect state-transfer protocol among distant emitters.


I. INTRODUCTION
Photons are weakly coupled to matter, so they rarely interact, making them perfect information carriers. Weak coupling constitutes, however, a double-edged sword, as it also hinders the readout process when the time comes to access the information being carried. The trade off is optimized in waveguide QED, where the photons are confined in one dimensional waveguides to enhance the light-matter coupling [1,2]. So far, different experimental platforms have been used to implement this coupling between quantum emitters (typically two level systems or qubits) and a one dimensional quantized electromagnetic field. Examples are superconducting circuits [3][4][5], optical waveguides [6] among others [7,8]. Waveguide QED can serve to control light-matter emission, to induce photon-photon interactions or to route the photons in quantum networks. Besides, by engineering the guides, more exotic interfaces can be implemented for quantum simulation [9], topological photonics [10], chirality [11,12] or quantum computing [13]. Consequently, waveguide QED may be a quantum technological solution.
Trying to optimize the light-matter coupling, several experiments have reached the so-called ultrastrong coupling regime (USC) between light and a single quantum emitter, both in cavity [14,15] and waveguide QED [16][17][18]. The USC is the regime where higher order processes, than the creation (annihilation) of one photon by annihilating (creating) one matter excitation play a role. Two main phenomena are paradigmatic of USC. The Rotating Wave Approximation (RWA) for the interaction breaks down and the atomic bare parameters get renormalized, either the Bloch-Siegert shift [19] in cavity QED or the renormalization due to the coupling to the continuum electromagnetic (EM) field in waveguide QED [20]. Besides, the ground state becomes nontrivial [21]. This has interesting consequences. Some of them are the possibility of transforming virtual onto real photons by perturbing the ground state [22][23][24][25], doing nonlinear optics with zero photons [26]. Further phenomenology in cavity QED can be found in recent reviews [27,28]. In this work we are interested in the USC regime in waveguide QED. Apart from the qubit frequency renormalization, there exist the localization-delocalization transition [29,30], particle production [31], non-linear optics at the single photon limit [32,33] or vacuum light emission [34].
In this work, we discuss the physics of bound states in the USC regime of waveguide QED. We focus on the lowest energy ones, discussing their existence and role in the spontaneous emission and thermalization. We also discuss the effective spin-spin models emerging when several emitters are ultrastrongly coupled to the EM field and envision protocols for perfect state transfer between distant atoms. To do this, we face a technical difficulty. The light-matter coupling is modelled via spin(s)-boson type Hamiltonians, a paradigmatic example of a non exactly solvable model [65]. Different techniques are available in the literature to deal with it. Matrix-product states (MPS) [29,32,33] , density matrix renormalization group (DMRG) [66] or path integral approaches [67,68], comprise the toolbox of numerical techniques. Analytical treatments are also used. They are based on different varational anstatzs: Polaron-like [30,34,[69][70][71][72] or Gaussian ones [73]. In this manuscript we will use a Polarontype approach that has been shown to be accurate in a wide range of parameters, including couplings well inside the USC.
The rest of the manuscript is organized as follows. In the next section, Sect. II, we will introduce the system, its model and the Polaron picture. In Section III, we treat the single emitter case. We discuss the ground state properties and the lowest bound state, discussing its existence conditions, energy, and localization length. Section IV develops the multiqubit case with emphasis in the effective tight-binding model and in protocols for perfect state transfer. We finish with some conclusions. Several technical issues are sent to the appendices. Finally, the link to the python codes used in the numerical calculations is given in App. D.

A. Model
In this manuscript we study the system sketched in Fig. 1(a). Several qubits are coupled to a cavity array forming the photonic medium. In the dipole gauge [74] and assuming that each qubit is coupled to a single cavity, the model is ( = 1 is set through the paper) Here, N q is the total number of qubits with level splitting ∆ (let us assume that all the atoms are identical). N is the number of sites; we will consider the thermodynamic limit N → ∞ in our analytical treatment. x j is the site to which the jth-qubit is coupled. Operators b † n and b n correspond to the bosonic creation and annihilation operators at site n, and σ z and σ x are the z and x Pauli matrices. To avoid extra parameters, we will consider that the qubit-resonator coupling is the same for all the qubits, g. The photonic medium (second and third term in Eq. (1)), which is a cavity array, is diagonalized introducing the bosonic operators in momentum space, which are the Fourier transform of their spatial counterparts: The dispersion relation, sketched in Fig. 1(b), is and the coupling per spin in momentum space is given by The dispersion relation is a finite band of width 4λ centred around ω 0 . Besides, the coupling constant is independent of the photonic mode and proportional to g, which justifies why, through out this work, we refer to both g and c k indistinctly as the coupling constant. Finally, it is convenient to define the spectral density, plotted in Fig. 1(c), conveniently rewritten in terms of the density of states [75] It is clear now that the state 0 1 , 0 2 , . . . , 0 Nq ; 0 with σ z j |0 j = −|0 j and b k |0 = 0 is the (trivial) ground state (GS) of the system and that the Hamiltonian preserves the number of excitations N , Consequently, within the RWA, the dynamics are split in subspaces with a fixed number of excitations which makes the low-energy dynamics amenable, at least numerically.

C. Polaron picture
It has been shown that in the low-energy sector of a spin(s)-boson model [(2)] is well approximated by an effective, excitation-number-conserving Hamiltonian derived from a Polaron transformation [70,71]. The basic idea is to construct a unitary transformation that disentangles the TLS from the bath. This unitary transformation depends on some parameters that are found with the variational principle. In this case, the ansatz is, Here |0 is the photon vacuum state (b k |0 = 0 for all k) and the spin state is arbitrary. The varational parameters are the c s -coefficients and {f k }, the N parameters Figure 1. (a) Schematic depiction of two qubits coupled to specific sites of a linear cavity array. Where g is the coupling constant, ∆ is the energy difference between the two states of the qubits, ω0 is the resonance frequency for photons in the cavity (omitted in the coupled cavities for aesthetic purposes) and λ is the hopping constant for photons travelling between cavities. The yellow shades represent localized-photon clouds. (b) Finite-band dispersion relation of the model. (c) Spectral density function for the model.
in the unitary U P . Building up on previous work from McCutcheon et al. [76] and Zheng et al. [77] we used a natural extension of the single-qubit Polaron transform valid for arbitrarily distant qubits Provided there is no privileged direction of travel, we can assume that for each boson with wavenumber k there will be another with −k, so that |f k | = |f −k |. From that, and the fact that the sine is odd, Eq. (9) factors as with It turns out that minimizing the energy of the spinsboson [(2)], is done by finding the ground state of the effective spin model with and the renormalized qubits frequency In the next section we will work explicit expressions in the case of one and two qubits. A generalized Polaron transformation is discussed in App. C, where it is shown that the much less cumbersome Eq. (9) is sufficiently good.

III. SINGLE QUBIT CASE
In the single qubit case, N q = 1, Hamiltonian [(2)] is nothing but the spin-boson model [20] [65,Chap. 3]. In this section, we tackle the ground-state properties, the single-qubit bound states, and the spontaneous emission within the USC regime for the cavity array model [Eq. (6)].

A. Ground state
Setting N q = 1, Eqs. (11) and (12) yield that the minimum of the energy is reached when σ z |s 1 = −|s 1 [Cf. (12)] with which is minimum when [69], Putting together Eqs. (14) and (16), we realize that the qubit frequency renormalizes to zero as the coupling strength increases. This is a well known result [20]. Besides, this renormalization is the responsible for the localization-delocalization phase transition that corresponds to the ferromagnetic-antiferromagnetic phase transition in the Kondo model [78]. The delocalized phase corresponds to ∆ r → 0, then the qubit state can be in either the symmetric or antisymmetric superpositions of the eigenstates of σ z . On the other hand, if c k = 0, the spin is at an eigenstate of σ z , which corresponds to the localized sector [79]. Using Eqs. (5) and (16) we can rewrite Eq. (14) as, Having a phase transition depends on J(ω) [80]. In our system it is not expected to have critical behaviour [81]. It is not within the aspirations of this work to study (the absence of) this phase transition, partly because it is not clear that the PT is valid in these ranges, so we will restrict our study to the so-called ultra-strong coupling region, g ∈ (0, ∼ 0.5), where we are confident that the Polaron ansatz works [34,72]. As we can see from Fig. 2(a), this region is characterized by a significant, albeit not complete, shrinkage of the tunneling frequency (∆), and as such we expect predictions from RWA to fail. We can further characterize the GS by computing spin observables as is P e = σ + σ − , the probability of having the spin excited. This is an insightful observable because it relates a measurable quantity, P e , to the renormalized frequency, ∆ r : Here we used σ + σ − = 1 2 (σ z + I) together with GS|σ z |GS = 0|U † P σ z U P |0 = − ∆r ∆ . In Fig. 2(b) we show the dependence of P e with ∆ and g, alongside is the GS energy plotted in the same parameter range, in Fig. 2(c). These are signatures of RWA failure, since within the RWA both P e and E GS are zero.
Another interesting observable is the spatial distribution of the photons, b † n b n . Some algebra (fully done in with f n = 1 √ N k e ik(n−N/2) f k being the Fourier transform of f k . We center the transformation at the qubit position that it is understood to be at the middle of the chain. Notice that f n has the clear interpretation of being the real-space variational amplitudes for the Polaron transformation. Figure 4 shows the spatial distribution of photons as calculated in Eq. (19). We observe that they are well localized around the impurity, exhibiting exponential decay f n ∼ exp{−κ GS (n − N/2)} with localization lenght (See App. A 3 for a proof) The photons dressing the impurity are commonly named virtual photons in reference to their special properties of being non-propagating and exponentially localized, that distinguish them from real photons [82,Chap. 1.3]. Again, this is an effect of being in the USC regime and is in stark contrast with the GS found with the RWA which is trivially |0 1 ; 0 .

B. Bound states
We discuss now the single excitation bound states (SEBS), which are the basis for creating effective interactions between the qubits. Before moving to the USC regime, let us summarize the existence of bound states within the RWA approximation where the number of excitations is conserved, see Sect. II B. In this case, the lowest energy bound states are localized eigenstates in the single excitation subspace. Its energy must be outside of the single-photon band. Given a general photonic model, its existence is not guaranteed; i.e., the eigenvalue equation may not have solutions for energies outside of the dispersion relation [47,59]. Notice that photons in these states cannot propagate. They can be thought of as particles trying to enter a potential barrier greater than their energy, and as such, their wavefunction must be exponentially decaying with the distance from the qubit. It turns out that, within the RWA, Hamiltonian [(1)] or [(2)], always accepts two exponentially localized eigenstates: one with energies above and other below the photonic band. In the full model [(2)] the number of excitations is not conserved and we cannot work in the single-excitation subspace. On the other hand, in the Polaron picture the effective Hamiltonian H P = U † P HU P is approximately number-conserving (see App. A 4 for details on the derivation), Here, h.o.t. stands for higher-order terms of order O(f 3 ) with two and more excitations.
Thus, in the Polaron picture, H P conserves the number of excitations and becomes tractable with the same techniques as RWA models; in particular we can compute the singleexcitation eigenstates. It is interesting to note that the GS obtained from the variational method is an eigenstate of H P with eigenvalue equal to the GS energy. This gives us a sense of consistency that confirms the effective RWA model is accurate: If the GS is well caught, one expects that the first excitations are single particle (quasiparticles) excitations over it.   In App. A 6 we show that Eq. (21) admits a bound state below the band, with energy E 1 , and that its localization length is given by, . See App. A 7 for the proof. In Fig. 3(b) we observe the exponential tails. We observe that, the higher the qubit bare-frequency is the peaks become shorter and they also get broader. Qualitatively, we can understand this as follows. The Polaron transformation is local in space [34], thus, for discussing the tails we can argue in the Polaron picture. Because the total number of excitations is 1 in this subspace, the sum of the values of the number of photons in each site must add up to 1, minus the amount taken up by the qubit (which is expected to decrease as the bare-∆ increases). That is why, as the peaks become smaller, they also get broader, in order to preserve the number of excitations. Figure  3(c) shows the de-excitation energy for several values of ∆, referred to the lower band limit, which means that there exists a bound state below the band for all values of ∆.
We can now compare the difference between the results provided by the Polaron transform to those obtained using the RWA. Figure 3(a) shows the relative energy difference between the SEBS calculated with each method.
The difference increases with g, becoming significant for g ∼ 0.1. We have shown only one value of ∆ for clarity, as all values behaved similarly, being the difference greater the smaller the value of ∆. The fact that the PT predicts different bound state energies is not sufficient to declare it superior to the RWA, it could be the case that these results were worse than those provided by the RWA. The definitive confirmation comes from Fig. 5, where we compare the bosonic spatial distributions from the PT and the RWA with those generated by exact diagonalisation of the Hamiltonian for different coupling strengths (see App. A 5 for details on the calculation). The results from the PT are in agreement with the numerical results, both in the ground and excited states. In addition, we see how the RWA and PT coincide for low values of g, but the RWA immediately begins to underestimate the number of photons when the value of g increases beyond g = 0.1. Exact diagonalisation is very limited because the state-space grows exponentially with the number of elements. In addition, exact might be an overstatement considering that one must limit the number of excitations per site in order to have a finite size Hamiltonian. That is why only 12 sites were used in the benchmark for g = 0.05, a number that had to be reduced for greater values of g in order to accommodate more excitations per site while maintaining the state-space size allowed by our  2)] has, at least, an even bound state |E 2 [32]. One can define a band of one-photon states over |E 2 : |k, E 2 [85]. The parity of these states is odd. The hypothetical bound state |E u 1 would also be odd, since it has one excitation in the RWA limit. This implies that, in order for this state to exist, it cannot be embedded in the band formed by |k, E 2 , since otherwise they would hybridize. A necessary condition is: Otherwise, |E u 1 does not exist. To demonstrate the latter, we note that the bound state energy is such that the two photon band). The overlap occurs (and thus the non-existence) if Putting it all together we arrive to the condition for existence given by Eq. (23). It seems a paradox, since this state does exist in the RWA for all ω 0 and λ. The puzzle is solved by noting that in the full model this state becomes a resonance with a lifetime that diverges in the RWA limit.

C. Spontaneous emission
To end our analysis of the single-qubit model we discuss the behaviour of the system during spontaneous emission. We assume the atom-waveguide at the GS, then the qubit is driven within a π-pulse. After the πpulse, the wavefunction is given by |Ψ(0) = σ + |GS . Since [σ x , U P ] = 0, we may work in the single excitation manifold in the Polaron picture. Employing the single excitation ansatz |ψ P = (βσ + + k β k a † k )|0; 0 , the solution is obtained as the inverse Laplace transform (24) The properties of the (inverse) Laplace transform determine the spontanteous emission. In particular, since 0; 0|a k H P σ + |0; 0 | 2 = 2∆ 2 r |f k | 2 , in the continuum limit the sum in Eq. (24) can be converted to an integral over the spectral density J(ω). Let us discuss the two main contributions to this integral. Far from the band limits, J(ω) is sufficiently smooth and the main contribution comes from the poles in the sum, yielding the exponential decay exp[−J(∆ r )t]. Notice, that this is analogous to the RWA result (where the spontaneous emission is J(∆)) but now it is renormalized [72]. The other important feature is the long time dynamics of β(t) which accounts for the qubit thermalization process. The final value theorem, lim s→0 sβ(s) = lim t→∞ β(t), tells that β = 0 if some divergence occurs in that integral. This occurs if bound states exist. Physically, this means that the initially excited state overlaps with the bound state [42,44]. This is conveniently calculated by chosing as a basis in the single excitation manifold, Where |E 1 is the bound state and |E 1 ⊥ p are all other eigenstates orthogonal to it. We recall that the bound state can be written in terms of the original states spanning the one-excitation subspace The first term corresponds to the initial state of the system |ψ 0 = |1 |0 , which indicates that the initial state has some projection on to the bound state, The projection onto the orthogonal basis states will contribute to the continuum and as such, it will not contribute the long time dynamics. The projection onto the bound state is responsible for the divergence and thus for the nonzero value for β(t → ∞). Doing the algebra and computing the observable (notice our return to the lab frame) we obtain that: In Fig. 6 we confirm this expression. The evolution converges to the stationary value predicted by Eq. (28). Also shown in Fig. 6 is the difference with σ z GS , which becomes significant as the ratio g/∆ increases, that is, as the system progresses into the USC regime. Let us emphasize that these results show that our theory is able to deal with the dynamics in USC confirming the peculiarities of the thermalisation process when both the light-matter coupling is non-perturbative and there exist excited bound states.

IV. TWO-QUBIT CASE
We tackle the case of two qubits coupled to the cavity array. Much like in the single qubit case, we first report the results for the ground state continuing with the bound states properties. We put emphasis in the qubitqubit interactions mediated by the cavity array. As an application, we devise a simple state transfer between two distant qubits that uses those interactions.  Figure 6. Evolution of σ z (magnetization) for an initially excited qubit as a function of ∆ for a fixed g = 0.3. Solid coloured lines represent the simulated evolution while dashed coloured lines mark the stationary value predicted analytically, Eq. (28). For contrast, the solid black line corresponds to a Markovian evolution calculated by applying the FGR to the excited and ground states. Coloured shaded boxes have been used to showcase the difference between the stationary magnetisation for each ∆ and the corresponding ground state magnetization, σ z GS .

A. Ground state
Setting N q = 2 and x = x 1 − x 2 in Eqs. (12) and (13) yields (see Apps. B 1-B 3) a spin model which can be diagonalized to yield a ferromagnetic GS of the form where |00 ≡ |s 1 = 0, s 2 = 0 and the coefficients are By Eq. (11), the GS mean energy is which is minimum for, see also Refs. [76,77] and App. B 4 for a detailed derivation, We have introduced the constant E = ∆ 2 r + J 2 to ease notation. It is immediate to check that, should the interaction constant (J ) vanish, we would recuperate the expression of f k that we found in the single-qubit case. This indeed happens when we set the qubits infinitely apart, as will be shown shortly. It is also evident that f k is even with respect to k, which matches the restriction we imposed so that the PT could be factored, see Eq. (10). Figure 8 shows that the dependence of J with x is exponential. This implies that the GS is a ferromagnetic state in a short-range Ising model. As such, in a multi-qubit scenario, only the interaction with firstnearest neighbours would have to be taken into account. Following our analysis of the single-qubit case, it is useful to study the renormalization of the bare frequency ∆ with g. We have used a distance of n = 2 sites to illustrate the deviation from the results obtained in the one-qubit scenario. Figure 7(b) shows that the influence of the neighbouring qubit sharpens the renormalization process, making the system go into full renormalization at lower values of g. Granted, this effect vanishes if one places the qubits further apart. Due to the exponential decay of J , we have found that at distances of around 20 sites the results obtained for one and two qubits are indistinguishable. Back when we studied the single-qubit system we showed that the probability of having an excited spin state was an observable directly related to the renormalization of the bare frequency. The extension to two qubits is straightforward Where |GS S = cos θ |00 + sin θ |11 now. Notice that at large distances, as J → 0 then cos θ → 1 and sin θ → 0, so Eq. (35) reduces to twice the probability found for a single qubit [Eq. (18)]. The effect of the Ising interaction is revealed at short distances where P e deviates from the single-qubit result, no longer equating to the sum of two non-interacting spins. We can again probe the spatial localisation of the bosonic cloud. Following the scheme presented in the single-qubit case, we obtain b † n b n = |f n,1 | 2 + |f n,2 | 2 + 2 cos θ sin θ Re f * n,1 f n,2 , See App. B 5 for details on this calculation. It is interesting to see the overlap between the two bosonic clouds surrounding each qubit. Figure 9 shows this phenomenon for a value of n = 3 where the overlap is significant. In the same figure, the inset monitors the effect of the coalescence of the clouds as the two qubits approach each other.

B. Bound states
Analogously to the single-qubit case, we seek an effective Hamiltonian for the two-qubit model that is a good approximation of the full Hamiltonian but conserves the number of excitations, allowing us to restrict our search for bound states to the one-excitation subspace. We have discussed how f k for two-qubits converges to the expression of f k for a single qubit, so if we assume J to be small, we can write This allows us to reach the effective Hamiltonian By construction, the GS obtained by applying the variational method is an eigenstate of H P and the eigenvalue also coincides with the variational energy. Once again, we can diagonalise the restriction of H P in search of states whose energy lies below the band limit and are, thus, bound. In this case, we expect to find two bound states, corresponding to the symmetric and antisymmetric combination of the wavefunctions corresponding to each single-qubit bound state (See Fig. 10). A general eigenstate of H P has the form The reader might recall that in the single-qubit case the one-excitation subspace was spanned by |1 |0 and |0 |1 k and wonder why we cannot substitute |S |1 k by |00 |1 k in the two-qubit basis. In that sense, it must be clarified that we seek to work in a subspace that is one excitation above the GS, regardless of however many excitations the GS contains. The proposed state, Eq.  (40), on the other hand, comes with a problem since the subspace spanned by Eq. (40) is not closed under the action of the Hamiltonian, e.g.
Fortunately, this second contribution is of second order in f k . Besides, the terms containing σ − 1 a † k in Hamiltonian [(39)] are of the order of f k . Thus they are h.o.t that, consistently with Eq. (21), are discarded. Figure 10 shows β n for the two lowest energy eigenstates. Where λ n is obtained by Fourier transforming λ k in Eq. (40). In the single-qubit case, we showed that there exists a bound state in the form of a cloud of virtual photons localised around the qubit. We have also shown that two sufficiently distant qubits do no interact and, as such, their wavefunctions do not overlap, contributing two bound states of equal energy to the spectrum. As the two qubits approach, we expect the increasing overlap to break the degeneracy, and split the two bound states into different energies. If that is the case, the interaction can cause the energy of the antisymmetric state to rise above the lower band limit, forcing it to no longer be bound (nor antisymmetric), as the corresponding photons have an allowed frequency in the waveguide, and as such they no longer exhibit exponential decay. These oscillating eigenstates are referred to as scattering states [64]. Figure 11 shows the aforementioned effect. In the figure we also compare our results wich the ones obtained within the RWA. We conclude that the latter understimates the interaction between the two bound states (see below). In Fig. 12 a comparison between the spatial distribution of photons for two different distances is drawn. As the two qubits approach, the difference in profiles becomes significant, and, should they reach n = 2, the antisymmetric bound state would cease to exist as it enters the allowed frequency band. See App. B 6 for a calculation of b † n b n .

C. State transfer
Inspired by Fig. 11 and restricting ourselves to the two bound states we can define a tight-binding Hamiltonian where |L represents the bound state of the left-most qubit and |R represents the bound state of the rightmost one. The eigenstates of H T B are the symmetric and antisymmetric combinations of |L and |R , provided τ = 0, with respective energies − τ and + τ . This simplified model is the basis for the study of effective interactions between bound states, which, as we introduced, provides a means to engineer lossless state transfer protocols through virtual photons, one of the main objectives of our work. Real, propagating photons can be used to transport information between distant qubits, but, even in one-dimensional arrays where an emitted photon travels non-dissipatively, there exist undesired losses intrinsic to the emission process. The reason is simple, in the absence of anisotropies, it is equally likely that a radiated photon will travel in the direction of the neighbouring qubit as it is for it to travel in the opposite direction, and thus be lost. The transmission of information via virtual photons bypasses this limitation. By virtue of them being non-radiative, there is no loss, information is shared between close qubits through the overlap of their photonic clouds.
To exemplify perfect lossless state-transfer using bound states, we propose the protocol shown in Fig. 13. Its purpose is to transfer the excited state from one qubit to the other deterministically. First, the left-most qubit is initialized in its excited state and the other is kept in its ground state while both are un-coupled from the waveguide. We assume that g 1 and g 2 , the coupling constant of each qubit to the waveguide, can be tuned independently, and so at t 0 , g 1 is increased adiabatically, so that the excited qubit entangles with the waveguide progres-sively, to become a bound-state, by means of the Adiabatic Theorem. Then, at t 1 , g 2 is increased diabatically to match g 1 . This sudden change in the Hamiltonian does not allow the state to evolve quasi-statically into the new eigenstate and instead gives rise to Rabi oscillations between the left and right bound states, whose symmetric and antysimmetric combinations are actually the eigenstates of the new Hamiltonian. Knowing the hopping frequency (τ ), we can interrupt the dynamics, by diabatically zeroing g 1 at t 2 , at the precise moment where the system is fully in the right bound state. Finally, g 2 is lowered adiabatically, so the right bound state transforms into an excited right-most qubit, succesfully completing the state-transfer protocol at t f . This protocol can be applied sequentially to a succession of qubits, to effectively transport a state along the waveguide.
It is important to note that this method is limited by the fact that the interaction decays exponentially, and this limitation is twofold. Firstly, an exponential decay means that, in order to transport a state between distant qubits, many ancilla quibts are required, placed in close formation, so that there exists an effective interaction amongst every pair of consecutive qubits. For every qubit added to the chain, the system becomes more susceptible to decoherence and losses. In addition, the hopping frequency (τ ), which is proportional to the coupling, is what determines the speed at which a single iteration of the protocol can be performed, so it is against our interest that the coupling decays so rapidly. One may even doubt if an exponentially decaying effective interaction would be, at any range, intense enough to not be overpowered by spurious dipole-dipole interactions between the qubits, which decay more slowly, with a power-law. Fortunately, we can assure that the effective interaction is orders of magnitude greater than dipole-dipole interactions, since the latter is of the order of 10 −4 eV for nearest neighbours withing a crystal lattice (d ∼ 1Å) [86]. In our set up, there is a non-negligible interaction up to distances of around 3 sites, which in experimental realisations of quantum circuits have sizes of millimetres.

V. CONCLUSIONS
We have discussed the main properties of bound states in waveguide QED beyond the RWA paradigm. In other words, we have quantified the corrections to the standard calculations where the qubits-photons interaction is assumed to be number conserving based on the perturbative character of the latter. We have shown that the Polaron technique is useful. It provides a unitary transformation that disentangles qubits and waveguide and the interaction, within this picture, is effectively number conserving. Therefore, it allows to export techniques as the Weisskopf-Wigner theory and intuitions to a broader range of light-matter coupling strengths where the RWA fails. The main results discussed in the paper are as follows. Figure 13. State transfer protocol between two qubits coupled to the same linear cavity array with distinct tunable coupling constants g1 and g2. Bound states are depicted as parabolic instead of exponential for aesthetic purposes.
We have extended the calculations for the spontaneous emission up to moderate light-matter couplings obtaining a renormalization of the rate (due to the qubit-frequency renormalization). The existence criteria for bound states has been generalised and its role in the thermalization of the qubits has been discussed. Finally, we have computed the effective spin-spin interactions both through vacuum fluctuations and bound states. We sketched a perfect state transfer protocol among bound states.
If we express the Polaron transform as U P = exp[A], we can apply the properties we just proved to arrive at the basic commutation relations.
In order to continue, we must first calculate U † p b † k b p U p . Since we have already taken f k as real in previous calculations, we assume it to be real here as well.
We are now equipped with the necessary ingredients to compute 00|U † p b † k b p U p |00 . Considering that the state |00 does not connect through the second and third terms, the mean value is just f k f p .
With that, we simply have

Exponential localisation of the GS
We recall that Noticing that . Therefore, which yields Eq (20) in the main text.

Single qubit effective Hamiltonian
The strict application of the Polaron transform to the original Hamiltonian, H P = U † P HU P , yields the transformed Hamiltonian We can further simplify it by expanding the exponential term. Making use of the Baker-Campbell-Hausdorff formula and taking {f k } real we get A power series expansions of the non constant terms gives Ignoring higher order terms, the right hand side of Eq. (A14) becomes Reintroducing this result in H P yields Considering that σ x σ z + σ x = 2σ − and −σ x σ z + σ x = 2σ + we can combine the second and second-to-last terms to arrive at the final expression for H P 5. Calculation of b † n bn for the SEBS of a single qubit The SEBS will be a state of the form The first term is constant so |v connects entirely yielding f k f p . The last term only connects the bosonic part of |v to give λ k λ p . Finally, the second term cross-connects the two components of |v resulting in λ 0 f p λ * k + λ 0 f k λ p . Reintroducing these partial results into Eq.

Existence of bound states in USC
We work in the Polaron picture. A non-normalized single excitation is, It is an eigenstate iff The solution for E is found by searching the zeros of the function F (E) [Cf. with the RWA case in Ref. 47] If E < min[ω k ], the state is a SEBS. Notice that the term in brackets is is a monotonically decreasing function with g. Therefore if a bound state exists for g → 0 + then it will exists for any finite value of g. For our model in the limit g → 0 + a bound state below the band exists [47], thus the existence of bound states in the USC is guaranteed.

Localization lenght
Apart from their existence the key property of bound states is their localization lenght. From, Eq. (A24b) we obtain that: In the log g-regime 2∆ r f k ∼ g and we can neglect the term 2∆ r k f k f k , therefore by simple inspection we see that Looking at Eq. (A22), Section A 3 and Eq. (20) the localization is given by κ −1 SEBS = max(κ −1 GS , κ −1 ) as given by Eq. (22).
We must now calculate U † p b † k b p U p in the two-qubit case.
The last two terms give, Putting everything together one has The ground state has no photons, so it only connects with itself through the first and second-to-last terms of The first term connects the GS with itself completely, while the other cross-connects the spin terms |00 and |11 . This yields Completing the Fourier transform one finally arrives at b † n b n = |f n,1 | 2 + |f n,2 | 2 + 4αβ Re f n,1 f * n,2 .
Where f n,1 is the fourier transform of f k, 1 , defined as 6. Calculation of b † n bn for the bound states of the two-qubit scenario The SEBS will be states of the form We saw in App. B 5 that Contrary to what happened with the GS, all terms must now be considered because the SEBS connect through them all in one way or another. Thus, we must study each term individually.
The first two are trivial, as they connect SEBS completely with them selves, so they will not be discussed. The second term is more interesting. Through the term λ 0 |01 |0 in |v becomes which connects with λ k (α |00 + β |11 ) |1 k to yield − λ 0 λ k (βf p,1 + αf p,2 ) .
Naturally, the term λ 1 |10 |0 in |v becomes which connects with λ 0 |01 |0 to yield the complex conjugate of its counterpart. Lastly, the term (α |00 + β |11 ) k λ k |1 k becomes and connects with (α |00 + β |11 ) k λ k |1 k to yield Finally, the term b † k b p connects the p th and k th photonic terms to yield λ * k λ p . Summarizing, we have Completing the Fourier transform, one arrives at b † n b n = |f n,1 | 2 + |f n,2 | 2 − 2λ 0 Re{λ * n (βf n,1 + αf n,2 )} − 2λ 1 Re{λ * n (αf n,1 + βf n,2 )} + 4λ 0 λ 1 Re f n,1 f * n,2 + 4αβ 1 − λ 2 0 − λ 2 1 Re f n,1 f * n,2 + |λ n | 2 (B34) Where f n,1 is the fourier transform of f k,1 , defined as with,α = k (f k b † k e ikxj −f * k b k e −ikxj ) andβ = k (l k b † k e ikxj −l * k b k e −ikxj ). When the new variational parameters {l k } vanish, Eq. (C2) reduces to Eq. (9). In fact, provided there is no privileged direction of travel, so that |f k | = |f −k | and |l k | = |l −k |, and the fact that the sine is odd, it can be seen that the transform factors as with Setting x 1 − x 2 = x, the minimization of the ground state energy using this new transform yields the following spin model with, This spin model is not exactly solvable, but performing perturbation theory on the term 2 (σ x 1 + σ x 2 ), we find a ground state energy of the form where η = √ 2 which is minimum when