Nonequilibrium, spatio-temporal formation of the Kondo screening-cloud on a lattice

We study the nonequilibrium formation of a spin screening cloud that accompanies the quenching of a local magnetic moment immersed in a Fermi sea at zero temperature. Based on high precision density matrix renormalization group results for the interacting single impurity Anderson model we discuss the real time evolution after a quantum quench in the impurity-reservoir hybridization using time evolving block decimation. We report emergent length and time scales in the spatiotemporal structure of non-local correlation functions in the spin- and the charge density channel. At equilibrium, our data for the correlation functions and the extracted length-scales show good agreement with existing results, as do local time dependent observables at the impurity. In the time-dependent data, we identify a major signal which defines a"light cone"moving at the Fermi velocity and a ferromagnetic component in its wake. Inside the light cone we find that the structure of the nonequilibrium correlation functions emerges on two time scales. Initially, the qualitative structure of the correlation functions develops rapidly at the lattice Fermi velocity. Subsequently the spin correlations converge to the equilibrium results on a much larger time scale. This process sets a dynamic energy scale, which we identify to be proportional to the Kondo temperature. Outside the light cone we observe two different power law decays of the correlation functions in space, with time and interaction strength independent exponents.


I. INTRODUCTION
Quantum impurities are among the most fundamental paradigms of strongly correlated quantum systems. Equilibrium properties of such systems have been subject to intense investigation and are nowadays well understood. A famous example is the Kondo effect, where a local spin- 1 2 degree of freedom interacts with the spins of a sea of free electrons [1]. The ground state of this system is a delocalized spin singlet, formed by the local moment and the spin of the free electrons, also called a screening cloud. The present work investigates how such a screening cloud develops over time when a local moment comes into contact with a free electron reservoir.
Quantum impurity systems, quite generally, feature an emergent screening length scale at low temperatures which provides the basis for their complex physics. In the 1950s, magnetic impurities have already been identified as the cause for a large resistivity anomaly at low temperatures when immersed in metallic hosts in dilute quantities [2,3]. It was found theoretically that the impurity's local magnetic moment becomes quenched below a certain temperature, known as the Kondo temperature [1,4], T K , to form a local Fermi liquid [5]. Increased spin-flip scattering between pairs of degenerate spin- 1 2 states then leads to an increase in resistivity below T K . Meanwhile, the Kondo effect has been observed also in nanoscopic devices like quantum dots [6][7][8][9][10][11][12], carbon nanotubes [13], and molecular junctions [14]. Here, the narrow, zero-energy resonance in the local density of states of the impurity, the Kondo-Abrikosov-Suhl resonance, leads * martin.nuss@tugraz.at to a well-defined unitary conductance in linear response. The Kondo effect has also proven essential to understanding tunneling into single magnetic atoms [15], adsorption of molecules onto surfaces [16], or defects in materials such as graphene [17]. On the theoretical side the Kondo effect lies at the heart of our current understanding of correlated materials, notably within the very successful dynamical mean-field theory [18][19][20].
Insight into the details of the screening cloud is important not only for the understanding of the physics of a single impurity but also for the understanding of the interplay of many magnetic impurities. Many impurities result in competing effects among conduction electrons and local moments, which form the basis for spin exhaustion scenarios [21,22] as well as for the Doniach phase diagram [23,24], which describes the relationship between Kondo [1] and RKKY interaction [25][26][27].
Experimental characterization of the structure of the singlet ground state, which is a bound state of the impurity spin and the reservoir electron "screening cloud," has proven difficult so far. Several proposals exist for how to measure the spatial extent of the spin screening cloud or its antiferromagnetic correlation with the impurity spin [28,29]. In principle, the real-space structure could be probed by performing nuclear magnetic resonance/Knight shift [30][31][32] measurements on bulk metals hosting dilute magnetic impurities, but the approach remains challenging [28]. Indirect observation by measurement of the Kondo resonance, for example, by photo emission, also remains elusive due to the too narrow resonance at the Fermi energy [33]. Other proposals suggest the use of scanning tunneling microscopy [34] and scanning tunneling spectroscopy to analyze adatoms or surface defects with Kondo behavior [35,36]. In the realm of nanodevices, proposals include experiments based on persistent currents [37] or in confined geometries [38,39]. Some progress has been made recently using single magnetic atoms [15], quantum corrals [40], or impurities beneath surfaces [35].  o/e that is proportional to the Kondo temperature T K . Region 2 (red) lies inside the light cone but outside the Kondo screening cloud. Here the spin correlations decay as a power law in space [53]. In region 3 (blue), which lies outside the light cone and outside the Kondo screening cloud, the correlation function at odd/even distances decays as a power law ∝ r −γ S/C o/e in space with exponents that are independent of time and interaction strength. in the Kondo model, only spin interactions survive and charge fluctuations are treated on an effective level [64], we take them into account explicitly. To our knowledge, our study is the first one analyzing the nonequilibrium properties of the screening length in the interacting SIAM.
Our results are summarized in Fig. 2, which also serves as a guiding map for this work. We identify a major signal following the quench, which propagates with the lattice Fermi velocity v F and defines a light cone for the propagation of information [55-57, [65][66][67][68][69]. Inside the light cone the timeevolved correlation functions converge to their equilibrium counterparts which exhibit the Kondo length scale. We find that Kondo correlations develop on two characteristic time scales. The main structure of the Kondo singlet is formed rapidly at v F . These correlations approach their equilibrium values exponentially in time, defining a dynamic energy scale α S o/e , which is proportional to the Kondo temperature T K . Outside the light cone, we find that correlation functions at odd/even distances decay as a power law ∝ r −γ S/C o/e in space, with exponents which are independent of time and interaction strength.
The structure of this paper is as follows: We summarize the specific model used in Sec. II. We define the Kondo singlet in Sec. III, present our numerical approach in Sec. IV, and provide an overview of the equilibrium situation in Sec. V. We start our presentation of nonequilibrium phenomena in Sec. VI, where we discuss the evolution of local observables. The main findings of this work are available in Sec. VII. There we discuss the nonequilibrium formation of the Kondo screening cloud in Sec. VII A. The situation outside the light cone is presented in Sec. VII B. The quality of our numerical data is assessed in Appendix A.

II. MODEL
We study a lattice realization of the SIAM [58], which consists of a single fermionic spin-1 2 impurity coupled via a standard hopping term to a reservoir of noninteracting tight-binding fermions (see Fig. 1). In particular, we consider a particle-hole symmetric impurity with on-site interaction U : The electronic annihilation (creation) operators f σ (f † σ ) obey the usual anticommutation relations with spin σ = {↑,↓}, and n f σ = f † σ f σ is the particle number operator [70]. The impurity is coupled via a tunneling term, to a one-dimensional tight-binding chain, such that the overall system, including the impurity, is of length L. We always take the reservoir FSĤ res half-filled. For large L, the reservoir mimics a semi-infinite one-dimensional tight-binding reservoir [71] with a semicircular density of states at the first site and bandwidth D = 4 t [72]. Studies of finite-size effects are available in Refs. [33,37,45], and [73][74][75][76]. The hopping parameter of the reservoir t is taken to be unity, and its coupling to the impurity t = 0.3162t combine to an equilibrium Anderson width where ρ reservoir (ω) denotes the density of states of the reservoir. At equilibrium, many characteristics of the SIAM are known, although it poses a difficult interacting problem. Seminal results for the ground-state and thermodynamic properties of the SIAM at equilibrium are available from perturbation theory [77][78][79][80], renormalization group [81][82][83][84], and the Bethe ansatz (BA) [85][86][87]. The Hirsch-Fye QMC [48,88] and the continuous-time QMC [89] accurately describe the imaginary time dynamics. Further, some physical results can be inferred from the Kondo Hamiltonian, which is related to the SIAM by the Schrieffer-Wolff transformation, to obtain its low-energy realization, in which charge fluctuations are integrated out [64,81].

III. THE KONDO SINGLET
at equilibrium, the SIAM features a characteristic length scale which, for finite interaction strength, is the Kondo length scale and is expected to correspond to the size of the singlet screening cloud. This length scale is defined as ξ K ≡ v F T K [28,41,49,[90][91][92]; i.e., it is proportional to the Fermi velocity v F ≈ 2t and the inverse Kondo temperature 1 T K [1,86,87]. T K can be extracted from many observables; most intuitive is the definition as the temperature at which the local moment becomes quenched, i.e., when the impurity entropy goes from ln(2), indicating the local moment regime, to ln(1), indicating the singlet state [93]. A scale proportional to T K is also available from the zero-temperature self-energy [94] or from the width of the Kondo resonance in the spectral function [95]. An analytic expression for T K , as obtained from the spin susceptibility, is available for the SIAM at particle-hole symmetry in the wide band limit with a linear dispersion [96] by the BA [85][86][87]: U . The Kondo singlet, therefore, is exponentially large in the interaction strength U : For typical Kondo materials, like dilute magnetic impurities in free electron metals [97], one finds v F ≈ 10 6 m/s and T K ≈ 1 K valid, for example, in gold with dilute iron impurities [98]. Thus, the screening length becomes macroscopic, ξ K ≈ 1 μm [49].
Here, we extract the screening length scale ξ K directly from correlation functions, and not via the Kondo temperature. The spin correlation function is defined as whereŜ r = (Ŝ x r ,Ŝ y r ,Ŝ z r ) [99] and r denotes the distance from the impurity in units of the lattice spacing (see Fig. 1). Due to the oscillations of S, it is convenient to distinguish between the spin correlation function for odd [S o (r,τ )] and that for even [S e (r,τ )] distances.
Length scales can be extracted from the crossover in the functional dependence of S o (r,τ ) or via determining zeros or minima in S e (r,τ ) [49,50,53]. Its charge analog is defined as [100] C(r,τ ) = σ σ n 0σnrσ (τ ).
Correlation functions without a time argument, S(r) and C(r), refer to the ground state of the equilibrium system, Eq. (1), i.e., an impurity coupled to the free electron reservoir. Steady-state correlation functions are indicated with τ → ∞. Later we show that in this limit the time-dependent correlation functions converge to the equilibrium correlations, S(r,τ → ∞) = S(r), as expected from the fact that the quench is intensive. An intuitive measure which quantifies how much of the singlet correlations is contained inside a distance r is the integrated spin correlation function, As discussed below and in Refs. [48] and [52], the screening length ξ k can be extracted from (r,τ ) by defining it as the length scale at which a certain fraction f (here we use f = 95%) of the correlation lies inside a given distance; i.e.,

IV. METHOD
Here we outline how the correlation functions, Eq. (4) and Eq. (5), are evaluated. We start with a short discussion of the noninteracting system in equilibrium. In this case we find where m r = 1 2 n r↑ −n r↓ , and the last result holds for the unpolarized case. Here, c † r /c r denote operators for any one of the spin directions σ = {↑ , ↓}. The opposite spin direction is denotedσ = −σ . For U = 0 at equilibrium [48] In the particle-hole symmetric and non-spin-polarized case the asymptotic limits can be analytically evaluated, using results of Ghosh et al. in Ref. [53], to be for odd r, with γ ≈ 0.577 216 the Euler-Mascheroni constant. The correlation function becomes 0 for even distances r. The behavior of the spin correlation function exhibits a crossover at distance ξ U =0 ≈ v F , which defines a length scale in the noninteracting system. We obtain both S(r,τ ) and C(r,τ ) for zero temperature from computer simulations using matrix product state [63] techniques. First, to study ground-state correlations, we employ the DMRG [59-61] on a system of length L, which is typically 500 sites. Second, to study the dynamic formation of the Kondo singlet, we start from a decoupled system in the state | = |↑ impurity ⊗ |FS reservoir , with a non-spin-polarized half-filled FS, at time τ = 0, and then switch on the tunneling term t = 0.3162t for times τ > 0. The evolution in real time is obtained from TEBD [62].
Matrix-product-state-based time evolution has proven to be a highly accurate method to evaluate the properties of one-dimensional strongly interacting quantum systems out of equilibrium [102][103][104][105][106][107][108][109][110][111][112][113][114]. The combination DMRG and TEBD is quasiexact as long as the quantum entanglement stays tractable. It has been shown that the main limitation arises due to the growth of entanglement after the quench [103,115], which ultimately restricts the available simulation time. Furthermore, since we are interested in the physics resulting from an infinite bath, the maximum available simulation time is restricted by reflections at the lattice boundary and therefore by the finite spatial extent of the system. We have been able to reliably evolve the system long enough to reach a local steady state for all presented data sets. We have checked the convergence of our correlation functions carefully by (i) making comparisons to exact data in the U = 0 system, (ii) systematically studying the dependence on the TEBD matrix dimension χ for finite U , and (iii) carefully analyzing the entanglement entropy. These analyses as well as details of the numerical approach and parameters are provided in Appendix A.

V. EQUILIBRIUM
We start our discussion by presenting the equilibrium spin [S(r)] and charge [C(r)] correlation functions. S(r) was first 085127-4 studied by Iishi [41], and C(r) by Grüner et al. [42,43], who determined the basic spatial dependence and properties. Seminal QMC data [48] have been extended with the use of the NRG [49,50] and recently also the DMRG [52]. Here we summarize the most important findings, relevant for the subsequent discussion, and we provide details specific to the finite-size model and numerical method used. In particular, we identify a length scale in the equilibrium spin correlation function and show later that our nonequilibrium correlation functions converge to the equilibrium correlation function for long times τ .
As shown in Fig. 3, both S(r) and C(r) exhibit an oscillating behavior, ∝ sin (k F r). Since the system is half-filled, the Fermi wave vector is k F = π 2 and the oscillation period is r = 2 sites. We first discuss the spin correlations for U = 0 using Eq. (9). In this case we find S U =0 (0) = 3 2 n 0↑ (1 − n 0↑ ) = 3 8 . Furthermore, S U =0 (r) vanishes for even distances r, which follows from general properties of tight-binding fermions [116]. The odd-r correlations S o (r) are negative and therefore antiferromagnetic with respect to the impurity. For U > 0, S o (r) stays negative and increases in magnitude [117]. At the same time, the spin correlation function for even distances S e (r) gradually develops ferromagnetic correlations at short distances, while it is antiferromagnetic at longer distances. On the one hand, it is the antiferromagnetic component which reflects the screening cloud and signals the formation of the singlet ground state. On the other hand, the ferromagnetic component can be attributed to Coulomb repulsion of opposite spins [48]. Neither the period nor the phase of the oscillations is changed by the presence of interactions [48].
The charge correlation for U = 0 is linked to the spin correlation via Eq. (10). There is oscillatory behavior between even and odd sites. For even sites the correlation function is unity, while for odd sites it increases monotonically towards unity. For finite interaction strengths we observe a suppression of these Friedel-like oscillations [118] with increasing U [119]. At even distances the charge correlations show behavior similar to that of the odd channel, however, of a smaller magnitude. The suppression due to the interaction can be traced back to the change in the impurity spectral weight, which develops a narrow Kondo resonance with a width proportional to T K at the Fermi energy [42,43].
While at U = 0 the characteristic length scale is ξ U =0 ∝ v F , for finite U , long-range correlations develop, which change the behavior at a distance ξ K ∝ v F T K . This crossover, characterizing the size of the Kondo spin compensation cloud, is visible in the spin correlation function S(r). Figure 4 (top) shows that the antiferromagnetic spin compensation is visible in the spin correlation function at odd distances, S o (r). S o (r) changes from a logarithmic dependence at small r v F to a power-law behavior at large r v F [see Eq. (11)] [53,120]. We note that this is different from the Kondo model, where the behavior is S(r) ∝ r −d for r < ξ K to S(r) ∝ r −(d+1) for r > ξ K , with d being the dimensionality of the conduction electron reservoir [49,50].
The crossover is difficult to extract directly from numerical data for S o (r) since very large system sizes and small are required to reach the low-r v F limit. We nevertheless found two ways to obtain an estimate for the crossover scale. First, , which displays a crossover between two behaviours at small and large r. This is particularly obvious in the U = 0 results, shown in the inset. There, S U =0 (r) displays the asymptotic behavior given in Eq. (11). The large-r behavior is shown by the dashed black line. Our interacting matrix product state calculations are determined for = 0.1, which corresponds to the thickest (green) line in this plot for U = 0. Middle: Integrated correlation function (r) of Eq. (6). Dashed vertical lines indicate the distances ξ K inside which 95% of the singlet cloud is contained, which we use to estimate the screening length. Bottom: Spin correlation function for even distances, S e (r). The position ξ e K of the minima (circles and dashed vertical lines) is used as the alternative definition for ξ K . Inset: These ξ e K [(green) squares] and ξ K [(orange) triangles]. As reference data we show the BA result and data obtained from an NRG calculation [101]; see text.
a screening length scale can be extracted from the integrated correlation function (r) [see Fig. 4 (middle)]. Similarly to Refs. [48] and [52], here we denote ξ K the distance at which 95% of the singlet correlations are covered, i.e., by Eq. (7). Second, we extract a length scale ξ e K from the spin correlation function at even distances S e (r), which, for finite U , contains both a ferromagnetic component at short distances and the decaying antiferromagnetic one at large distances. As shown in Fig. 4 (bottom) the function S e (r) displays a 0 and a minimum and is fit well by a Morse potential [121]. We take the position of the minimum as a measure for the crossover scale ξ e K . The numerical results obtained with these two crossover scales agree very well and they also agree qualitatively with that obtained by locating the crossover length between r −1.4 and r −(1+1.4) behavior in the S o (r) data, which can be estimated from Fig. 4 (top).
In the inset in Fig. 4 (bottom) we show that our two estimates, ξ K and ξ e K , agree well with established results for the equilibrium screening length. An analytical result, ξ BA K [Eq. (3)], for the screening length is available via its relation to the Kondo temperature, which can be obtained from the BA in the wide-band limit by calculating the static spin susceptibility, Eq. (3). A second benchmark is provided by accurate numerical data from the NRG [101,122], where T NRG,S K is defined as the temperature at which the impurity entropy reaches S = ln (2) 2 [123]. However, while the large-U behavior is universal for all these definitions, the small-U expression, as well as the overall coefficient, depends on the specific observable from which it is extracted (spin susceptibility, entropy, etc.). Our data, ξ K and ξ e K , agree well with the NRG result, ξ NRG,S K ; they are all compatible with a simple exponential growth in U , For U > 2 this agrees with the BA prediction, Eq. (3), which features an additional factor of √ U . We note that for U 2 no well-defined local moment has formed [52]; i.e., U is too small for the system to develop a pronounced local moment regime in between the low-and the high-temperature limit. Our data also compare very well with those presented in an extensive study of length scales in the SIAM on a finite lattice in equilibrium in Ref. [52].
These results indicate that the method presented here is reliable in producing unbiased correlation functions in equilibrium which exhibit the characteristic features of a Kondo screening cloud. The cloud is well contained within the numerically tractable lattice size L 500 for U 6 (see Appendix A). Therefore, we focus our calculations on U 6 when discussing the time-dependent correlation functions.

VI. TIME EVOLUTION OF LOCAL OBSERVABLES
Before beginning the discussion of the temporal evolution of spatial correlations we present the time evolution of the local observables with a focus on the impurity site. At time τ = 0 we start with a spin-↑ particle at the impurity and a nonspin-polarized half-filled FS: | = |↑ impurity ⊗ |FS reservoir . For the connected equilibrium system in the thermodynamic limit we expect a uniform and non-spin-polarized density that is n 0↑ (τ → ∞) = 0.5, n 0↓ (∞) = 0.5, n 0 (∞) = 1, and S z 0 (∞) = 0. The impurity double occupation in a noninteracting system or in the high-temperature limit is n 0↑ n 0↓ U =0 (∞) = 0.25 [124]. For finite interaction strength the evolution is nontrivial. Figure 5(a) shows the expectation values of the spin-z projection S z r (τ ) = 1 2 ( n r↑ (τ ) − n r↓ (τ )). Due to particlehole symmetry, the total charge density n r (τ ) is unity.
Indeed we find that, following the hybridization quench, the excess spin-↑ on the impurity is transported away. This happens essentially at the Fermi velocity v F ≈ 2t, as shown by the major signal in Fig. 5(a).
The resulting missing spin-↑ density is exactly compensated by the spin-↓ density due to particle-hole symmetry. This compensation takes place simultaneously and completely symmetrically in both spin channels as is evident from the spin-↑ and spin-↓ currents shown in the inset in Fig. 5(a). The time integral over the spin current reveals that half a particle is transferred in or out of the impurity in a time of the order of ≈3 −1 for U = 3 . Figure 5(b) shows the local evolution of expectation values as a function of the time and interaction strength. All expectation values converge to their respective, exactly known, equilibrium values as noted above. The time-evolved double occupancy also converges to the equilibrium results obtained by DMRG. This indicates that our time evolution is accurate and unbiased, at least for large times. For more convergence checks and uncertainty estimates we refer the interested reader to Appendix A.
At a certain distance r from the impurity, a resulting signal arrives at τ ≈ r v F . This signal is oscillating and strongly damped in time [see Fig. 5(b)]. With increasing interaction strength U , the initial spike becomes dampened in amplitude, but the oscillating tail gains weight. The signal at r = 40 in the double occupancy has the same structure on a scale of 10 −3 around its equilibrium value.
In the following we consider the temporal decay of the spin-z density at the impurity in detail. Previous studies using the time-dependent NRG for the SIAM [125] and analytical calculations at the Toulouse point of the anisotropic Kondo model [126] found that the initial dynamics of S z o (τ ) is governed by a fast time scale, ∝ 1 , while the eventual relaxation exhibits a long time scale, ∝ 1 T K , governed by Kondo physics. These results were confirmed by bold-line QMC simulations [127] on the SIAM, which showed these two time scales collapsing into one for an applied bias voltage.
From our data we find that, as expected, the spin-z density at U = 0 decays in a single-exponential manner, hence it features the fast hopping time scale T U =0 ∝ 1 . For finite U , a double-exponential decay develops: similar to T U =0 . The corresponding coefficient c 2 decreases in magnitude with increasing U . In contrast, the more interesting slow exponential decay G 1 has a coefficient c 1 which becomes more and more prominent with increasing U . In particular, the coefficient c 1 exhibits a linear behavior in U : The slow decay rate G 1 is exponentially small in U : This implies that the Kondo physics manifests itself in the local dynamic observable S z 0 (τ ) in the form of a slow time constant, T slow ∝ e 2(0.19±0.02) U , which shows the same U behavior as the Kondo temperature [cf. Eq. (12)].
The double occupancy n 0↑ n 0↓ (τ ) converges to its equilibrium value with the same dominant slow decay (within numerical uncertainty) as observed in the spin-z density for finite U . At U = 0 the fast decay rate is twice the rate observed in the spin-z density at U = 0.
Performing the same analysis for distances r away from the impurity that considers S z r (τ ), we again observe the same decay as at the impurity site within the fit uncertainty [see Fig. 5(b)]. This supports the quasiparticle picture introduced in Ref. [56], which translates the physics at the impurity via emission of spin-dependent quasiparticles to a given distance r.

VII. TIME EVOLUTION OF THE SCREENING CLOUD
A very interesting question is how the spatial structure of the Kondo screening cloud develops, i.e., whether and how it is approached in a nonequilibrium time evolution starting from an initial state without Kondo physics. The question was recently first studied in pioneering work on the case of an exactly solvable model, namely, the anisotropic Kondo model at the Toulouse point [55, 56,128]. A complementary numeric study using the time-dependent NRG [57] was performed shortly afterward on the isotropic Kondo model, extending and confirming the analytical results from the Toulouse limit.
After the quench in the hybridization t , we observe a strong signal in S(r,τ ), traveling at the Fermi velocity v F ≈ 2t, which defines a light cone. It has been attributed to the propagation of quasiparticles in Ref. [56]. The propagating signal front divides the space time into two regions: (i) a region at large times and small distances, where the correlation function is directly affected by the impurity and Kondo correlations develop; and (ii) a region at small times and large distances, where small structures from the quench are observed. In Sec. VII A and Sec. VII B we discuss the detailed behavior of the correlation functions within these two regions. The signal front itself carries a large chaotic disturbance in its wake and a small monotonic perturbation ahead of it.
As discussed below in detail, the time-evolved data S(r,τ ) converge to the equilibrium correlation functions S(r) within the light cone. Already a look at the almost-vertical structures in Figs. 6(c) and 6(d) for times τ 8 −1 and a comparison of the line plots for τ = 6.5 −1 and τ = 9 −1 for small distances r hint at the convergence to a local steady state within the light cone. Figure 6(d) reveals that, as expected from the equilibrium state, a ferromagnetic correlation develops for even distances r in S e (r,τ ) within the light cone for finite U , while outside the light cone this correlation function is always antiferromagnetic. As shown in Figs. 6(a) and 6(c), the wake behind the light cone carries a ferromagnetic signal also at odd distances r that is in the otherwise antiferromagnetic S o (r,τ ) for all U . We interpret this signal as remnant information of the spin which occupied the impurity at τ = 0 before the quench. Following the signal wake, all characteristic features of the equilibrium correlation function S(r) develop quickly on a qualitative level. Far behind the signal wake the antiferromagnetic component coincides with S o (r,τ ), and the ferromagnetic component with S e (r,τ ).
A closer look, as provided in Fig. 7, reveals that the nonequilibrium correlation functions gradually develop the characteristic features of the equilibrium correlation functions S(r) and C(r) for r < v F τ . As a precursor of the equilibrium structure, the spin correlation function S(r,τ ) develops the oscillatory behavior of its equilibrium counterpart inside the light cone. That is, it oscillates from an antiferromagnetic correlation at odd distances r to a ferromagnetic correlation at even r for finite U or to 0 at U = 0. This structure of the phase and period of these oscillations in space is fixed over time inside the light cone. However, the light cone induces a phase shift of π in the nodal structure of the correlation function. We 085127-9 , and even, S e (r,τ ) (right), distances are depicted as a function of the distance r for three times-τ = 2.5 −1 , τ = 5 −1 , and τ = 7.5 −1 -in a log-log fashion (dashed lines). We plot −S o (r,τ ) since it is almost entirely negative, while S e (r,τ ) is positive inside the light cone and negative outside (see Fig. 6). Blue lines represent |S e (r,τ )| in regions where S e (r,τ ) is negative. The key depicted in the left panel is valid for both panels. Green arrows mark the direction of increasing time τ . Data from the equilibrium simulation are plotted in solid black and referred to as τ = ∞ in the key. The vertical cyan lines in the left panel mark those distances at which cuts through the data as a function of τ are presented in Fig. 9. All data shown are for U = 3 . attribute this phase shift to the initial state of the FS. It takes place across the broad signal behind the light cone visible at around r ≈ 30 in Fig. 7. As a function of U the same behavior is present inside the light cone as in equilibrium, apart from the chaotic disturbance at the light cone. The qualitative functional form of the correlation functions develops quickly in the wake of the light cone. However, its amplitude overshoots the expected equilibrium value slightly and then gradually decays to it at a much slower time scale (see discussion in Sec. VII A).
The charge correlation function C(r,τ ) gradually develops reduced Friedel-like oscillations with increasing U , as observed at equilibrium. We find C(r,τ ) < 1 except at distances r < 3 and in the vicinity of the signal front.
In the following we investigate in detail the interplay of characteristic time and length scales and their dependence on the interaction strength.

A. Inside the light cone
Next we discuss the spin correlation function S(r,τ ) inside the light cone. Figure 8 shows the convergence of S o (r,τ ) and S e (r,τ ) to their equilibrium S o (r) and S e (r) values for large times in detail. For large times the odd component is antiferromagnetic, while the even component exhibits a sign change from ferromagnetic at small distances to antiferromagnetic at large distances (blue curves) as discussed in the equilibrium results. The vanishing ferromagnetic component represents a related measure for the extent of a screening cloud [48].
In the following we identify a time scale at which large correlations with the impurity develop inside the light cone, i.e., for distances r v F τ (see Fig. 2). In Fig. 9 (left) we show the difference between the time-dependent spin correlation function and the equilibrium result: S o/e (r,τ ) = |S o/e (r,τ ) − S o/e (r)|. This quantity exhibits contributions from the traveling signal, which show up in the form of large spikes at times τ ≈ r v F . We first focus on the convergence in time at fixed distances r. For times beyond the signal wake τ ∝ r v F , the qualitative structure of correlations has established itself; i.e., Kondo correlations have reached the given distance r. We find that soon after the signal wake S(r,τ ) converges to the equilibrium result exponentially in time, [see Fig. 9(a), inset]. Note that this implies that the curves move "as a whole." We determine S o/e by a single-exponential fit in time of S o/e (r,τ ), successively for distances r ∈ [40,120]. We observe that S o/e (r,U ) is only weakly dependent on r, with odd distances r being especially stable [see Fig. 9(b)], while S e has larger uncertainties and some drift at large r. The uncertainty increases slightly with distance r, which is also due to the smaller available fit intervals in τ . A two-exponential decay as in Sec. VI, featuring also a fast time scale ∝ 1 and independent of U , might be present in the data but cannot be identified due to the presence of the signal at the light cone, which overshadows this fast decay. In general, the fit quality improves with increasing U . Details on the data analysis and uncertainty estimates are provided in Appendix B.
In order to condense this information we consider a mean value,  o/e from α BA may be due to the fact that it is particularly difficult to reach the common asymptotic limit in space and in time for large U . Note that S(0,τ ) S(r,τ ) S(r,τ ), thus the connected correlation function displays essentially the same behavior as S(r,τ ).
We conclude that the formation of Kondo correlations inside the light cone is a process which involves two major time scales. The first time scale is fast and determined by the lattice Fermi velocity v F , which defines the light cone and develops qualitatively correct correlations in S(r,τ ) and C(r,τ ). The second time scale is slow and depends exponentially on U . This process sets in after the qualitatively correct correlations have built up with v F and renormalizes the correlation functions, which then converge at an exponential rate, Eq. (14), α S o/e ∝ T K , to the equilibrium result.
The SIAM is related to its low-energy realization, the antiferromagnetic, symmetric Kondo model via the Schrieffer Wolf transformation [64], which effectively integrates out charge fluctuations. The two models share common features in their low-energy behavior, most prominently the Kondo scale T K . Note, however, that the correlation functions of the two models have very different spatial structures in general. It is therefore interesting to compare our results to recently obtained ones for the Kondo model. In Ref. [57] Lechtenberg et al. studied a coupling quench in the symmetric Kondo model using the time-dependent NRG as well as second-order perturbation theory. Similarly to our results for the SIAM, they found that in the Kondo model spin correlations develop rather rapidly, on the scale of the Fermi velocity. In the linear response to a magnetic field, at equilibrium they observed another, slower time scale similar to 1 T K . Our results unambiguously and quantitatively identify this common slower scale 1 T K beyond linear response, directly from the nonequilibrium time evolution of correlation functions.

085127-11
Charge correlations in equilibrium do not exhibit Kondo physics. We observe that the charge-time-dependent correlation functions C(r,τ ) do exhibit qualitatively the same convergence to equilibrium as the spin correlations S(r,τ ), that is, with a time constant exponentially large in U (not shown). The same analysis as for the spin using Eq. (14) yields respective coefficients for the charge correlation function α C o/e ≈ (0.3 ± 0.1). That is, the exponent is the same as for the spin, albeit with a larger uncertainty. We attribute this to the resolution of the spin in the correlators present in C(r,τ ). Note that this is true neither for the local density, which does not show such a scale, nor for the mean-field result, C mf (r,τ ) ∝ 1.

B. Outside the light cone
For distances r > v F τ , i.e., outside the light cone (see Fig. 2), we find decaying correlation functions S(r,τ ) and C(r,τ ) as a function of r (see Fig. 10). As before, both spin and charge correlation functions show alternating behavior from site to site. The overall magnitude of both correlation functions decreases over time and the charge correlation function is of a smaller magnitude than the spin correlation function for all except very early times. To identify the correlations generated by the quench, we subtract the initial correlation S(r,τ = 0) and C(r,τ = 0) from the time-dependent data.
The second main result of this work is that correlations outside the light cone are power law suppressed, with slightly time-dependent exponents γ S o/e and γ C e . Due to the finite size of the system, we only have a limited set of data available to extract the asymptotic decay outside the light cone. We start the extraction of power-law exponents at distances r s = v F τ + 35 to avoid spurious contributions from the light cone and end it at r e = L − 70 to avoid a bias originating from the boundary at L = 450. From the separate fits for odd/even distances we obtain γ S o ≈ 1.9 ± 0.3 and γ S e ≈ 4.8 ± 0.9. The charge correlation function exhibits a power-law decay γ C e ≈ 1.7 ± 0.3 for the odd component, while the even component's behavior cannot be identified within our numerical accuracy due to the small magnitude of the correlations. The fit has been performed in the same fashion as presented in Appendix B but here we estimate the uncertainty in the γ 's from the fluctuations of the respective γ upon changing the start (r s ) and end point (r e ) of the fit. Within this uncertainty, the values are independent of U and τ .
Considering the fact that extracting exponents from numerical data is challenging, our results agree quite well with two recent studies of similar models exhibiting comparable low-energy physics. First, in Ref. [56] Medvedyeva et al. In an analytic calculation in several limits, neglecting Friedel oscillations, they showed that outside the light cone the commutator spin-z correlation function [Ŝ z 0Ŝ z r (τ )] − , which is related to the linear response to a perturbation, vanishes. For the anticommutator, which is proportional to our S(r,τ ) [see Eq. (4)], however, they obtained a power-law decay r −2 at zero temperature (see Eq. 27 in their work). They found the initial entanglement in the reservoir FS to be responsible for the power-law decay of the anticommutator correlation function.
Moreover, second-order perturbation theory results [57] suggest that initial correlations of the FS transfer to the timedependent correlations outside the light cone. Here again a r −2 power-law decay outside the light cone was obtained, this time for the isotropic Kondo model with antiferromagnetic coupling J . Our study of the symmetric SIAM finds an r −2 decay for S o (r,τ ) outside the light cone, which we attribute to the same structures of the electronic reservoir in the three studies. We are not aware of any previous reports of even-distance decay exponents γ S e ∝ r −5 .

VIII. CONCLUSIONS
We studied the time-dependent formation of the spin screening cloud in the SIAM. Starting from an unentangled state we switched on the impurity-reservoir hybridization and followed the subsequent dynamics in real time. We used the DMRG to obtain ground states and TEBD to obtain spin and charge correlation functions. From these correlation functions we obtained characteristic time and length scales. Our results agree with previous calculations at equilibrium and, for local observables, out of equilibrium. We found that the nonequilibrium correlation functions converge to the equilibrium results for long times.
In the time-dependent data, we identified a linear spreading of signals traveling at the lattice Fermi velocity, which has been referred to as a light cone in recent literature on the buildup of a screening cloud at the Toulouse point of the anisotropic Kondo model [55,56]. We observed a ferromagnetic response in the wake of the signal at the light cone. While Refs.
[55] and [56] identified the Kondo temperature as an inverse time scale in the anisotropic Kondo model outside the light cone, for the symmetric Kondo model it was observed as a time scale in an equilibrium linear response calculation to a magnetic perturbation following an initial fast decay [57]. We observe directly from the nonequilibrium time evolution of correlation functions that, in the SIAM too, the structure of the correlation functions inside the light cone emerges on two time scales. The qualitative core of the correlation functions develops rapidly, at the lattice Fermi velocity. This includes the phase and period of oscillations as well as fixed ferromagnetic and antiferromagnetic domains. These correlations then reach their equilibrium values exponentially slowly in time, defining a dynamic rate which has the same exponential U dependence as the Kondo temperature.
Outside the light cone, we find a power-law decay of the correlation functions ∝ r −γ S/C o/e , with essentially interactionstrength-and time-independent exponents, Eq. (15). In addition to the r −2 decay also observed in the Kondo model [55-57], we find a decay ∝ r −5 .
Our results could be experimentally verified in onedimensional optical lattices featuring two fermionic species. By monitoring the evolution of the spin correlations in time, our findings provide the basis for extracting information about the dynamic scale and, therefore, indirectly about the Kondo screening cloud dynamics as well as the system parameters.
Possible future extensions to this work include the study of the inverse process. Starting from a coupled impurity-reservoir system and investigating the Kondo destruction after switching the hybridization to 0 would allow study of the time-reversed situation. It would also be very interesting to study the effects of a bias voltage on the Kondo screening process using a twoterminal setup as in Ref. [115]. Further interesting extensions involve the study of conduction bands with singularities or testing of implications of the nonequilibrium fluctuationdissipation theorem. Also, calculations away from particlehole symmetry or with applied magnetic fields are feasible.

ACKNOWLEDGMENTS
We gratefully acknowledge fruitful discussion with Sabine Andergassen, Masud Haque, Fabian Heidrich-Meisner, Kerstin T. Oppelt, and Shreyoshi Ghosh. We thank RokŽitko for providing the NRG LJUBLJANA code [101]. This work was partly supported by Austrian Science Fund (FWF) Grant No. P24081-N16 and SFB-ViCoM projects F04103 and F04104 as well as NaWi Graz. M.N. thanks the Forschungszentrum Jülich, in particular, the Autumn School for Correlated Electrons, for hospitality.

APPENDIX A: NUMERICAL DETAILS
In this Appendix we specify details about our numerical analysis carried out via the DMRG [59,60] and TEBD, [62] and we present the DMRG and TEBD parameters used. In addition, we discuss finite-size effects and the convergence as a function of auxiliary parameters specific to the numerical method applied, as well as the stability of the real-time evolution. Our numerical implementation of the DMRG and TEBD is flexible, is parallelized, and exploits two Abelian symmetries: particle numberN and spin projectionŜ z . To find ground states we use the two-site DMRG algorithm with successive single-site DMRG steps. The time evolution is based on a second-order Suzuki-Trotter decomposition of the propagator [61,63].
After extensive studies of the dependence of our results on auxiliary system parameters we found converged results for a Trotter time step of δτ = 0.05t −1 . We used DMRG and TEBD matrix dimensions of χ = 2000 states and always started the DMRG optimization from a half-filled system in the canonical ensemble where alternating up and down spins are chosen as the seed. A detailed discussion is available in Ref. [115] in the context of previous work. Figure 11 shows the equilibrium DMRG calculation of the correlation functions. The influence of the finiteness of the lattice is twofold: (i) The equilibrium spin correlation function S(r) displays an even-odd effect as a function of the total system size L: While for even L, S o (r) converges from above to its L → ∞ value, for odd L it converges from below. S e (r) converges in the opposite way. For odd L an extra spin-↑ gives a spurious total magnetization. For the equilibrium simulations, in the main part of the paper, we have chosen L = 450, since it supports a half-filled and non-spin-polarized system. The spin correlation function at r 150 is converged, as can be seen in Fig. 11 by comparing the L = 450 and L = 300 results. Larger distances are influenced by L because S(r) is a nonlocal quantity. Nevertheless, even for larger distances, no qualitative differences are observed between the L = 450 and the L = 300 data. When performing the time evolution we use L equilibrium + 1 lattice sites, including the impurity, so 085127-13 that the reservoir is nonmagnetized and half-filled. With this choice the correlation functions of the equilibrium and the nonequilibrium simulations become comparable.
(ii) The size of the Kondo screening cloud becomes exponentially large in U . It is therefore important to identify the characteristics of finite-size effects with increasing U . In Fig. 11 (right) we plot data with increasing U for fixed L and study the behavior of S o (r). From U = 0 to U = 6 the correlation function follows a monotonic trend and qualitatively the same behavior. However, the curves for U = 10 and U = 20 are qualitatively different. This indicates that these values of U are too large for the given L, as expected from the size of ξ BA K , which becomes of the order of L ≈ 200 sites here [see Eq. (3)]. In the present work we therefore restrict ourselves to values of U 6 .
Next we show that our time evolution yields a controlled accuracy using a DMRG/TEBD matrix dimension of χ = 2000.
The bipartite entanglement ω(i,τ ) = −tr[ρ L/R (τ )ln(ρ L/R (τ ))] [63] provides an estimate of the time when TEBD becomes unreliable for a fixed χ . This is signaled by a sharp increase in ω. Hereρ L/R denotes the reduced density matrix to the left (L) or to the right (R) of a lattice bipartition at bond i. Figure 12 (left) shows the entanglement increase, which turns out to mostly affect the region next to the impurity and the major propagating signal at r = v F τ . In our simulations we find that χ = 2000 is sufficient to account for the additionally generated entanglement, which is not much larger than in the equilibrium case. In addition, we investigate the direct influence of increasing χ on the interacting spin correlation function S χ (r,τ ) by comparing results using χ = 2000 with results obtained at a smaller χ . Figure 12 ω(i,τ ). We subtracted the ω(i,τ = 0) data to highlight changes caused by the time evolution. Inset: Cuts through the ω(i,τ ) raw data at constant times. The black line is the result of a corresponding equilibrium simulation. The area hidden by the inset is homogeneously dark blue, which corresponds to ω(i,τ ) − ω(i,τ = 0) ≡ 0. Data shown is for U = 3 . Middle: Convergence of the interacting spin correlation function with increasing TEBD matrix dimension χ . Modulus of the residuals |S 2000 (r,τ ) − S χ (r,τ )|, benchmarking the quality of the time evolution with increasing TEBD matrix dimension χ . We show results comparing χ = 2000 with χ = 500 (blue lines) and χ = 2000 with χ = 1000 (orange lines). We show the residuals averaged over distance and interaction strength as a function of time τ . Inset: Spatially resolved residuals plotted at time τ = 2 −1 and for U = 3 . Right: Comparison of the noninteracting spin correlation function as obtained by TEBD, S(r,τ ), and the noninteracting spin correlation as obtained exactly, S exact (r,τ ). Spatially averaged absolute distance. Inset: Spatial resolution for two times, τ = {2,6} −1 . Note that each blue (orange) line belongs to one data set only, which is alternating.
with systematic signatures at the light cone and beyond it, while the interior of the light cone looks chaotic. The results are almost independent of U . We find that the space r and interaction U averaged deviation grows over time and becomes of the order of O(5 × 10 −4 ) for χ = 500 and O(1 × 10 −4 ) for χ = 1000 within the reachable simulation time. Furthermore, for U = 0 we compare the correlation functions obtained via TEBD with the numerically exact ones [Eq. (8)] in Fig. 12 (right). As one can see, the maximum deviation occurs at the boundary far from the impurity, with a maximum error of ≈1 × 10 −5 .
We conclude that for simulations of nonlocal correlation functions within the DMRG and TEBD in the SIAM the major limiting factor is the computation time T ∝ L(χ ) 3 . This is due to the large matrix dimensions χ needed to achieve accurate results and is, furthermore, complicated by the fact that the SIAM exhibits exponentially long correlation lengths, which require large lattice sizes L. The accuracy can be controlled by benchmarking against exactly known U = 0 data and, for finite U , by increasing the TEBD matrix dimension χ . Furthermore, all the scales extracted in the text, α S o/e and γ C/S o/e are retrieved from two subtracted correlation functions, in which we expect errors to further compensate.

APPENDIX B: EXTRACTION OF THE DYNAMIC ENERGY SCALE
In the following we provide details of the data analysis of the dynamic scale α o/e as discussed in Sec. VII A, which is valid for both even and odd distances. First, we obtain the time dependence of the spin correlation function by performing a nonlinear fit in time τ to the spin correlation function for fixed distances r and given interaction U : S(τ |r,U ) (see Sec. VII A), using f (φ = (c 1 , (r,U )),τ ) = c 1 e − (r,U )τ with two fit parameters φ. The data are single exponential plus oscillations and exhibit an eventual systematic bias close to the lattice border and due to the signal front at the light cone. For each r we manually choose intervals [τ s (r,U ),τ e (r,U )] for the fit in time in order to minimize these influences. Typically we choose fit intervals which start r s ≈ 10 sites behind the light cone and extend up to r e ≈ 250 for large U . For small U the data become noise before this r e is reached and we adjust the end points accordingly. We estimate the fit uncertainty (r,U ) by φ i ≈ √ C ii , where C = (J † J)η 2 is the estimated covariance, J = ∂f (φ,τ i ) ∂α j is the fit Jacobian, and η 2 = r T r N τ (r,U )−p is the mean square error defined by the residuals r i = S(τ i |r,U ) − f (φ,τ i ) on N τ (r,U ) data points in time S(τ i |r,U ). These estimates are consistent with those obtained from fluctuations upon changing τ s (r,U ) and τ e (r,U ). Second, we condense the r dependence by averaging (r,U ) over distances r. We make use of a Bayesian approach with Gaussian error statistics. We obtain the weighted mean value (U ) = 1 P r 1 (r,U ) 2 (r i ,U ) with P = r 1 (r,U ) 2 and a weighted error estimate (U ) = 1 √ P , where the weights are obtained from (r,U ). Third, we obtain the U dependence of the exponent considering data for (U ) for N (U ) = 6 data points at U = {1,2,3,4,5,6} . The data (U ) can be fitted very well by a single exponential in U : (U ) = c 2 e −αU . The same scheme as in the first step is used to estimate the uncertainty . We assume correlated data, i.e., η 2 = r T r N eff , with N eff ≈ N(U )−p 2N corr ≈ 6−2 2×6 , which enlarges the uncertainty by a factor of √ 3 compared to the naive value.