Non-Markovian Dynamics of a Qubit Due to Single-Photon Scattering in a Waveguide

We investigate the open dynamics of a qubit due to scattering of a single photon in an infinite or semi-infinite waveguide. Through an exact solution of the time-dependent multi-photon scattering problem, we find the qubit's dynamical map. Tools of open quantum systems theory allow us then to discuss the general features of this map, find the corresponding non-Linbladian master equation, and assess in a rigorous way its non-Markovian nature. The qubit dynamics has distinctive features that, in particular, do not occur in emission processes. Two fundamental sources of non-Markovianity are present: the finite width of the photon wavepacket and the time delay for propagation between the qubit and the end of the semi-infinite waveguide.


I. INTRODUCTION
Waveguide quantum electrodynamics (QED) is an emerging area of quantum optics that investigates coherent coupling between one or more emitters (qubits) and a one-dimensional (1D) photonic waveguide [1][2][3][4][5]. Novel correlations among injected near-resonant photons result from the nonlinearity of the qubits, and intriguing interference effects occur because of the 1D confinement of the light. The field has focused on qubits in a local region for which these correlation and interference effects can be used for local quantum information purposes such as single-photon routing [6], rectification of photonic signals [7][8][9][10], and quantum gates [11][12][13]. This regime of waveguide QED involves neglecting delay times: the time taken by photons to travel between qubits is far shorter than all other characteristic times. However, an important goal for photonic waveguides is to carry out long-distance quantum information tasks such as quantum state transfer between remote quantum memories [14,15]. As these necessarily involve distant qubits, delay times cannot be neglected, leading to different kinds of photon correlation and interference effects through the non-Markovian (NM) nature of the system. Here, we study a model waveguide-QED system with large delay time. We apply recent developments in the theory of open quantum systems (OQS) in order to quantitatively assess the qubit's degree of non-Markovianity.
A large variety of waveguide-QED setups have been experimentally demonstrated in recent years [4,5,16,17]. Because of high photon group velocities and small systems, these experiments are mostly described FIG. 1. A qubit coupled to a waveguide scatters a single-photon wavepacket injected from the left. The qubit may start in an arbitrary state. The semi-infinite waveguide sketched here is terminated by an effective mirror, thus introducing a delay time. by a Markovian approach in which delay times are neglected. In contrast, recent experiments have started entering the regime of non-negligible delay times [18][19][20][21][22], an area that is expected to grow rapidly due to interest in extended systems and long-distance quantum information. Accounting for photon delay times is, however, a challenging theoretical task: only recently have the dynamical effects of long delay times started being investigated [23][24][25][26][27][28][29][30][31][32][33]. A major consequence of delay times is that NM effects are introduced that affect the physics profoundly, as predicted for e.g. qubit-qubit entanglement in emission [24,32] and second-order correlation functions in photon scattering [25,28,29,31].
Clarifying the importance of NM effects and the mechanisms behind their onset is thus pivotal in waveguide QED. At the same time, the theory of OQS [34,35] is currently making major advances, yielding a more accurate understanding of NM effects [36][37][38][39]. Through an approach often inspired by quantum information concepts [40], a number of physical properties such as information back-flow [41] and divisibility [42] have been spotlighted as distinctive manifestations of quantum NM behavior and then used to formulate corresponding quantum non-Markovianity measures. These tools have been effectively applied to dynamics in various scenarios [37,38], including in waveguide QED with regard to emission processes [26,32,43] such as a single atom emitting into a semi-infinite waveguide [26].
Motivated by the need to quantify NM effects in photon scattering from qubits, we present a case study of a qubit undergoing single-photon scattering in an infinite or semi-infinite waveguide (see Fig. 1), the latter of which is the basis of the proposed controlled-Z [12] and controlled-NOT [13] gates. We aim at answering two main questions: What are the essential features of the qubit open dynamics during scattering? Is such dynamics NM? The key task is to find the dynamical map (DM) of the qubit in the scattering process, which fully describes the open dynamics and is needed in order to apply OQS tools [38]. A distinctive feature of our open dynamics is that the bath (the waveguide field) is initially in a well-defined single-photon state [44,45]. Toward this task, we tackle in full the time evolution of multiple excitations (in contrast to those limited to the one-excitation sector [44,[46][47][48][49]), a problem that has become important recently [30,31,33,[50][51][52][53][54][55][56][57][58].
Intuitively, one may expect that the dynamics is fully Markovian in the infinite-waveguide case and NM in the semi-infinite case due to the atom-mirror delay time. We show that this expectation is inaccurate in general, mostly because it does not account for a fundamental source of NM behavior namely the wavepacket bandwidth. This NM mechanism is present in an infinite waveguide, while in a semi-infinite waveguide it augments the natural NM behavior coming from the photon delay time. Recently, NM effects in infinite-waveguide scattering were addressed in Ref. [45]. There, however, the qubit is always initially in the ground state, while a fair application of non-Markovianity measures should be based on the entire DM thus requiring consideration of an arbitrary initial state of the qubit.
The paper is organized as follows. We first define the system under consideration in Sec. II. Next, in Sec. III, we find the general form of the the qubit's DM in a single-photon scattering process and discuss its main features. In Sec. IV, we present the time-dependent master equation (ME), which is fulfilled exactly by the qubit state at any time. In Sec. V, we discuss the explicit computation of the dynamical map in the infinite-and semi-infinite-waveguide case (most of the details regarding the former are given in the Appendix). Since this task requires the time evolution of the scattering process, we in particular find a closed delay partial differential equation that holds in the two-excitation sector of the Hilbert space for a semi-infinite waveguide. In Sec. VI, we assess the non-Markovian nature of the scattering DM by making use of non-Markovianity measures. In this way, we identify two fundamental sources of NM behavior: the finiteness of the wavepacket width and the time-delayed feedback due to the mirror. We finally draw our conclusions in Sec. VII. Some technical details are given in the Appendices.

II. SYSTEM
Consider a qubit with ground (excited) state |g (|e ) and frequency ω 0 , which is coupled at x = x 0 to a photonic waveguide (along the x-axis) with linear dispersion. We model the system via the standard [29,59,60] real-space Hamiltonian under the rotating-wave approximation (we set = c = 1 throughout) where the bosonic operatorâ R (x) [â L (x)] annihilates a right-going (left-going) photon at x,σ − =σ † + = |g e|, and V is the qubit-field coupling strength such that the qubit decay rate into the waveguide is Γ = 2V 2 . For an infinite waveguide, the upper integration limit is ξ = +∞ and x 0 = 0, while for a semi-infinite waveguide ξ = 0 + and x 0 = −a (see Fig. 1).

III. DYNAMICAL MAP
By definition the qubit's DM, Φ t , is the superoperator that when applied on any qubit state at t = 0, ρ 0 , returns its state at time t > 0 [34,35], The DM fully specifies the open dynamics of the qubit coupled to the waveguide field, with the latter serving as the reservoir. We now find the DM for a single-photon wavepacket. LetÛ t = e −iĤt be the unitary evolution operator of the joint qubit-field system. The initial state for single-photon scattering is σ 0 = ρ 0 |ϕ ϕ| (tensor product symbols are omitted), where ρ 0 = ρ gg |g g| + ρ ee |e e| + (ρ ge |g e| + h.c.) and |ϕ = dx ϕ(x)â † R (x)|0 is the incoming single-photon (normalized) wavepacket (|0 is the waveguide vacuum state). At time t, the atom-field state is σ(t) =Û t σ 0Û † t =Û t ρ 0 |ϕ ϕ|Û † t . By plugging ρ 0 into σ(t), we get hence for any ρ 0 the knowledge of the pair of elementary unitary processesÛ t |gϕ andÛ t |eϕ fully specifies the time evolution of σ(t). Due to the conservation of the total number of excitations [see (1)], the joint evolved state in the two processes has the form where |Ψ n (t) is the joint wavefunction at time t for n excitations. Here, |φ 1 (t) and |ψ 1 (t) are unnormalized single-photon states, and |χ 2 (t) is an unnormalized two-photon state. Note that (4) [(5)] describes the joint dynamics of a single photon scattering off a qubit initially in the ground [excited] state, which takes place entirely in the one-excitation [two-excitation] sector of the Hilbert space. In particular, Eq. (5) is a stimulated emission process [61].
The qubit state at time t is the marginal ρ(t) = Tr field σ(t). This partial trace can be performed by placing Eqs. (4) and (5) into Eq. (3), which yields ρ(t) = ρ gg φ 1 (t) 2 + ρ ee χ 2 (t) 2 |g g| + ρ gg |e(t)| 2 + ρ ee ψ 1 (t) 2 |e e| where we took advantage of orthogonality between one-photon and two-photon states. Since {|gϕ , |eϕ } are normalized, so are (4) and (5) due to unitarity ofÛ t . Thus, φ 1 (t) 2 + |e(t)| 2 = χ 2 (t) 2 + ψ 1 (t) 2 = 1. By defining three time functions we have φ 1 (t) 2 = 1 − p g (t) and χ 2 (t) 2 = 1 − p e (t). Therefore, changing to the matrix representation, Eq. (6) takes the form (with ρ ee = 1 − ρ gg ) where we defined Note that both p g (t) and p e (t) are the qubit excited-state populations but in the two different processes (4) and (5), respectively. We refer to the qubit DM (8) as the "scattering DM." Since the atom-field initial state σ 0 is a product state andÛ t is unitary, the map Φ t is necessarily completely positive (CP) and trace preserving [40]. In contrast to pure emission processes at zero temperature [34], here both the one-and two-excitation sectors are involved. We stress that the DM is fully independent of the qubit's initial state ρ 0 , being dependent solely on the Hamiltonian (1) and the field initial state |ϕ . This dependance occurs through the functions of time p g/e (t) and c(t) in (8), where p g/e (t) determine the qubit populations while c(t) governs the coherence.
The DM's form is best understood in the Bloch-sphere picture [40] in which a qubit state ρ is represented by the Bloch vector r = (2 Re ρ ge , 2 Im ρ ge , ρ gg − ρ ee ) with r ≤ 1. In this picture, the map Φ t is defined by the vector identity where v(t) = (0, 0, 1 − p e (t) − p g (t)) T and M t is the 3 × 3 matrix Here, 0 = (0, 0) while R θ(t) is a standard 2 × 2 rotation matrix of angle θ(t) = arg[c(t)]. Thus, apart from the rigid displacement v(t) and rotation around the Z-axis, the scattering process shrinks the magnitude of the XY -and Z-components of r(t) by the factors |c(t)| 2 and |∆(t)|, respectively. Since these two factors are generally unequal, the DM transforms the Bloch sphere into an ellipsoid. Such a lack of spherical symmetry does not occur in emission processes [62], thus providing a hallmark of scattering open dynamics. In addition, a careful look at Eqs. (10) and (11) shows that the Bloch vector undergoes a reflection across the XY -plane whenever ∆(t) < 0. This is a further distinctive trait of scattering dynamics, not occurring in emission processes [62], which is shown below to be relevant to the onset of NM behavior.

IV. NON-MARKOVIAN MASTER EQUATION
The most paradigmatic Markovian dynamics is the one described by the celebrated Lindblad ME [34], whereĤ is self-adjoint, and all the rates γ ν 's are positive constants. In our case, the DM (8) is not described by a Lindblad ME; instead, we show that it is described by a time-dependent ME [63]: whereĤ(t) is a time-dependent Hamiltonian, and the jump operatorsL ν =σ ν describe three non-unitary channels with the time-dependent rates γ + (t) for absorption, γ − (t) for emission, and γ z (t) for pure dephasing [explicit forms are given below in Eqs. (16)]. We note that the dephasing term, which reflects the lack of spherical symmetry of the evolved Bloch sphere discussed above, does not occur in spontaneous emission. The general form for a time-dependent ME isρ where L t is a time-dependent linear (and traceless) map, which is fulfilled by ρ(t) as given by Eq. (8). The standard recipe for carrying this out starting from the DM is to first take the derivative of Eq.
The task now reduces to explicitly calculating L t and expressing it in a Lindblad form so as to end up with Eq. (13). This task is efficiently accomplished in the generalized 4-dimensional Bloch space. Recall that the set of four Hermitian operators , respectively -is a basis in the qubit operators's space and fulfills Tr{Ĝ iĜj } = δ ij . We express both Eqs. (13) and (15) in this basis and equate them; some details are presented in Appendix A. The resulting expressions for the time-dependent Hamitonian as well as the three time-dependent rates in ME (13) are given bŷ It can be checked that whenĤ(t) and these rates are placed in it, ME (13) is exactly fulfilled by Eq. (8) at all times t.
Before concluding this section, we note that an exact, differential system (DS) governing the same open dynamics that applies in the case of an infinite waveguide was worked out in Ref. [64] and more recently further investigated and generalized in Refs. [65,66]. For the present case of a single-photon wavepacket and a qubit, this DS has overall three unknowns: two density matrices, one of which is ρ(t), and a traceless non-Hermitian matrix. In contrast to ME (13) here, the DS has the advantage that its time-dependent coefficients are known functions of the wavepacket functional shape (in the time domain). However, since A qubit coupled to a semi-infinite waveguide with qubit-mirror distance a (left) and its equivalent model featuring a qubit coupled to a chiral infinite waveguide at the two points x = ±a (right). such a DS is not closed with respect to ρ(t) it is less suitable for analyzing the general properties of the qubit's dynamical map, which is a major goal of the present work. Finally, while ME (13) holds for a qubit coupled to a generic bosonic bath under the rotating-wave approximation, the DS relies on the further hypothesis of a white-noise bosonic bath (hence, in particular, it does not hold in the semi-infinite-waveguide case).

V. EXPLICIT COMPUTATION OF DM
For the initial state of the waveguide, throughout this paper we consider an exponential incoming wavepacket of the form where k is an arbitrary central frequency, α ≡ δk/Γ is the wavepacket bandwidth in units of Γ, and θ(x) is the step function. This choice of the wavepacket shape is often made (see e.g. [61]) as it has at least three advantages. An exponential shape allows for closed-form solutions in the Laplace domain in some cases. In addition, in a numerical approach, its well-defined wavefront leads to a significant reduction in computational time, which is important in the two-photon sector when there is a long time delay. Finally, such wavepackets can be generated experimentally [67][68][69][70][71][72] by either spontaneous emission of a qubit or tunable, on-demand sources. Our general approach is to plug the ansatz for |Ψ n (t) , Eqs. (4) and (5), into the Schrödinger equation i∂ t |Ψ n (t) =Ĥ|Ψ n (t) to obtain a system of differential equations for the amplitudes that we solve for the three functions p g (t), p e (t), and c(t) in (7) and hence for the DM (8). For an infinite waveguide, this can be accomplished analytically in both the one-and two-excitation sectors as shown in Appendix B. Here, we focus on the far more involved case of a semi-infinite waveguide.
In this case, it is convenient to "unfold" the waveguide semi-axis at the mirror (x = 0) by introducing a chiral field defined on the entire real axis byâ(x) =â R (x)θ(−x) −â L (x)θ(x) (the minus sign encodes the π-phase shift due to reflection from the mirror); see Compared to Eq. (1), the form of the free-field term shows that only one propagation direction is allowed (chirality) while the term ∝ V shows the bi-local coupling of the field to the qubit at points x = ±a (these can be seen as the locations of the real qubit and its mirror image, respectively).

A. Semi-infinite waveguide: one-excitation sector
The wavefunction ansatz in the one-excitation sector is [cf. Eq. (4)] which once inserted into the Schrödinger equation yields the pair of coupled differential equations Integration of the photonic part yields the formal solution [cf. Eq. (17) of Ref. [23]], which once plugged back into the equation for e(t) [Eq. (20b)] yields the delay (ordinary) differential equation (DDE) Equation (22) is the same as the well-known DDE describing the spontaneous emission process in Refs. [23,73,74] but with the presence of the extra source term ∝ Γ/2 due to the incoming photon wavepacket. For the initial conditions e(0) = 0 and φ(x, 0) = ϕ(x) [cf. Eq. (4)], the solution for e(t) obtained by Laplace transform reads where p = k − ω 0 + iΓ/2(1 − α) and γ(n, z) is the incomplete Gamma function [75]. The corresponding solution for φ(x, t) follows straightforwardly by using (23) in Eq. (21).

B. Semi-infinite waveguide: two-excitation sector
The ansatz for the time-dependent wavefunction [cf. Eq. (5)] reads The Schrödinger equation then yields the system of coupled differential equations (25b) The formal solution for χ(x 1 , x 2 , t) is thus where note that χ is symmetrized under the exchange x 1 ↔ x 2 . By placing Eq. (26) into Eq. (25a) we find a spatially non-local delay partial differential equation (PDE) for ψ(x, t): (27) is the two-excitation-sector counterpart of the DDE (22). Mathematically, such a spatially non-local delay PDE is far more involved than the DDE (22) or conventional delay PDEs [76] that are local in space. A spacetime diagram is shown in Fig. 3.
In our case [cf. Eq. (5)], the initial conditions are ψ(x, 0) = ϕ(x) and χ(x 1 , x 2 , 0) = 0. Due to the latter condition, the terms on the last line of Eq. (27) are identically zero. Hence, overall, the differential equation features four source terms that are non-local in x and t and are non-zero for x > −a. In the region x ≤ −a, the equation takes the simple form By taking the Fourier (Laplace) transform with respect to variable x (t), this equation is turned into an algebraic equation whose solution is given bȳ whereφ(q) = αΓ 2π e iqa /(k − q − i αΓ 2 ) is the Fourier transform of ϕ(x). Performing the inverse Fourier transform with respect to q then yields where we used that, since x < −a, only the pole q = k − iαΓ/2 contributes to the integral. Upon inverse Laplace transform with respect to s term by term, we finally find for x < −a. This solution can be expressed compactly as ψ(x, t) = ϕ(x − t)e sm (t), where e sm (t) is the qubit excited-state amplitude in the spontaneous emission process [23,73,74] namely the solution of Eq. (22) for the initial conditions e(0) = 1 and φ(x, 0) = 0. This is physically clear: since the qubit starts in the excited state [cf. Eq. (5)] so long as the photon has not reached its location x = −a the system's evolution consists of the free propagation of the input single-photon wavepacket and the spontaneous emission as if the field were initially in the vacuum state. The next natural step would be finding the wavefunction for −a ≤ x ≤ a. However, a look at Eq. (27) shows that such a task is non-trivial. Specifically, two of the source terms, ψ(−x − 2a, t − x − a) and −ψ(−x, t − x − a), enter the differential equation which forces one to find the solution "tile by tile" as discussed in the Supplementary Material [77], a challenging and in the end impractical task. We choose instead to solve the delay PDE numerically by adapting the finite-difference-time-domain (FDTD) method [78][79][80]; our approach is described in Ref. [80]. Note that the effectiveness of our code is crucially underpinned by the knowledge of the exact solution for x < −a discussed above [80].
Finally, once the solutions in both number sectors are found, the three functions p e/g (t) and c(t) can be obtained explicitly as where e(t) is given by Eq. (23), φ(x, t) is given by Eq. (21), and ψ(x, t) is obtained from FDTD.

VI. NON-MARKOVIANITY
Despite having a Lindblad structure, the time-dependent ME (13) is not in general a Lindblad ME, not even one whose Lindblad generator is time-dependent, because the rates {γ ν (t)} are not necessarily all positive at all times [34,35,38,81]. The condition γ ν (t) ≥ 0 for any ν and t is indeed violated if ∆(t) < 0 at some time t [recall the definition of ∆(t) in Eq. (9)].
Indeed, since ∆(0) = 1, if ∆ becomes negative during the time evolution then there exists an instant at which both∆ < 0 and ∆ < 0. Then, since γ + + γ − = −∆/∆ [see Eq. (13) and related text], at least one of the rates {γ ± (t)} must be negative at some time. When this happens, the DM is not "CP-divisible": the dynamics cannot be decomposed into a sequence of infinitesimal CP maps [34] each associated with a time 0 ≤ t ≤ t and fulfilling (13) with positive rates {γ ν (t )}. Equivalently, it is not governed by a Lindblad ME even locally in time [82]. Thus, according to the criteria in Refs. [42,81], the dynamics is NM. Note that the negativity of ∆(t) is also sufficient to break P-divisibility (a weaker property than CP-divisibility) since it ensures that the sum of at least a pair of time-dependent rates in ME (13) is negative [38]. Negativity of ∆ can occur already for an infinite waveguide if the wavepacket width is in an optimal range. Indeed, from the analytic expressions for p e/g (t) in the infinite-waveguide case given in Appendix B 3, we get (time in units of Γ −1 and k = ω 0 ) Based on elementary analysis, two behaviors are possible (Appendix B 4): for α > 5, ∆(t) ≥ 0 always, while for 0 < α ≤ 5, ∆(t) has a minimum at a negative value at some time. To illustrate this, we plot ∆'s negativity, N ∆ (t) ≡ −min[0, ∆(t)], in Fig. 4. For 0 < α ≤ 5, N ∆ (t) is initially zero, then exhibits a maximum at a time of the order of Γ −1 and eventually decays to zero. For α 10 −2 , this maximum is in fact negligible reaching at most N ∆ ∼ 10 −6 [ Fig. 4(b)]. For practical purposes, then, ∆(t) becomes negative for an optimal range of wavepacket widths δk around α 1, i.e., δk Γ, which excludes small α and hence in particular quasi-plane waves.
Rigorously speaking, it should be noted that -as typically happens with time-convolutionless MEs [38] -in general there may be singular times at which ME (13) is not defined and correspondingly the DM not invertible. In the present case, these are the times at which ∆(t) and/or c(t) vanish [cf. Eqs. (9), (16), and (B13)]. The above sufficient condition should thus in general be complemented with the additional requirement that c(t) does not simultaneously vanish. This is always the case in Fig. 4, which is easily checked with the help of the analytical expression for c(t) [Eq. (B13c)].
Further light on the onset of NM effects can be shed by studying in detail a non-Markovianity measure, which by definition is a function of the entire DM (i.e., at all times) [38]. Out of the many proposed [38], we select the geometric measure (GM) [62] for its ease of computation and because it facilitates a comparison with the spontaneous emission dynamics in the semi-infinite waveguide where the GM was already used [26]. The GM is defined in terms of the DM's determinant as [62] where the integral is over all the time intervals in which | det M t | grows in time, and [cf. Eq. (8)]. Note that N depends on the modulus of the determinant, which is the volume of the ellipsoid into which the Bloch sphere is transformed by the DM [see Eq. (10)]. Hence, a non-zero N means this volume increases at some time, in contrast to dynamics described by the Lindblad ME in which such an increase cannot occur [62]. It is known [38] that a non-zero GM implies that the dynamics is NM also according to the BLP measure [41], which in turn entails NM behavior according to the RHP measure [42]. A remarkable property following from Eqs. (34) and (35) is that if there exists a time such that ∆(t) < 0 and c(t) = 0 then det M t must grow at some time. This then brings about that the dynamics is NM according to the GM (34) and hence NM even according to the BLP and RHP measures. We thus in particular retrieve the sufficient condition for breaking P-and CP-divisibility discussed at the beginning of this section since non-zero BLP (RHP) measure ensures violation of P-divisibility (CP-divisibility) [38]. Figure 5(a) shows N for the infinite-waveguide case for a wavepacket carrier frequency resonant with the qubit (solid line). Similarly to the negativity of ∆, N takes significant values only around α 1 (i.e., δk Γ), being negligible in particular for quasi-monochromatic wavepackets. The values of α yielding N = 0 are also such that the dynamics is NM according to the BLP measure [41] and even the RHP measure [42], the latter meaning that rates γ ±,z (t) in the ME (13) break the condition of being positive at all time.
The behavior of N changes substantially for a semi-infinite waveguide, as shown in Fig. 5 for several values of the qubit-mirror distance a. First, non-Markovianity is generally larger, even by an order of magnitude in some cases [note the difference in scale between panels (a) and (b)]. Second, N can be significant even at α 0 (the plane-wave limit), in which case it matches its value in the corresponding spontaneous-emission process [26]. For our parameters, the maximum non-Markovianity in this limit occurs near k 0 a = 4π. Third, N exhibits a more structured behavior as a function of α, the shape of N (α) being dependent on k 0 a (recall k 0 = ω 0 ).
The feedback due to the mirror, evident in Eqs. (22) and (27), generally introduces memory effects in the qubit dynamics that are expected to cause NM behavior. These add to the finite-wavepacket effect already occurring with no mirror, leading generally to enhanced non-Markovianity-note that the semi-infinitewaveguide curves in Fig. 5 typically lie above the mirrorless one. The non-Markovianity can be especially large when either 2k 0 a, the phase corresponding to a qubit-mirror round trip, is an integer multiple of 2π and/or the corresponding photon delay time τ = 2a is large compared to the atom decay time Γ −1 . In the former case, enhanced NM behavior occurs because a standing wave can form between the mirror and the qubit under these conditions; indeed, it has been shown that a bound state in the continuum exists in this system [23]. In the latter case, the fact that the qubit decays completely before the photon returns causes a periodic re-excitation of the qubit-a kind of revival.
We found numerically that the scattering DM (8) reduces to that for spontaneous emission [26] in the limits of very large and very small α, thus in particular explaining the behavior of N at α 0. In the infinitewaveguide case, this property can be shown analytically (see Appendix B 3). Physically, these limits can be viewed as follows. When α 0, the wavepacket is so spread out spatially that the photon density at the qubit is negligible: the qubit effectively sees a vacuum, hence behaving as in spontaneous emission. This clarifies why NM effects cannot occur without the mirror for a quasi-plane-wave (the emission DM in an infinite waveguide is clearly Markovian). When α 1 in contrast, the photon is very localized at the qubit position. The energy-time uncertainly principle then implies that the photon passes too fast for the qubit to sense, hence the qubit again behaves as if the field were in the vacuum state.
Since non-Markovianity measures are generally not equivalent [38] it is natural to wonder whether the outcomes of our analysis in Fig. 5 for the GM hold qualitatively if a different measure is used, for instance the widely adopted BLP measure [41]. While a comprehensive comparative study of different measures is beyond the scope of the present paper, we computed the BLP measure N BLP for some representative values of the parameters. For an infinite waveguide, the behavior of N BLP as a function of α is analogous to that of the GM [see Fig. 5(a)]. In the semi-infinite-waveguide case, N BLP overall behaves similarly to the GM but exhibits a less structured shape: for instance, the inflection point for k 0 a = 0.5π in Fig. 5(a) is absent.

VII. CONCLUSIONS
We studied the open dynamics of a qubit coupled to a 1D waveguide during single-photon scattering, presenting results for its DM, the corresponding time-dependent ME, and rigorous non-Markovianity measures developed in OQS theory. The qubit dynamics was shown to have distinctive features that, in particular, do not occur in emission processes. To compute the DM for a semi-infinite waveguide, we solved the scattering time evolution by deriving a spatially non-local delay PDE for the one-photon wavefunction when the qubit is excited. For an infinite waveguide, NM behavior occurs when the photon-wavepacket bandwidth δk is of order the qubit decay rate Γ. For a semi-infinite waveguide (mirror), time delay effects are an additional source of non-Markovianity, resulting in generally stronger NM effects.
The system we studied here, a semi-infinite waveguide plus a qubit, is the simplest waveguide QED system with a time delay. Yet the nature and effects of the time delay should be completely generic as there is no fine tuning in our system. We thus expect these main conclusions to also hold in, for instance, the case of two distant qubits coupled to a waveguide which is relevant for long-distance quantum information.
It is interesting to note [see Eq. (35)] that ∆(t) has the same sign as det M t , hence times can exist at which det M t < 0. Among qubit CP maps, those with negative determinant are the only ones that break the property of being "infinitesimally divisible" [83]. This class does not include spontaneous-emission DMs-in particular vacuum Rabi oscillations-where the determinant is always non-negative [62]. In sharp contrast, the scattering DMs studied here do belong to this class.
Finally, we note that some results here rely solely on the DM structure (8) that in turn stems solely from having an initial Fock state for the field and the rotating wave approximation. Further investigation of this class of open dynamics is under way [84]. In this Appendix, we present some details of the derivation of the time-dependent ME, Eq. (13). In particular, we express both Eqs. (13) and (15) using as a basis the four Hermitian operators 2} with i = 0, 1, 2, 3, respectively; recall that Tr{Ĝ iĜj } = δ ij . An operator ρ (and so in particular a density operator) can be decomposed as ρ = 3 i=0 r iĜi with r i = Tr{Ĝ i ρ}, hence the 4-dimensional real vector r is a representation of the density operator ρ. A map is analogously represented by a 4×4 transformation matrix.
For L t we start by noting that the dynamical map can be expressed as where we used the linearity of Φ t and defined the entries of the 4×4 matrix F as The matrix F thus represents the map Φ t (we drop the time dependance for simplicity). The composition of two maps is correspondingly turned into the matrix product of the associated matrices [each defined analogously to Eq. (A2)]. Hence, if L is the 4×4 matrix associated with map L t [see Eq. (15)], it is given by We are thus led to compute the (time-dependent) matrix F, calculate its derivativeḞ and inverse F −1 , and finally take the matrix product (A3). To calculate F we use Eqs. (8) and (A2), where the matrix elements ρ ij entering Eq. (8) are now the matrix entries of operators {Ĝ j } (for instance,Ĝ 1 =σ x / √ 2 has entries (G 1 ) ee = (G 1 ) gg = 0 and (G 1 ) eg = (G 1 ) ge = 1/ √ 2). By proceeding in this way, matrix F reads Using this and Eq. (A3), we find that matrix L is This shows that Eq. (14) holds with the generator L t whose 4×4-matrix representation is given by Eq. (A5). The remaining step is to show that the generator can indeed be expressed as the right-hand side of (13). To this aim, we consider Eq. (13) without specifyingĤ(t), γ ± (t) and γ z (t), work out its 4×4-matrix representation, impose that it yieldsṙ = Lr with L given by Eq. (A5) and solve forĤ(t), γ ± (t) and γ z (t). Thus let us definẽ and callL the associated 4×4 matrix. To computeL, in Eq. (A6) we replace ρ = i r iĜi and calculate Tr{ÂĜ iB } withÂ,B =Ĝ 0 , ...,Ĝ 3 , obtaining We next require that Eq. (A7) equals Eq. (A5). Upon comparison of these two equations, we immediately get µ(t) = 0, while S(t) = −2 Im [ċ(t)/c(t)]. Moreover, by requiring the entries L 22 , L 44 and L 41 of matrix (A7) to match the corresponding ones of (A5), we find the three rates given in Eq. (16).

Appendix B: Calculations for the infinite-waveguide case
Here, we present details of the calculation of the time-dependent wavefunctions in both the one-and twoexcitation sectors that are needed for the explicit calculation of the dynamical map Eq. (8) in the infinitewaveguide case. Following the main text, we refer to ϕ(x) as a single-photon exponential wavepacket of the form Eq. (17). Further technical details, including the study of other possible initial conditions, are given in the Supplementary Material [77].

One-excitation sector
This is the scattering process corresponding to Eq. (4) in the infinite-waveguide case, based on which the ansatz for the time-dependent wavefunction reads where φ R/L (x, t) is the wavefunction of the right-/left-going photon. Imposing the Schrödinger equation, i∂ t |Ψ 1 (t) =Ĥ|Ψ 1 (t) , yields the three coupled equations The equations for φ R/L (x, t) can be formally integrated by Fourier transform, yielding where we set θ(0) ≡ 1/2. The first term on each righthand side describes the free-field behavior, while the second one can be interpreted as a source term originating from qubit emission at an earlier time. Note that causality is preserved as it should be. Eqs. (B3) immediately entail φ R (0, t) + φ L (0, t) = φ R (−t, 0) + φ L (t, 0) − iV e(t), which once substituted in Eq. (B2c) yields a time-local first-order differential equation for e(t) Imposing the initial conditions φ R (x, 0) = ϕ(x), φ L (x, 0) = e(0) = 0, we obtain, By using Eq. (B5) in Eqs. (B3), one then obtains the photon wavefunctions φ R/L (x, t).

Two-excitation sector
Based on Eq. (5), the ansatz for the time-dependent wavefunction reads where ψ R/L (x, t) is the probability amplitude to have a right-/left-propagating photon at position x with the qubit in the excited state, while χ αβ (x 1 , x 2 , t) is the probability amplitude to have an α-propagating photon at position x 1 and a β-propagating photon at position x 2 (with the qubit unexcited). Terms ∝ χ LR have been incorporated in those ∝ χ RL by exploiting the symmetrization property The Schrödinger equation then yields five coupled differential equations that read Note that the equations for χ RR and χ LL are symmetrized because of the bosonic statistics. Similarly to the previous subsection, we first formally solve for the purely photonic wavefunctions and find Next, we plug these solutions back into Eqs. (B7a) and (B7b), which are those featuring the qubit degree of freedom, under the initial condition that χ αβ (x 1 , x 2 , 0) = 0 for any α, β = L, R. The resulting pair of (B9b) These coupled differential equations are non-local with respect to both x and t, the non-locality being due to the rightmost "source terms" that feature the double step functions. Based on the arguments of the step functions, it is natural to partition space-time into the three regions , and x > t (R 3 ) in the case of Eq. (B9a) and x ≤ −t (L 3 ), −t < x ≤ 0 (L 2 ), and x > 0 (L 1 ) in the case of Eq. (B9b), as shown in Fig. 6. Then, the differential equations (B9) can be analytically solved in four steps as follows: (i) Solve Eq. (B9a) for ψ R (x, t) in region R 1 under the initial (i.e., boundary) condition ψ R (x, 0) = ϕ(x). In this region, the source term is identically zero.
(iii) Solve Eq. (B9a) for ψ R (x, t) in region R 2 under the boundary condition ψ R (0, t) [this being fully specified by the solution found at step (i)]. In this region, the source term is non-zero but is fully specified by the solutions ψ R (x < 0, t) and ψ L (x > 0, t) worked out at the previous steps (i) and (ii), respectively. Note that the initial condition automatically guarantees that the wavefunction is continuous at x = 0.
(iv) Solve Eq. (B9b) analogously for ψ L (x, t) in region L 2 under the boundary condition ψ L (0, t) = 0. In this region, the source term is again fully specified by the solutions ψ R (x < 0, t) and ψ L (x > 0, t) obtained in the previous steps.
Finally, ψ R (x, t) vanishes identically in region R 3 and, likewise, so does ψ L (x, t) in region L 3 . This is because causality prevents the wavefunction outside the light cone from being affected by the qubit or input wave. Since initially the wavefunction is zero in this region, it remains so at all times. Hence, the wavefunction is non-zero only in regions R 1 , R 2 and L 2 .
With the help of Mathematica (see Supplementary Material [77]), the above procedure straightforwardly yields analytical expressions for the wavefunctions. We checked that, in the steady-state limit t→∞, the above solution for the wavefunctions in the stimulated-emission problem yields results in full agreement with those obtained via a time-independent approach [61]. In particular, the two-photon scattering outcome probabilities P RR , P RL and P LL of Ref. [61] are recovered as with α, β ∈ {R, L}. We finally mention that, in the case of an incoming two-photon wavepacket (not addressed in the main text), one or more terms χ αβ (x 1 , x 2 , 0) are non-zero and Eqs. (B9) feature additional terms. For instance, in the case of a left-incoming two-photon wavepacket, the additional term −i Γ/4 χ RR (x − t, −t, 0) + χ RR (−t, x − t, 0) must be added to the right-hand side of Eq. (B9a). In this case, in the steady-state limit t → ∞ known results for two-photon scattering (in particular second-order correlation functions) [25,29,85,86] are recovered, which confirms the effectiveness of our real-space time-dependent approach.
The time derivative of function (33) For α =1, at a stationary point of ∆(t), thereby, curves f (t) and g(t) cross. Note that the positive functions f (t) and g(t) both monotonically increase with time and so do all their derivatives. Thus there exist either zero or only one crossing point, whose occurrence depends on whether f (t) is above or below g(t) at t = 0 and t → ∞. A simple calculation yields .
By noting that f (0)/g(0) > 1 for α < 1 and f (0)/g(0) < 1 for α > 1, we see that three cases occur. For α < 1, f (t) is above g(t) at t = 0 and below it at t → ∞, hence a single crossing point occurs. For 1 < α ≤ 5, f (t) lies below g(t) at t = 0 and above it at t → ∞, hence a single crossing point occurs in this case as well. Finally, for α > 5, f (t) lies below g(t) both at t = 0 and at t → ∞, hence no crossing points occur. Function ∆(t) thereby has a single stationary point for 0 ≤ α ≤ 5 and none for α > 5. One can show that this stationary point is indeed minimum, concluding the proof. Note that we have now included the case α = 1 since this yields ∆(t) = e −2t e t (3 − 2t) − 2 , which exhibits a single stationary point that is a minimum.