Markovianity of an emitter coupled to a structured spin chain bath

We analyze the dynamics of a spin 1/2 subsystem coupled to a spin chain. We simulate numerically the full quantum many-body system for various sets of parameters and initial states of the chain, and characterize the divisibility of the subsystem dynamics, i.e. whether it is Markovian and can be described by a (time dependent) master equation. We identify regimes in which the subsystem admits such Markovian description, despite the many-body setting, and provide insight about why the same is not possible in other regimes. Interestingly, coupling the subsystem at the edge, instead of the center, of the chain gives rise to qualitatively distinct behavior.


I. INTRODUCTION
The subject of open quantum systems (OQS) focuses on the description of quantum systems coupled to a (typically much larger) environment. In general, solving the time evolution of the total system is out of reach due to the macroscopic number of environmental degrees of freedom and the exponentially large Hilbert space. Instead one tries to obtain an effective reduced description which involves the degrees of freedom of the OQS only. In this context, the distinction between Markovian and non-Markovian dynamics is a central theme. [1][2][3] Originally, the former denoted situations that allowed for the derivation of a 'Markovian' master equation. 4 This is a specific differential equation generating the dynamics of the OQS that has the property of being 'memoryless' (Markovian), in the sense that the evolution of the OQS at a given time depends only on its state at that time.
The same OQS perspective can be applied to study the evolution of a subsystem in a closed many-body quantum system. For instance, although the full system is in a pure state, the computation of local observables requires only to know the state of a (small) subsystem, for which the rest of the system would play the role of environment. Again, solving the long time evolution of the full system is in general out of reach due to build-up of entanglement, such that it might be desirable to find an effective reduced description for the subsystem (OQS) only. However, the standard 'Markovian' master equation derivation is based on weak coupling and a separation of time scales between open system and environment; 5 two conditions that are generally not fulfilled in the quantum many-body setup. It is thus interesting to analyze whether there are cases in which the dynamics still admits a reduced description in this setting, and, in case such a description exists, whether it is Markovian.
In recent years a variety of non-Markovianity measures were put forward that go beyond the original Markovianity conditions mentioned above and attempt to quantify deviations from Markovian dynamics from a quantum information theoretical perspective.  Two of the most widely used are the one introduced by Breuer, Laine, Piilo 7,31 (BLP measure), which detects the non-monotonicity of the trace distance between pairs of states evolving in time, and the more stringent one introduced by Rivas, Huelga, Plenio 8 (RHP measure), which detects non-divisibility of the quantum channel mapping initial OQS states to their time-evolved states. These measures are not equivalent, since there exist cases which are characterized as BLP-Markovian but RHPnon-Markovian. 4,12,[32][33][34] For cases with Markovian dynamics according to RHP, a Markovian reduced description in the above sense exists, i.e. the (time dependent) equation governing the reduced dynamics does not explicitly depend on past system states and is called 'time dependent Markovian' master equation. 8, 35 We quantify non-Markovianity by its robustness, as originally introduced for 'snapshots' of quantum evolution in reference 6 (and here generalized to continuous evolution), i.e. how much noise can be added to the dynamics before RHP Markovianity is recovered, thus providing a physical interpretation.
In this paper we explore the above questions in the particular case of a spin (the OQS) coupled to a XY spin chain, which plays the role of environment. In particular, we identify regimes that allow for a description via a 'time dependent Markovian' master equation and provide insight into what prevents such a description. We consider two scenarios in OBC: (i) spin coupled to the center; (ii) spin coupled to the first site; and we analyze different initial states of the chain. In some particular cases (namely in scenario (ii), 36 and for (i) in the thermodynamic limit 37 when the chain is initially in the vacuum) an exact solution is possible. In more general cases, we use tensor network methods (MPS and MPO, as we are in 1D) to simulate the evolution of the full system.
The initial state of the chain can be empty (i.e. in the fermionic vacuum of the XY chain) or contain 'excitations' (if some of such modes are occupied). We find that while in the first case, BLP and RHP Markovianity is equivalent, in the second case, only the divisibility (RHP) measure detects all the non-Markovianity appearing in regions of the parameter space. One possibility we explore to populate fermionic modes in the initial environment state is using thermal states. Whilst small temperature induces additional non-Markovianity, we find arXiv:1912.09151v1 [quant-ph] 19 Dec 2019 for scenario (i) that increasing the temperature gradually removes the non-Markovianity until at high temperature the dynamics is captured by a 'time-dependent Markovian' master equation. This applies even at the band edges where the spectral density diverges, a scenario that is often associated with strong non-Markovian behavior. [38][39][40][41] In contrast, in scenario (ii) we find that any RHP non-Markovianity of the vacuum case survives at all temperatures. We show that this remarkable difference between the two cases can be anticipated from the different decay of the environment correlation functions at high temperature in both cases.
The paper is structured as follows: we start in II by introducing the non-Markovianity measures and the 'time dependent Markovian' master equation. In III we present the details of the model and the specific form the non-Markovianity measure adopts in this model. We also discuss the conditions that allow for the derivation of a standard 'Markovian' master equation. The section closes with a summary of our numerical methods. In IV we review the analytically solvable case in which the spin is coupled to the center of the chain initialized in the vacuum. In V we introduce initial environmental excitations in this setup which yields the main results of this paper. In VI we explore the qualitatively different nature of the dynamics obtained by coupling the spin to the first site in the chain and provide insight into why the two scenarios differ so much with respect to non-Markovianity. Finally, in VII we conclude and summarize the main results.

