Analysis of the buildup of spatiotemporal correlations and their bounds outside of the light cone

In non-relativistic quantum theories the Lieb-Robinson bound defines an effective light cone with exponentially small tails outside of it. In this work we use it to derive a bound for the time evolution of the correlation function of two local disjoint observables if the initial state has a power-law decay. We show that the exponent of the power-law of the bound is identical to the initial (equilibrium) decay. We explicitly verify this result by studying the full dynamics of the susceptibilities and correlations in the exactly solvable Luttinger model after a sudden quench from the non-interacting to the interacting model.


Introduction
The rising interest in the out-of-equilibrium dynamics of physical systems has posed the question of the propagation of information in many-body quantum systems. While relativistic models intrinsically contain the notion of a maximal speed, i.e. the speed of light, it was historically unclear whether there is an equivalent bound in generic non-relativistic models. A milestone in this context was achieved by Lieb and Robinson in 1972. [1] They demonstrated that the dynamics of a lattice system with finite-range interactions is effectively constrained inside a region of the spatiotemporal plane, with exponentially small "leaks" outside of it. As a consequence, any perturbation effect requires a finite amount of time to spread from one point of the system to another. Mathematically, the theorem provides a bound for the operator norm of the commutator of two local observables O A and O B defined over two disjoints sets A and B and reads where O A (t) (O B (t)) denotes the time-evolved observable with norm ||O A || (||O B ||), N min = min(|A|, |B|) and c and ξ are additional constants that depend on the interaction, the local Hilbert space dimension and the lattice structure. L is the distance between the sets A and B (smallest number of vertices connecting the two regions) defined over the lattice. The parameter v, which has the dimension of a velocity, is called Lieb-Robinson velocity. The value of this important parameter is anyway not fixed by the theorem, but depends on the lattice structure and the interaction. [2] The previous expression which is also called "Lieb-Robinson bound" (LR bound) tells us that the commutator between the two operators is exponentially suppressed for L > v|t|. Keeping the distance L fixed, the value of the bound grows in time until t = L/v where it becomes O (1). It thereby divides the full spatiotemporal plane into a part inside and outside of a light cone which is defined by the Lieb-Robinson velocity v.
Since the "leaks" outside of the light cone practically vanish, the two observables effectively commute in this region (also called "space-like" region in analogy to relativistic theories). Therefore, the finite velocity v defines the fastest possible information propagation. Historically, the theorem first only applied to translationally-invariant quantum spin lattices with short-range interactions until some of the requirements were loosened by several improvements and generalizations. [3,4,5,6] In its most recent version, it does not require a finite local Hilbert space, since it was extended to systems with an unfinite local Hilbert space, for example lattice bosons. [4,5,7,6,8] Moreover, it can be generalized to long-range interacting systems [9,10,11,12,2], where a faster-than-ballistic propagation,with algebraic or even logarithmic light cone, is allowed. However, the propagation observed in some specific models appears to be significantly slower than the one predicted by the bounds [13,14,15,16]. Despite the interest for the latter, we will focus in this paper on short-range interacting systems, where a ballistic light cone is present and postpone generalizations with lesser requirements to further studies.
In the discussion about the LR-bound it is important to notice that the commutator of two observables as it appears in Eq. (1) is closely related to linear response theory (see for example Ref. [17]). It measures the response of a system to a local perturbation: the system is shaken at one spatial point and the corresponding effect is observed at some other point later in time. Therefore, the commutator is also known as the susceptibility. In a semi-classical argument one models the spreading of the signals to be governed by quasiparticles that move through the system. [18] The light cone velocity, which was a non-fixed parameter in the LR result, can then be determined by the microscopic theory to be equal to the maximal group velocity of the excitations defined in the system. Experiments, however, do not only focus on susceptibilities, but they can also access other quantities like correlation functions of the form | O A (t)O B (t 0 ) |. These are directly connected to the anticommutator {O A (t) , O B (t 0 )}. Notable developments which focus on the study of the dynamics of the correlation functions have been achieved in cold-atomic gas experiments. [19,20] Additionally, the spreading of correlations has been studied using several advanced numerical techniques like tDMRG [21,22,23,24,25], time-dependent variational Monte-Carlo [26,13], and artificial neural network techniques [27]. Finally, exact analytic results on the dynamics in solvable models [28,29,30,31,32] have helped to understand the mechanisms and to build up a physical picture where one can interpret correctly the observed phenomena.
The original LR-bound as in Eq. (1), however, does not provide information about the behavior of such quantities which strongly depend on the initial state. Nevertheless, there are few mathematical bounds for the correlation function which rely on assumptions about the initial state. They require either the exponential [33,8,34] or algebraic [35] decay of initial correlations or a spectral gap above the ground state. All these results are derived using the LR-bound in different ways. The main difference between the commutator and the anti-commutator, i.e. the susceptibility and the correlation function, is the behavior outside the light cone. The correlation function in fact exhibits "leaks" in that region, meaning that correlations are non-zero from the very beginning of the time evolution. This does not spoil the main result of the LR bound, correlations between two space-like separated regions do not imply superluminal propagation of information, like in the EPR-paradoxon [36]. This difference between commutator and anticommutator has been shown explicitly in the Kondo model [17]. There, one can see a light cone in the susceptibility (with exponentially small "leaks"), while the correlation function has algebraic "leaks" outside of it. The extension of the LR bound to non-equal-time correlations and its relation to the initial state are the main results we want to present.
This paper is organised as follows: In the next section we derive a bound for the non-equaltime correlation function in the case where the initial correlation function is algebraically decaying. We will demonstrate that also for this case the correlations have a light cone with a velocity determined by the Lieb-Robinson velocity. This generalizes other results present in the literature for equal-time correlations [33,35]. Moreover, we connect the leaks outside of the light cone with the initial (equilibrium) correlation function. Thereafter, we derive exact expressions for the dynamics of a system described by a Luttinger liquid. We calculate and discuss the correlation buildup and decay regimes for the commutator and anticommutator of the density-density correlation for both zero temperature and a thermal initial state.

