Time Evolution within a Comoving Window: Scaling of signal fronts and magnetization plateaus after a local quench in quantum spin chains

We present a modification of Matrix Product State time evolution to simulate the propagation of signal fronts on infinite one-dimensional systems. We restrict the calculation to a window moving along with a signal, which by the Lieb-Robinson bound is contained within a light cone. Signal fronts can be studied unperturbed and with high precision for much longer times than on finite systems. Entanglement inside the window is naturally small, greatly lowering computational effort. We investigate the time evolution of the transverse field Ising (TFI) model and of the S=1/2 XXZ antiferromagnet in their symmetry broken phases after several different local quantum quenches. In both models, we observe distinct magnetization plateaus at the signal front for very large times, resembling those previously observed for the particle density of tight binding (TB) fermions. We show that the normalized difference to the magnetization of the ground state exhibits similar scaling behaviour as the density of TB fermions. In the XXZ model there is an additional internal structure of the signal front due to pairing, and wider plateaus with tight binding scaling exponents for the normalized excess magnetization. We also observe parameter dependent interaction effects between individual plateaus, resulting in a slight spatial compression of the plateau widths. In the TFI model, we additionally find that for an initial Jordan-Wigner domain wall state, the complete time evolution of the normalized excess longitudinal magnetization agrees exactly with the particle density of TB fermions.