II. QUANTIFYING NON-MARKOVIANITY
Over the past decade multiple inequivalent characterisations of non-Markovianity based on quantum information theory have been introduced. 6-30 In the next paragraphs we review the measures that are relevant for the rest of our paper.
A. Non-Markovianity robustness and 'time-dependent Markovian' master equation We consider the set T of finite dimensional, completely positive (CP), trace-preserving, linear maps (quantum channels) T : The time evolution of finite dimensional quantum systems is given by a one parameter family of quantum channels, known as dynamical map, that maps the initial state of the system, described by the density matrix ρ(0), to the time evolved state ρ(t) = T (t)[ρ(0)]. Denoting by T (t 2 , t 1 ) the map for the evolution from time t 1 to t 2 , we have, by continuity, T (t + ∆t) = T (t + ∆t, t)T (t) and thus T (t + ∆t, t) = T (t + ∆t)T (t) −1 . (1) If for all > 0 there is a finite ∆t < such that these maps are CP for all t, T (t) is infinitesimal divisible (T (t + ∆t, t) → 1 as ∆t → 0). Note that this condition is more restrictive than mere infinitesimal divisibility of a 'snapshot' of T (t) at a given time, since our decomposition needs to follow the dynamics at all times.
To be consistent with the recent literature we drop the word 'infinitesimal' and call the time evolution divisible or (RHP 8 ) Markovian if the above is true. Notice that for T (t + ∆t, t) to be defined unambiguously, T (t) −1 needs to exist. Since T (t → 0) → 1, for t small enough T (t) will be invertible. For later times, T (t) −1 may not exist, in which case one lacks essential information as a consequence of being blind to the environment part of the whole system. In that case one can resort to pseudoinverse techniques. 8,42 For the remainder of this paper ∆t denotes a small (but finite) time step. We can represent the map by a matrix 35 where {F α } α=1,...,d 2 is an orthonormal basis in M d . We will use the canonical basis {|i j|} i,j=1,...,d . In the rest of the paper by T and dT we denote the matrix representations of T (t) and T (t + ∆t, t) respectively and omit their time dependence for convenience. The density matrix can be written as a linear combination of this basis with components ij|ρ , where A, B = tr A † B is the Hilbert-Schmidt scalar product and we identify |ij ↔ |i j|. Given a map T on a d-dimensional space, its Choi state 43,44 is where ω is a maximally entangled state ω = |ω ω|, |ω = Divisibility is equivalent to ρ(t) being a solution of a time-dependent Lindblad ('time-dependent Markovian') master equation: 8,35 where γ i ≥ 0 and L i are called rates and Lindblad operators of the (time-dependent) Lindbladian L(t), and H(t) = H † (t). If the dynamics is described by a timedependent Lindblad master equation, the dynamical map T can be decomposed into infinitesimal 'pieces' dT = e Ldt , where L is the matrix representation of L(t). Note that the corresponding Choi state L Γ is hermitian. We can now quantify non-divisibility of the reduced dynamics by calculating how much log dT deviates from a valid Lindbladian. Following reference 6 this amounts to checking if (a) (log dT ) Γ is hermitian and (b) if there exists a branch of the logarithm such that ω ⊥ (log dT ) Γ ω ⊥ ≥ 0, where ω ⊥ = 1 − ω is the projector onto the orthogonal complement of the maximally entangled state. If dT is hermiticity-preserving (a necessary condition for it to be CP), its eigenvalues are either real or come in complex conjugate pairs. The set of hermitian logarithms of dT is then parametrized by a set of integers m c ∈ Z, L m = L 0 + 2πi c m c (P c − P c ), where L 0 denotes the principal branch and P c and P c are projectors onto the eigenspaces associated to a complex conjugate pair of eigenvalues λ c and λ c of dT . We define If L Γ 0 and thus A 0 is hermitian, the dynamics is Markovian iff for any time there exists {m c } such that If (4) is not satisfied during a given time step, adding noise may remove the non-Markovianity. In 6 the non-Markovianity is measured by its robustness, i.e. by the minimum amount of isotropic noise µ that achieves this: is then a valid Lindbladian. Note that if dT has some real negative eigenvalue, L Γ 0 is non-hermitian and log dT cannot be made a valid Lindbladian by adding a finite amount of noise. If dT is hermiticity-preserving and does not have real negative eigenvalues, we assign robustness according to Eq. (5), otherwise µ = ∞. µ > 0 at some time implies that the evolution cannot be described with a valid 'time-dependent Markovian' master equation, which is equivalent to non-divisibility. It is thus obvious that this provides a necessary and sufficient criterion to decide about Markovianity.
In practice we compute T at discrete times separated by time steps ∆t until a final time t fin = K∆t. We then compute the minimum noise µ (n) required to make the n−th time step Markovian. In order to compute a non-Markovianity measure for the whole evolution interval, we choose to average: µ = 1 K K n=1 µ (n) . We then use the normalized degree of non-Markovianity 6 where we have N ∈ [0, 1] and the dynamics is Markovian until time t fin iff N = 0. We choose ∆t sufficiently small such that N is converged with respect to the time step. Notice that this non-Markovian measure depends on the final time t fin .
Note from Eq. (5) that µ (n) only depends on the most negative eigenvalue of A (n) 0 are chosen in a way to minimize the magnitude of the most negative eigenvalue. We can gain further insight into the dynamics by looking at the full spectrum. We are interested only in the non-zero eigenvalues λ (n) i ∈ R and corresponding eigenvectors v If γ i . 45 In the limit ∆t → 0 the time dependent equation can be recovered. This is a 'time-dependent Markovian' master equation if the rates are non-negative at all times. If we denote the minimum rate and its corresponding operator by γ min and L min respectively (omitting their time dependence for convenience), we have and thus the measure is non-zero if any rate ever becomes negative. The rates contain the required robustness information whilst providing additional insight into what kind of process is responsible for making a given time step non-Markovian.