Bound for the correlation decay outside the light cone
In this part we prove that the correlation decay outside of the light cone is defined by the spatial decay of the initial (equilibrium) correlation function.Bounds for the correlation functions at equal time have already been found for both exponentially [33] and algebraically [35] decaying correlations in the initial state. We extend these results to the non-equal time case, i. e. for The c subscript indicates that we consider the connected correlation function defined as Moreover, optimizing our bound analytically we show how the spacelike region decay is closely connected to the initial correlations. Physically, the theorem states that any leaks outside the outer light cone defined by the LR velocity are bounded by the initial (equilibrium) correlations which have been "dragged" during the time evolution, see Fig. 1. The following results hold for all models for which the famous Lieb-Robinson (LR) theorem in its form displayed in Eq. (1) applies. In the derivation we consider a quantum spin lattice with short-range interaction (also called "local Hamiltonian"). As previously mentioned, it is possible to loosen this condition to a more generic setting which anyway does not change the main message of our result [4,5,6]. One example is to consider a continuous system (not on a lattice). In section 5 we will explicitly demonstrate how the bound naturally arises in the case of the continuous Luttinger model. For the moment, we postpone the generalization to systems where more general LR-bounds have been found (for example for a system with long-range interaction) to further studies. In order to bound the correlation function we assume that the time evolution is governed by a local Hamiltonian and that the initial state has the property that all connected correlation functions decay with a power law (algebraically) given by Here, O A and O B are two normalized (||O A ||, ||O B || ≤ 1) operators which initially act on the compact disjoint supports A and B, which are apart by the distance L. The additional parameterã defines a positive cut-off if L = 0. Following Ref. [33], we now study how fast the effective part of the time-evolved operator O A (t) (O B (t)) grows in space. A different way to bound the equal time correlation function using the LR bound can be found in Ref. [35]. Mathematically, it is clear that O A (t) requires the full support over the entire system, because the unitary evolution contains all powers of the Hamiltonian. Yet, we will exploit the LRbound that tells us that effectively it only acts on a finite region which grows linearly with time. Therefore, we define an operator, O l A A (t), that acts on the region A extended by l A sites in all directions like O A (t) and outside of this region as a trivial unity. The outside region of the ball is denoted by S. It is defined by The normalization makes sure that the norm is given by the inside part. Additionally, One can now use the Lieb-Robinson bound to prove that the operator norm of the difference of the full observable and the one only acting inside the ball is bounded by One notices that the difference is basically unbounded (exponentially large) if l A < vt, i. e. inside the light cone and exponentially decreasing, if it is outside. Hence, this states that effectively O A (t) only acts in the region A enlarged by l A in all directions (if l A is chosen appropriately). With this result one finds an upper limit for the connected correlation function at different times: It is important to note that we have two light cones with an identical velocity v, but different final times t and t 0 . In the next step the last term which by construction only acts inside the two compact balls is bounded using the assumption Eq. (2). The initial support of A has grown by l A after time t and of B by l B after time t 0 in all directions. As a consequence the distance between the two regions shrinks by l A + l B such that one assumes that for all times Plugging this into Eq. (5) one gets For t = t 0 , this result is the same as the one found by Kastner in Ref. [35] for algebraically decaying initial states up to different constants. However, Eq. (5) depends on l A and l B in a non-trivial way. We now optimize it analytically to get a bound depending just on L, t, and t 0 which are the physically relevant variables. It is important to note that l A and l B need not to be constant, but must fulfill l A + l B < L. In this way the supports of O l A A and O l B B do not overlap. We then hit the light cone if l A + l B = L. By setting l A = (L + v(3t − t 0 ))/4 and l B = (L − v(t − 3t 0 ))/4 we bound both terms for t > 0 and t 0 > 0 to with ξ = 4ξ,Ã = 2ã andC = 2 βc . It is worth to notice how the two optimized parameters l A and l B depend non-trivially on the two times t and t 0 . l B decreases with t and it increases with t 0 while l A increases with t and decreases with t 0 . In summary, one can bound the correlation function by a sum of a subleading exponential and a dominant algebraic bound. Outside the light cone, where L > v(t + t 0 ) ã, the exponential bound is negligible and the total bound is a power-law where the exponent is given by the initial (equilibrium) correlation decay (at t = t 0 = 0) where we replacedÃ withã, because both parameters are small compared to L − v(t + t 0 ). Inside the light cone, where L < v(t + t 0 ), the correlation function is quasi-unbounded as  8)) over the spatiotemporal plane outside of the light cone. One notices that its scale is different from the initial correlation depicted in blue (Eq. (2)), however, the decay exponent is identical. The inside of the light cone is marked by the grey area and is basically unbounded.
expected from the Lieb-Robinson bound. Close to the light cone, where L ∼ v(t + t 0 ), the light cone is dominated by the reciprocal cut-off valueÃ −β with the exponent β defined by the equilibrium spatial decay of the correlation function.
It is possible to formally extend the previous results in a more rigorous mathematical language where one would speak, for example, of a set of vertices that create a lattice and a general metric. In this publication, however, we choose to stick to a simple physical picture. In this picture the locality of the Hamiltonian means that the interaction only couples nearby degrees of freedom which are sitting on a lattice that can be infinite. It is, of course, also possible to apply the results to other spin systems with a more general short-range interaction. As stated above, the requirement is that the considered model fulfills the conditions for the Lieb-Robinson bound as given in Eq. (1). The mathematical result in Eq. (8) will be explicitly verified in section 5 for the Luttinger model presented in section 3. The calculation is an example how the theorem also holds in the case of a continuous theory with point-like interaction. Analogously to our derivation, it is also possible to extend the result [33] by Brayvi et al. to non-equal times. The assumption there is that the initial correlation is given by

