Dynamical signatures of bound states in waveguide QED

We study the spontaneous decay of an impurity coupled to a linear array of bosonic cavities forming a single-band photonic waveguide. The average frequency of the emitted photon is different from the frequency for single-photon resonant scattering, which perfectly matches the bare frequency of the excited state of the impurity. We study how the energy of the excited state of the impurity influences the spatial profile of the emitted photon. The farther the energy is from the middle of the photonic band, the farther the wave packet is from the causal limit. In particular, if the energy lies in the middle of the band, the wave packet is localized around the causal limit. Besides, the occupation of the excited state of the impurity presents a rich dynamics: it shows an exponential decay up to intermediate times, this is followed by a power-law tail in the long-time regime, and it finally reaches an oscillatory stationary regime. Finally, we show that this phenomenology is robust under the presence of losses, both in the impurity and the cavities.

In this work, we study the signatures of bound states in the spontaneous decay in a waveguide-QED scenario. We choose a prototypical model where an impurity is coupled to a bosonic medium: a tight-binding model, which gives a cosine-shaped band. Due to the finite width of the band, two states appear bound to the impurity. This problem has already been treated in the literature when the energy of excited state of the impurity is in the middle of the band [14] or when it is close to its inferior limit, so the superior limit of the band can be neglected [12]. Here, we solve it for general values of the parameters, which are the energy of the excited state of the impurity with respect to the band and the ratio between the impurity-photon coupling and the bandwidth.
We find an energy shift of the emitted photon with respect to the energy required to excite the impurity, provided the * edusb@unizar.es latter is not in the middle of the photonic band. Naively, one could argue that spectral features in the spontaneous emission should also appear in the scattering, since a scattering process comprises both absorption and emission of the photon by the impurity. However, there is no shift in the single-photon scattering [39][40][41]. Second, we study the spatial profile of the emitted photon and discuss the differences with respect to a photonic medium with a linear dispersion relation. Finally, we find a rich dynamics in the excited state of the impurity. First, it decays exponentially at a decay rate different from that given by Fermi's golden rule, oscillating with a phase which is shifted with respect to the bare energy of the impurity. This shift, which corresponds to the Lamb effect, is different from the shift found in the energy of the emitted photon. After the initial exponential decay, the dynamics presents an algebraic decay due to the presence of singularities in the density of photonic states. These power tails are robust under the presence of losses, both in the impurity and in the cavities. Eventually, it reaches a stationary oscillating regime.
The paper is organized as follows. In Sec. II, we introduce the Hamiltonian and summarize both its spectrum in the singleexcitation subspace and its one-photon-scattering properties. In Sec. III, we discuss the main results of the paper. First, we present the already mentioned frequency shift of the emitted photon as a function of the coupling constant and the energy of the excited state of the impurity. Then we study the spatial distribution of the emitted photon. We next characterize the spontaneous emission when the impurity is initially excited and discuss the effect of the losses. We end with the conclusions in Sec. IV. Some technical details are described in the appendixes.

A. Hamiltonian and bound states
The photonic medium is an infinite chain of discrete bosonic sites coupled to an impurity located at site x 0 = 0. The Hamiltonian of the combined system is (h = 1) where a x and a † x annihilate and create, respectively, a photon at position x, and b and b † annihilate and create excitations at the impurity. This impurity can be a two-level atom or qubit, another resonator, a spin, or any system equivalent to a qubit in the single-particle subspace. The energy of the excited state of the impurity is . Henceforth, we borrow the nomenclature from molecular physics and refer to as the exciton energy; in the same way, b and b † will annihilate and create an exciton. The band of free photons is defined by a dispersion relation which depends on both the on-site photon energy and the hopping parameter J : ω k = − 2J cos k, k being the dimensionless momentum and ω k the corresponding energy. In consequence, the bandwidth is 4J . The momentum k lies in [−π,π). The group velocity is v k ≡ dω k /dk = 2J sin k. The interaction Hamiltonian [last term in Eq. (1)] is the dipole-field Hamiltonian in the rotating-wave approximation, which is given by the celebrated Jaynes-Cummings model, where g is the coupling constant. A scheme of the system and the dispersion relation ω k are shown in Figs. 1(a) and 1(b), respectively. This model can be realized with the instances of quantum technologies enumerated in Sec. I.
Due to the rotating-wave approximation, Hamiltonian (1) commutes with the number operator N ≡ x a † Thanks to this symmetry, this model is analytically solvable in the single-excitation subspace. A complete basis is formed by the scattering eigenstates | k [39] and the bound states | ± [42,43]. We introduce the bound states now and we leave the scattering ones for the next subsection. They read  State |0 represents the vacuum state of the system (a x |0 = b|0 = 0). The factor N ± is a normalization constant, 1/|κ ± | is the localization length, and d ± is the exciton amplitude. The energy of | ± is ω ± = − J (e −κ ± + e κ ± ). The expressions of d ± and N ± , as well as the computation of κ ± , are given in Appendix A. The quantities κ ± fix the properties of the bound states; namely, their energies ω ± , exciton amplitudes d ± , and normalization factors N ± .
We plot the bound-state energies ω ± as a function of the coupling constant g, as well as the band limits in Fig. 2. Two cases are shown: (i) the exciton energy in the middle of the band ( − = 0; solid lines) and (ii) closer to the band bottom ( − = −J ; dotted-dashed lines). The energies of the bound states lie outside the band, thus they are localized (not propagating). As g → 0, ω −(+) approaches the bottom (top) of the band. If the exciton energy coincides with the band center, the energies of the bound states are symmetrically located. Otherwise, if the exciton energy is below the center, − < 0, the energy of the lower bound state ω − moves away from the exciton energy more rapidly than the energy of the upper bound state ω + does, and vice versa. Therefore, the position of the exciton energy with respect to the band center originates an asymmetry between ω + and ω − .