B. BLP measure
Another broadly used measure of non-Markovianity is the BLP measure, 7,31 which is based on the study of the time behavior of the trace distance D between pairs of density matrices evolving in time. For a pair ρ 1,2 it is defined as D ρ 1 , ρ 2 = 1 2 tr ρ 1 − ρ 2 . The idea is that since the trace distance is contractive under CP maps, Markovian processes, described by 'Markovian' or 'timedependent Markovian' master equations, cannot increase it during time-evolution.
The rate of change of the trace distance is defined as σ t, ρ 1,2 (0) = d dt D ρ 1 (t), ρ 2 (t) . The BLP measure is defined as Notice that the definition involves a maximisation over initial pairs of states ρ 1,2 (0). A discretized version, in line with our construction in the previous section, uses a summation over time steps ∆t, in which the trace distance has increased: D ρ 1 (t + ∆t), ρ 2 (t + ∆t) − D ρ 1 (t), ρ 2 (t) > 0. Whilst the BLP measure does not require knowledge of the dynamical map, the maximisation over initial states in general cannot be performed exactly. The trace distance allows for an immediate information theoretic interpretation in the sense that D ρ 1 , ρ 2 = 1 if the states are perfectly distinguishable while it is zero if they are identical. 46,47 The BLP measure exploits this by interpreting an increase of the trace distance as information backflow from the environment into the system making the states more distinguishable. Such a backflow is then identified as non-Markovian as the time-evolution of the states at that time depends on information about the states that flowed into the environment at previous times ('environment keeps memory'). We call the dynamics BLP-Markovian if such a back flow never occurs N BLP = 0. RHP (divisibility) implies BLP Markovianity, but not the other way round. 4,12,[32][33][34] In particular, the family of P-divisible dynamics, for which T (t + ∆t, t) is positive, but not necessarily CP, is BLP-Markovian. 7,31

III. SETUP, MODEL AND MEASURE
As open quantum system we consider a single spin 1/2 (S), coupled to a spin chain environment (E) of length N governed by a Hamiltonian of XY type: where σ µ m (µ ∈ {x, y, z}) are Pauli operators acting on site m and h is an external magnetic field in the zdirection. The system is coupled to the m 0 -th spin of the chain via an exchange interaction of strength Ω where τ µ are simply the Pauli operators acting on the spin. Finally, the system Hamiltonian is where τ ± = 1 2 (τ x ± iτ y ). The environment is exactly solved in terms of diagonal fermionic modes d k with energies E k = 2J cos πk N +1 + 2h (k = 1, . . . , N ). 48 In the continuum limit this gives an energy band from 2h − 2J to 2h + 2J with diverging density of states at the edges. The field h acts like a chemical potential that allows us to move the band up and down in energy, hence, the detuning is defined as ∆ h = ∆ − 2h. We consider different initial states for the environment, corresponding to either ground states at different magnetic fields, or to thermal equilibrium states.

A. Spin dynamics
Under the assumption of no initial system-environment correlations, the dynamics of an OQS is given by the dynamical map where ρ(0) and ρ E are the initial states of system and environment respectively and H is the total Hamiltonian. T (t) can be obtained using quantum process tomography, 49 which requires knowledge of ρ(t) for a number of different initial (pure) states ρ(0).
The total Hamiltonian H = H E + H S + H SE conserves the total spin along the z-direction or, in terms of the fermion operators of the chain, In a slight abuse of notation, we will refer to this conserved quantity as the number of 'excitations' in the full system. Furthermore, we consider initial states of the environment which commute with the total number of fermions. Thus, using Eq. (14), the channel takes the following form where a and c are the excited state populations at time t for the open system initialized in the excited and ground state respectively and b evolves the coherences: e|ρ(t)|g = b e|ρ(0)|g . We omit the time dependence of the channel elements for notational convenience. Here we labeled the basis elements i = e, g, corresponding to excited and ground states of H S . The block structure of T translates in a straightforward way to dT .
Following the procedure explained in section II and exploiting the structure of dT , we find the following. First, A c = 0, which simplifies Eqs. (4) and (5). Also, the eigenvectors v i are independent of time and the reduced system dynamics is described by the following time-dependent differential equation: where γ i (t) 50 is real (a and c take values between 0 and 1) and in Eq. 17 we take the principal branch. The coherent part of the equation with the Lamb shift energy E LS is obtained after subtracting the dissipative part from log dT .
Since we identify non-Markovianity with the occurrence of negative rates, in the following we focus on these expressions. In particular, for a − c ≥ 0, the dynamics is non-Markovian if the time-derivative of at least one of the fractions inside the logarithms is negative.

B. Conditions for deriving a 'Markovian' master equation for the spin
It is interesting to compare the previous constructions with the usual steps followed in the quantum optics formalism in order to derive a 'Markovian' master equation, which is a (time-independent) Lindblad equation.
The standard derivation assumes a separation of time scales between system and environment. In the simplest case, where the system is coupled to a single environment operator R, the environment time scale is characterized by the correlation time τ c after which the environment correlation function α(τ ) = tr ρ ER (τ )R(0) , with interaction picture operatorsR(t), has decayed. 5 As an aside, by system time scale it is typically meant the time scale on which the system changes due to its interaction with the environment. Thus, it depends on the intensity of the system-environment coupling Ω. The standard derivation is valid if the following condition is satisfied: 5 This means that the coupling has a weak effect during the correlation time of the environment fluctuations or, in other words, the reduced system dynamics at a given time only weakly depends on previous system states (Markovianity) since the environment memory of those persists only on the τ c time scale. 5 We now discuss how the validity of the 'Markovian' master equation can be checked in our model. The typical starting point for its derivation is to express the reduced system dynamics by a specific integro-differential equation, which is valid to second order in the coupling (Born approximation). 1 For our spin it reads: 38 To simplify this equation we use the Jordan Wigner transformation to map the spin chain to free fermions: and c i are fermionic operators. In this language the interaction Hamiltonian reads Eq. (22) can then be rewritten as follows: 1,51 where are the environment correlation functions in our model. From Eq. (24) we see that if the kernels Re(α + (t)e −i∆t ) and Re(α − (t)e i∆t ) decay sufficiently fast, i.e. Ωτ c 1, the integrand has decayed before the coupling has had a significant effect on the evolution of the spin such that we can replace ρ(t − s) with its unperturbed evolution. The equation is then time-local in the sense that it does not depend on the state of the spin at previous times ('memoryless') and can be further manipulated until the 'Markovian' master equation is obtained.
In Appendix D we show how to compute the correlation functions for thermal and ground states of the chain. They can be written α ± (t) = N k=1 e ±iE k t α ± k (t) such that the kernels of Eq. (24) are: with α ± k (t) given in Eqs. (D4) and (D5). In certain regimes of our model, namely, when the environment is initially in the vacuum or when the spin is coupled to the first site m 0 = 1, the string operator is the identity u m0 = 1, and the coefficients adopt a time-independent form: α ns+ 48 and f k = 1 1+e βE k is the Fermi-Dirac distribution. We labeled such cases 'ns' for 'no string operator'.
In such case, we can check condition (21) using a self-consistency argument. In the thermodynamic limit and going to energy space W m0,k → W m0 (E), we introduce the spectral density , and an analogous expression for α ns+ (t). Now the argument proceeds as follows. We assume that α ns± (E) is flat around ∆ over an energy range set by the characteristic frequency of the spin Γ (e.g. dissipation rate). Within the integral over past times s, Eq. (24), it can then be replaced by its value at ∆, and, using the relation where P denotes the Cauchy principal value, we can replace the respective kernel by 38 Hence, under that assumption, the spin dynamics is captured by the 'Markovian' master where we set Γ = max(Γ ns± ). 38