The Luttinger model
In this section we briefly review the Luttinger model in the bosonization picture. This topic is now standard textbook knowledge and the readers already familiar with that can skip this section.
An increasing interest in one-dimensional electron models has initially motivated Luttinger to present an interacting fermion model in 1963. [37] This model resembled the one that Tomonaga had studied 13 years before and which featured the splitting into left-and right-moving excitations for the first time. [38] However, Luttinger assumed a simpler interaction among the fermions and treated the Dirac sea as being infinite. [39] This was corrected two years later by Mattis and Lieb who presented an exact solution for a generalized model by introducing bosonic field operators. [40] They also demonstrated the breakdown of the Fermi liquid theory in one dimension. Further important contributions contained the calculation of one-and two-particle correlation functions [41] and the emergence of a nonuniversal power-law decay and the finding that the low-energy behavior suggests that it belongs to his newly defined "Tomonaga-Luttinger universality class". [42,43,44] In the past decade the interest in the Luttinger model has revived in the context of nonequilibrium dynamics induced by quantum quenches. Cazalilla, for example, has studied the behavior of the model after different quenches [45,46,47], while Dora and Moessner investigated the out-of-time-ordered density correlators. [48] This work aims at extending the findings of Cazalilla in his extensive study of correlation functions. We focus on investigating the correlation buildup and their behavior in the full space-time plane.
The initial Hamiltonian of interest describes a one-dimensional system with interacting (spinless) electrons and is given in terms of fermion creation and annihilation operators. Since only the low-energy excitations are considered in the following, one assumes that the dispersion relation takes on a linear form. To be able to apply the bosonization machinery, two additional requirements are needed. [49] First, one needs a discrete momentum number k which is achieved by imposing antiperiodic boundary conditions. Second, one extends the dispersion relation to infinity by adding "positron states" which do not have a physical interpretation. In the low-energy regime these states, however, do not contribute. After a complete bosonization procedure the Tomonaga-Luttinger model can be written in terms of bosonic collective excitations which represent a particle-hole pair given by the density operator ρ α . Hereby, α = L, R denotes the left-and right-moving excitations with opposing sign s L/R = ∓1, respectively. Our system of length L is initially prepared as a thermal state with respect to the free theory which only contains a kinetic part with v F the Fermi velocity. The double dots denote fermionic normal-ordering with respect to a reference state which is chosen to be the Fermi sea. [44] The thermal state is set by a bath at inverse temperature β (β = ∞ describes the ground state). At time t = 0 the bath is removed and an interaction quench is performed, i. e. we turn on the interaction parameters g 2 and g 4 adding the interaction part H int to the previous Hamiltonian. It is given as The interaction is entirely point-like and describes scattering processes within one branch (g 4 ) and the forward scattering between the two branches (g 2 ). Transforming the full Hamiltonian into momentum space using yields where several constants which do not affect the dynamics were omitted. To be able to sum over all q we introduced a regularization exp(−a|q|/2) which is discussed in the next section. One remark to the choice of interactions is necessary: The simplified model with only pointlike repulsion bears implications which, however, become irrelevant in the low-energy regime (long-wavelength or q ∼ 0). [47] The bosonic operators b † q (b q ) create (annihilate) a particle-hole excitation with momentum q. This Hamiltonian can now be diagonalized with a Bogoliubov rotation given by with B q denoting the new bosonic operators in the diagonal basis. The Bogoliubov angle φ is given by tanh(2φ) = g 2 /(v F + g 4 ). The diagonalized Hamiltonian defines the dispersion relation ω(q) = ν|q| with the renormalized velocity ν = ((v F +g 4 ) 2 +g 2 2 ) 1/2 . In the diagonalized form the time-dependence is trivial such that one finds a resulting transformation with f (q, t) = cos(ν|q|t) − i sin(ν|q|t) cosh(2φ) and g(q, t) = i sin(ν|q|t) sinh(2φ).