B. One-photon scattering
Let us now review the form of the single-particle scattering eigenstates of (1) and their physical implications. They read [39] The coefficients t k and r k are the transmission and reflection amplitudes for an incident plane wave, respectively. They are given by A well-known feature in this system [39][40][41]44,45] is that it presents perfect reflection, R k ≡ |r k | 2 = 1, if the energy of the input photon is equal to ; see Eqs. (4) and (5). This is illustrated in Figs. 3(a) and 3(b), where R k is plotted as a function of (ω k − )/J for several values of and as a function of (ω k − )/J and ( − )/J , respectively. Considering the input as a single-photon-spectroscopy probe, we could be tempted to argue that, as in scattering, the impurity emission is also maximum at resonance. We show that, due to the presence of bound states, this is not the case.

III. SPONTANEOUS DECAY
We now discuss the spontaneous emission of the exciton. For that, we consider that the impurity is excited at t = 0, | (0) = b † |0 , and compute the time evolution of the system. Spanning this state in bound and scattering eigenstates, Eqs. (2) and (3), respectively, the state at time t is with In the following, we exploit these formulas to obtain our results. First, we discuss the behavior of the mean energy of the emitted wave packet. Then we describe the spatial profile of that photon depending on . Finally, we study the dynamics of the exciton.