If an initial state for such a study is prepared within a finite system, boundary effects such as Friedel oscillations interfere with a passing signal. System boundaries also limit the time span for signal tracing before non-trivial reflections occur at the boundaries. The maximum time is even more severely restricted by entanglement which develops across the system and which requires a computational effort that can drastically increase with time [45,46]. This has greatly hampered the analysis of large time asymptotic behaviour [19,28,34]. Boundary effects do not appear in infinite systems, for which the ground state and its time evolution can be efficiently calculated with MPS methods [6,[47][48][49]. However, these methods require complete translation invariance and can therefore not be applied to studying signal propagation.
In this paper we present a simple method to simulate the propagation of local signals on an infinite chain using MPS time evolution, without any finite size effects distorting the signal front. For related approaches to boundary effects, see [50][51][52] 5 . We study the time evolution of the transverse field Ising (TFI) model and of the spin-1/2 XXZ chain after local quantum quenches up to large times, which were not accessible before using conventional MPS techniques. In both models we observe distinct magnetisation plateaus developing over time close to the signal front similar to the case of TB fermions [41,42,44], and which also exhibit similar asymptotic scaling. Surprisingly we find an exact agreement at all times and positions between the magnetisation in the TFI model and the density of TB fermions for a particular type of signal. For the XXZ chain we observe interaction effects between individual plateaus, which can be tuned via the model parameters.
For our method we consider a spin chain of infinite size with nearest neighbour interactions, initially prepared in a state, such as the ground state, which is translation invariant for sites n > n 0 to the right of some site n 0 . At time zero, the system is excited by a quantum quench like one or more spin flips at sites ⩽ n n 0 or a modification of the Hamiltonian at ⩽ n n 0 . For local interactions it is known from the Lieb-Robinson bound [53,54] that wave fronts generated by such quenches can at most propagate with a characteristic maximum velocity v max , i.e. within a 'light cone' even in a non-relativistic system, as recently also seen experimentally [1]. Any correlations beyond the light cone are exponentially suppressed. In the following we will consider right moving signals for the sake of concreteness.

Method
Our approach is to introduce a division of the system into three parts, namely a comoving window (CMW), which moves towards the right with the wave front, and two half-infinite parts, a uniform one in front (i.e. to the right) of the window, and an arbitrary one to the rear. The window is chosen wide enough to contain the complete signal front, including the exponentially damped part to the right of the main front, to high precision. The signal therefore does not affect the uniform system to the right of the window. Likewise, when the window moves with v max , modifications in the rear part do not affect the CMW and need not be calculated. The method is therefore fit for studying fronts of propagating signals, in particular those generated by local quenches. Since bipartite entanglement [14,26,55,56] spreads at most with v max , the bipartite entanglement entropy is significantly lower around the wave front than in the bulk, allowing for reduced computational effort when using the CMW.
We mark the left and right boundary of the CMW with indices ℓ and r, respectively, and divide the system into left part ⩽ ℓ j , CMW ℓ ⩽ ⩽ + j r 1 , and right part ⩾ + j r 1 Low energy states of the overall system are well approximated by Matrix Product States (MPS) [4,5] and we write the wave function as an MPS in the so-called mixed canonical form as are left-and right-orthogonal matrices, respectively, defined inside the CMW, and ℓ⩽ ⩽ λ k r are diagonal matrices containing the Schmidt values of a bipartition at bond (k, k + 1). For a finite system, the left and right ends of (2) are terminated by contractions with boundary vectors; we however consider the infinite size limit. For a graphical representation of this MPS see figure 1.
The matrices R sj describe the uniform half-infinite system in the front and are therefore constrained to be translation invariant. We use a 2-site unit cell, i.e. = + R R s s j j 2 . The matrices A sj and B sj describe the CMW and are site dependent. For the matrices L sj , which describe the left part, we impose no uniformity restrictions. They represent initial conditions for the left boundary of the CMW and remain unchanged throughout the simulation. Additional matrices are added to this collection of L sj whenever the CMW is moved.
Let us consider one step of unitary time evolution for the entire system. Inside the CMW, between sites ℓ + 1 and r, we employ the time-dependent Density Matrix Renormalisation Group (tDMRG [8]), using a second-order even-odd Suzuki-Trotter decomposition [57] with local operators u e j j h , 1 i j j , 1 ( )τ = τ + − + and small time steps τ. In order to connect time evolution inside and outside the CMW we introduce two different approaches, which we now sketch for the case of the right (front) and the left (rear) boundary, respectively. Details can be found in appendix A.
In Method I (Uniform Update), applied to the right boundary, the matrices R sj of the right part are first updated by infinite system time evolving block decimation (iTEBD 5 The original preprints of [50][51][52] (arXiv:1207.0652, 1207.0678, 1207.0691) and of the present work (arXiv:1207.0862) appeared at the same time. [47]). We then evolve the junction bond (r, r + 1) by applying ˆ+ u r r , 1 and we exploit the right-orthogonality of + R sr 1 to update B sr and to ensure the gauge consistency of the MPS matrices around the junction bond.
For 1, such that all changes in the left part are compressed into the boundary matrix ℓ+ A s 1 , and the matrices ⩽ℓ L sj remain unchanged. The Uniform Update has some immediate advantages. It is easier to implement and it is also applicable in case of a timedependent H R . It does, however, require translation invariance of the right part. The Renormalised Update does not preserve the structure of the Suzuki-Trotter decomposition at the boundary and therefore continually introduces small perturbations there. In appendix D we compare both methods to analytical results and to a reference simulation on a very large stationary lattice and show that both methods work well. As errors in our new Uniform Update, when applied to the right boundary, are only of order ( ) − O 10 8 and thus smaller by several orders of magnitude than for the renormalised update, we use the uniform update for the right boundary.
For the left boundary, the simplest approach is to disconnect the left part by setting h 0 , 1ℓ ℓ = + , which already works quite well (see appendix D) when the window moves with v max , as then any perturbations are confined to the neighbourhood of the rear boundary. Since the perturbations there are, however, smallest with the Renormalised Update, we use this method for the left boundary in the present paper For further details on the boundary updates and how to move the CMW along with a propagating signal see appendix A.

Transverse field Ising (TFI) model
The spin-1/2 TFI model [9][10][11][12][13][14][15][16][17][18][24][25][26][27] can be solved exactly [59,60] (see also appendix B), and the time evolution of local observables can in principle be calculated [10,11]. For the longitudinal magnetisation S x (n, t) (order parameter), analytical calculations are, however, difficult and some results have become available in the literature only recently [10,12], but to our knowledge not for local quenches on infinite systems. In the ferromagnetic phase h < h c = 0.5 the ground state is twofold degenerate and there is long-range order in S x . We prepare the system in the maximally symmetry broken ground state ⟩ |⇓ (appendix A.1) with ⟨ˆ⟩ = < S S : 0 x n x GS using iDMRG [6,48] and study the time evolution of several initial states excited from ⟩ |⇓ . In figure 2 we show the results for a Jordan-Wigner (JW) excitation on site n 0 inside the window, where † c c , are JW fermion operators (see [61] and appendix B). This corresponds to a spin flip in the z-direction at site n 0 and a domain wall in the x-direction between sites n 0 − 1 and n 0 . Window movement is triggered by bipartite entanglement entropy, resulting in window velocities consistent with exact maximum velocities (appendix B). We use a second order Suzuki-Trotter decomposition with a step size of τ = 0.002 and maximum matrix dimension = m 120 max during time evolution. The time evolution inside the CMW (figure 2) shows that boundary effects are indeed removed at both ends of the CMW. In appendix D we show that results inside the CMW are unperturbed to very high accuracy (about 10 −8 ) at all times.
When the window is not moved (figure 2, inset), the signal is absorbed by both boundaries temporarily, but reflections emerge eventually with both methods. This remains true also for additional models studied in appendix F, in all cases. We also investigate a pure domain wall (DW) excitation (ˆ) ⟩ ∏ |⇓ < S 2 n n n z 0 between sites n 0 − 1 and n 0 and a spin flip in Step structure. Despite different global shapes (see appendix E) for the different excitations, we find that a step structure always develops in S x (n, t) at the signal front at large times (figure 3), similar to the time evolution from an initial DW state for TB fermions [41,42,44]. The step structure takes much longer to develop for FlipX and DW excitations than for the JW case. The transverse magnetisation S z (n, t) does not show such a step structure.
The step structure is expected to be related to the ballistic nature of propagation at the signal front [36,37,44], like for TB fermions, where the steps are now fully understood as individual propagating particles [44]. For the TFI model, in different quench scenarios where two initially separate chains are joined, the beginnings of steps were previously visible in the results of [27], but were not investigated further. We are not aware of other occurrences for the symmetry broken phase. In the paramagnetic phase at large 2h = 10, TB-like scaling was observed in [25] for the transverse magnetisation S z (n, t) after joining two initially separate chains at different temperatures. No steps occurred for the longitudinal magnetisation. Due to their quantum origin these steps appear not to be accessible [44,62] by semi-classical approaches such as in [13].
We find that the proper quantity to analyse our results is the normalised excess longitudinal magnetisation Figure 3 shows that at large times this quantity indeed obeys the same scaling behaviour as the particle density of TB fermions [44] at the signal front. For the DW and FlipX cases, there is an additional proportionality factor ≠ C 1. The  and H(y) are the density and entropy scaling functions for TB fermions [44]. The lines are successively offset by 0.25 in the vertical direction. asymptotic scaling function G(y) for TB fermions [44] is approached from different directions for different excitations. For DW and FlipX excitations, the exponent α with the best data collapse depends on h, whereas for the JW case it is independent of h. 3 where v = h is the TFI signal velocity (appendix B.2) and ( ) N n vt , TB is the particle density of TB fermions at time vt after a DW excitation (steplike initial density as in [44]). We find this identity to hold up to the numerical precision of our data  for all sites n and times t for < h h c , i.e. in the ferromagnetic phase, but for the longitudinal magnetisation only.
The steps in ( ) N n t , TB have been shown to correspond to individual propagating particles [42,44] and we note that in the case of the TFI model a similar interpretation in terms of individual quasi-particles can only be given to the scaled excess longitudinal magnetisation M(n, t) after a JW excitation in the ferromagnetic phase. Due to the twofold degeneracy of the ground state in this phase the application of a local perturbation in the fermion picture generates a topologically non-trivial excitation by creating a domain wall (plus spin flip) in the spin picture, which then decays like a domain wall of TB fermions with time scale vt. In the paramagnetic phase the same excitation would create a local excitation also in the spin picture, i.e. no domain wall.
Other observables, however, are different between the TFI model and TB fermions. The transverse magnetisation ⟨ˆ⟩ S z is finite in the TFI model (see appendix B) while the corresponding quantity ⟨ ⟩ † + c c vanishes for TB fermions. The bipartite entanglement ( ) S n t , ent in the TFI model also develops a step structure, but it is at all times smaller than for TB fermions (see appendix E) and it exhibits different scaling behaviour (see figure 3). This fact only becomes fully apparent at large enough times, which our approach can provide. It would be interesting if the above identity between TB fermions and the TFI model could be understood in more detail analytically.

XXZ model
Inspired by the above observations in the TFI model in the symmetry broken ferromagnetic phase, we also investigate the XXZ antiferromagnet [19-23, 28-38, 39], in the gapped symmetry broken phase for several ∆ < −1, where the ground state is also twofold degenerate. We prepare the system in the maximally symmetry broken ground state ⟩ |⇓ with staggered magnetisation ˜( ) ⟨ˆ⟩ = − < S S 1 0 z n n z GS using iDMRG and again study the evolution of a JW excitation Notice that due to = S 0 x GS a JW excitation is locally indistinguishable from a simple domain wall according to the magnetisation and that the roles of x and z are interchanged with respect to the TFI results. Additionally, we also study a spin flip in the z-direction at site n 0 (FlipZ). Window movement is triggered by bipartite entanglement entropy, resulting in window velocities consistent with exact results (see appendix C). During the time evolution we use a second order Suzuki-Trotter decomposition with a step size of τ = 0.01 and maximum matrix dimensions of The signal front again develops a step structure. To our knowledge this had not been realised before our study, however it was recently confirmed [36,39] after the preprint version of our study, but not further investigated. We also observe a pairing effect between neighbouring spins, leading to an additional internal step structure, which stems from the spinon like nature of elementary excitations created by the quench (figure 4 inset). Due to the dynamics generated by (7), elementary spinons can only hop by two lattice sites at a time.
We find that at very large times, which are virtually impossible to access with conventional MPS techniques [28,34], the staggered normalised excess magnetisation at the signal front shows the same scaling behaviour as TB fermions, albeit with an additional horizontal scaling constant a, which is parameter dependent and increases with |∆| (figure 5 and inset). We therefore again interpret magnetisation steps as due to individual propagating quasi-particles, which, however, show interaction effects by getting squeezed together more and more around the signal front with increasing |∆|. This behaviour can be explained by the fact that particles repel each other more with increasing interaction, but at the same time they are confined within the light cone dictated by the Lieb-Robinson bound.
Since the particle density is much lower around the signal front, more and more particles are pushed towards the signal front and get squeezed together there. Our data, however, suggest that this effect saturates around |∆| ≈ 5 (see inset of figure 5). It would be very interesting to understand these interaction effects between individual steps in more detail analytically. The asymptotic scaling function G( y) is approached differently for different ∆, but the scaling exponents appear to be independent of ∆ for all the quenches investigated. For M(n, t) they are equal to the TB case with value 1/3, whereas we again find a different effective exponent of ≈ 1/4 for the bipartite entanglement entropy (figure 5).

Conclusions
We have introduced an easy-to-implement method combining finite and infinite system MPS techniques that can follow the propagation of a signal front on an infinite spin chain unimpeded and free from finite size effects for very long simulation times and with very high precision, considerably improved over other approaches. We note that even when the window is not moved, local signals can be simulated on the background of an infinite system, without perturbations emanating from the boundary. In this scenario the signal can be temporarily absorbed by the boundary, although it is always reflected eventually.
Furthermore, the method is not restricted to the evolution of excitations under uniform Hamiltonians. For example, the AKLT model [63] with inhomogeneous bond interactions or 1D quantum systems under exponential or hyperbolic deformation [64,65] have uniform ground states, whereas the Hamiltonians are not uniform.
To simulate the time evolution of a signal front of width L propagating with velocity v up to some time t, our method requires a numerical effort of the order ( ) Lt O , whereas for the same calculation using standard finite size MPS techniques the numerical effort would scale as ( ) with an additional v-dependent factor which scales quadratically in simulation time. We want to emphasise that additionally, standard finite size MPS techniques would also suffer from finite size effects such as boundary effects or the absence of exact ground state degeneracies in symmetry broken phases.
We have found that for all local quenches investigated in the symmetry broken phases of the TFI and the XXZ model, distinct magnetisation plateaus develop at the emerging signal front at very large times, where the scaled excess magnetisations in both models show the same long time limit scaling behaviour as the particle density of TB fermions after an initial domain wall excitation. For TB fermions these plateaus have recently been understood as being due to individual propagating particles [44]. Because of their quantum origin these plateaus cannot be studied [44,62] by means of semiclassical approaches such as in [13].
Our method has enabled us to calculate the time evolution of the order parameters of both models around the signal fronts generated by local quenches and investigate their features, which to our knowledge are available neither analytically nor semi-classically. In all cases it is important to reach very large simulation times, which are easily accessible through our approach, in order to reach the proper scaling regimes.
In the XXZ model we have observed an additional internal step structure due to the spinon nature of the involved elementary excitations, as well as parameter-dependent interaction effects between individual plateaus in the form of increasing spatial compression of the plateau width close to the signal front. This effect appears to saturate for |∆| 1. For the TFI model we have additionally found a surprising exact agreement of the normalised excess longitudinal magnetisation after a JW excitation with the density of TB fermions after a domain wall excitation. This exact mapping, however, does not apply to other observables such as, e.g. bipartite entanglement.
It would be interesting to understand both the interaction effects between plateaus in the XXZ model and the exact agreement between the TFI model and TB fermions in more detail analytically.

Acknowledgments
We would like to thank Th Barthel, V Eisler, F Maislinger, M M Rams, D Schuricht, U Schollwöck, and F Verstraete for valuable discussions. This work was supported by the Austrian Science Fund (FWF): F4104 SFB ViCoM and by the EP-SRC under grant EP/I032487/1. TN acknowledges the support of Grant-in-Aid for Scientific Research (C) No. 22540388.

Appendix A. CMW time evolution and boundary update methods
In this appendix we illustrate one time evolution step for the entire system when following a right moving signal. We describe the procedure in the following order. We first evolve the part of the system contained within the CMW (appendix A.2) before updating the right part using Method I (appendix A.3) and updating the left part using the more involved Method II (appendices A.4 and A.5). Note that this is the setup used in the main text; however, in principle any of the two methods can be used at any boundary. A detailed assessment of different setups is given in appendix D. We also describe the process of moving the CMW along with a propagating signal (appendix A.6). A short sketch of both boundary update methods, illustrating their advantages and restrictions, along with a motivation of the above choice is given in the main text.

A.1. System initialisation
In the main text in particular we use a setup dividing the system into a semi-infinite, initially translation invariant left part, a finite-size CMW (inside of which a signal will be created) and a semi-infinite, at all times translation invariant right part. We initialise the system by first determining a uniform MPS representation of the respective model's ground state on an infinite chain using iDMRG [6,48]. We then set all MPS matrices inside the CMW (matrices σ A j and σ B j ), the semiinfinite right part (matrices σ R A and σ R B forming this part's two-site unit cell) and the semi-infinite left part ( all matrices σ L j ) to this uniform MPS ground state representation after appropriate (left or right) orthonormalisation [5,48], i.e. we initialise the entire system to be in the infinite system's translation invariant ground state. Subsequently, we locally excite the system out of its ground state to generate several different kinds of local signals by applying suitable operators to one or more MPS matrices inside the CMW.
For other purposes the generalisation to different initial conditions is straightforward.