Density-density correlation function
In the following we study the density-density correlation function n(x, t)n(0, t 0 ) where the average is taken with respect to an initial thermal state (at T = 0 the ground state). Separating the density operator in left and right movers n(x, t) = ρ L (x, t) + ρ R (x, t), the full density-density operator contains the four contributions where we omitted terms that oscillate fast with x. This corresponds to averaging over ∆x ∼ (2k F ) −1 . [50] This averaging procedure erases any non-universal power-law decay that is non-Fermi liquid like. Similar quantities like Eq. (16) have already been studied in the literature, however the full calculation at non-equal time has not been done to give a closedform expression. [46,47,48] In Ref. [47] the authors discussed the sum which appears in the calculation in the infinite-time limit without performing it. The calculation of the full density-density correlation function requires plugging Eq. (12) and (15) in (16) which results in four terms. In order to be able to sum over q, a regularization exp(−|q|a/2) defined by the "effective bandwidth" a −1 is introduced.
[43] One usually considers the case where a → 0 + , however a small but finite value is closer to the real world, since it constrains possible scattering events. It will become clear in the following that this basically only affects the region around the light cone, where a finite a limits the peak height. Away from it one can safely let a → 0 + .

Zero temperature
One finds at T = 0 for a finite system size L that with γ = g 2 /ν denotes the quench parameter which causes the time-dependence (since the g 4 term commutes with the non-interacting Hamiltonian). We will now sometimes use the notation τ for any combination of x ± ν(t ± t 0 ). To discuss the functions U (τ, a) and V (τ, a) we naively let a → 0 + . This limiting produce gives the known result where U (τ ) = d(τ |L) −2 , where d(τ |L) = L sin(πx/L)/π is the conformal distance which becomes d(τ ) = τ in the thermodynamic limit. [46] The imaginary part of Eq. (17) converges to a sum of derivatives of δ-distributions on the left-and right-moving light cone (where τ = 0), because V (τ ) = −π∂ τ δ(τ ) and consequently only contributes there. However, considering a to be very small, but non-zero, changes the behavior close to the light cone. In this case one needs to approximate U (τ, a) and V (τ, a) for small τ and a and finds U (τ ≈ 0) ≈ (τ 2 − a 2 )/(τ 2 + a 2 ) 2 with U (τ = 0, a) = −1/a 2 on the light cone and V (τ ≈ 0) ≈ 2aτ /(τ 2 + a 2 ) 2 with S(τ = 0, a) = 0. The arising divergences create two light cones in both directions (we will consider x > 0 from now on): an outer light cone at x − ν(t + t 0 ) = 0 and an inner light cone at x − ν(t − t 0 ) = 0. The first thing to note is that the behavior close to a light cone is independent of the system size (it is identical to the result in the thermodynamic limit). While U (τ, a) is symmetrically algebraically peaked with a maximum value −1/a 2 , the width of the V (τ, a)-double peak depends strongly on a. ¿From its functional form one can see that V (τ, a) is pointsymmetric with respect to τ = 0: Considering the peak around a right-moving light cone (τ = x − ν(t + t 0 ) = 0) the time-like region is defined by x < ν(t + t 0 ) or τ < 0, while the space-like region is given by τ > 0. Thus, it has a negative peak inside the light cone, while the outside peak is positive. For a left-moving light cone it is τ > 0 (τ < 0) for the time-like (space-like) region, but the negative sign in front of the corresponding function V (x + ν(t + t 0 )) assures the same behavior holds for both light cones (see also Ref. [ [48]]). Finally, one can recover the result by Cazalilla [[46]], if one considers only the right-moving density-density correlation at equal times (t = t 0 ). This requires neglecting the terms proportional to γ 1 + γ 2 , since these terms derive from L − R and R − L contributions. Moreover, one needs an additional factor of 1/2, since both the L−L and R−R are identical in this case.
The function that describes the initial correlations can be found by setting t = t 0 = 0 in Eq. (17). These correlations are those that were present before the time-evolution. They are given by For x a it is approximately n(x, 0)n(0, 0) ini T =0 ≈ −1/(2π 2 x 2 ). Generally, we call the correlations that are not γ-dependent pre-quench correlations (also indicated by "pq"). For equal-time correlators (t = t 0 ) these contributions are identical to the initial correlations, however, if t = t 0 the different measurement times change the correlation function. In the space-like region (outside of the light cone) where x ν(t + t 0 ) > ν(t − t 0 ) a, however, the difference is negligible, because The third part of Eq. (17) are the quench-induced contributions (indicated by "qi") which depend on γ. They are thoroughly discussed in section 6.