A MPS for an open boundary system of N sites with physical dimension d and local basis
m are D × D matrices, except for the first and the last, which are 1×D and D ×1 vectors respectively. The bond dimension D sets the number of free parameters in the ansatz. 52-54 MPS yield good approximations to ground states of gapped, local Hamiltonians. 55,56 Efficient numerical algorithms exist to find MPS approximations to ground states of much more general situations, and also to simulate real time evolution. 57-61 On the other hand, thermal states of local Hamiltonians can be efficiently approximated by an analogous ansatz in the space of operators, 62-64 referred to as matrix product operators (MPOs). 65,66 We write the state of the full system as a MPS (if the environment is initially in the ground state) or a MPO (in the thermal case), and apply standard MPS methods 59,62 to simulate real time evolution. For convenience, we include the system in the m 0 -th site of the chain (which then has physical dimension d 2 ). The initial state is built as the tensor product of the desired initial spin state and the MPS (resp. MPO) approximation of the spin chain state found with standard MPS algorithms. In the thermal case, we evolve a purification MPS, 61,62,67 which ensures positivity of the evolved state. Bond dimension, system size and trotter step δ were chosen such that for the evolution times reported in the text the results are converged. In particular, we used D ≤ 200, N = 200 and δ = 0.01.

IV. VACUUM INITIAL STATE
Let us consider that the spin is coupled to the center of the chain, m 0 = N/2. The simplest scenario for this setup is when the chain is initialized in the fermionic vacuum ρ E = |0 0|, where d k |0 = 0 for all k. Note that |0 is equivalent to the ground state if h > J (h < −J yields the same physics due to particle-hole symmetry). Then the number of excitations in the total system N exc is set by the initial system state and the only sectors involved are those of zero and one excitations, which are not mixed under the dynamics. We can write the total system-environment state at any time as where C k (0) = 0 and C g (t) = C g (0) does not evolve. The dynamical map Eq. (15) in this particular case has elements a vac = |C e | 2 , b vac = C e , c vac = 0. The expressions in Eqs. (18) to (20) thus simplify to γ vac In this case, with at most one excitation present in the whole system, the fermionic chain we consider is completely equivalent to the bosonic one studied in reference, 37 and the analytical calculation of C e in the thermodynamic limit presented in that work is also applicable to our setup. We used this result to obtain the non-Markovianity degree N , which we plot for a wide range of Hamiltonian parameters in Fig. 1. We find a non-Markovian region for detunings ∆ h around the band edges with a width that increases with the coupling strength. Note that for detunings in this region condition (28) for deriving a 'Markovian' master equation is obviously violated because the spectral density D(E) diverges at the band edges. If we detune far outside the band |∆ h ± 2J| Ω the measure vanishes, which is what we expect because the system effectively decouples from the environment and we have a closed quantum system with coherent dynamics and thus γ i = 0 in Eq. (16). In contrast, strong coupling induces strong non-Markovian behavior because the model effectively reduces to the one of the system coupled to a single spin. 68 We observe that the largest non-Markovianity is not at the band edge, but slightly shifted inside. This is due to the fact that at short times |C e | 2 reaches values close to zero, smaller in this case than exactly at the band edge, leading to larger magnitude γ vac 3 , as shown in Fig. 2  (black and blue lines).
For ∆ h further inside the band, the system is more Markovian, in line with a classification based on condition (28), which is closer to being satisfied due to smaller values of the spectral density and its derivative. 69 In the middle of the band (∆ h = 0), until intermediate times, the time dependent |C e | 2 exhibits a monotonic decay, almost exponential, but modulated by oscillations (clearly appreciated in γ vac 3 ) at a frequency approximately equal to 2J (Fig. 2, orange lines). The dynamics is captured by a 'time-dependent Markovian' master equation. Only when most of the population has decayed (tJ ≈ 25), monotonicity of |C e | 2 is broken and the rate becomes negative. We focus on characterising the behavior at times before that happens and thus those late times do not enter our calculation of N . When the coupling is increased while staying at the center of the band, the oscillations become stronger (see Fig. 3) until they break the monotonicity of |C e | 2 and the description of the dynamics via a 'time dependent Markovian' master equation is no longer possible.
As shown in 37 (and reviewed in Appendix A for completeness), in the calculation of C e one identifies a number of terms which we call resonant, edge and bound state contributions, since intuitively they can be connected to the overlap of the initial state |e, 0 with continuum eigenstates of H close to ∆ and close to 2h ± 2J and with its two bound states 70 respectively. C e is the result of summing up these five contributions. The absolute square of each of these contributions is either a monotonically decaying or a constant function of time. Hence, the nonmonotonicity of |C e | 2 results from cross terms between pairs of those contributions, which oscillate with the difference of their corresponding frequencies ν. The different non-Markovianity behavior obeys to which contributions are important for a given choice of parameters. At short times, and for the detunings discussed above, |C e | 2 is in general dominated by the exponentially decaying resonant contribution. 37 If ∆ h is close to the band edge, also the corresponding edge (power law decaying) and bound state (constant magnitude) contributions become important, 37 and give rise to non-Markovianity via cross terms and to the incomplete decay of |C e | 2 observed in Fig. 2a (black and blue lines).
If ∆ h is close to the middle of the band, the cross terms between the resonant contribution (ν r ≈ 2h) and edge and bound state contributions (ν e± ≈ ν b± ≈ 2h ± 2J) oscillate with a frequency approximately equal to 2J (|ν r − ν e± |, |ν r − ν b± |). Their magnitude depends on the coupling strength, so that only if this is strong enough (see Fig. 3, pink lines) non-Markovianity appears at short times. However, at long times, non-Markovian behavior can appear even for weak coupling. That at long times the dynamics cannot be described by a 'Markovian' master equation, which is a stricter criterion, is well-known. There is always a transition from exponential to power law decay behavior. 5 In here we identify this non-Markovianity at long times with the relevant cross terms decaying exponentially only with half the rate with which the absolute square of the resonant contribution decays. The corresponding oscillations thus dominate at long times (Fig. 2, orange lines). This is in fact a general feature of this model: no choice of parameters results in perfect Markovianity at all times, as for large enough times the constant contributions from bound states at both edges always give a non-monotonic behavior, after the other (non-constant) contributions have decayed. At sufficiently long times, when the cross terms involving the resonant contribution have decayed, cross terms involving edge and bound state contributions from opposite sides that oscillate with a frequency approximately equal to 4J become visible (|ν e+ − ν e− |, |ν e+ − ν b− |, . . . ).
We may ask how much of the non-Markovian behavior described here is detected by the BLP measure. It is in fact easy to see (see Appendix B) that in the vacuum case, since we have a−c = |b| 2 = |C e | 2 , information flows back in the sense of BLP iff d dt |C e | 2 > 0, and thus all (RHP) non-Markovian behavior is detected by the BLP measure.

