{'fraction_non_alphanumeric': 0.040700502675531054, 'fraction_numerical': 0.015749148694665153, 'mean_word_length': 4.249175444196191, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Real-time energy dynamics in spin-1/2 Heisenberg chains

We study the real-time dynamics of the local energy density in the spin-1/2 XXZ chain starting from initial states with an inhomogeneous profile of bond energies. Numerical simulations of the dynamics of the initial states are carried out using the adaptive time-dependent density matrix renormalization group method. We analyze the time dependence of the spatial variance associated with the local energy density to classify the dynamics as either ballistic or diffusive. Our results are consistent with ballistic behavior both in the massless and the massive phase. We also study the same problem within Luttinger Liquid theory and obtain that energy wave-packets propagate with the sound velocity. We recover this behavior in our numerical simulations in the limit of very weakly perturbed initial states.


I. INTRODUCTION
The understanding of transport properties of lowdimensional systems with strong correlations still poses viable challenges to theorists. These include, on the one hand, the fundamental problem of calculating transport coefficients for generic models such as the Heisenberg chain, 1,2 and on the other hand, the theoretical modeling of experiments that typically require the treatment of spin or electronic degrees of freedom coupled to phonons, in particular, in the case of the thermal conductivity. 3,4 Most theoretical work has focussed on the linear-response regime, in which the properties of current-current autocorrelation functions determine transport properties (see Refs. 1,2 for a review).
More recently, the out-of-equilibrium properties of onedimensional systems have evolved into an active field of research, one reason being recent advances in experiments with ultracold atoms. 5 These have paved the way for studying the dynamics of quantum many-body systems that are driven far away from equilibrium in a controlled manner, with little or no coupling to external degrees of freedom. Much attention has been paid to the question of thermalization, typically studied in so-called quantum quenches (see Ref. 6  usually do not induce any finite net currents (of either spin, energy, or particles), we will be particularly interested in set-ups that feature finite net-currents. Such situations are realized in, for instance, the sudden expansion of particles in optical lattices after the removal of trapping potentials. 7 Further examples are spin and/or particle currents induced by connecting two regions with opposite magnetizations or by letting two particle clouds collide (see, for instance, Refs. [8][9][10][11].
Theoretical work in this context ranges from the expansion dynamics of bosons and fermions in optical lattices [12][13][14][15][16][17][18] over the dynamics of wave-packets in spin chains, [19][20][21][22][23][24][25][26][27][28] to the demonstration of signatures of spincharge separation in such set-ups. 29,30 In the aforementioned examples, non-equilibrium situations were studied with either finite spin or particle currents. In our work, we address the energy dynamics for a model that is prototypical for systems with strong correlations, namely the spin-1/2 XXZ chain: where S µ i and µ = x, y, z are the components of a spin-1/2 operator acting on site i and S ± i are the corresponding lowering/raising operators. The global energy scale is set by the exchange coupling J, ∆ is the exchange anisotropy in the z-direction, and L denotes the number of sites. Equation (1) describes either interacting quantum spins or, via the Jordan-Wigner transformation, 31 spinless fermions.
Specifically, we follow the time-evolution of the local energy density h i starting from initial states that are far from the ground state of Eq. (1) and that feature an inhomogeneous profile in the local energy density (see Fig. 1 for a sketch). We emphasize that, in the main part of our work, we choose the initial conditions such that only finite energy currents exist, whereas the spin (particle) density is constant during the time evolution, hence all spin (particle) currents vanish. Obviously, an initial state with an inhomogeneous spin density profile leads to both finite spin and energy currents, and we revisit this case, previously studied in Refs. 20 and 28, as well.
Our work is motivated by and closely related to a specific experiment on a spin-ladder material. Many lowdimensional quantum magnets are known to be very good thermal conductors with heat predominantly carried by magnetic excitations at elevated temperatures. 32,33 Examples for materials that exhibit particularly large thermal conductivities are (Sr,La,Ca) 14 Cu 24 O 41 (Refs. 34 and 35) and SrCuO 2 (Ref. 36). While these experiments are carried out under steady-state conditions and in the regime of small external perturbations, more recently, time-resolved measurements have been performed on La 9 Ca 5 Cu 24 O 41 (Ref. 37). For this spin ladder material, two approaches have been implemented: A timeof-flight measurement, in which one side of the sample is heated up with a laser pulse and the time-dependent response is recorded on the other side. Second, a nonequilibrium local heat distribution was generated in the surface of the material by shining laser light on it. It is possible to record the heat dynamics via thermal imaging that uses the response of an excited thin fluorescent layer placed on top of the spin ladder material.
It is the latter case that we mimic in our work: The time evolution of local energy densities induced by inhomogeneous initial distributions. We utilize the time-dependent density matrix renormalization group (tDMRG) [38][39][40][41][42] technique. It allows us to simulate the dynamics of pure states whereas in the experiment, temperature likely plays a role. Our work thus addresses qualitative aspects in the first place, while a direct comparison with experimental results is beyond the scope of this study: The goal is to demonstrate that in a spin-1/2 chain described by Eq. (1), the energy dynamics is ballistic, irrespective of how far from equilibrium the system is and also irrespective of the presence or absence of excitation gaps. To this end, we use the same approach as in Ref. 20: We classify the dynamics based on the behavior of the spatial variance σ 2 E (t) of the local energy density: The ballistic case is σ 2 Our main result for the XXZ chain, based on numerical tDMRG simulations, is that energy propagates ballistically at sufficiently long times, independently of model parameters (such as ∆). One can then interpret the prefactor V E in σ 2 E (t) = V 2 E t 2 as a measure of the average velocity of excitations contributing to the expansion. The velocity V E can be calculated analytically and exactly in non-interacting models, which (in the absence of impurities or disorder) typically have ballistic dynamics, and we consider two examples: (i) the noninteracting limit of the XXZ-Hamiltonian (∆ = 0), i.e., spinless fermions and (ii) the Luttinger liquid, which is the universal low-energy theory in the continuum limit of Eq. (1) for |∆| < 1. We show that our tDMRG results agree with the exactly known expansion velocity V E in these two examples.
Our main result, namely the numerical observation of σ 2 E (t) ∼ t 2 independently of initial conditions or model parameters such as the exchange anisotropy ∆, is consistent with the qualitative picture derived from linear-response theory. Within that theory transport properties of the XXZ chain have intensely been studied in recent years, both the energy 43-47 and the spin transport. 1,2,[47][48][49][50][51][52][53][54][55][56][57][58][59][60][61][62] Ballistic dynamics is associated with the existence of non-zero Drude weights. Since the total energy current of the anisotropic spin-1/2 chain is a conserved quantity for all ∆, the thermal conductivity κ(ω) diverges in the zero-frequency limit and is given by Re κ(ω) = D E δ(ω), where D E is the thermal Drude weight. [43][44][45][46] This behavior is different from the spin conductivity σ(ω). This quantity takes the form Re σ(ω) = D s δ(ω) only at the noninteracting point ∆ = 0, whereas for 0 < ∆ ≤ 1, many numerical studies 1,2,60,61 indicate D s (T > 0) > 0, with a finite weight at finite frequencies, though. Therefore, for 0 < ∆ ≤ 1, Re σ(ω) = D s δ(ω) + σ reg (ω). Recent fieldtheoretical and numerical work suggests that the regular part σ reg (ω) of σ(ω) in massless phases is consistent with diffusive behavior. 54,56,62 A finite value of the current-current correlation function in the long time limit is associated with a finite Drude weight. Finite Drude weights can be traced back to the existence of conservation laws, 43,60 and in consequence, a potential relation between integrability 63 and ballistic behavior -in the sense of non-zero Drude weights -has been intensely discussed (see, e.g., Refs. 1,2,47,60-62 and further references cited therein). Very recently, Prosen has presented results that provide a lower bound to the spin Drude weight that is non-zero for ∆ < 1. 60 This is in qualitative agreement with earlier exact diagonalization studies. 2,47,48 The particular point ∆ = 1 is still discussed controversially: 47,49,52,54,[59][60][61] First, no finite lower bound to the Drude weight is known, 60 and second, the qualitative results of exact diagonalization studies seem to depend on details of the extrapolation of finite-size data to the thermodynamic limit and the statistical ensemble that is considered. 47,48,59 Our approach that analyzes the time-dependence of spatial variances, albeit restricted to the analysis of densities, is numerically easily tractable and is an alternative to the numerically cumbersome evaluation of current correlation functions. tDMRG has, for instance, been applied to evaluate current-current autocorrelation functions in the thermodynamic limit. 54 However, the accessible time scales are quite limited (t ∼ 10/J), making an unambiguous interpretation of the results difficult and the approach is not applicable to non equilibrium. Our approach allows us, at least in principle, to study the entire regime of weakly perturbed states to maximally excited ones. An earlier analysis of spin-density wave packets in various spin models has yielded the following picture (all based on the time-dependence of the spatial variance): 20 In massless phases, ballistic dynamics is seen, whereas in massive ones, examples of diffusive dynamics have been identified. It is important to stress that the observation of a variance that increases linear in time is a necessary condition for the validity of the diffusion equation.
Finally, to complete the survey of related literature, recent studies have addressed steady-state spin and energy transport in open systems coupled to baths, with no restriction to the linear-response regime. 28,61,[64][65][66] These studies suggest spin transport to be ballistic in the gapless phase of the XXZ spin chain and to be diffusive in the gapped phase with a negative differential conductance at large driving strengths. The heat current has been addressed in Ref. 64 where Fourier's law has been validated for the Ising model in a tilted field.
A by-product of any tDMRG simulation is information on the time-evolution of the entanglement entropy. While this is not directly related to this article's chief case, it nevertheless provides valuable information on the numerical costs of tDMRG simulations. Qualitatively, speaking (see the discussion in Ref. 42 and references therein), the faster the entanglement growth is, the shorter are the time scales that can be reached with tDMRG. We here show that the quenches studied in this work generate a mild logarithmic increase of entanglement, which is why this problem is very well suited for tDMRG. Such a behavior is typical for so-called local quenches. 67 This result might be useful for tDMRG practitioners.
This paper is organized as follows: First, we introduce the model and the quantities used in our analysis in Sec. II. Section III A reviews the framework of bosonization, which is applied in Sec. III B to give an analytical derivation of ballistic spin and energy dynamics in the low-energy case, valid in the massless phase of Eq. (1). Sections IV and V contain our numerical results. First, we study the energy dynamics in the absence of spin currents in Sec IV. To this end, we generate an initial state consisting of a variable number of ferromagnetic bonds in the center of an antiferromagnetic chain. We calculate the time evolution of these states under Eq. (1) finding ballistic energy dynamics independent of the phase and the strength of the perturbation. To supplement these findings we derive an observable, which depends on the local currents, and whose expectation value is timeindependent whenever σ 2 E (t) ∼ t 2 . The numerical calculation of this quantity indicates ballistic dynamics as well. Section V revisits the scenario of Ref. 20 where local spin and energy currents are present during the dynamics as we start from states with an inhomogeneous spin density. In that case the energy density shows ballistic dynamics in the massless phase with a velocity matching the bosonization result in the limit of small perturbations. In the massive phase we observe a different behavior of the two transport channels, i.e., ballistic energy dynamics while the spin dynamics looks diffusive. 20 Finally, we summarize our findings in Sec. VI. Additionally, we discuss the entanglement growth induced by coupling two regions with an opposite sign of the exchange coupling in the Appendix.

A. Preparation of initial states and definition of spatial variance
In this work, we focus on spin-1/2 XXZ chains of a finite length L given by Eq. (1) where our goal is to study the dynamics of an inhomogeneous distribution of the local energy density originating from a local quench of system parameters. The inhomogeneous distributions are generated by preparing the system in the respective ground states of the following Hamiltonians that are perturbations of H XXZ from Eq. (1). First, where h i is defined in Eq. (1), and second, where In the first case, we quench site-dependent exchange couplings. In this scenario we obtain initial states with large local energy densities. Typical initial states that are ground states of Eq. (2) are shown in Fig. 2: These states have b bonds with ferromagnetic J i < 0 in the center while the rest has antiferromagnetic J i > 0. We refer to this setup as the J i quench.
In the second case, the dynamics is driven by an inhomogeneous spin density, enforced by an external magnetic field applied in the initial state. This allows us to generate smooth spatial perturbations of h i with small differences in energy compared to the ground state of Eq. (1). We refer to this setup as the B 0 quench. A more detailed discussion of the initial states generated by a J i quench will be given in Sec. IV A. The B 0 quench was introduced in detail in Ref. 20.
The definition of the local energy density from the Hamiltonian Eq. (1) is not unambiguous. For instance, it is always possible to add local terms to the Hamiltonian whose total contribution by summation over all lattice sites vanishes. However, this seeming ambiguity can be resolved up to constants by requiring that any block of adjacent lattice sites m i=l h i is Hermitian and yet to have the same structure as the total Hamiltonian H. These details seem to be rather specific, yet for the definition of the appropriate local energy density within the Luttinger liquid description, see below, these formal considerations are important. For the XXZ chain the local energy density is therefore determined by the bond energies h i . To classify the dynamics of a density e i we study its spatial variance.
where µ is the first moment of e i . The e i are the normalized distribution linked to the energy density via is the energy difference between the initial state |ψ init [i.e., the ground state of either H J init or H B init ] and the ground state |ψ 0 of Eq. (1), both energies measured with respect to the unperturbed Hamiltonian from Eq. (1): On physical grounds, the energy density should be normalized by the amount of energy transported by the propagating perturbation. This is well approximated by the energy difference δE between the initial state and the ground state of Eq. (1), as we have verified in many examples. In some cases, though, the propagating energy is, on a quantitatively level, better described by estimating the area under the perturbations, as δE may also contain contributions from static deviations from the ground state bond-energies in the background. Nevertheless, δE does not depend on the overall zero of energy and is an obvious measure of how far the system is driven away from the ground state. This, all together, justifies our definition of the e i .
To remove static contributions depending only on the initial distribution e i (t = 0), we subtract σ 2 E (t = 0) and study δσ 2 is expected to grow quadratically in time in the case of ballistic behavior, where V E has the dimensions of a velocity.
For diffusive behavior, we expect, from the fundamental solution of the diffusion equation, 69 that δσ 2 E (t) ∼ Dt grows linearly in time, where D is the diffusion constant (see, e.g, the discussion in Ref. 20). Within linear response theory the diffusion constant can be related to transport coefficients via Einstein relations, see e.g. Ref. 70. To be clear the observation of δσ 2 E ∼ t 2 or δσ 2 E ∼ t is a necessary conditions for the respective type of dynamics and time-dependent crossovers are possible.

B. Spatial variance in the non-interacting case
For pedagogical reasons and to guide the ensuing discussion we next calculate the spatial variance in the noninteracting limit of Eq. (1), i.e., at ∆ = 0. Using the Jordan-Wigner transformation, we can write the Hamiltonian as where c † i creates a spinless fermion on site i. A subsequent Fourier transformation diagonalizes the Hamiltonian: Since we will compare with numerical results on systems with open boundary conditions, we obtain Next, we compute k , we finally obtain, after straightforward calculations: i.e., ballistic dynamics independently of the initial state. Terms linear in t will be absent if in the initial state, the density is symmetric with respect to its first moment, i.e., e µ+δ = e µ−δ and if the wave packet has no finite centerof-mass momentum at t = 0 already. In the remainder of the paper we will work under these two additional assumptions that are valid for all initial states considered in our work. The prefactor V 2 E is given by: where v k = ∂ǫ k /∂k and is the difference between the momentum distribution function (MDF) in the initial state and the one in the ground state of Eq (1).
Since we use open boundary conditions, we compute n k from We can also express δE via δn k : The expression Eq. (13) suggests that V E is the average velocity of excitations contributing to the propagation of the wave packet. Characteristic for ballistic dynamics, V 2 E is fully determined by the initial conditions through δn k .
For completeness, we mention that an analogous calculation can be done for the spatial variance σ S of the spin density. This quantity is defined as The normalization constant N measures the number of propagating particles. The spin density is, in terms of spinless fermions: The result for the spatial variance of the spin density is with Although we started from the Hamiltonian for ∆ = 0, we stress that Eqs. (12), (16), (13), and (17) are valid for any dispersion relation ǫ k , irrespective of the presence of a gap, provided that k has the meaning of a momentum.

C. Energy current
Another aspect worth noting is that the time-evolving state carries a nonzero energy current, a situation that usually does not appear in the case of a global quench. From the equation of continuity for the energy density one can derive the well-known expression for the local energy current operator 43 where˜ S = (S x , S y , ∆S z ). With periodic boundary conditions, the total current J E = i j E i is a conserved quantity, i.e., [H, J E ] = 0 (see Ref. 43). On a system with open boundary conditions such as the ones that are well-suited for DMRG, this property is lost, yet the dynamical conductivity still has a quasi-Drude peak at very low frequencies, reminiscent of the true Drude peak Reκ(ω) = D E δ(ω) of a system with periodic boundary conditions. 68 The latter form is recovered on a system with open boundary conditions as L → ∞ (Ref. 68), showing that ballistic dynamics due to the existence of globally conserved currents can still be probed on systems with open boundary conditions.
To connect the local energy currents to the spatial variance of the time-dependent density one can rewrite the time derivative of σ 2 E (t) using the equation of continuity, assuming no current flow to sites at the boundary (this assumption is justified in our examples as long as we restrict ourselves to times before reflections occur at the boundary in our simulations): If σ 2 E (t) = V 2 E t 2 + b and µ = µ(t), then using J E = 0 leads to: If we interpret this equation as an operator equation, then we see that we can define a quantity J * E via: If for a given initial state and over a certain time window, J * E (t) = const, then we have identified a regime with ballistic dynamics, δσ 2 E (t) ∼ t 2 . If J * E (t) =const holds for all times and initial states, then J * E is a conserved quantity, [H, J * E ] = 0. This is the case at ∆ = 0, the non-interacting limit of Eq.
We emphasize that we have here identified a operator that connects the phenomenological observation of a quadratic increase of σ 2 E (t) to the local energy currents. In ballistic regimes, its expectation value becomes stationary.
For completeness, we mention an analogous result in the diffusive regime where σ 2 E ∼ t. Then, expectation values of the operator are time independent. Obviously, similar expressions can be written down for the spatial variance associated with the spin density.

III. PROPAGATING ENERGY AND SPIN WAVE-PACKETS IN A LUTTINGER LIQUID
In the gapless phase, i.e., for |∆| < 1, the low-energy and low-momentum properties of the XXZ chain can be described by an effective Luttinger liquid theory. 71 In the following we want to analyze the energy density and the spin dynamics of the XXZ chain in this exactly solvable hydrodynamic limit. Specifically, we show that at least asymptotically for large times, the spatial variance always grows quadratically both in the case of spin and energy dynamics. In addition, we work out the precise dependence of the prefactor in front of the t 2 increase of the spatial variance on system parameters. Since our DMRG results to be presented in Secs. IV and V show that σ 2 E (t) ∼ t 2 at any ∆, we did not investigate the influence of marginally relevant perturbations at ∆ = 1 on the wave-packet dynamics. In passing, we mention that in the massive phase, where the appropriate low-energy theory is the sine-Gordon model, the expansion velocity could also be derived at the Luther-Emery point (this case was studied, in, e.g. Refs. 25,26).

B. Ballistic dynamics in the gapless phase
Within the Luttinger liquid description for ∆ < 1, an initially inhomogeneous local energy density profile always propagates ballistically independently of the details of the perturbation as can be seen from general arguments. For the effective low-energy Hamiltonian the probability distribution e(x, t) associated with the local energy density is given by where |ψ init is the initial state, For the exact definition of the fields ϕ Assuming a L ↔ R symmetry in the initial state, i.e., a state with zero total momentum, one obtains for the variance from Eq. (5) for all times t with V E = v g = u. This results can also be obtained from evaluating Eq. (13) in the continuum limit.
In the case of an initial L ↔ R asymmetry in the initial state we get δσ 2 E (t) → (v g t) 2 for t → ∞, but the short time behavior may differ. Thus, within the validity of a Luttinger liquid description the energy transport is always ballistic for all initial conditions. This is evident from a physical point of view as all excitations propagate with exactly the same velocity v g , the left-movers to the left and the right-movers to the right. Note that the applicability of a Luttinger liquid description is manifestly restricted to cases in which the initial energy density profile is a smooth one in the sense that the associated excitations do not feel the nonlinearity of the fermionic dispersion relation. Thus, the time-evolution starting from initial profiles such as the ones shown in Fig. 2 are beyond the scope of this low-energy theory.
In analogy to the above arguments, the dynamics of spin-density wave-packets is also ballistic in the XXZ chain for ∆ < 1 in the Luttinger liquid limit. In the bosonic theory, the spin density is proportional to ρ L (x)+ ρ R (x) up to a constant, see Sec. III A. The associated probability distribution ρ(x, t) = Q −1 ρ L + ρ R /2π, with Q = dx ρ L + ρ R /2π, can again be separated into a left-and a right-moving contribution, i.e., Thus, similar to the case of the energy dynamics, one finds ballistic behavior for |∆| < 1 consistent with the numerical results of Ref. 20.

IV. DMRG RESULTS FOR THE Ji QUENCH
Now we turn to the numerical simulations. Using the adaptive time-dependent DMRG 38-42 method we can access the real-time dynamics of initial bond energy distributions. Within this approach we can probe the microscopic dynamics including the time dependence of bond energies or the entanglement entropy starting from various initial states in an essentially exact manner without limitations in the range of parameters. We discuss the pure energy dynamics in the absence of spin currents induced by the J i quench in this section. We detail the construction of initial states and their specific features, then move on to the analysis of the time evolution of the bond energies. We calculate the spatial variance and the related quantity J * E and discuss the emergent velocities of the energy dynamics. Within the numerical accuracy of our simulations we find a quadratic increase of σ 2 E (t) in all cases studied. However, it seems that for a J i quench a large number of different velocities contribute as opposed to the Luttinger Liquid theory result, the latter valid at low energies. Our study of the energy current during the time evolution and the time evolution of the expectation value J * E (t) , defined in Eq. (21), gives additional insights into short-time dynamics and further validates the conclusion of ballistic energy dynamics.

A. Initial states
Let us first describe the typical shape of initial states induced by a J i quench on a few bonds in the middle of the spin chain. To be specific, in the Hamiltonian Eq. (2) we set which provides us with initial states with an inhomogeneous energy density profile with a width of 2b of the ferromagnetic region. Outside this ferromagnetic region we obtain antiferromagnetic nearest-neighbor correlations.  Fig. 2, the system forms a region with ferromagnetic nearest-neighbor correlations in the middle of the chain. Note that for ∆ = 1, h i is the sum of the nearest-neighbor transverse and longitudinal spin correlations, the latter weighted with ∆. In the regions with antiferromagnetic J i > 0, the local energy density oscillates, reflecting the antiferromagnetic nearest-neighbor correlations. Figure 3 shows the energy difference δE. As a function of b, the energy difference δE increases linearly once the smallest possible ferromagnetic region has been established. The minimum energy difference δE = E init − E 0 is of the order of 2J, i.e., initial states that are only weak perturbations of the respective ground state cannot be generated using a J i quench.
At the isotropic point ∆ = 1, we can explain the dependence of the initial state on the width b in a transparent manner. The ground state energy per site for the antiferromagnetic ground state is known from the Bethe Ansatz to be lim L→∞ E 0 (L)/L = −ln(2) + 1/4, 74 while for the ferromagnetic ground state E 0 /(L − 1) = 1/4 excluding the boundary sites which gives rise to a very small system-size dependence. By growing the ferromagnetic region symmetrically with respect to the center of the chain and taking E 0 (L) from the unperturbed ground state with open boundaries, we obtain states with an en- ergy that increases as for our finite system size (δE 0 is simply an off-set). Equation (31) exactly reproduces the data for ∆ = 1 shown in Fig. 3.
Furthermore, at ∆ = 1, the total spin is a conserved quantity. Since the ground state calculation only respects the conservation of magnetization (S z tot = 0) we obtain a hierarchy of states with S > 0. This can be easily understood by considering the block structure of the initial state. Taking, e.g., a total of L = 100 spins and assuming a ferromagnetic region of only two spins (i.e., b = 1), the two ferromagnetic spins are fully polarized with a total spin of S = 1 while each of the antiferromagnetic blocks have 49 spins and therefore a total spin of S = 1/2. Thus, the total spin of the whole chain is S tot = 2. Increasing the width of the ferromagnetic region by one, i.e., to b = 2, we have S = 2 in the middle and the antiferromagnetic blocks are of even length, which both have S = 0 in their ground state. This pattern repeats itself upon increasing the length 2b of the ferromagnetic region.

B. Time evolution of bond-energies after a Ji quench
Now we focus on the time-evolution of the local energy density induced by the aforementioned perturbation. At time t = 0 + , we set all J i = J and then evolve under the dynamics of Eq. (1). The DMRG simulations are carried out using a Krylov-space based algorithm 75,76 with a time-step of typically 0.25J and by enforcing a fixed discarded weight. We restrict the discussion to times smaller than the time needed for the fastest excitation to reach the boundary.  1). Besides the outer rays that define a light cone structrue, Fig. 4 unveils the presence of more such rays inside the light-cone. Since our particular initial states have a sharp edge in real space, there ought to be many exciations with different momenta k contributing to the expansion.

Ji quench: Spatial variance
Our main evidence for ballistic dynamics in both phases is based on the analysis of the spatial variance, shown in Fig. 5. Fitting a power-law (straight lines) to the data, i.e., σ 2 E (t) − σ 2 E (0) = αt β yields a quadratic increase with β ≈ 2, classifying the dynamics as ballistic.
In order to estimate uncertainties in the fitting parameter α, we compare this to the results of fitting a pure parabola σ 2 E (t) − σ 2 E (0) = V 2 E t 2 to the data. Typically, V 2 E deviates from α by about 10% while the exponent of the power-law fit is usually different from 2 by 5%. As an example, for ∆ = 0.5, b = 1 we obtain β = 2.03 and α = 0.53 vs. V 2 E = 0.6J 2 . The main reason for the deviation of β from 2 is, in fact, that the short time dynamics is not well described by a power law at all over a bdependent time-window. We shall see later, in Sec. IV C, that the ballistic dynamics sets in only after the block of ferromagnetically correlated bonds has fully 'melted'. Indeed, by excluding several time steps at the beginning Fitting a power-law (straight lines) to σ 2 E (t) − σ 2 E (0) = αt β yields a quadratic increase with sufficient accuracy, classifying the dynamics as ballistic. For instance we find α = 0.53, β = 2.03 for ∆ = 0.5 and b = 1 [black circles in (a)]. We do not find any qualitative difference between the massless (|∆| ≤ 1) and the massive (∆ > 1) phase. The deviations between the fit and the tDMRG data in the ∆ = 1.5 curves at the largest times simulated are due to boundaries. of the evolution from the power-law fit, we observe that β → 2 and α → V 2 E . Therefore, we will present results for V 2 E , obtained by fitting σ 2 E (t) − σ 2 E (0) = V 2 E t 2 to our tDMRG data.

Exploiting SU (2) symmetry at ∆ = 1 for the Ji quench
Before proceeding to the discussion of the expansion velocity V 2 E , we wish to discuss the long-time limit, which can be accessed in the case of ∆ = 1. Since our perturbation is proportional to the operators for the local energy density, global symmetries of the unperturbed Hamiltonian are respected by the initial states of the type Eq. (30). Therefore, at ∆ = 1, we can exploit the conservation of total spin S, a non-Abelian symmetry. This can be used to push the simulations to much longer times, since we can perform the time evolution in an SU (2) invariant basis. 77 The number of states needed to ensure a given accuracy is reduced substantially compared to a simulation that only respects U (1) symmetry. Therefore, we can work with larger system sizes and study the longtime dynamics of the energy density. As we can reach longer times, we can also analyze and discuss finite-size effects for ∆ = 1 here. Figure 6 shows our result for the . We still find a quadratic increase of σ 2 E (t) and thus ballistic dynamics for times up to t ∼ 60/J and in addition, the prefactor does not depend on the system size. Both simulations were carried out keeping the discarded weight below 10 −4 which requires at most m = 900 states using only U (1) symmetry on L = 100 sites versus a maximum of m = 400 using SU (2) for L = 200 sites.

Expansion velocity
The results for V 2 E are collected in Fig. 7 and plotted as a function of δE for ∆ = 0, 0.5, 1, 1.5. In the noninteracting case, ∆ = 0, V 2 E is constant for b ≥ 2, while at b = 1 (the smallest possible δE), V 2 E = 0.5J 2 . For all ∆ > 0, V 2 E slightly decreases with δE and V 2 E is much smaller than v 2 g given by Eq. (24), suggesting that indeed, many velocities contribute during the expansion of the energy wave-packet.
Intuitively, one might associate the decrease of V 2 E , which is a measure of the average velocity of propagating excitations contributing to the expansion, to band curvature: The higher δE, the more excitations with velocities smaller than v g are expected to factor in.
It is instructive to consider the non-interacting limit first by comparing the numerical results obtained from a time-evolution with exact diagonalization to the analytical (and also exact result) from Eq. (13). To that end we need to compute the MDF [see Eq. (14)] of the initial state. Our results for ∆ = 0, which are shown in Fig. 8, unveil a peculiar property: The J i quench always induces changes at all k, i.e., the system is not just weakly perturbed in the vicinity of k F . This is not surprising since our initial states have sharp edges in real-space (compare E t 2 as a function of δE for ∆ = 0, 0.5, 1, 1.5 and Ji quenches with b = 1, 2, 3, 4, 5 (for ∆ = 1.5, we show b = 1, 2, 3 only). Fig. 2). Moreover, the J i quench changes the MDF in such a way that δn k (b) = n init k (b)−n k is point-symmetric with respect to k F = π/2 , where k F is the Fermi wavevector. As Fig. 7 shows, V 2 E as extracted from fits to δσ 2 E (solid symbols) and V 2 E from Eq. (13) [open symbols] perfectly agree with each other, as expected.
The MDF of initial states for the interacting systems are also such that δn k = 0 at all momenta and we may therefore conclude that the observation V E < v g is due to the fact that the J i quench induces many excitations with velocities smaller than v g (compare the data shown for ∆ = 0.5 shown in Fig. 8). Of course, Eq. (13) is not directly applicable to the interacting case since, first, it does not account for the correct eigenstates at ∆ = 0 and second, in general, h i = J(S + i S − i+1 + h.c.)/2 . Nevertheless, by numerically calculating δn k for the interacting system and by using the renormalized velocity in Eq. (13) instead of J [i.e., J → v g (∆)], we obtain an estimate for This reproduces the qualitative trend of the tDMRG results for V 2 E as we exemplify for ∆ = 0.5 in Fig. 7. To summarize, the overall picture for the time evolution of the bond energies after a J i quench is: Energy propagates ballistically with an expansion velocity V E that is approximately given by Eq. (33). Combined with the observation that on a finite system, a J i quench induces changes in the MDF at all momenta k, we conclude that many excitations contribute to the wave-packet dynamics, resulting in V E < v g , both in the non-interacting and in the interacting case.

C. Energy currents
To conclude the discussion of the J i quenches we present our results for the local energy currents at ∆ = 1 in Fig. 9. By comparison with Fig. 4(b), we see that the local current is the strongest in the vicinity of the wave packet. The energy current in each half of the system becomes a constant after a few time steps, i.e., J E L/2 := j E i reaches a constant value. We plot the absolute value of J E L/2 for ∆ = 0.5, 1, 1.5 for b = 1 in Fig. 10(a). The qualitative behavior is independent of ∆: As soon as the initial perturbation has split up into two wave packets, we have prepared each half of the chain in a state with a constant, global current J E L/2 =const. For a system with periodic boundary conditions, the total current J E = i j E i is a conserved quantity. 43 Since the effect of boundaries only factors in once these are reached by the fastest excitations, we directly probe the conservation of a global current with our set-up, after some initial transient dynamics. Therefore, we can link the phenomenological observation of ballistic wave-packet dynamics to the existence of a conservation law in the system.
While the currents J E L/2 clearly undergo some transient dynamics [see Fig. 10(a)], we have derived a quantity in Sec. II, called J * E , whose expectation value is stationary if σ 2 E ∼ t 2 . We now numerically evaluate J * E (t) from Eq. (21) which provides an independent probe of ballistic dynamics. Figure 10(b) shows our results for ∆ = 1 and J i quenches with b = 1, 2, 3, 4, 5. It turns out that J * E (t) is indeed constant at sufficiently large times, consistent with the observation of δσ 2 E ∼ t 2 . In Sec. IV B, we have noted that δσ 2 E ∼ t 2 at short times t b/J. This renders J * E (t) a time-dependent quantity over the same time window: Clearly, the time window over which J * E (t) = const, depends on b [see Fig. 10(b)], which suggests that the deviation of ballistic dynamics is associated to the 'melting' process of the region with ferromagnetic correlations. We have carefully checked that these observations are robust against errors in the calculation of time derivatives in Eq. (21) induced by the finite time step. Since J * E (t) is time-dependent (at least at short times), we conclude that J * E is not a conserved quantity in the interacting case. Finally, within our numerical ac- curacy and as an additional consistency check, we find that J * E /δE = α in the stationary state as expected from the discussion in Sec. II C.
To summarize, J * E (t) =const whenever δσ 2 E ∼ t 2 but J * E is very sensitive to the initial transient dynamics in the energy dynamics and becomes constant after a time ≈ bJ. Furthermore, our set-up serves to prepare each half of the system in a state with a finite global energy current J E L/2 that, after some transient dynamics, does not decay since the global energy current operator is a conserved quantity.

V. COUPLED SPIN AND ENERGY DYNAMICS
After focussing on the energy dynamics in the absence of spin-/particle currents we now revisit the case of spin dynamics starting from states with S z i (t = 0) = 0. Thus, during the time evolution, the local spin and energy currents are both non-zero. In Ref. 20, the dynamics of the magnetization was studied, where the inhomogeneous spin density profile was induced by a Gaussian magnetic field in the initial state. We take the initial state to be the ground state of Eq. (3) in the sector with zero global magnetization, i.e., S z tot = i S z i = 0. Such a perturbation naturally also results in an inhomogeneous energy density in the initial state, which is coupled to the spin dynamics during the time evolution. 28

A. Massless phase
In Fig. 11(a) we compare the initial magnetization (black solid line) and the local bond energies (dashed red line) induced by a Gaussian magnetic field with B 0 = J and σ 0 = 5 at ∆ = 0.5 finding qualitatively the same pattern: Both the spin and the energy density follow the shape of the magnetic field, resulting in a smooth perturbation with small oscillations in the background away from the wave packet.
For the time evolution of the bond energies at 0 < ∆ ≤ 1, we perform an analysis of their spatial variance analo- gous to the discussion of the J i quench, finding ballistic dynamics in the massless phase. Since with a B 0 quench, initial states with very small δE can be produced, we next connect our numerical results to the predictions of LL theory, valid in the limit δE ≪ J (compare Sec. III).
Since we enforce zero global magnetization, we draw magnetization from the background into the peak. 14 Therefore, one has to carefully estimate the contributions to δE that do not contribute to the time-dependence of bond energies yet change the background density n bg . The latter, in turn, affects the expected group velocity and we thus expect to recover the LL result derived for the half-filled case, i.e., propagation with v g from Eq. (24), in the limit of large systems where n bg → 1/2. Furthermore, B 0 quenches induce 2k F -oscillations in the spin and energy-density. 20 To account for this we use coarse graining, i.e., averaging the energy density over neighboring sites, and we take the sum only over the area of the peak when estimating δE: We obtain where h i 0 denotes the ground state expectation value. From this quantity we calculate the velocity via V 2 E → V 2 E · δE/δE peak , which is shown in Fig. 11(b). Note that while δE peak is the correct normalization to obtain the correct velocities, we label our initial states via δE. At ∆ = 0 (blue circles), V 2 E decreases linearly as a function of δE. Next we compare the result from the low-energy theory from Sec. III (solid symbols at δE = 0) to our tDMRG data. For both ∆ = 0 and ∆ = 0.5, V 2 E for L = 200 sites is approximately 5% smaller than v 2 g from Eq. (24), which is mainly due to the deviation of the background density from half filling. While it is hard to get results for larger systems than L ∼ 200 in the interacting case, we can solve the ∆ = 0 case numerically exactly in terms of free spinless fermions, allowing us to go to sufficiently large L to observe V 2 E (L) → v 2 g as L → ∞. The inset of Fig. 11(b) shows the finite-size scaling of V 2 E (L) for ∆ = 0 using L = 100, 200, ..., 800 which yields V 2 E → v 2 g in the limit L → ∞, taking first δE → 0 for each system size. We thus, in principle, have numerical access to the dynamics in the low energy limit well described by Luttinger Liquid theory using a B 0 quench.  16), were found in the massive phase, which were interpreted as an indication of diffusive dynamics. We now demonstrate that while the spin dynamics may behave diffu- sively, i.e., δσ 2 S ∼ t over a certain time window, the energy dynamics in the same quench is still ballistic, i.e., δσ 2 E ∼ t 2 .
In Fig. 12, we show the full time evolution of the bond energies for a Gaussian magnetic field with B 0 = 1.5J and σ 0 = 5 on a chain of L = 200 sites at ∆ = 1.5. It consists of two rays propagating with opposite velocities. In Fig. 13, we compare the spatial variance of the magnetization σ 2 S (t) to the one of the bond energies σ 2 E (t) calculated in the same time evolved state. The main panel of Fig. 13 shows σ 2 E (t) − σ 2 E (0) which is very well described by a power-law fit with an exponent β = 1.98 on the accessible time scales. The inset of Fig. 13 displays the data for δσ 2 S (t) = σ 2 S (t) − σ 2 S (0) taken from Ref. 20. The spatial variance of the energy density is quadratic in time, even at times t 12/J where the spatial variance of the magnetization increases only linearly. This example reflects the qualitative difference between spin and energy transport in the massive phase of the XXZ model at zero global magnetization: The conservation of the global energy current is consistent with the observation of ballistically propagating energy wave-packets while spin clearly does not propagate ballistically. Our result, obtained in the non-equilibrium case with a zerotemperature background density, is consistent with the picture established from both linear-response theory 51,55 and steady-state simulations. 28,61,64 Very recently, Jesenko andŽnidarič have also studied the time evolution of spin and energy densities induced by a B 0 quench. 28 They concentrate their analysis on the velocity of the fastest wave-fronts, contrasting energy against spin dynamics. Based on the presence of these rays of fast propagating particles, they claim that the wave-packet dynamics still has ballistic features. However, their analysis neglects the influence of slower excitations that also contribute to the dynamics of the wave packet, which is captured by the variance, and it ignores the decay of the intensity in the outer rays that we typically observe whenever δσ 2 S ∼ t. 20 The latter is, if at all, weak in a ballistic expansion characterized by δσ 2 S ∼ t 2 . Therefore, while the analysis of Ref. 28 unveils interesting details of the time evolution of densities during a B 0 quench, we maintain that the variance is a useful quantity to identify candidate parameter sets for spin diffusion in, e.g., the non-equilibrium regime. Final proof of diffusive behavior then needs to be established by either demonstrating the validity of the diffusion equation or by computing correlation functions, see, e.g., Ref. 55,61. For instance, in Ref. 28, Jesenko andŽnidarič analyze the steady-state currents in the ∆ > 1 regime at finite temperature and obtain diffusive behavior.

VI. SUMMARY
We studied the real-time energy dynamics in XXZ spin-1/2 chains at zero temperature in two different scenarios. First, we investigated the energy dynamics in the absence of spin currents induced by a local sign change in the exchange interactions. The spatial variance behaves as δσ 2 E (t) ∝ t 2 for all ∆, consistent with ballistic dynamics. In the gapless regime, the velocity of the fastest excitation present in the dynamics is the group velocity v g of spinons, yet our particular quench also involves excitations with much smaller velocities resulting in expansion velocities V E < v g . Furthermore, the ballistic dynamics can be related to properties of energy currents. While the total current vanishes in our set-up, i.e., J E := i j E i = 0, the current in each half of the chain J E L/2 > 0 takes a constant value, after some transient dynamics. Therefore, in each half of the system, we prepared a state with a conserved global current, allowing us to make a direct connection to the predictions of linear-response theory where the existence of ballistic dynamics is directly linked to conservation laws that prohibit currents from decaying. 43 Moreover, we identified an observable J * E built from local currents whose expectation value J * E (t) is time independent if δσ 2 E ∝ t 2 and vice versa. This carries over to other types of transport as well and, in fact, the analysis of the time-dependence of J * E can be used as an independent means to identify ballistic regimes, or to unveil the absence thereof.
In the second part, we studied the energy dynamics induced by quenching a Gaussian magnetic field, with two main results. These quenches allow us to access the regime of weakly perturbed initial states and in that limit, we recover the predictions from Luttinger Liquid theory for the wave-packet dynamics: Their variance simply grows as δσ 2 E = v 2 g t 2 .
In the massive phase a very interesting phenomenon occurs, since the energy dynamics is ballistic on time scales over which the spin dynamics behaves diffusively although both are driven by the same perturbation. This resembles the picture established from linear-response theory, 43 there applied to the finite-temperature case, in the non-equilibrium setup studied here. While our numerical results cover spin chains on real space lattices and initial states far from equilibrium, the extension of our work to a finite temperature of the background will be crucial to tackle the most important open questions.
of the reduced density matrix ρ A = T r B ρ, where ρ = |ψ(t) ψ(t) and |ψ(t) is the time-evolved wave function, for a bipartition in which we cut the chain into two halves of length L/2 across the central link. Our results are plotted in Fig. 14. We observe that the von-Neumann entropy grows at most logarithmically (purple dashed line), in agreement with Ref. 21. The overall largest values of S vN (t) are found at the critical point ∆ = 1 (red squares). This behavior is very similar to the observations made in Ref. 19 for spin dynamics starting from a state with all spins pointing up(down) in the left(right) half.