Non-zero temperature
In the following we extend the previous results to non-zero temperature where the initial state is a thermal state at inverse temperature β.
At T > 0 additional terms add to the zero temperature result. The full result for a finite system of size L can be found in the appendix. We now focus on the result in the thermodynamic limit, where one can execute the infinite sum and gets: with the derivative of the DiGamma-function DiG 1 (z) = d dz ln Γ(z). This general result is new and has not been calculated before. The behavior close to the light cone is dominated by the sharp peaks of the T = 0 case, since the DiGamma-functions only give a constant contribution there (2 DiG 1 (1) = π 2 /3). In the space-like region from the light cone the limit a → 0 + can be taken and the additional terms can be summarized to DiG 1 (1 − iτ /(βv F )) + DiG 1 (1 + iτ /(βv F )) = (βv F ) 2 /τ 2 − π 2 sinh −2 (πτ /(βv F )). The algebraic decay terms now exactly cancel the T = 0 terms such that it finally becomes n(x, t)n(0, t 0 ) therm = 1 where the temperature induced correlation length ξ = βv F /π has been introduced. Another form of writing the expression above is taking the T = 0 result and substituting the conformal distance d(τ ) with the newly defined (thermal) distance d(τ, β) := ξ sinh(τ /ξ) as hinted at in Ref. [[46]]. This is, of course, only true in the thermodynamic limit or if the correlation length scale is way shorter than the system size (ξ = βv F /π L). The thermal correlation length ξ now separates two different regimes: if x ± ν(t ± t 0 ) ξ we observe again the algebraic decay regime as in the T = 0 case which is now followed by an exponentially decaying regime for x ± ν(t ± t 0 ) ξ. Both regimes will be discussed in the form of the anticommutators in section 6.

Explicit bound for the Luttinger model
In this section we explicitly verify the bound presented in section 2 for the two observables O A = n(x) and O B = n(0) at non-equal times t and t 0 . The compact support contains only the spatial point B = 0 and A = x, respectively. To do this we define a function C ini (x) which describes the absolute value of the initial correlations At t = t 0 = 0 and zero temperature we know the exact expression for the Luttinger model from Eq. (17): Here, we used that n(x) = 0. As discussed earlier C ini (x) describes a power-law decay with an exponent β = 2, whilec = 1/2π 2 . From the theorem Eq. 8 we thus expect a bound which has the same exponent. We prove that this is indeed fulfilled by bounding the full density-density correlation function for t > t 0 > 0. Although we have calculated the exact correlation function also for the inside of the two light cones in the previous sections, Eq. (8) states that it is quasi-unbounded inside of the outer light cone. We therefore only focus on the decay outside of the outer light cone i. e. where x ≥ ν(t + t 0 ). In this region we can estimate that The bounds are tight if t = t 0 = 0. However, bounding each term one finds that | n(x, t)n(0, t 0 ) c | = 1 Analogously to the initial case, there is again the upper limit of Eq. (26) with C ini (x − ν(t + t 0 )) ≤ 1/(2π 2 )/((x − ν(t + t 0 )) 2 + a 2 ) that takes on the form used in the theorem Eq. (8).
Thus, we have explicitly demonstrated that in the Luttinger model the correlations outside of the outer light cone can be bounded by the initial correlations shifted by the current position of the light cone. In a physical picture one can think of the correlation as being "dragged" with the quasiparticles which move with the light cone speed. As soon as any quasiparticle has passed a region the correlation function is basically unbounded. In a real system the bound does not need to be tight, but the real expression can fall off even quicker (see discussion of the exact decay in section 6).

The spatiotemporal correlation build-up
In this section we want to study the dynamics of the anti-commutator and commutator of the density-density correlation function that have been calculated in part 4 at zero and non-zero temperatures. Hence, we define the functions where [·, ·] ± denotes the anticommutator or commutator. Since the anticommutator is just twice the real part of the density-density correlation and the commutator the corresponding doubled imaginary part of Eqs. (17) and (22), one immediately finds the resulting expressions. The average is taken with respect to the ground state or the thermal state defined by an initial distribution exp(−βH kin )/Z at inverse temperature β. Hereby, we consider the thermodynamic limit, since we focus on the decay regimes and the finite size effects have been discussed elsewhere. [46] In the following, when we study the anticommutator we often subtract the initial contributions C ini + (x) := C + (x, 0, 0) such that we actually analyze This is done, because we want to study how the correlations buildup without having the background of the equilibrium (if we do include it, it is mentioned).