V. EXCITATIONS IN THE ENVIRONMENT
In order to investigate how the presence of excitations in the environment affects the non-Markovianity analysis of the previous section, we study two scenarios. On the one hand, we consider the environment in a thermal state, for a chain that has the fermionic vacuum as ground state, h = J. This allows us to recover the previous case in the limit of low temperature βJ → ∞. On the other hand, by tuning the parameter h, the ground state of the chain can be chosen to contain the desired number of occupied modes. For these two types of thermal states, the dynamical map is still of the form of Eq. (15), but it is no longer determined by a single parameter. Instead, γ 1 and γ 2 , given by the general expressions in Eqs. (18) and (19) do not vanish, and can now become negative and give rise to non-Markovianity. In this case, no analytical (exact) solution is available, and the dynamical map is computed using tensor network techniques.

A. few excitations induce (new) non-Markovianity
Thermal and ground states with few excitations correspond to populating the lower edge of the band, which can be achieved, respectively, by a low temperature (βJ 1) or suitable chemical potential (0 < 1 − h/J 1). In the rest of this section, we consider the cases h = J at βJ = 10 and h = 0.95J at βJ → ∞ (ground state).
For ∆ h = 0, i.e. in the middle of the band, we observe that the map parameter a varies very little with respect to the vacuum case (see Fig. 4), which is not surprising, as the only excitations in the system are far off-resonant, close to the lower band edge. However, we obtain a contribution to non-Markovianity at short times from γ 2 . This is originated from the monotonicity breaking oscillations of c, clearly observed in Fig. 4a/b, which exhibit approximately the same frequency as the ones of a. As we argue in Appendix A, for the case of a single initial excitation in the environment, we may expect that the time dependence of the c component of the dynamical map is determined by the same frequencies that appear in the vacuum case 71 , and that its oscillations may break monotonicity already at early times since the resonance contribution doesn't have the same dominating effect it has on a vac . Numerically, we find that this signature of the cross terms between resonant, edge and bound state contributions seems to explain qualitatively also the case of few excitations. Next we consider setups, in which the detuning is chosen close to the lower band edge, for which the non-Markovianity was large in the vacuum case. We observe that in this case, the crossing of a and c results in vanishing denominators in Eqs. (18) to (20), such that the rates diverge and change sign (see Fig. 5). The non-Markovianity in this case is thus more dramatic. Looking at Eq. (18), we notice that the early time non-Markovianity (γ 1 < 0) is due to |b| 2 decaying similarly whilst a − c decaying faster than their respective values in the vacuum case.
It is interesting to notice that this non-divisibility at early times, corresponding to the divergence of γ 1 , is not witnessed by the BLP measure. As shown in Appendix B, there is no information backflow in the sense of BLP when a − c ≥ 0, d dt (a − c) ≤ 0 and d dt |b| 2 ≤ 0. These conditions are in fact satisfied in this setup, until the time when a and c cross. They are also satisfied in the setup discussed in the previous paragraph, for ∆ h = 0, so that most of the new non-Markovian phenomena induced by a few excitations in the environment are not detected by BLP.

