Distant emitters in ultrastrong waveguide QED: Ground-state properties and non-Markovian dynamics

PHYSICAL REVIEW A 104, 053701 (2021) Distant emitters in ultrastrong waveguide QED: Ground-state properties and non-Markovian dynamics Carlos A. González-Gutiérrez ,*Juan Román-Roche , and David Zueco Instituto de Nanociencia y Materiales de Aragón (INMA), CSIC–Universidad de Zaragoza, Zaragoza 50009, Spain (Received 15 June 2021; accepted 12 October 2021; published 1 November 2021) Starting from the paradigmatic spin-boson model (SBM), we investigate the static and dynamical properties of a system of two distant two-level emitters coupled to a one-dimensional Ohmic waveguide beyond the rotating wave approximation. Employing static and dynamical polaron Ansätze we study the effects of finite separation distance on the behavior of the photon-mediated Ising-like interaction, qubit frequency renormalization, ground-state magnetization, and entanglement entropy of the two-qubit system. Based on previous works we derive an effective approximate Hamiltonian for the two-impurity SBM that preserves the excitation-number and thus facilitates the analytical treatment. In particular, it allows us to introduce non-Markovianity arising from delay-feedback effects in two distant emitters in the so-called ultrastrong coupling (USC) regime. We test our results with numerical simulations performed over a discretized circuit-QED model, finding perfect agreement with previous results, and showing interesting dynamical effects arising in ultrastrong waveguide QED with distant emitters. In particular, we revisit the Fermi two-atom problem showing that, in the USC regime, initial correlations yield two different evolutions for symmetric and antisymmetric states even before the emitters become causally connected. Finally, we demonstrate that the collective dynamics, e.g., superradiance or subradiance, are affected not only by the distance between emitters, but also by the coupling, due to significant frequency renormalization. This constitutes another dynamical consequence of the USC regime. DOI: 10.1103/PhysRevA.104.053701 I. INTRODUCTION In most scenarios, photons are weakly coupled to matter, and they travel fast relative to the matter dimensions. In waveguide QED (wQED), however, they are confined in onedimensional waveguides to enhance the light-matter coupling [1–3]. Also, well-separated emitters can be coupled to the guide, such that photons take an appreciable amount of time traveling between them. In these cases and depending on the experimental conditions and architectures [3–11], neglecting strong light-matter correlations and/or retardation effects is not always justified. When taken into account, emitter dynamics becomes non-Markovian, i.e., the system exhibits memory. This is interpreted as a back-flow of information between the traveling photons and the emitters. Non-Markovianity occurs already at the level of a single emitter interacting with structured environments or reservoirs [12–19], in the presence of bound states [20], or in the ultrastrong coupling (USC) regime [10,11]. In this particular regime, higher order processes beyond the creation (annihilation) of one photon by annihilating (creating) one matter excitation, play a fundamental role. Entering the USC regime implies that the rotating-wave-approximation (RWA) for the interaction breaks down and the atomic bare parameters get renormalized, by either the Bloch-Siegert shift [21], the effective qubit-qubit couplings [22,23] in cavity QED, or renormalization due to the coupling to the continuum electromagnetic (EM) field in wQED [24]. The ground state also becomes nontrivial [25]. *[email protected] Apart from non-Markovianity, USC in wQED has interesting consequences, such as the localization-delocalization transition [26,27], particle production [28], nonlinear optics at the single photon limit [29,30], or vacuum light emission [31]. In this work, we will consider this regime of coupling. When several emitters are present, radiation retardation effects act as an alternative source of non-Markovianity. This has been recently studied outside of the USC regime showing that even in a simple system of two separated two-level emitters, the collective dynamics can be significantly affected by non-Markovian interference caused by radiation-delayed feedback between them [32–38]. In this case retardation effects become important when the atomic lifetime γ−1∼x/vg, for a separation distance x. In the USC regime this condition depends on the coupling and frequency renormalization, as we will see below. The case of three distant emitters has also been reported recently in [39]. This retarded back-action can lead to new collective states in which the emission rate of photons by the emitters can be enhanced or inhibited beyond the usual Markovian limit. This effect, termed superduperradiance [36], allows decay rates larger than 2γ, for two emitters with individual decay rate γ, which has been predicted to scale linearly with the number of qubits and numerically confirmed up to 100 emitters coupled to a one-dimensional waveguide [37]. Non-Markovian features in the emission spectrum of a driven two-qubit system were also investigated recently, together with an experimental proposal in circuit QED using transmon qubits and Josephson-Junction arrays [40]. Another interesting proposal for observing modification of collective phenomena has been laid out in [41], by emulating the emitter dynamics in optical waveguide arrays. 2469-9926/2021/104(5)/053701(14) 053701-1 ©2021 American Physical Society CARLOS A. GONZÁLEZ-GUTIÉRREZ et al. PHYSICAL REVIEW A 104, 053701 (2021) In general, memory effects can be viewed as a problem of quantum feedback, where the system is fed with a quantum signal after some time delay. A general theory of quantum feedback using tensor networks was developed in [42] and later generalized for multiple delays in [43]. Efficient numerical methods to treat more complex systems have also been proposed in [44,45], and the capabilities of feedback to generate universal states for quantum computing was demonstrated in [46]. So far, all previous works rely on neglecting counterrotating terms in the light-matter coupling. In this work we present a first step in generalizing the effect of distance to the USC coupling regime, combining both sources of non-Markovianity. This is done by considering an Ohmic environment with nonflat spectral density, and by taking into account delay-memory effects due to finite distance separation between emitters. Using the polaron transformation, a method that has been shown to be useful for discussing the dynamics of wQED in the USC regime, we first review the ground-state properties of two distant emitters as done in [47]. We introduce a discrete model to support our numerical calculations, matching the Ohmic spectral density used throughout the paper. We discuss and characterize the localization-delocalization quantum phase transition, finding closed formulas for the entanglement entropy. Then, we discuss the dynamics of the emitters in the low-energy sector. A detailed discussion of the Fermi two-atom problem in the USC regime is also given. As well, we compute the effective decay rates for initially correlated symmetric (antisymmetric) initial states. In short, we present a formalism for the study of the different sources of non-Markovianity with various emitters in the nonperturbative regime of wQED. The rest of the paper is organized as follows. The next section deals with the model, the polaron formalism, and introduces the discrete model for the waveguide. Then the ground state is discussed. Section IV is the main section of this work. We discuss the low-energy sector of the dynamics. The paper ends with the conclusions, while several technical aspects are assigned to the Appendixes. II. THEORETICAL MODEL We start reviewing the spin-boson model (SBM) for Nq spin-1/2 particles (or qubits) located at positions xjin one spatial dimension. This model is also known as the Nq-impurity SBM, and it is considered a paradigmatic Hamiltonian for the understanding of decoherence, dissipation, and the physics of open quantum systems [24,48]. Applications of this model include problems from fundamental physics such as the quantum to classical transition [48], electronic transport in biological complexes [49], and quantum simulators [50]. The general SBM is described by the following Hamiltonian (¯h=1): H= 2 Nq  j=1 σz j+ k ωka† kak + Nq  j=1 σx j k gk(ake−ikxj+H.c.),(1) where σz j=|ee|j−|gg|j,σx j=|eg|j+|ge|j.The SBM was originally introduced for a single two-state system as a fundamental physical model for decoherence [24]. It describes the interaction of Nqspin-1/2 particles with the surrounding environment considered as a collection of harmonic oscillators at finite temperature (heat bath). In wQED, it is the standard model for either actual spin-1/2 systems or quantum emitters with anharmonic spectrum after performing the two-level approximation. In the latter case, the two-level approximation can be questioned at enough strong couplings [51,52]. In this work, we fix our attention on natural or artificial two-level emitters, for which the SBM provides an excellent description in the USC regime, our main interest here. This ultrastrong coupling regime of light and matter has been recently achieved in several experimental systems involving superconducting circuits [53,54], semiconductors [55], organic aggregates [56], optomechanical systems [57], and others [11,58]. The system-environment coupling can be completely encoded in the so-called spectral function of the environment, defined as [48] J(ω)=2π k|gk|2δ(ω−ωk),(2) where gkis the coupling strength to the kth mode of frequency ωk. The explicit form of this spectral function depends on the physical realization of the corresponding bath. In this paper we will be interested in modeling the environment as an Ohmic waveguide for which gk∝√ωk,giving rise to an Ohmic spectral function J(ω)=παω in the continuum limit, being αthe coupling strength parameter (see Appendix A1b). For the single SBM (Nq=1) and zero temperature it is well known that, as the coupling increases, there exists a critical value where the system suffers a quantum phase transition called localization transition with vanishing spin magnetization σz=0[59]. Experimental realizations of the SBM through the coupling of a flux qubit to an open 1D transmission line using a shared Josephson junction have been achieved in Ref. [54]. In the case of strong coupling, the spontaneous emission rate from the qubit to the waveguide was measured to be γ∼2π×88MHz for a qubit frequency of ∼2π×3.99GHz. From the microscopic description of the SBM we know that α=γ/π (see Appendix A1d), which results in a value of α≈7×10−3.In this regime results are well described within the RWA. On the other hand, to explore the USC regime we require that γ∼, which was also achieved in the same experiment (with a different device), measuring values of γ∼2π×9.24 GHz, and ∼ 2π×7.68GHz, from which we extract the spin boson coupling α≈0.38. This is a clear manifestation that the system has entered the nonperturbative USC regime [54]. From a general view, it seems that experiments with open transmission lines (propagating photons) can explore a wide range of coupling regimes in regards to the SBM, ranging from the underdamped (α<0.5) to the localized regime (α>1.0). Another recent experiment in the context of the driven SBM [60] has shown that values of α∼0.8 are within the current reach. In general, the SBM can not be diagonalized, but a useful approach based on the introduction of a variational displaced 053701-2 DISTANT EMITTERS IN ULTRASTRONG WAVEGUIDE … PHYSICAL REVIEW A 104, 053701 (2021) oscillator basis can be employed for the study of static and dynamical properties of such systems [19,31,47,61–63]. The accuracy of the polaron approach has been tested in several works for one-qubit [27,64–66] and multiqubit systems [67,68]. The variational approach is based on the following multiqubit polaron transformation: UP[{fk}]=exp − Nq  j=1 σx j k (fka† keikxj−f∗ kake−ikxj), = Nq  j=1 Uj,(3) where Uj=exp[−σx jk(fka† keikxj−H.c.)]. The factorizationinEq.(3) in local operators acting on the qubits can be done once isotropic propagation of bosons (photons) is assumed, i.e., |fk|=|f−k|. The ground-state Ansatz can be defined as the application of the unitary in Eq. (3) to a nonentangled state of the qubits and photons |GS[fk,ζ s]=UP sj ζs|s1,...,sNq⊗|0,(4) where |s1,...,sNqis an arbitrary state of the qubits and |0= |0,...,0is the multivacuum state of the photonic waveguide. Within this approach the task of finding the ground-state energy of the SB system is completely equivalent to minimizing the energy over the following effective Nq-qubit Hamiltonian: HS=r 2 Nq  j=1 σz j− i<j Jijσx iσx j+2 k fk(ωkfk−2gk),(5) where the qubit frequencies become renormalized according to the rule r=exp −2 k|fk|2,(6) and the photon-mediated Ising interaction Jij =2 k fk(2gk−ωkfk)cos [k(xi−xj)].(7) The Hamiltonian in Eq. (5) resembles the one describing the Ising model in the presence of a transverse magnetic field. In the case of a single-qubit, the ground-state energy in the polaron basis is EGS =−r/2+kfk[ωkfk−2gk], which can be minimized with respect to the variational parameters {fk}. The corresponding free-energy minimization results in a self-consistent relation for the variational parameters [61] fk=gk ωk+r .(8) In the case of a single qubit, the renormalization frequency can be computed explicitly, obtaining the well-known formula in the scaling limit /ωc1[24,48], r≈ ωcα/(1−α) .(9) It is worth recalling that this renormalization is responsible for the localization-delocalization phase transition that FIG. 1. (a) Spectral function of the Ohmic waveguide J(ω). Open circles are obtained from the discrete model simulations with parameters: {N,ω c,α}={1001,1.0,0.1}and the solid line with the formula J(ω)=παω. (b) Pictorial illustration of two two-level emitters separated a distance xinteracting with the guided modes of an Ohmic one-dimensional waveguide. (c) Circuit-QED implementation of the system (see Appendix A1 for further details). corresponds to the ferromagnetic-antiferromagnetic phase transition in the Kondo model [69]. Here rdrops to zero as the coupling increases, triggering the transition at the critical coupling αc=1. Discrete microscopic model Below we will face situations, mainly when we discuss the dynamics of the system, for which an analytical solution is not possible. Then we approximate the environment using a finite number of modes Nby taking advantage of the discrete model for a transmission line resonator exposed in Appendix A1b. In the continuum limit (N→∞) this model reproduces the spectral linear function of an Ohmic environment with an appropriate cutoff frequency: see Fig. 1(a), where the points are the spectral function constructed from the discrete model, recovering the linear dispersion relation ωk=vg|k|. This finite N-model has already been successfully used to test the dynamics of few emitters interacting with waveguides within [35], and beyond [62] the RWA, with numerical methods based on matrix-product-state simulations (MPS), showing remarkable agreement with the advantage of low computational cost. In Appendix A1, for completeness, we sketch the main characteristics of the discrete model in the context of circuit QED. Of course, when possible, the convergence of this model to the Ohmic (continuum) model has been verified. III. TWO-QUBIT SPIN-BOSON MODEL: STATIC PROPERTIES The ground state of two qubits in an Ohmic environment was discussed in [47]. They used a variational Ansatz equivalent to the variational polaron treatment. In this section, we review ground-state properties that will be needed when discussing the dynamics of the system. This also is useful to benchmark the discrete model for numerical simulations. For the two-qubit case, the relative distance between emitters 053701-3 CARLOS A. GONZÁLEZ-GUTIÉRREZ et al. PHYSICAL REVIEW A 104, 053701 (2021) FIG. 2. Ising-like coupling as a function of the normalized relative distance between the emitters, for different coupling parameters. Solid lines are plotted according to Eq. (12), and empty dots correspond to numerical results based on the discrete Ohmic model with parameters {N,/ω c}={1001,0.2}. plays an important role in the static and dynamical properties. In this situation the induced Ising interaction of the effective Hamiltonian modifies the single-qubit behavior, giving rise to collective effects. For Nq=2 the Hamiltonian (5) reads HS=r 2 2  j=1 σz j−Jσx 1σx 2+2 k fk(ωkfk−2gk),(10) where J(x)=2kfk(2gk−ωkfk) cos(kx),being x=x2− x1the distance between the qubits. In this case the groundstate energy is EGS =− 2 r+J2+2kfk(ωkfk−2gk), and its minimization results in fk=gk ωk E+Jcos(kx) E+Jcos(kx)+2 r/ωk ,(11) where E=2 r+J2. The effective Ising model is characterized by the ratio J/r. In what follows, we compute both parameters numerically by using the finite-Nmodel. However, it is always appealing to have analytical expressions even if approximates. Taking the limit J→0inEq.(11) we recover the one-qubit relation in Eq. (8). Interestingly, the induced Ising coupling can be computed in the large coupling approximation for which fk−→ gk/ωkas r→0, which gives the result [47] J≈1 πωc 0 dωJ(ω) ωcos(ωx/vg) =αωcsinc(ωcx/vg),(12) where we have used the Ohmic spectral function J(ω)= παω, and assumed a linear dispersion relation ωk=vg|k| in the continuum limit. This Ising coupling shows longrange damped oscillations as a function of the distance. This means that, depending on the distance, ferromagnetic (J> 0) or antiferromagnetic (J<0) interaction can be induced. In Fig. 2we show the behavior of the Ising-like coupling as a function of the emitters relative distance for different coupling strengths values. Continuous lines are calculated using the approximation in Eq. (12), and open circles indicate results obtained from numerical simulations based on the microscopic discrete model presented in Appendix A1b.We observe a very good agreement between the continuous and FIG. 3. Renormalized frequency as function of the coupling strength for different relative distances between emitters. Other parameters are the same as in Fig. 2. the discretized microscopic Ohmic model. A more accurate but cumbersome expression for the Ising coupling can also be obtained by using the exact expression for fkin Eq. (8), which includes small correction terms to the formula given in Eq. (12). We show this expression in Appendix A2. Obtaining an expression for the renormalized qubit frequency rfor two qubits is not as simple as in the single-qubit case. In general, it has to be found by numerical means. Assuming that JisgivenbyEq.(12) and therefore does not play any role on the energy minimization, the optimal variational parameters can be approximated by fk≈gkE/(ωkE+2 r), which after replacement into Eq. (6)gives r=exp −2 k g2 kE2 ωkE+2 r2.(13) This equation shows that ris now distance-dependent through the Ising-like coupling. For infinite separation distance (J=0), the renormalized frequency in the scaling limit is given by r≈(e/ωc)α/(1−α), which shows the quantum phase transition at α=1, in agreement with the single impurity case in Eq. (9). For slightly separated qubits in the limit J/ 1, we can take E→J, obtaining that r≈(e2/ωcJ)α/(1−2α),(14) predicting the critical coupling at α=1/2 in stark contrast with a single qubit [47]. It is worth mentioning that this shifted critical coupling can also be inferred by studying the properties of the dynamical Kernel in the general non-Markovian case [62]. In Fig. 3we plot the renormalized frequency varying the relative distance x. We show the lower and upper bounds, given by the expressions for zero and infinite separation distances, respectively. It is clear that the effective interaction between qubits modifies the localization transition point depending of their relative position. In particular, at intermediate relative positions (e.g., x=8xin Fig. 3), r drops abruptly to zero in a discontinuous way, and the localization transition occurs for values of the coupling strength α∈[0.5,1.0]. 053701-4 DISTANT EMITTERS IN ULTRASTRONG WAVEGUIDE … PHYSICAL REVIEW A 104, 053701 (2021) Ground-state and entanglement entropy The ground state of the effective two-spin Hamiltonian (10) can be found by straightforward diagonalization. It reads |GSS=cos θ|gg+sin θ|ee,(15) where cos θ=r+E (r+E)2+J2(16) and sin θ=J (r+E)2+J2.(17) The ground state of the complete system of two-qubits plus the bosonic bath is then obtained by applying the polaron transformation to the product of the ground state |GSSand the multimode vacuum state of the bosonic environment according to the Ansatz (4) |GS=UP|GSS|0 =(cos θσx 1+sin θσx 2)UP|e|g|0.(18) After applying the transformation we get the following entangled state in terms of multimode coherent states of the photonic bath |GS=1 2(sin θ+cos θ)|−−|χ−− +1 2(sin θ−cos θ)|−+|χ−+ +1 2(sin θ−cos θ)|+−|χ+− +1 2(sin θ+cos θ)|++|χ++,(19) where we chose the eigenstates of σx jfor representing the state of the qubits, i.e., |±j=(|ej±|gj)/√2, and the multimode coherent states are given by |χ±±=exp −i k f2 ksin kx k|±(fk+fkeikx ),(20) |χ±∓=exp i k f2 ksin kx k|±fk∓fkeikx,(21) with |αk=D(αk)|0, being D(αk) the displacement operator. Notice that in the deep-strong coupling (localized regime) when r→0, the ground state becomes |GS= (|−−|χ−−+|++|χ++)/√2, which can be recognized as a Schrödinger cat state of light and matter analogous to the one-qubit SBM [70]. Having obtained the ground state, one can also compute relevant physical observables, for instance, the qubit magnetization and the entanglement entropy. In the first case we find that σz jGS =GS|σz j|GS=−r cos 2θcos [φ(x)],(22) where we have introduced the auxiliary function φ(x)= 4kf2 ksin kx. It is also interesting to compute the linear entanglement entropy, defined as SL=1−Tr2 q, where the FIG. 4. Entanglement entropy as a function of the coupling strength for increasing distance between the emitters. Parameters are the same as in Fig. 3. reduce density matrix of the two-qubit system can be obtained by tracing out the bath degrees of freedom, i.e., q= TrenvGS =Trenv|GSGS|. This measure quantifies the degree of mixedness of the two-qubit ground state with the bosonic bath. Remarkably, a rather general expression of the linear entropy for arbitrary emitter separations can be obtained, but the expression is too lengthy to be shown here (see Appendix A3). However, we show the limiting formula for infinite distance SL=1 43−2r 2 −r 4.(23) It is well known that entanglement entropy reflects the appearance of a quantum phase transition, so we expect to observe signals of the localization transition in the entanglement as the coupling is increased, in a similar fashion as observed for rin Fig. 3. Quantum phase transitions manifest as a nonanalyticity in the entanglement contained in the total state of the system and its environment [71]. We confirm this by plotting the entanglement entropy as a function of the coupling in Fig. 4for increasing relative distance. The dashed line in Fig. 4is given by Eq. (23), indicating a smooth behavior of the entropy with no discontinuities for large couplings. This is expected from the one-qubit results for the renormalized frequency discussed in Sec. II for which αc∼1. The asymptotic entropy in this case is a complete mixed state of the qubits with SL→3/4, which can be seen directly from Eq. (23) in the limit r→0. For finite separation we have a richer behavior, as sudden decay of the entropy occurs at specific couplings, reflecting the influence of the effective qubit-qubit interaction on the localization transition. Notice also that in this case the steady-state entropy as a function of the coupling is SL=1/2. This is because the Ising interaction dominates over the magnetic field produced by rin the two-spin effective model, and the ground state becomes degenerate with reduced density matrix q=1/2. This is also clear from the general expression for the linear entropy in Eq. (A20)inthe case of sin θ=cos θ=1/√2. IV. DYNAMICS Not only the ground state but the first excited states are well captured by the polaron Ansatz [31,62,63,72]. It is convenient to work in the polaron picture (1): HP=U† PHUP≈Heff, 053701-5 CARLOS A. GONZÁLEZ-GUTIÉRREZ et al. PHYSICAL REVIEW A 104, 053701 (2021) which can be approximated by [62] Heff =r 2 2  j=1 σz j−2r 2  j=1 σz j kk fkfkei(k−k)xja† kak + 2  j=1 k2rf(0) k+(r−ωk)σ− ja† keikxj+H.c. + k ωka† kak−Jσx 1σx 2+2 k fk(ωkfk−2gk).(24) Here f(0) kis the variational parameter of a single qubit given by Eq. (8), i.e., when the ferromagnetic interaction is negligible (J=0), and is a small correction. The key fact here is that Heff preserves the excitation number. Therefore, this effective Hamiltonian can be now treated using standard RWA methods. In the following we will study the generated dynamics in the polaron frame in order to extract relevant information of the system in the laboratory picture. As the Hamiltonian Heff preserves the number of total excitations, we can project the dynamics in the single excitation sector spanned by the basis: {|e|g|0,|g|e|0,|GSS|1k}. This means that photonic excitations are created from the polaron ground state of the two-spin system. The state of the system at any time is postulated via the following dynamical polaron Ansatz [62]: |(t)P=c1|e|g+c2|g|e+ k ψka† k|GSS|0,(25) where the time-dependent probability amplitudes require to fulfill the normalization condition i|ci|2+k|ψk|2=1. These coefficients satisfy the coupled set of equations, i˙c1=−Jc2+ k Gke−ikx1ψkcos θ, (26) i˙c2=−Jc1+ k Gke−ikx2ψkcos θ, (27) i˙ ψk=−(˜ r−ωk)ψk+Gk j eikxjcjcos θ +2rfk j,k fkei(k−k)xjψkcos 2θ. (28) Here we have defined the functions Gk=2rf(0) k+(r− ωk) and ˜ r=rcos 2θ+Jsin 2θ. Using the discrete model, this set of equations are integrated numerically and, then, transformed back to the original laboratory picture. A. Dynamics in the polaron frame For future convenience and understanding, it is useful to start analyzing the dynamics in the polaron picture. The reason is twofold. It is easier, so some explicit formulas can be obtained, and it will be convenient to understand the dynamics in the laboratory frame. As well, if the coupling is small enough, both polaron and polaron frame are equivalent. Moving to the rotating frame for the photon amplitudes, ˜ ψk=e−i(˜ r−ωk)tψk, and performing the integration, i.e., tracing out the photonic modes, and setting ˜ ψk(0) =0, we arrive to i˙ci=−Jcj−2i2 rcos2θ πωc 0 dωJ(ω) (ω+r)2 ×t 0 dτei(˜ r−ω)τci(t−τ) −2i2 rcos2θ πei˜ rx/vgωc 0 dωJ(ω) (ω+r)2 ×t−x/vg 0 dτei(˜ r−ω)τcj(t−x/vg−τ),(29) where it was assumed that |Gk|2≈42 rf(0) k 2for small . Here the presence of non-Markovian processes is evident. The second and third contributions in the r.h.s both depend on the spectral function J(ω) of the waveguide (environment). So there is an unavoidable source of non-Markovianity as consequence of the memory of the environment, which cannot be considered as flat (Markovian) in general [73]. On the other hand, the third term contains an additional source of non-Markovianity due to the finite distance effect between emitters, which also enters in the nonlocal parts of the Kernel and therefore invalidates the Markov approximation. These two non-Markovian timescales are different in general and need to be distinguished in the dynamical evolution. Unfortunately, this is not so simple, and we have to get rid of one of such memory time scales in order to study the behavior of Eq. (29). At this point we follow the lines of the WeisskopfWigner theory of spontaneous emission [74], and neglect the memory effects of the environment in order to evaluate the contribution of the last two terms in Eq. (29). We can then make the replacement ci(t−τ)≈ci(t) and extend the time integral to infinity. By using the identity [75] lim t→∞t 0 dτei(˜ r−ω)τ=πδ(ω−˜ r)+iP1 ˜ r−ω,(30) we arrive to the coupled delay-differential equations ˙ci(t)=iJcj(t)−γ 2ci(t) −γ 2ei˜ rx/vgcj(t−x/vg)(t−x/vg),(31) with the spontaneous emission rate defined by γ=42 rJ(˜ r) (r+˜ r)2cos2θ, (32) and the function (τ) being the Heaviside step function. Obviously, once the approximation ci(t−τ)≈ci(t) is made, the only source of non-Markovianity is due to the finite distance of the emitters, which induces back-action of the field, leading to a radiation feedback phenomenon in the system. As a consequence, these equations resemble the delay dynamics obtained for two distant emitters coupled to a waveguide with flat spectral density of the field modes [36]. However, in this case a coherent coupling term between the emitters appears as a consequence of the Ising coupling Jinduced by the polaron transformation. As this coherent interaction depends on the distance between emitters [see Eq. (12) and Fig. 2], it contributes to the dynamics at small separation distances. 053701-6 DISTANT EMITTERS IN ULTRASTRONG WAVEGUIDE … PHYSICAL REVIEW A 104, 053701 (2021) Notice also that the spontaneous emission rate and the relative phase gained by the retarded radiation gets renormalized as they are given in terms of the polaron qubit frequency r. For infinitely distant emitters, the spontaneous emission rate reduces to the result of a single emitter with γ=J(r)= παr[19,62]. It is worth mentioning that in deriving the delay differential equations (31), we have neglected the Lamb-shift correction to the qubit frequency, which is formally given by δL=2 π2 rPωc 0 dωJ(ω) (ω+r)2(˜ r−ω).(33) A detailed analysis of the consequences of collective spontaneous emission and Lamb shift is given in [35]. The magnetization dynamics (or population inversion) in the polaron image can be expressed in terms of the qubit probability amplitudes, σz i(t)P=cos2θ[2|ci(t)|2−1] −sin2θ[2|cj(t)|2−1].(34) Let us now analyze the emergence of non-Markovian dynamics due to the back-action induced by the delay term in Eq. (31) for different initial states of the emitters. Naturally, we expect these results to agree with those obtained within the RWA with a flat environment for very small couplings [36]. We first explore the case of large but not infinite separation distance between the emitters for which J≈0. The steadystate solutions in the polaron picture can be computed from the standard final value theorem [76], by taking the Laplace transform in Eqs. (31). We first take initial symmetric (antisymmetric) states in the polaron Ansatz, i.e., c1(0) =c2(0) = 1/√2orc1(0) =−c2(0) =1/√2. For these initial conditions the Laplace transformed amplitudes read ˜c±(s)=1 √2[s+γ/2±γe(ir−s)x/vg/2],(35) where ˜c+=˜c1=˜c2, and ˜c−=˜c1=−˜c2. For a renormalized qubit frequency, and separation distance satisfying the condition rx/vg=2nπ, being n∈Z, the symmetric initial state will decay to the polaron ground state, i.e., limt→∞ c+(t)= 0. However, if the initial state is antisymmetric we get a bound state with finite excitation probability amplitude given by limt→∞ c−(t)=(1 +γx/2vg)−1/√2[36]. This gives the steady-state magnetization for the antisymmetric state σz −P=(1 +γx/2vg)−2−1, which is different from −1 for finite-distance emitters [see Fig. 5(a)]. This bound state in the continuum (BIC) is the result of a trapped stationary excitation between the two emitters acting as an effective cavity formed by two perfectly reflective mirrors [77,78]. B. Laboratory frame The previous analysis will help us discuss the actual dynamics. They can be obtained by back-transforming with the polaron unitary. It turns out this can be done exactly, as we have restricted our analysis to the single-excitation polaron manifold. In particular, the magnetization in the laboratory FIG. 5. Dynamics of the magnetization σzfor initial symmetric |+ 0and antisymmetric |− 0states. Dash-dotted lines indicate the evolution in the polaron picture. To enable comparison with the true dynamics, the polaron lines are renormalized by r/.From(a)to (d) the increasing couplings are α=0.01, α=0.1, α=0.5, and α=1.0. Other parameters are L=40π,N=1001, x=2πvg/. picture reads σz i(t)=P(t)|U† Pσz iUP|(t)P =r σz iP+4 cos θReci k fkψ∗ k −4sinθRecj k fkψ∗ k +4 cos 2θ kk fkfkψkψ∗ k,(36) which shows the dynamical contributions due to each emitter and field correlations. In (36) we see that the differences between both frames are O(fk); therefore in the weak coupling both frames give the same results, as expected by construction. In the same region, one can perform the RWA approximation already in the original spin boson model (1). This is illustrated in Fig. 5(a) for which we recover RWA results. As we increase the system-bath coupling α, nonMarkovian dynamics emerges as consequence of the memory of the bath and also from interference caused by retardation effects [36]. At this point the collective two-qubit dynamics is no longer captured by the effective delay dynamics dictated by Eqs. (31). We show these results in Figs. 5(b),5(c), and 5(d) for a fixed separation distance between the emitters. We also expect that for very strong couplings, the localization transition takes place (r→0), and the two-qubit system freezes in a state of zero magnetization. In this deep-strong coupling regime retardation effects are no longer present in the magnetization. 1. The Fermi problem in the USC regime Another interesting fact is the apparent causality violation that can be observed in Fig. 5(c). For moderate couplings (α=0.01,0.1), where RWA still applies, the dynamics of symmetric and antisymmetric initial states are exactly the 053701-7 CARLOS A. GONZÁLEZ-GUTIÉRREZ et al. PHYSICAL REVIEW A 104, 053701 (2021) same before the retardation time, i.e., independent exponential decay of both emitters at the same rate (Fermi golden rule) occurs. As we increase the coupling we can see from Fig. 5(c) that this is no longer the case. Symmetric and antisymmetric initial states start to have different dynamics long before they reach the vertical line indicating the light-cone overlap between the two emitters. This apparent paradox was first studied by Fermi in 1932, and it is known as the Fermi problem [79]. Fortunately this paradox can be explained in terms of correlations between spacelike distant events arising in the ultrastrong coupling regime and that are absent in usual RWA models [80]. The Fermi paradox can be formulated explicitly as follows: in a system of two distant emitters, where initially one of them is in the excited state and the other is in its ground state, the following interesting question arises: is it possible to excite the second emitter through the spontaneous emission of the first one at a time t<x/vg? It turns out that this is not possible, and there is no causality violation in this problem, as was pointed out by Fermi [79], and better justified later by several authors [81–83]. However, as was shown by Sabin et al. [80], if nonlocal correlations are shared initially by the two emitters, the probability of finding the second emitter in the excited state and the first one in the ground state can be different from zero even at times t<x/vg. We will show here that the existence of these correlations explains the different dynamics experienced by symmetric and antisymmetric initial states at times before the retardation time observed in Fig. 5. We now move for convenience to the Heisenberg picture, and write the equation of motion for the magnetization using the original laboratory Hamiltonian (1) σz j(t)=σz j(0) +2t 0 dτσy j(τ)V(xj,τ),(37) where we have defined V(xj,t)=V0(xj,t)+V1(xj,t)+V2(xj,t) (38) with the field operator V0(xj,t)= k gk(ake−ikxje−iωkt+H.c.) (39) and the qubit operator Vi(xj,t)=− kt 0 dτ|gk|2σx i(τ)eik(xj−xi)e−iωk(t−τ)+H.c. (40) The expectation value of the qubit magnetization is then the sum of different contributions σz j(t)=σz j(0)+2t 0 dτσy j(τ)V0(xj,τ) +2t 0 dτσy j(τ)V1(xj,τ) +2t 0 dτσy j(τ)V2(xj,τ),(41) where the last two correlation terms suggest some influence or back-action of the respective distant emitter. However, it can be shown that both contributions are independent FIG. 6. Dynamics of σz 1,2starting from the initial state |eg. From (a) to (d) the parameters are the same as in Fig. 5. of such distant emitter [80]. To make this clear let us focus in the first emitter magnetization, i.e., σz 1(t).Inthis case one can show that the correlation σy 1(t)V2(x1,t)∝ d dxσy 1(t)σx 2(t−x/vg)(t−x/vg), showing that this function is zero for times t<x/vg.An analogous calculation for the autocorrelation function of the first emitter shows that it is in fact independent of the second emitter, so the only contribution to the qubit magnetization is due to the second term in Eq. (41). We can examine this contribution for initial symmetric or antisymmetric states of the emitters, which are initial entangled states in the polaron picture. Then σy 1(τ)V0(x1,τ)=± 0|σy 1(τ)V0(x1,τ)|± 0,(42) where |± 0= 1 √2(σx 1±σx 2)UP|g|g|0.An explicit calculation of these qubit-field correlations is complicated as consequence of the nontrivial action of the polaron transform, but one can see that, in general, they are different from zero for all times and have different value for symmetric and antisymmetric initial states. This is in strong contrast with the RWA case, where these correlations are totally absent and this contribution is exactly zero, as the polaron transform does not enter in the initial state. This analysis explains the different dynamics experienced by both initial states when the system enters in the ultrastrong coupling regime as shown in Fig. 5(c). It is also worth mentioning that although we were interested in the difference in time evolution of symmetric and antisymmetric initial states, the above analysis is in fact independent of the initial state of the two-qubit system. These correlations emerge as the coupling of the two distant emitters to the Ohmic waveguide increases, breaking down the usual RWA. We have additionally explored the dynamics for the initial state |e|g|0in the polaron image. This is shown in Fig. 6, where the magnetization of each emitter is plotted for increasing coupling and for the same parameters used in Fig. 5. For weak coupling we observe exponential decay (Fermi golden rule) of the excited emitter and the subsequent excitation of the noninitially excited emitter at time t=x/vg. For α=0.1 [Fig. 6(b)] the spontaneous decay rate gets renormalized, and we are able to see the damped coherent exchange of the excitation between the two emitters until it reaches the expected equilibrium magnetization given by −r/. Excitation of the second emitter before retardation time is shown in Fig. 6(c) for α=0.5 with notorious effects on the dynamics 053701-8 DISTANT EMITTERS IN ULTRASTRONG WAVEGUIDE … PHYSICAL REVIEW A 104, 053701 (2021) FIG. 7. Dynamics of the instantaneous decay for increasing coupling α. (a) Symmetric initial state |+ 0, and (b) antisymmetric initial state |− 0. Other parameters are the same as in Fig. 5. at x/vgand the expected localization transition for α=1.0 in Fig. 6(d), where localized and small amplitude (2r/) Rabi oscillations between the emitters are observed due to the coherent Ising interaction moving the excitation between the two emitters [see inset in Fig. 6(d)]. 2. Time-dependent decay rate The analysis of non-Markovian dynamics allows a description in terms of a master equation that is local in time [73]. It turns out that the time-dependent generator for this master equation can be written in the standard Lindblad form with time-dependent decay rates. Following this approach, and from our analysis in the polaron picture, we can define a time-dependent decay rate for each emitter as γ(i)(t)=−2Re˙ci(t) ci(t).(43) In general, as in the case of probability amplitudes in the strong coupling regime, an analytical evaluation of this formula is not possible, and we must resort to numerical calculations based on the discrete model. The time evolution of this quantity will contain both sources of non-Markovianity, the one due to the memory of the bath, and the influence of retardation effects. Of course, as Eq. (43) is defined in terms of the polaron probability amplitudes, we need to transform it back to the laboratory frame. In Fig. 7we show the numerical calculation for the collective dynamics of the time-dependent decay rate for symmetric (γ+) and antisymmetric (γ−) initial states. As expected, for weak coupling, the decay rate starts to oscillate around the Fermi golden rule value (almost without renormalization) γ±=γ=παr∼ =πα until the feedback from the other emitter suddenly accelerates the decay process, reaching a maximum value of ∼4.4γfor the symmetric state in Fig. 7(a). At this point, the emitters see each other, forming a collective state and emitting in a superradiant manner beyond the Markovian Dicke limit of 2γfor two emitters located at the same position. This non-Markovian collective enhancement of spontaneous emission has been reported recently within the RWA limit in Ref. [36], where it was shown that larger values can be reached at a particular critical distance. We must emphasize that even for weak coupling the Ohmic environment model used here is not flat, so it is reasonable to expect that our calculation of the decay rate does not match exactly the results for a flat environment. Also, it is clear that non-Markovian processes involving retardation are also present beyond the RWA limits. On the other hand, for the antisymmetric initial state in Fig. 7(b) we have similar behavior at times t<x/vg. Once they are causally connected, in this case, the antisymmetric state is a dark state and no decay is observed. For weak coupling this is compatible with the appearance of bound states in the continuum for subradiant initial states in RWA [36,77]. We also have verified that for slightly distant emitters in the weak coupling we recover the usual Markovian dynamics for initial symmetric (antisymmetric) states. As the coupling is increased we observe the suppression of emission caused by the qubit renormalization frequency, i.e., γ=παraccompanied by an oscillating behavior accounting for emission and reabsorption of radiation, which explains the negative values of the decay rate. In the deep-strong coupling regime, where α→1, the system undergoes a localization transition with zero magnetization and the two-qubit system does not radiate anymore. It is worth recalling that what determines the nature of the collective state generated after both qubits become causally connected is the actual phase gained by the radiation field during propagation inside the waveguide. This phase is responsible of the constructive (destructive) interference of the delay dynamics [36,40]. Within the RWA this phase is given by x/vg, which is independent of the qubit-field coupling. In the USC regime the accumulated phase is now coupling and distance-dependent for finite distant emitters as we saw in Sec. III, through the renormalized qubit frequency, dropping to zero in the localized phase for a particular critical coupling. This is a purely USC phenomenon, in which non-Markovian interference effects are also affected by the coupling to the waveguide as we observe in Fig. 7. In particular superradiant and subradiant behaviors arise as appropriate combinations of the symmetry of the initial state, the coupling strength, and the relative distance between the emitters, reflecting the complex relation between the Hamiltonian parameters in this regime, which has to be taken into account in real experiments aiming to investigate non-Markovian delay dynamics in waveguide QED [84]. As complementary information we show in Fig. 8 the individual emission (absorption) rates, where the initial state is |e|g|0in the polaron picture in correspondence with Fig. 6, showing the anticorrelated dynamics of both emitters. In this case one can clearly see the influence of Rabi oscillations in the two-qubit system over the spontaneous emission (absorption) rates at strong coupling in the localized regime. V. SUMMARY AND CONCLUSIONS Employing theoretical tools based on the polaron transformation, we have addressed the problem of two distant emitters interacting with a one-dimensional Ohmic waveguide beyond RWA. We have tested early results for the static properties of the two-impurity SBM in the continuum limit, comparing those with a discretized finite N-model based on circuit QED, 053701-9