Analysis of the commutator C − (x, t, t 0 )
As has been discussed before, [48,47] the commutator is just a complex number (independent of the inverse temperature β) and as such it is independent of the initial state (also compare Eq. (17) and (23)). Due to the very sharp δ-distribution form of the commutator, it becomes clear that the only contribution results from the left-and right-moving light cones where x ± ν(t − t 0 ) = 0. Physically, this can be understood by regarding the Luttinger model as being effectively a relativistic theory with a well-defined speed of light ν. The linear dispersion creates an identical ballistic propagation for all quasiparticles (particle-hole excitations) with exactly this speed (no curvature in the dispersion). Therefore, there is also no signal at equal times when t = t 0 . Generally, it follows that we do not observe any "leak" outside the light cone and one sees that the LR-bounds are trivially fulfilled. Not taking the limit a → 0 + broadens the peak to a Lorentzian form with FWHM= 2a. However, the arising contributions outside the light cone do not violate the LR-bounds.

Analysis of the anticommutator C + (x, t, t 0 )
This simple behavior can, of course, not be seen when studying the anticommutator which is not bounded by the Lieb-Robinson bounds without further restrictions (see derivation of a bound in part 2). While the commutator is independent of the initial state (it is an operator norm), the anticommutator strongly depends on it. In this sense, the anticommutator has to be understood as a statistical property, since one needs to perform many weak measurements to determine its exact shape. The equal-time (t = t 0 ) anticommutator has already been calculated in Ref.
[ [46]], but with a focus on the long-time convergence to a generalized Gibbs ensemble and not on the spatiotemporal correlation buildup.
The analysis starts with the T = 0 case, where the anticommutator is defined by twice the real part of Eq. (17). Moreover, we simplify by considering a → 0 + and taking the thermodynamic limit L → ∞. This can be done, because a finite a is only needed in the vicinity of the light cone where τ ≈ 0. If a finite a is required, we will from now on explicitly discuss the change.
In order to study the correlation buildup induced by the interaction quench we subtracted the initial correlations (if t = t 0 , it is initial = pre-quench correlations) which can be found by setting t = 0 in Eq. (17): They show the usual Fermi gas or free fermion behavior, since initially the model is noninteracting. Having subtracted this contribution the quench-induced part of the anticommutator reads In the last expression we have expressed the anticommutator in terms of the initial correlations. From now on, we will always use this description to be able to relate the quench-induced correlations to the initial ones. The absolute value of this function |C + (x, t)| is depicted in Fig. 2a for the right-moving light cone where the initial correlations have been subtracted. In a physical picture the figure shows the correlations which are induced by the movement of the quasiparticles after the quench. [18] In the figure one sees that for a fixed x the anticommutator is zero initially, but then growing to the sharply peaked light cone on the diagonal with time. This corresponds to a vertical cut through the plot. This behavior can be understood as follows: At t = 0 particle-hole excitations get excited by the global quench which then spread in both directions with velocity 2ν and possibly correlate all regions that they have passed. This is called the correlation buildup. Since all quasiparticles move at the same speed, they will continue their coherent movement forever. Thus, there is a huge correlation peak (named light cone) that indicates that all excitations come from the same point x = 0 (at t = 0) to which they are therefore highly correlated. In order to discuss the growth rate and the decay regime outside of the light cone one needs to approximate Eq. (30) in the space-like region, i. e. where x 2νt. This results in In this region, one finds that the growth of correlations is proportional to t 2 /x 4 resulting in algebraically increasing "leaks" outside the light cone if one considers a fixed x (vertical cut of Fig. 2a). Closer to the light cone, x ∼ 2νt, the full algebraic growth ∝ 1/(x − 2νt) 2 takes over before the light cone peaks to −1/a 2 (or infinitely high if there is an infinite number of excitations, i. e. a → 0 + ). While the correlation function grows with time for a fixed x, the behavior is different if you focus on the decay outside of the light cone for a fixed time t: As has been shown in section 2 the decay outside the light cone can be effectively bounded by the initial (equilibrium) decay that moves with the speed of the quasiparticles. Only the amplitudes might change. In the Luttinger liquid we notice in section 5 that the decay is only tightly bounded in the close vicinity of the light cone. In this sense, the bound we derived in section 2 can be interpreted as the initial correlations being "dragged" with the quasiparticles during the time evolution. Mathematically, this is given by the simple substitution x → x − v (t + t 0 ) in the correlation function. In the very deep space-like region where x 2νt, however, the correlation function decays faster with 1/x 4 . This is depicted in Fig. 2b which cuts Fig. 2a horizontally for a fixed time. Finally, we discuss the structure that is visible inside the light cone. There we can identify a second signal where correlations vanish (see Fig. 2b). There, the time-independent stationary, or time-like, (x 2νt) correlations C tl +,alg (x) = −(γ 1 + γ 2 − γ 2 )C ini +,alg (x) cancel the term coming from the light cone. In other words, it is exactly at x = (2νt)/ √ 3 where C + (x, t) = 0. Since the time-like decay is stationary, it is also the remaining contribution in the infinitetime limit. Adding the initial correlation again one notices that the initial decay exponent of the non-interacting Fermi gas prevails, but with a reduced strength −(γ 1 + γ 2 − γ 2 ) + 1 which is bounded by 1/2 from below. In this sense, we have an interacting Fermi liquid. This also confirms our expectation to not see any non-universal power-law, since we averaged over ∆x ∼ (2k F ) −1 . [50] There is an interesting difference to results in the literature where only the R − R densitydensity correlation has been studied. [46] There, the infinite-time limit is always larger than the initial correlations (for a repulsive interaction). However, taking the full physical fermion density into account shows the physically expected behavior: the interaction limits the movement of the quasiparticles and thus reduces the initial correlations. The exponent of the decay, however, remains the same.
The previous results drastically change when one considers a thermal initial state. In that case, one needs to analyze the real-part of Eq. (23) (with a → 0 + ). The initial correlation decay is found to be C ini + (x) = −1/(πξ sinh(x/ξ)) 2 . We depict the full spatiotemporal behavior without the initial contributions as defined by Eq. (28) in Fig. 3a. Analyzing the defining equation one now observes two different regimes which are governed by the temperature via the new correlation length ξ. If x ξ one finds the algebraic regime as in Eq. (29). Physically, in this regime we see length scales much shorter than the one induced by the temperature. However, for x ξ the decay becomes exponential with In summary, the temperature introduces a new correlation length which separates two different decay regimes. This behavior of the time-independent initial correlation has already been described in Ref. [[50]]. We now show that it extends to the quench-induced correlations by studying the full dynamics. We start with the correlation buildup far away from any causal interference, i. e. the space-like region (x 2νt) where we see that the two initial decay regimes persist in the full spatiotemporal plane. Closer to the light cone, where |x ± 2νt| < ξ, we observe an algebraic decay region with similar expressions for the space-like and time-like case as the ones found for T = 0. Outside of this region, however, we find an exponential decay of correlations which has been observed in many other systems. The space-like expression for this region is   30)). The signal on the diagonal is the light cone with an algebraic (power-law-like) decay. A horizontal cut through the plot for time t/2ν = 10 is shown on the right (Fig. 2b) with an algebraic decay indicated by the dashed line. The grey shaded area marks the inside of the light cone with a stationary algebraic correlation decay starting at x = 0. Due to different signs this contribution cancels the inner decay coming from the light cone leading to an inner structure with vanishing correlation at x root = t/(2ν √ 3).
where the initial correlation decay is now given by Eq. (32). It looks similar to the algebraic version (Eq. (29)), but with a different prefactor, a different decay of the initial correlations (exponential instead of algebraic) and ξ takes over the role of x in the buildup process which is again proportional to t 2 .
In conclusion, one finds that the algebraic regime present around the peaks on the light cones at T = 0 shrinks with increasing temperature. It is followed by a new, temperature induced exponential decay. Both regimes reflect the initial correlations and can be seen in the log − log-plot in Fig. 3a and 3b.
We now switch to the non-equal-time anticommutator where t = t 0 . We can write C + (x, t, t 0 ) = C(x, t − t 0 , t + t 0 ), i. e. it only depends on the difference and sum of the two times. We can therefore write the expressions in terms of time differences t−t 0 and find that t+t 0 = t−t 0 +2t 0 where t 0 is fixed. Fig. 4a shows the entire spatiotemporal plane of the anticommutator C + (x, t − t 0 ). We show the positive time axis which corresponds to values t > t 0 . Physically, this means that we   (23)). The Fermi velocity is chosen such that ξ = β = 30. The signal on the diagonal is the light cone with an algebraic (power-law-like) decay for x < √ 3ξ followed by an exponential decay. The algebraic regime is found between the white lines. A horizontal cut through the plot for time t/2ν = 10 is shown on the right (Fig. 3b). The grey shaded area marks the inside of the light cone with a stationary algebraic correlation decay starting at x = 0. Due to different signs this contribution cancels the inner decay coming from the light cone leading to an inner structure with vanishing correlation at x root (see text for formula).
Fermi velocity 2ν in both directions. At the first measurement at time t 0 they have propagated by 2νt 0 . Now we measure again at a different position x after an additional time t − t 0 . We find two distinct signals: One describes the outer light cone at x = ν(t + t 0 ) = ν(t − t 0 + 2t 0 ) which represents the spread of quasiparticles due to the quench until the first measurement after t 0 . The other one is the light cone on the diagonal, i. e. at x = ν(t − t 0 ) and shows the correlation of the first measurement to the second measurement. Since both light cones have different signs, we again see a light cone of vanishing correlations in between the peaks. This double light cone structure has also been observed in other models and appears, because the Luttinger model with its well-defined propagation speed is effectively a relativistic theory. [17] In the following we want to express the corresponding decay regimes in terms of the prequench correlation decay. The pre-quench correlations are those that we would see without quenching (γ = 0) for x ν(t − t 0 ). They read The leading (slowest) decay term is identical to the initial correlation decay at equal times in Eq. (29). The next order expansion term reflects the change of the correlation due to different times and is negligible if one considers a point far away from the light cone. Thus, we will approximate C pq + (x, t, t 0 ) ≈ C ini + (x). The decay outside of the outer light cone (space-like region) shows a behavior closely related to the equal-time correlation function. At T > 0 the algebraic decay which is the only decay at T = 0 is enclosed by an exponential decay   (Fig. 4b). Due to different signs the two light cone contributions at (t − t 0 )ν = 30 and (t + t 0 )/ν = 110 cancel exactly in the middle of the peaks. We again see an algebraic decay regime outside of the outer light cone. This is directly related to the initial correlations (see text).
similar to the equal-time case. In the time-like region one finds that the quench induced correlations decay proportional to 1/t 2 , since Hence, C tl +,alg (x, t, t 0 ) vanishes in the infinite-time limit and only the initial correlations remain. However, for a finite time difference t − t 0 Eq. (37) demonstrates that the spatial decay inside of the light cone is independent of x in the first order approximation (see Fig. 4b).