B. high temperature leads to Markovian dynamics
We may ask how the picture changes with an increasing number of excitations in the environment, either due to a higher temperature, or to a lower chemical potential. First of all, we observe that for thermal states at a higher temperature (see Fig. 6a/c) and for ground states at a lower chemical potential (see Fig. 6b/d) the crossing of a and c described above does not take place. On the other hand, when, for the ground state case, we set the detuning at the Fermi level, we get a crossing (see insets in Fig. 6b/d). Together, these observations suggest that the occurrence of the crossing phenomenon, which is always accompanied by diverging non-Markovianity, is linked to the Fermi-Dirac distribution f (E) changing sharply across ∆ h (compare insets in Fig. 5 and Fig. 6).
In Fig. 6c we can also observe how the increased tem- perature smoothes out non-Markovian effects either introduced by few excitations (Fig. 5) or already present in the vacuum scenario (dotted line). It turns out that it is possible to obtain entirely Markovian dynamics if one chooses a high enough temperature. This is illustrated in Fig. 7 where we plot the transients of the channel elements and rates for system energies close to the band edge (∆ h = −1.8J, left panels) and at the band center (∆ h = 0, right panels) at high temperature. We observe that the time dependence of the channel elements is monotonic and that the rates are positive for the times we can access with our simulations. One might ask if this Markovian behavior may be anticipated in the sense that condition (21) for deriving the 'Markovian' master equation is satisfied at high temperature. Using the explicit expression Eq. (D4), we compute one of the kernels Eq. (27) for βJ = 0.05 (see Fig. 8), and find that it decays rapidly to zero (at such high temperature the other one behaves qualitatively the same). Remarkably, the environment correlation time τ c is essentially the same at the center and the edge of the band, although the corresponding spectral densities are very different (diverges for the latter), and, in the literature, the flatness of the spectral density is often associated to short correlation times. [38][39][40][41] At infinite temperature we can write a closed form for the correlation functions: α ± (t) = 1 2 e ±i2ht e −J 2 t 2 . 72 Their superexponential decay confirms that at high temperatures and small enough (but still intermediate) coupling strength (which sets the system time scale), the 'Markovian' master equation becomes a valid description at all detunings. 73 In summary we found in this section that while few initial excitations, e.g. the ones present at small temperature, introduce a number of new non-Markovian features (early time negative γ 2 for ∆ h deep within the band and early time diverging negative γ 1 for ∆ h close to the lower band edge), these phenomena get smoothed out together with any non-Markovian features already present in the vacuum case as one increases the temperature to large values such that, at high temperature, we obtain divisible dynamics. We found that this high temperature Markovianity could already be anticipated from the ob- servation that the conditions for deriving the 'Markovian' master equation are satisfied at intermediate couplings.
The Markovianity at high temperatures is nevertheless not a completely general effect, as we discuss in the next section.

VI. EXACTLY SOLVABLE SETUP
If the system is coupled to the first site of the chain (m 0 = 1), the full model can be solved analytically, as it can be mapped to a quadratic fermionic Hamiltonian. Using the exact solution, we analyze here the divisibility properties of this setup, and compare them with the case discussed in the previous sections. Notice that the non-Markovianity in this scenario was already studied in reference, 36 with the chain initially in a ground state, but with a focus on the (less strict) BLP measure, which in some cases gives a qualitatively different picture.
The dynamical map is still of the form of Eq. (15) and we compute explicit expressions for its elements using Gaussian methods for free fermions (see Appendix C). We find that in this case a − c is the same for any environment initial state that is a ground state with different h or any other thermal state, i.e., a − c = a vac . Thus, different from the setup considered in the previous section, introducing a few excitations cannot affect the non-Markovianity dramatically (see Eqs. (18) to (20)), in particular, the crossing between a and c cannot happen. Also, a − c = |b| 2 for the environment initial states we consider, which implies γ 1 = 0 Eq. (18).
Under these conditions, information backflow in the sense of BLP is equivalent to having d dt (a − c) > 0 (see Appendix B). Thus, BLP-Markovianity is independent of such environment states and it is enough to check the simplest one, i.e. the vacuum (ground state for h > J). For this particular case BLP-Markovianity and divisibility turn out to be equivalent, and reduce to the condition d dt a vac ≤ 0 ∀t, in accordance with what we found in section IV. Thus, for fixed parameters ∆ h /J and Ω/J, non-divisibility in the vacuum case implies BLP-non- Markovianity (and hence non-divisibility) independent of the considered environment states. This is in contrast to the setup discussed in the previous section, in particular, it is not possible to obtain Markovian dynamics at high temperatures if the corresponding vacuum case is non-divisible.
It is interesting to compare the predictions based on the behavior of the high temperature environment correlation functions for both setups. In Appendix E we show that (at infinite temperature, where they can be obtained in closed form) whilst in the previous case, with its superexponentially decaying correlation functions, the derivation of the 'Markovian' master equation is valid at all detunings if Ω J 1, in the present case (∼ t − 3 2 ) this is only true at detunings away from the band edges | ± 2 − ∆ h J | Ω J . If we move the coupling site into the chain, our numerics show the power law decay getting steeper as ∼ t m 2 0 + 1 2 (Appendix D). The simple expression for the BLP measure in this setup allows us to reproduce one of the main results of: 36 for each value of the coupling Ω/J ≤ 1 there is a specific detuning ∆ h = 2J − Ω 2 J for which the BLP measure is identically zero for any of the environment initial states we consider. We refer to this situation as the BLP-Markovian point. For any other set of parameters the measure is different from zero, with the largest non-Markovianity occurring at the center of the band (∆ h = 0). 36 The divisibility of the channel, on the contrary, depends on the initial state of the environment. In general, c does not vanish, and the evolution may be non-divisible, even for the detuning set at the BLP-Markovian point. This is shown explicitly in Fig. 9: for ∆ h /J = 1, the BLP-Markovian point at coupling strength Ω/J = 1, and an initial state different from the vacuum, the rate γ 2 becomes negative at some times.
On the other hand, for any other parameters and any of the considered initial states, since the BLP measure is always non-zero, the channel is not divisible. However, it is interesting to analyze the time dependence of this non-Markovianity. For the most non-Markovian setup (∆ h = 0), shown in Fig. 10, we find that the dynamics (in the vacuum case) is indeed divisible according to our measure until intermediate times. The BLP-non-Markovianity arises only from the late times (tJ 40), when the population has decayed so much that the oscillations break monotonicity of a vac . This is analogous to our observation in section IV (see Fig. 2, orange lines), in which we obtain large non-Markovianity contributions from late times because a vanishing population leads to diverging γ vac 3 . A significant difference between this setup and the one discussed in the previous sections is that, as shown in reference, 36 in the case of the system coupled at the beginning of the chain, there is a single bound state of the interacting Hamiltonian for Ω/J < 1.5, and only if |∆ h | ≥ 2J − Ω 2 J . Although the analytical calculation of section IV (see also Appendix A) is not directly applicable to this setup, we observe that the oscillations of a vac at early times (during the exponential transient) and at late times still have a frequency approximately equal to 2J and 4J respectively (see Fig. 10), suggesting that they still originate in cross terms involving resonant and edge contributions.
As the detuning is shifted closer to the edge of the band, we observe that a vac decays as a power law modulated by damped oscillations, until, at the BLP-Markovian point, they do not break monotonicity anymore (Fig. 11). This suggests a competition between the cross term involving both edge contributions on the one hand and the monotonic power law decay of the absolute square of the relevant edge contribution on the other hand (for detunings ∆ h ≤ 2J − Ω 2 J ). Generalising from the previous setup, the decreasing relevance of the oscil- − Ω 2 J ) we have a bound state, and cross terms will always break monotonicity of a vac at some times (see e.g. Fig. 11, green line). Finally, for Ω/J > 1, even at ∆ h = 2J − Ω 2 J the cross term involving both edge contributions is still strong enough to break monotonicity of a vac (see Fig. 12). Thus, at a fixed coupling strength, a region (in ∆ h ) without bound states is necessary (but not sufficient) for the existence of BLP Markovian points, and such a region can only exist if the spectral density does not diverge at the band edges (Appendix A).
In the previous section and in Appendix A we argued that signatures of single excitation physics survive in the characteristic frequencies that modulate the channel elements, also in the setups with few excitations. The m 0 = 1 model provides an extreme example of this where these frequencies are present in setups with an arbitrary number of initially populated fermionic modes since a − c does not depend on the initial state of the environment, as far as it commutes with the total number of excitations.
Notice that the BLP-Markovian points are located quite close to the band edge, where the standard derivation of the 'Markovian' master equation is not valid as discussed in Appendix E. 74 Still, in the vacuum case, they are captured by a 'time-dependent Markovian' master equation at all times.