A. Energy shift
The state given by Eq. (7) can be used to obtain the average value of Hamiltonian (1). As it is a conserved quantity, it must be equal to the value at t = 0, which is The average energy for the propagating field is with P lig ≡ |c + | 2 + |c − | 2 . Using Eq. (10), ω ph can be written in a more convenient way, which shows that the energy of the emitted photon is typically different from because of the presence of the bound states. In short, the amount of energy going to the propagating states must compensate that going to the bound ones so that the total energy is conserved. This is the physical origin of the energy shift. This confirms the nonequivalence between scattering and emission spectra, since the scattering resonance always occurs when the input energy is ; see Eq. (4) and Fig. 3. The energy of the emitted photon (ω ph − )/J is plotted as a function of ( − )/J in Fig. 4(a) for several values of g. The closer is to the band edges, the more ω ph departs from . In fact, if is close to the bottom of the band, where the frequency shift is larger, the effect of the upper bound state is negligible (|c + | 2 |c − | 2 ), and vice versa. In conclusion, the frequency shift survives in waveguides without an upper cutoff. The shift increases monotonically with g. Eventually, as g/J → ∞, the emitted energy coincides with the middle of the band for all . Note that, when the exciton energy is in the middle of the band, i.e., when = , the following relation holds: Inserting this in Eq. (12), we conclude that the emitted energy is equal to the exciton one, ω ph = . This is related to the symmetry of the energy of the bound states, discussed in Sec. II A (see Fig. 2).
We also study the energy distribution of the emitted photon |c k | 2 . We plot it as a function of (ω maximum for ω k = , the difference being larger the closer is to one of the band edges. This deviation of the maximum away from implies a frequency shift of the emitted photon, as shown in Eq. (12) and Fig. 4(a). The reason is simple. In spontaneous emission some energy is released into the bound states, with a mean energy that does not generally match the exciton energy. Therefore, the coupling into flying photons must compensate for this imbalance. However, because bound and scattering states are orthogonal, the former do not play any role in the latter. It is worth emphasizing that this mechanism is rather general. In any photonic system supporting singleparticle bound states, the frequency of the flying photon arising from spontaneous emission will present a shift with respect to that of the exciton.
We also characterize the emission probability into propagating modes, P emission ≡ 1 − P lig = 1 − |c + | 2 − |c − | 2 , in Fig. 5. Two effects are observed. First, the emission into bound states is negligible (P emission 1) in the range g/J 1. With an increase in this ratio, P emission decreases. In addition, the closer is to the band gap, the smaller P emission is. At any rate, the emission probability is appreciable for very large values of the ratio g/J : for instance, g/J 2.5 yields P emission 0.25 for the values of considered in Fig. 5.

B. Emitted field
We now study the spatial profile of the emitted field. We compute the amplitudes in position space, φ x (t) ≡ 0|a x | (t) The probability |φ x | 2 is mostly confined within the causal cone. For |x| > x max , it is not 0 but it decays exponentially, as expected for the free-field scalar propagator (Sec. 4.5 in [46], Sec. 2 in [47]). If is in the middle of the band, the emitted photon has a momentum distribution peaked around k = π/2, where v k = v max . If = , the velocity of the emitted photon is not peaked around v max so the maximum of |φ x | 2 is below x max (see the dashed red curve in Fig. 6, where − = −J ). Finally, note that the emitted photon would be well peaked around |x max | in position space, independently of the value of , if the dispersion relation were linear.

C. Impurity dynamics
We finish with a detailed study of the exciton dynamics. From Eq.
First, we focus on c s e (t), with The behavior of c s e (t) is determined by the kernel F (y), which is related to the density of photonic states as a function of the dimensionless energy y = cos k. This kernel is plotted in Fig. 7(a). At sufficiently long times, the oscillating term in the integral, (14), e i2yJ t , cancels out any smooth contribution of F (y). Therefore, the asymptotic relaxation dynamics is governed by the sharpest peaks and the singularities of F (y). There are three main contributions: (i) a Lorentzian peak, associated with a pole of F (y) in the complex plane; (ii) two peaks appearing at y * ± , with y * ± close to ±1 [see Fig. 7(b), where we zoom in on F (y) around y = −1]; and (iii) the singular points at y = ±1, where the first derivative of F (y) is discontinuous. All these features are clearly shown in Figs. 7(a) and 7(b).
The Lorentzian peak gives an exponential decay c s e (t) ∼ e −(iϕ+1/2τ 0 )t . This is equivalent to an excited atom emitting photons into free space. This contribution is the fastest and main one for short enough times, t < τ 0 , since it comes from the widest peak in F (y); see Fig. 7(a). We compare the (numerical) exact results for τ 0 and δϕ ≡ ϕ − , computed by integrating Eq. (14), with those obtained with the Lorentzian approximation of F (y) in Fig. 8. We also compare the results to those obtained with Fermi's golden rule, τ FGR 0 = J sin k /g 2 , with k such that ω k = , and ϕ FGR = . Fermi's golden rule describes accurately the exact results when is around the middle of the band, but corrections are necessary when gets closer to the band edges and when the coupling g increases. The exciton energy appears in the phase of the exponential up to a correction: ϕ = + δϕ. Thus, δϕ is the Lamb shift due to the coupling to the photonic bath. Note that this Lamb shift is different from the energy shift of the emitted photon [compare Fig. 4 to Figs. 8(c) and 8(d)], even though both converge to in the limit g/J → 0. In fact, as stated, there is another characteristic energy of the system with a different behavior: the single-photon reflection resonance, which occurs exactly at the bare excitation energy [see Fig. 3 and Eq. (5)].
At later times, t τ 0 , the singular parts of F (y) are relevant. Singularities give nonexponential decays [5,11]. In particular, the contribution of the peaks of F (y) at y * ± , with y * ± ±1, starts to dominate. Let us define the widths of these peaks at y * ± as y ± ≡ |y * ± ∓ 1| [see Fig. 7(b)]. For short enough times, when e i2Jyt can be considered to be constant for y ∈ (−1, − 1 + y − ) and y ∈ (1 − y + ,1), kernel F (y) can be approximated by setting g = 0 [dotted black curve in Fig. 7(a)]. At g = 0, the kernel diverges as 1/ 1 − y 2 when y → ±1. This kind of singularity gives an algebraic decay t −1/2 for c s e (t). For long enough times, when e i2J ty cannot be taken as a constant, we have to consider the full kernel, with the actual value of g. Therefore, the mentioned divergences are rounded off and the algebraic decay is modified by exponential factors. In other words, these peaks provide a contribution c s e (t) = t −1/2 (a − e −i2J t e −t/2τ 1,− + a + e i2J t e −t/2τ 1,+ ), with τ 1,± = (4J y ± ) −1 . The values of the constants a ± , as well as details on the computation, are reported in Appendix C.
Eventually, these exponential contributions vanish. The only surviving contribution comes from the singularities at the band edges. There, F (y) is not differentiable and gives a nonexponential (power-law) contribution for all times to c s e (t), which dominates for t τ 0 ,τ 1,± . We show in Appendix C that this contribution goes as c s e (t) ∼ t −3/2 cos(2J t − 3π/4). This transition between t −1/2 and t −3/2 decay is discussed in [12], but those authors did not see the oscillating factors, since they took the exciton energy very close to the lower part of the band, neglecting the contribution of the upper bound state. As mentioned, this decay with t −3/2 originates from a discontinuity in the derivative of the density of photonic states and is quite common in impurity decay problems [48], both for continuous systems [49,50] and for discrete ones [51][52][53].
We sum up all this information in Fig. 9, where we plot the impurity dynamics for − = 0 and g = J /5 (same parameters as in Fig. 7), using logarithmic scale. For the sake of clarity, we average the oscillations coming from the different contributions: a − e −i2J t + a + e i2J t (arising from the peaks around y * ± ), cos(2J t − 3π/4) (from the singularities at y = ±1), and cos[(ω + − ω − )t] [from c b e (t)]. The population P e (t) ≡ |c e (t)| 2 is shown as the dotted black curve. It first decays as e −t/τ 0 . In addition, the bound-state term dominates over the remaining contributions from the scattering states. Therefore, after a transient period, P e (t) achieves the stationary regime of P b e (t) (dashed purple curve; remember that we do not show the oscillations). We also show P s e (t) ≡ |c s e (t)| 2 with the solid red curve. After the initial exponential decay with e −t/τ 0 , where P s e (t) P e (t), it decays subexponentially. To show the different contributions to this subexponential decay more clearly, we plot it in the inset at log-log scale. After the mentioned exponential decay with e −t/τ 0 , it follows a decay with t −1 e −τ 1 for τ 0 t τ 1 (as − = 0, τ 1 ≡ τ 1,+ = τ 1,− ; in particular, τ 1 200τ 0 for the chosen parameters). Eventually, as t τ 1 , P s e (t) goes with t −3 . The agreement between the analytical predictions (dashed blue curve for e −t/τ 0 , dotted green curve for t −1 e −t/τ 1 , and dotted-dashed orange curve for t −3 ) and the exact (numerical) integration is clear in the figure.
Finally, even though we have focused on the case with in the middle of the band, the mathematical analysis shown in Appendix C is general, so another choice of parameters will give the same qualitative behavior.

Losses
Here we incorporate losses to the model. We add an imaginary part both to the exciton energy and to the cavity energy,˜ = − iγ e /2 and˜ = − iγ c /2.
The dynamics is still given by Eqs. (14) and (15), by changing and by˜ and˜ , respectively. We take γ e/c /g ∼ 0 − 0.15. Considering losses in the exciton, the integrand F (y) resembles that in the lossless case (see Figs. 7 and 10), apart from the fact that now it is a complex function; the same happens if we instead add losses to the cavities. Therefore, we can repeat the analysis of the lossless case.
We illustrate the modifications with γ e = 0 in Fig. 11. Initially, it still decays exponentially, but the decay rate is a sum of the previous one, 1/τ 0 , and γ e : the amplitude reads c sc (t) ∝ e −(iϕ+1/2τ 0 +γ e /2)t [see Fig. 11(a)]. The power law with t −1 , c s e (t) = t −1/2 (a − e −i2J t e −t/2τ 1,− + a + e i2J t e −t/2τ 1,+ ), is preserved. The coefficients a ± , whose expressions are shown in Appendix C, get modified 10 −5 % at most for the chosen values of γ e . Finally, the asymptotic decay with t −3 does not depend on (see Appendix C). The robustness of the power-law tails is shown in Fig. 11(b).
If we instead consider lossy cavities, γ c = 0, there is a global factor e −γ c t/2 multiplying c sc (t) [see Eq. (14)]. When integrating F (y), the imaginary part of adds an increasing exponential e γ c t/2 to c sc (t), contrarily to [see the denominator of F (y), Eq. (15); and have opposite signs]. This increasing exponential cancels out with the global factor e −γ c t/2 . Therefore, no modifications are seen in the initial exponential regime [see Fig. 12(a))]. The global factor e −γ c t/2 suppresses the power laws in the long-time limit. If the characteristic time of the losses 1/γ c is larger than τ 1,± , we can see the power-law tails for intermediate times [see Fig. 12(b)].