Conclusions
In this paper we discussed the spreading of local observables following a quantum quench in many-body quantum systems. In particular, we elucidated the difference between the commutator and anticommutator and the influence of the initial state on the dynamics. First, we derived a bound for the expectation value of local observables at non-equal times for systems with algebraically decaying spatial correlations in the initial state. This result is obtained using the Lieb-Robinson bound and generalizes the findings by Brayvi et al. [33] and Kastner [35] to non-equal time correlations. By optimizing the bound we demonstrated that it effectively describes a light cone which velocity is the Lieb-Robinson velocity. Outside the light cone, correlations are non-zero but algebraically suppressed with a power law determined by the equilibrium (initial) correlations. Moreover, the same power law determines the fast buildup of correlations close to the light cone, where the maximum occurs. The dependence on the initial state reflects the difference of the bounds: whereas the LR-bound is an operator norm and thus completely independent of the initial state, the bound for correlation functions depends non-trivially on the initial state. In a next step, we verified the bound by studying the time evolution of an exactly solvable model: the Luttinger liquid. We investigated the dynamics induced by a sudden switch on of the interactions starting from different initial states: the ground state, which has algebraically decaying correlations, and a thermal mixture, where correlations decay exponentially at long distances. Using bosonization techniques, we derived exact expressions for the density-density correlation at different positions and times. We discussed both the commutator and the anticommutator of the mentioned observables in the full spatiotemporal plane. While the commutator gives a simple delta function on the light cone and as such trivially fulfills the LR-bound, the anticommutator has a much richer behavior that can be compared to the bound we derived. First of all, it exhibits a clear light cone for both exponentially and algebraically decaying initial correlations, as expected. In the space-like regime (large distance and small time), we find a correlation buildup which is algebraic for the zero temperature case with a power law determined by the initial state. This result does not saturate the bound derived in section 2 but has a slower growth. In the case of non-zero temperature, we dicussed the appearance of the thermal correlation length also outside of the light cone where the correlation buildup is exponential, yet again closely related to the initial state. The situation close to the light cone is much more interesting. We find that the correlation function grows algebraically fast in that region independently on the initial state. The exponent of this growth is exactly the one that controls the spatial decay of the observable in the initial state at zero temperature. Comparing the exact calculation to the general bound we derived before, we find that the correlation buildup in the Luttinger liquid is very close to the maximal allowed buildup. Moreover, the connection between the initial state and the time evolution suggested by the bound is perfectly respected in the Luttinger liquid. Finally, in the time-like region, the bound diverges exponentially fast, meaning we cannot use it to describe the real-time evolution. However, our exact calculations can be interpreted via the quasiparticles approach developed by Cardy and Calabrese. The spectrum of the excitations is linear and consequently they all have the same velocity. From our calculations we find that another ballistic signal with vanishing correlations is present inside the time-like regime. This is due to the fact that the light cone and the relaxation value (infinite-time value) of the anti-commutator have opposite signs. This confirms that even in a model with an extremely simple spectrum, observables can show a non-trivial behavior that cannot be guessed by more general results like the Lieb-Robinson bound.

Acknowledgments
We are grateful for valuable discussions with M. Schmitt.
Funding information Authors are required to provide funding information, including relevant agencies and grant numbers with linked author's initials.Financial support from the Deutsche Forschungsgemeinschaft (DFG) through SFB/CRC 1073 (Project B03) is gratefully acknowledged.