VII. CONCLUSION
In this work we have explored the dynamics of a single spin coupled to a quantum spin chain, when considered as an open quantum system. We have used a simulation of the real time evolution of the whole system to compute the dynamical map that governs the evolution of the spin, and to characterize and measure its non-Markovianity. We have identified situations, determined by the parameters of the system (coupling and detuning) and the initial state of the chain, in which the dynamics of the spin, at least until intermediate times, admits a description in terms of a time-dependent Markovian master equation, i.e. the map is divisible. Some of these scenarios occur in regimes that do not allow a standard derivation of the master equation.
We studied two scenarios. In the first one, we couple the spin to the center, in the second, to one edge of the environment chain.
In the first case, and when the chain is initialized in the vacuum, we find a Markovian parameter region when the detuning of the spin is deep within the band of the singleparticle spectrum of the environment and the coupling is small to intermediate. This is in line with a characterisation based on the validity of the standard derivation of the master equation. 37 Setting the detuning close to the band edges produces strong non-Markovian effects, as does a strong coupling.
If the initial state of the chain contains a few excitations close to the lower band edge, the scenario changes: the Markovian regions disappear and the non-Markovianity close to the edge increases dramatically. The latter effect persists also beyond few excitations if one initializes the environment in a filled Fermi sea and sets the detuning close to the Fermi level. On the other hand, a high temperature initial state of the chain, which introduces a large number of excitations evenly spread across the spectrum, results in Markovian behavior, even for detunings close to the band edges (where the spectral density diverges).
In the exactly solvable case of the chain initialized in the vacuum, the Markovian or non-Markovian character of the map can be completely explained in terms of the eigenstates of the full Hamiltonian and how they contribute to the scattering amplitude. In particular, non-Markovianity obeys to the presence of sufficiently strong cross terms between different contributions. A qualita-tively similar picture holds in the case of few initial excitations.
In the second case, when the spin is coupled to the edge of the chain, the problem is exactly solvable. A remarkable difference to the first case is that, if the chain is initialized in the vacuum, there are points in parameter space that are Markovian at all times. We explained this phenomenon with the cross term argument above and found that a non-diverging spectral density at the band edges is paramount to the existence of such points. Another difference to the first case is that while any Markovianity still disappears on introducing few excitations, a dramatic increase of non-Markovianity does not occur. Finally, high temperature does not impose Markovian dynamics. Instead, any non-Markovianity of the vacuum case survives at all temperatures. At high temperature this is in stark contrast to the Markovianity of the first case, but is consistent with the behavior of the environment correlation functions that we have computed, which, whilst showing a superexponential behavior inducing Markovianity in the first case, are characterized by a power law behavior for the second case, ruling out the standard derivation for detunings close to the band edges. The decay of the correlation functions becomes steeper as the position of the coupling is moved away from the edge of the chain.
We define Markovianity as divisibility of the evolution, but we can also compute other non-equivalent non-Markovianity measures. In particular, we have compared the results to the widely used BLP measure, less restrictive, which does neither detect the non-Markovianity deep within the band, nor the early time onset of the dramatic non-Markovianity close to the lower band edge, when few excitations are present in the environment.

VIII. ACKNOWLEDGEMENTS
We gratefully acknowledge the insights of and discussions with T. Shi, A. González-Tudela, I. de Vega and R. Verresen. JR is particularly grateful to T. Shi for pointing the way to the method in Appendix D, and to A. González-Tudela for in-depth discussions on the vacuum case. JR acknowledges ExQM and IMPRS for support. This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC-2111 -390814868, and by the European Union through the ERC grant QUENOCOBA, ERC-2016-ADG (Grant no. 742102).