IV. CONCLUSIONS
We have discussed the differences between spontaneousdecay and scattering spectra. As we argue above, naively we could expect that the scattering resonance should coincide with the spontaneous-emission energy. However, whereas the scattering resonance is always equal to the exciton energy, we have shown that the emission frequency is shifted. In particular, this shift is more clear as the coupling increases and/or the exciton energy is closer to the band edges. We have also seen that the profile of the emitted photon strongly depends on the exciton energy with respect to the photonic band. Finally, the presence of bound states and a nontrivial density of states makes the impurity dynamics nontrivial, with three dynamical regimes: exponential decay, power law with a transition from t −1 to t −3 , and oscillatory asymptotic. This dynamics has proven to be robust under the presence of losses, both in the atom and in the cavities. Even though the population of the power-law regime is very small, it could be measured. In fact, such power laws have already been measured in the context of dissolved organic materials, where the fluorescence follows an algebraic decay at long times [54].
Some features, such as the spectroscopic shifts in the spontaneously emitted photons, are detectable by tuning up and down the frequency of the exciton with respect to the band edge. For probing the dynamics, we suggest using a more sophisticated protocol that (i) places the exciton energy at the right frequency, (ii) then excites it, and (iii) after a finite time t detunes the exciton and probes its excited-state population dispersively. All these ideas can be implemented in state-of-the-art setups with superconducting cavities and transmon qubits [25] and also with quantum dots in photonic crystals [32][33][34].