A.2. Time evolution within the CMW (CMW update)
Without loss of generality we consider a CMW with an even number of sites and first-order, even-odd, Suzuki-Trotter decomposition [57] with local operators ˆ( )τ = τ and finite time steps τ. The generalisation to higher order Suzuki-Trotter decompositions and windows containing an odd number of sites is straightforward. All the simulations in this work were performed using second-order Suzuki-Trotter decomposition and windows with an even number of sites. For one time step inside the CMW we use tDMRG [8] and 1 and (r, r + 1) are thus defined to be even bonds (see figure A1). By choosing this order we preserve the structure of the Suzuki-Trotter decomposition of the CMW and the right part, when Method I is used to update the right boundary.
At this stage all the even and odd bonds have been updated, except for the junction bonds (ℓ ℓ ) + , 1 and (r, r + 1), i.e. the boundary matrices ℓ+ A s 1 and B sr are not yet fully updated. Note that an implementation of this update using time evolving block decimation (TEBD [7]) is equivalent. For a graphical representation see figure A1.

A.3. Method I (uniform update).
We use this easy to implement procedure for the right boundary. Due to the assumed translation invariance over a 2-site unit cell, this part can be described by two right-orthogonal matrices R A s j and R B s j , such that the wavefunction in MPS representation around the right boundary reads The evolution of the matrices R A s j and R B s j is performed by iTEBD (or variations thereof) using local operators ˆ( ) τ u A and ˆ( ) τ u B [47,66], where ˆ( ) τ u A acts on odd bonds and ˆ( ) τ u B acts on even bonds.
In a first step, we apply an odd bond iTEBD update in the right part to get where ° denotes matrices having received an odd bond update.
Here the decomposition of the result of the right-hand side of (A.2) is implicitly assumed. It can be done by an SVD either involving a division by Schmidt values following [47] or avoiding the division by Schmidt values by using the approach of [66]. The wavefunction at this point reads . In parallel we perform an even bond iTEBD update in the right part to get where again the decomposition of the result of the right side is implicitly assumed.
All the bonds have now been updated. Since there is negligible influence of the signal around the right boundary by construction, the state of the right part should be the same as for a time evolved uniform system without signal up to high precision, i.e. we can also assume Φ =    Note: For a graphical representation see figure A2.
The procedure can also be easily translated to the left boundary exploiting left orthogonality, where the translation invariance of the left-orthogonal matrices L sj is then required.
Method I is also applicable when H R is time dependent, e.g. in case of a global quench.