Appendix A
In reference 37 the amplitude C e (t) was explicitly computed for the case of an emitter coupled to a bosonic tight-binding chain, which, in the thermodynamic limit, is equivalent to our model with m 0 = N/2 in the vacuum case. This is done by expressing the amplitudes as where α ∈ {e, k}, and using the structure of singularities in the complex energy plane of the retarded Green functions with the self energy The ±-sign depends on whether Re(z − 2h) ≷ 0.
In particular it was found that C e (t) can be decomposed into a sum of contributions (roughly) due to different parts of the interacting spectrum where we chose labels representing the parts close to the (upper,lower) band edge (UE,LE), close to resonance ∆ (RS) and the upper and lower bound states 70 (UBS,LBS) of the interacting spectrum.
If we consider that the environment in the initial state contains a single excitation, the channel element c from Eq. (15) can be computed as 75 such that, following reference, 37 its behavior will be determined by the singularities of Eq. (A3).
In the (large) finite case the resonance and edge contributions in Eq. (A5) correspond to terms that are dominated by a sum over terms | ñ|e, 0 | 2 e −iẼnt running over (scattering 70 ) states |ñ close to resonance ∆ and close to the band edges respectively (H|ñ = E n |ñ ).
The bound state contributions are simply | BS ± |e, 0 | 2 e −iẼ BS ± t . With this intuitive picture one can anticipate that these contributions oscillate with frequencies ν approximately given by ∆ (ν r ≈ ∆), the band edges (ν e± ≈ 2h ± 2J) and the bound state energies (ν b± =Ẽ BS ± ) respectively, which we confirm by explicit computation. The magnitudes of the contributions decay exponentially (RS), with a power law (UE,LE) or are constant (UBS,LBS). 37 Strictly, the resonance and bound state contributions stem from singularities in G e (z). In 70 it was proven that if Im Σ e (z) diverges at the band edge, where it is proportional to the spectral density D(E), 37 there always exists a pole associated to a bound state. This is the case when m 0 = N 2 , where the spectral density is proportional to the (diverging) density of states; otherwise it is not necessarily the case. 76 These singularities are still present 77 in G k (z): which has an additional pole corresponding to the free propagator of the mode k, leading to a constant magnitude free propagator contribution (oscillating with frequency ν f = E k ) to C k (t). If the initial environment excitation is close to the lower band edge (E k ≈ 2h − 2J), we thus expect to have terms in C k (t) that oscillate with similar frequencies as those in C e (t), but since |C k (0)| 2 = 0, in this case the resonant contribution cannot be dominating all the other contributions. 78 This is the mechanism behind the early time non-Markovianity (at ∆ h = 0) in V A. Fig. 13 illustrates this for a finite chain as studied in our paper. We set the initial state to d † N |g, 0 , i.e. the environment contains a single excitation in the lowest energy mode, and compute c with exact diagonalisation, which, in this case, can be done efficiently since the dynamics is restricted to the one excitation sector. At the band center (∆ h = 0), we find that the frequency of the oscillations is approximately equal to 2J (|ν r − ν e± |, |ν r − ν b± |, |ν r − ν f |) and that this feature is stable upon increasing the chain length (upper panel). Also, in the inset, we observe that there is a transition after which the frequency is approximately equal to 4J (|ν e+ − ν e− |, |ν e+ − ν b− |, |ν e+ − ν f |, . . . ). These two frequencies are characteristic for the interplay between these contributions. The frequency survives in the few excitations case (bottom panel), where an increasing number of excitations is obtained by initialising longer chains in their ground state at fixed chemical potential. Fig. 14 illustrates the importance of the free propagator contribution, since for the case that the environment in the initial state contains a single excitation somewhere within the band, the frequency of the oscillations of the component c is approximately given by |ν r − ν f |. Looking more carefully at the oscillations in Fig. 4b, we find that their frequency lies somewhere between the highest occupied mode (at the Fermi level E k = 2h − 1.9J) and the lowest occupied mode (at the lower band edge In the m 0 = 1 model of section VI we find that a − c = a vac is independent of the environment initial state and can thus show no signature of the free propagator contributions. This means that beyond the vacuum case, free propagator contributions must also contribute to a such that their effect cancels out in a − c. We remark that the presence of oscillations due to free propagator contributions is independent of the existence of bound states and band edges, and should induce non-Markovianity also in unstructured environments (constant spectral density) if few initial excitations are present.
To prove the other direction we choose the initial pair ρ(0) = |e e|,ρ(0) = |g g|. Taking into account the general structure of the channel Eq. (15), we have v 1 = v 2 =ṽ 1 =ṽ 2 = 0 at all times and thus: Hence, since a − c ≥ 0, we find thatȧ −ċ > 0 at some time implies that the trace distance is increasing (BLP measure nonzero) and the statement is proved.

Appendix C
Here we give the concrete expressions of the channel elements in the case of the spin coupled to the edge of the chain (section VI).
The total system can be mapped to the free fermion Hamiltonian: with real space fermionic operatorsc i = e iπ i−1 j=0 σ + j σ − j σ − i , where σ α 0 ≡ τ α 0 are the spin operators on the subsystem. The channel elements of Eq. (15) can be written in terms of the Heisenberg picture fermionic operators as follows: where I ∈ {e, g, x+, y+} denotes the global initial state |I I| ⊗ ρ E and we have defined |x+ = 1 √ 2 (|e + |g ), |y+ = 1 √ 2 (|e + i|g ). Note that in contrast to fundamental fermion models we can have c 0 = 0 since the 0-th mode corresponds to the spin subsystem. The parity (of N exc ) is nevertheless conserved, such that if one starts with a superposition in the spin, one can solve each sector independently and then add them up, where in each of the calculations c 0 = 0. These expectation values can be computed exactly.
We write down the Heisenberg equation of motion where H is the real, symmetric, N + 1 dimensional tridiagonal matrix in H = N ij=0c † i H ijcj . Taking the expectation value with respect to global initial state I on both sides and defining, for each initial state, a matrix M I ij (t) ≡ c i (t)c † j (t) I , we get the matrix equations: where in our convention X 00 denotes the first matrix element of a N +1 dimensional matrix X with matrix indices running from 0 to N . Similarly, and thus |b| 2 = a − c. This immediately implies that γ 1 = 0 (see Eq. (18)).
spin (e.g. decay rate), which in our model is set by the coupling strength Ω. Thus for small enough coupling, At the band edges, where the oscillations are too slow to kill the integral, the slow power law decay violates the Markovian approximation at any Γ > 0.
We thus find that whilst, deep within the band and at small enough coupling, the dynamics is captured by the 'Markovian' master equation in both cases, the underlying mechanism is completely different: superexponentially decaying correlation functions in (i); rapid oscillations of the correlation functions in (ii). The self consistency condition (28) is blind to what is the fundamental origin of its validity. Too close to the band edges (| ± 2 − ∆ h J | Γ J ) the Markovian master equation is only valid for (i) (at small enough coupling).
its PB form W PB m 0 ,k = 1