APPENDIX A: BOUND STATES
We provide the explicit expressions for d ± , κ ± , and N ± appearing in the text [Eq. (2)]. The excited-state amplitude of the impurity d ± is In order to compute κ ± , we define η ± ≡ e −κ ± and use the eigenvalue equation H | ± = ω ± | ± [43]: This equation has four solutions. However, we have two constraints: (i) Re(κ ± ) > 0, because the photonic cloud must be localized around the impurity and cannot explode at x → ±∞; and (ii) Im(κ ± ) = 0,π , since the energies ω ± = − J (e −κ ± + e κ ± ) are real. With these restrictions, there are only two solutions for η ± , which can be found numerically.
If we take the limit J → ∞, where the dispersion tends to be linear, the valid solutions for η ± are ±1, so Re(κ ± ) = 0. Therefore, | ± are not bound anymore. In fact, they converge to the scattering states | k with k = 0 and k = π , that is, those at the band edges.

APPENDIX B: EMITTED FIELD
The profile of the emitted field φ x (t) = 0|a x | (t) is given by where we have used Eq. (7). In order to compute the amplitude 0|a x | k we take the expression of | k , Eq. (3), for k > 0: 0|a x | k = e ikx + r k e −ikx , x < 0; t k e ikx , x 0.

Exponential decay
In order to extract the first exponential decay, we can approximate F (y) by L(y) = a p /(y − y p ), y p being the pole corresponding to the peak of F (y), with −1 < Re(y p ) < 1 and Im(y p ) > 0, and a p the residue of F (y) at y = y p . The value of y p is found numerically, equating the denominator of F (y) to 0 [see Eq. (15)]. The residue a p is computed by definition. We extend the integration domain to ±∞. Then, applying the residue theorem, c s e (t) = i8a p (g/J ) 2 e −i t e i2y p J t .
By computing this numerically, we obtain the decay rate τ 0 = [4J Im(y p )] −1 and the phase ϕ = − 2J Re(y p ), as shown in Fig. 8 in the text.

Subexponential regime: t −1/2
Kernel F (y) shows a sharp behavior around y * ± . In fact, it diverges when y → ±1 if g = 0. In order to take this contribution into account, we can approximate F (y) by F (y)| g=0 [see dashed blue curve in Fig. 7(a) If 2 y ± J t 1, with y ± = |y * ± ∓ 1|, the oscillatory term e i2yJ t will not be sensitive to the difference between F (y) and F (y)| g=0 when y is close to the edges. As we are concerned in the contribution around ±1, we can approximate the integral as c s e (t) These integrals are analytical c s e (t) In consequence, P s e (t) decays with (J t) −1 after the initial exponential decay if τ 0 t τ 1,± , with τ 1,± = (4J y ± ) −1 . We can rewrite the last expression by adding the decaying exponentials with τ 1,± : c s e (t) (C5)