A.4. Method II (renormalised update)
We use this procedure for the left boundary. For this Method we follow a similar approach as introduced by Cazalilla and Marston [58] (Method II is similar to the algorithm introduced in [50], where preprints of [50] and of the present paper appeared at the same time), such that matrices L sj in the left part remain unchanged at all times during time evolution.
The effect of the left part is encoded in a renormalised formulation of ˆˆ◃ ℓ , 1 , which is exactly the renormalised Hamiltonian used in standard DMRG formulations (see e.g. [4,5]). All changes in the left part are then solely encoded in an update of the boundary matrix ℓ+ A s 1 . Note that for this method the matrices L sj in the left part need not be translation invariant.
Since matrices L sj are not changed during this update, we rewrite the wavefunction in MPS form after the CMW update in terms of the auxiliary basis states where ℓ a is the left index of matrix A s 1 ℓ°+ . The right-hand side of (A.7) is formally just the semi-infinite product of all matrices to the right of site ℓ. The overall state vector after the CMW update can thus also be written .    Notice also that ◃ ℓ+ U , 1 eff breaks the structure of the evenodd Suzuki-Trotter decomposition in the left part. This introduces an additional error, which is of the same order as the Suzuki-Trotter error and can in principle be made arbitrarily small by using higher order Suzuki-Trotter decompositions and smaller time steps τ at the cost of increased computational time. The effect of this additional error is investigated in detail in appendix D. It could be avoided by using the renormalised imaginary-time transfer matrix, as used in finite temperature DMRG [67], to update ℓ°+ A s 1 . For an algorithmic summary see table A2, for a graphical representation of this update see figure A4.    accumulates the renormalised Hamiltonian containing all sites ⩽ k j (see e.g. [69]).
To determine (A.17) we need a way to calculate the semiinfinite matrix product [ℓ] E . For the moment we consider the case where both Ĥ L and the matrices L sj are translation invariant. In this case F [ j] is also translation invariant and [ℓ] E can be calculated by, e.g. finding the dominant left eigenvector of F [ j] , as explained in [72].
However, here we follow an approximate but sufficiently accurate approach for calculating [ℓ] E , which is inspired by standard DMRG formulations. For this we relax the condition of semi-infinity for the left Hamiltonian ˆ◃ ℓ+ H , 1 and approximate it with a finite size Hamiltonian, which we increase in size until we get a converged result. The finite size version of ˆ◃ ℓ+ H , 1 in MPO form is thus contracted also on the left side by exploiting the left-orthogonality of the matrices L sj . For a graphical representation see figure A3(b E . In case of a translation invariant left part, its calculation is very similar to the renormalisation steps of an iDMRG simulation [6,48] (no eigenvalue/SVD steps). The number of renormalisation steps is dependent on the effective correlation length induced by the uniform MPS matrices L sj .
In practice, it takes about 75 renormalisation steps for the TFI model at h = 0.45 (m 0 = 30) and about 100 steps for the XXZ model at J z = −2 (m 0 = 88) for convergence in energy up to an accuracy of 10 −15 , where m 0 is the bond dimension of the ground state MPS representation. The overall computational effort here is comparable to a few time evolution steps within the CMW.

A.6. Window movement
We describe the window movement by a single site. For a shift by a 2-site unit cell, the same procedure as for a single site is applied twice.
If can grow with successive window shifts. An impinging signal can therefore be partly absorbed such that immediate perturbations are considerably suppressed (see also appendix F).
We trigger the window shift when the relative change of the bipartite entanglement entropy at some site sufficiently far away from the right boundary rises above a certain threshold. The margin between this site and the right boundary should be large in comparison to the correlation length of the initial state such that the exponentially suppressed correlations reaching beyond the Lieb-Robinson light cone [53] are negligible. For all simulations in the main text we use a margin of 24 sites and a threshold of 1%. If known beforehand, the window can also be moved directly with v max .
Here we have already taken the thermodynamic limit while considering periodic boundary conditions (a boundary term arising from the JW transformation and periodic boundary conditions is neglected as it is of the order ( ) L O 1/ where L is the system size).
A subsequent Bogoliubov transformation [73] to fermionic operators η k , † η k in momentum space then diagonalises the Hamiltonian. The coefficients a k and b k are real and satisfy and can be determined as The Hamiltonian then readŝ † and the ground state corresponds to the vacuum state ⟩ |0 in terms of the fermionic operators η k and † η k .

B.2. Signal velocity in the TFI model
The propagation of a signal induced on top of the ground state ⟩ |0 of the TFI model can be understood as the excitation and propagation of a superposition of non-interacting particles with momenta k and corresponding energies ε k created by † η k . In this picture, the maximum velocity v max of the signal can be exactly calculated as the maximum of the group velocity = − + +∆ can be solved, e.g. by means of the coordinate Bethe ansatz [74]. We seek solutions for the ground state and elementary excitations of the XXZ antiferromagnet with ∆ < −1 in the thermodynamic limit, which can be found, e.g. in [75].
In the thermodynamic limit the roots of the Bethe equations become dense and their distribution for the ground state is characterised by a density function g 0 (x), which for ∆ < −1 satisfies the following integral equation . The solution to this integral equation is given by is a Jacobian elliptic function [76], K(m) the complete elliptic integral of the first kind The root density g 0 (x) can then be used to calculate various quantities such as the ground state energy and elementary excitations.

C.2. Signal velocity in the XXZ model
To calculate the maximum signal velocity v max as a function of ∆ we first determine the dispersion relation ε k for the elementary excitations. As for the TFI model in appendix B.2 we then obtain v max as the maximum of the group velocity =  Note: These values are obtained from numerically finding the maximum of (C.12) with a numerical precision of 10 −8 .
The dispersion of the elementary excitations is given by [75] is the finite energy gap present in this phase and x 0 (k) has to be determined by inverting for a given momentum k. From this we can calculate the group velocity where we need Using some properties of Jacobian elliptic functions [76] and defining ( ) are also Jacobian elliptic functions. Defining we can then write , , .

Appendix D. Test of precision of results
To assess the accuracy of the CMW approach, we compare it with a reference system on a very large lattice and with exact results obtained in appendix B.3. We investigate simulations of a JW excitation on top of the infinite system ground state in the TFI model at h = 0.45 for windows of different sizes, different boundary update methods and different margins between signal and right boundary for triggering the window shift. Note that the correlation length of this system is ξ ≈ 4.36 sites. It can be obtained from the second largest eigenvalue in magnitude λ 2 of the MPS transfer matrix [4,77]. For all simulations we use secondorder Suzuki-Trotter decomposition with time step τ = 0.002 and maximum bond dimension = m 120 max . These are the same simulation parameters as used for the investigation of a JW excitation in the TFI model in the main text.
The reference simulation is also performed using the CMW algorithm, but starting with the translation invariant initial state inside a non-moving window of very large size of N = 1000 sites. This means that the window is never shifted. Boundary effects are removed by using Method I for both boundaries. For the reference simulation we perform time evolution up to t = 800, such that the signal induced in the centre of the system at t = 0 does not reach the boundaries. For a plot of the reference simulation see figure D1. There we show the transverse magnetisation S z (n, t) and the bipartite entanglement entropy ( ) S n t , ent . It can be seen that boundary effects are indeed removed for the non-moving window with Method I (otherwise disturbances would constantly radiate from the boundaries) and that the signal is still about ≈150 sites away from the boundaries at t = 800. We compare results from setups with different CMW sizes N and different numbers of margin sites between signal and right boundary (in sites, see appendix A.6), as well as different setups for using Method I and II for the updates at the left boundary (LB) and right boundary (RB). We find that the accuracy of the simulation strongly depends on the boundary update method used at the right boundary and the margin between the signal and right boundary, whereas the window size N has virtually no impact on the accuracy. For a selection of compared setups see table D1.
For comparison we will consider the transverse magnetisation, since only this quantity is available analytically. We display the absolute value of the difference in transverse magnetisation,   figure D2. For other observables, analytic results are not available, but we can compare to the reference simulation. We find that comparison of the magnetisation in x direction S x (n, t) and of the bipartite entanglement entropy ( ) S n t , ent to the reference simulation yield results that look very similar to figure D2 and the obtained absolute differences are also of the same orders of magnitude. We note in addition that comparison between left and right column in figure D2 confirms the absence of boundary effects in the reference simulation to high precision.
In the following we discuss the comparisons of the 3 cases listed in table D1.
Case (1)  In case (2), where Method II is used at the right boundary, differences inside the CMW rise up to ( ) − O 10 5 for both comparisons, i.e. they are considerably higher by about 3-4 orders of magnitude in comparison to case (1), where Method I is used. This can be explained by the fact that Method II breaks the structure of the Suzuki-Trotter decomposition at the boundary, which introduces additional perturbations. These perturbations can in principle be reduced by using higher order Suzuki-Trotter decompositions and smaller time steps and thus increasing computational effort, but they are always present. Method I, however, is completely devoid of this kind of perturbations. Also, the renormalised Hamiltonian ▹ H r eff necessary for Method II is only calculated up to a finite precision. We, however, find the perturbations to be largely independent of the precision used to calculate ▹ H r eff as described in appendix A.5. We conclude that using Method I at the right boundary yields results which are better by about 3-4 orders of magnitude in precision than using Method II when employing second-order Suzuki-Trotter decomposition with a time step of τ = 0.002. In case (3) the left part has been disconnected from the CMW altogether by setting h 0 , 1ℓ ℓ = + ('cut') as described in the main text. Also the margin between signal and right boundary is reduced to 3 sites. Due to the cut, perturbations around the left boundary are now considerably higher and go up to ( ) − O 10 2 both for the comparison to analytic results and the reference simulation. These perturbations, however, again remain confined around the left boundary at all times. Differences inside the CMW are now ( ) − O 10 6 for both comparisons. This can be explained by the fact that the margin of 3 sites is now smaller than the correlation length ξ ≈ 4.36 and the exponentially suppressed correlations reaching beyond the Lieb-Robinson light cone [53] induce perturbations at the right boundary.
In conclusion, both Method I (Uniform Update) and Method II (renormalised update) work quite well. Furthermore, the easy to implement Method I yields results with a precision of about 10 −8 , still better by several orders of magnitude than Method II when used at the right boundary. For the methods to work, the margin between the signal and right boundary needs to be considerably larger than the correlation length. At the left boundary the easiest approach, a simple cut, already works well when the very rear of the CMW is not of interest.
Overall we have shown that the error produced by the CMW approach, especially with Method I, is very small and remains virtually constant for large times during the simulation when the margin between the signal and right boundary is kept sufficiently larger than the correlation length in the initial state.

Appendix E. Unscaled time evolution results
In this section we show time evolution results before scaling for the TFI model and the XXZ model, for the signals investigated in the main text.

E.1. TFI model
In our simulations we use a Trotter step size of τ = 0.002 and a maximum matrix dimension = m 120 max . The unscaled

Note:
We compare CMW results on N = 120 sites with analytic results (anal.) and with results from a reference simulation on N = 1000 sites (ref.). 'Margin' specifies the number of sites kept between the signal and the right boundary of the CMW as explained in appendix A.6. The precision P ref./anal. is the resulting maximum absolute difference in transverse magnetisations (D.1) inside the CMW away from the left boundary, between the CMW simulation and the reference simulation or analytic result (black dashed lines in figure D2). All simulations were performed using = m 120 max and second-order Suzuki-Trotter decomposition with time step τ = 0.002 up to t = 800. Case (1) corresponds to the setup used for data analysis in the main text. For (3), 'cut' means that the CMW is disconnected from the left part by setting ˆℓ ℓ = , corresponding to the simplest to implement setup, as described in the main text. A comparison between cases (1) and (2) shows that Method I yields very precise results, better by several orders of magnitude than Method II. magnetisation S x (n, t) in the TFI model for the three different quenches employed is shown in figure E1 for times t = 0 and t = 90. The global shapes are quite different, while developing plateaus are visible for all three quench types at t = 90. It can also be seen that around the signal front, the magnetisation of a single spin flip is always larger than of a domain wall, which in turn is always greater than the magnetisation of a JW excitation. This fact is reflected in the different values for the constant C in figure 3 of the main text.
The unscaled S x (n, t) at h = 0.45 after a JW excitation in the infinite system ground state versus absolute position n at large times 500 < t < 1000 is shown in figure E2. The ballistic propagation of the signal front as well as magnetisation steps near the front are clearly visible. No such steps appear in the transverse magnetisation. A scaling behaviour of the magnitude and distance to the signal front of the steps can be conjectured. This scaling behaviour is discussed in detail in the main text.
Other observables and signals, such as single spin flip and domain wall excitations qualitatively show the same propagation, shape and step structure. Their scaling behaviour, however, varies in scaling exponents and quality with varying field strength h.
We also show the bipartite entanglement entropy    The scaling behaviour of the larger step structure is investigated in detail in the main text. The overall shape of the unscaled staggered magnetisation ˜( ) S n t , z looks similar to the shape of the longitudinal magnetisation S x (n, t) of the TFI model with a JW excitation as shown in figure E1. Different signals such as single spin flips yield similar results.

Appendix F. Boundary reflections
In this appendix we consider the case of signals impacting the boundaries of a non-moving window for several different models. We study the time evolution beyond the time where a signal reaches the boundaries, both with Method I and Method II. In all cases we observe reflections from the boundary after some time. The nature of these reflections generally depends on the boundary update method as well as the initial uniform state and the type of the signal.
The models and signals that have been studied in particular are the TFI model with a JW excitation and a single spin flip in the x-direction, the XXZ model with a JW excitation and a single spin flip in the z-direction, the S = 1 Heisenberg model with a spin up excitation (this particular case is also studied in [50] with a method similar to Method II, but only for shorter times), and the S = 1 AKLT model [63] with a spin up excitation. We observe reflections from the boundary after some time in all cases.
In the following we show results for the two cases of the TFI model with a JW excitation and the S = 1 AKLT model with a spin up excitation, where we have used Method II for the left boundary and Method I for the right boundary to see their respective behaviour.

F.1. TFI model with a JW excitation
We again consider the TFI model at h = 0.45 after a JW excitation in the infinite system ground state. We use a nonmoving window with N = 50 and maximum bond dimension  can be seen in figure F1. The signal reaches the boundaries at ≈ t 40 and reflections start to emerge at ≈ t 90. We compare the magnetisation S x (n, t) of this simulation with the magnetisation ( ) S n t , x ref of the reference simulation of appendix D and show their absolute difference ( ) ( ) ( ) ∆ = − M n t S n t S n t , , , figure F2, where subplot (a) shows ( ) ∆M n t , at the left and right boundaries of the N = 50 non-moving window (n = 1 and n = 50, respectively) versus time t and subplot (b) shows ( ) ∆M n t , versus position n inside the non-moving window at various times t.
In figure F2(a) it can be seen that initially the deviations at the right side (Method I) are much lower than at the left side (Method II) until ≈ t 50. The deviation at both boundaries then increases exponentially further until ≈ t 100, where it becomes of the order ( ) O 1 . We notice that the deviations for the right boundary are always a bit lower than for the left boundary. We conclude that for the investigated case Method I performs slightly better than Method II in absorbing a signal for a limited time.

F.2. AKLT model with spin up excitation
We also consider the S = 1 bilinear, biquadratic chain at the AKLT point [63] defined by the Hamiltonian The ground state is a valence bond state and has an exact MPS representation with bond dimension m 0 = 2 (see, e.g. [5]). We induce a signal on top of the infinite system ground state by applying the spin ladder operator + S n 0 . We use a non-moving window with N = 60 sites and maximum bond dimension  entropy ( ) S n t , ent and the magnetisation S z (n, t) can be seen in figure F3.
Here the signal impacting at ≈ t 35 is reflected almost immediately. This stems from the fact that the MPS matrices at the boundary sites have to absorb all the information about excited states contained within the propagating signal. Here these matrices, however, have bond dimension m 0 = 2, which is much too small for the matrices to absorb this information for a long time span.