Time evolution within a comoving window: scaling of signal fronts and magnetization plateaus after a local quench in quantum spin chains

The spin-1/2 TFI model [9–18, 24–27] on an infinite chain defined by

$$\begin{eqnarray}&&\hat{H}=-\underset{n}{\mathop \sum}\hat{S}_{n}^{x}\hat{S}_{n+1}^{x}-h\underset{n}{\mathop \sum}\hat{S}_{n}^{z}\end{eqnarray} \tag{ 3 }$$

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 $|\Downarrow \rangle$ (appendix A.1) with $S_{\text{GS}}^{x}:=\langle \hat{S}_{n}^{x}\rangle <0$ using iDMRG [6, 48] and study the time evolution of several initial states excited from $|\Downarrow \rangle$. In figure 2 we show the results for a Jordan–Wigner (JW) excitation

$$\begin{eqnarray}&&(c_{{{n}_{0}}}^{\dagger}+{{c}_{{{n}_{0}}}})|\Downarrow \rangle =\underset{n<{{n}_{0}}}{\mathop \prod}(-2\hat{S}_{n}^{z})(2\hat{S}_{{{n}_{0}}}^{x})|\Downarrow \rangle \end{eqnarray} \tag{ 4 }$$

on site n₀ inside the window, where ${{c}^{\dagger}},c$ are JW fermion operators (see [61] and appendix B). This corresponds to a spin flip in the z-direction at site n₀ and a domain wall in the x-direction between sites n₀  −  1 and n₀. 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 $\tau =0.002$ and maximum matrix dimension ${{m}_{\text{max}}}=120$ 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⁻⁸) at all times.

Zoom In Zoom Out Reset image size

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 $\prod\nolimits_{n<{{n}_{0}}}(2\hat{S}_{n}^{z})|\Downarrow \rangle$ between sites n₀  −  1 and n₀ and a spin flip in the x-direction (FlipX) $(2\hat{S}_{{{n}_{0}}}^{z})|\Downarrow \rangle$ at site n₀.

3.1.1. 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.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&M(n,t)\equiv \left[{{S}^{x}}(n,t)-S_{\text{GS}}^{x}\right]/|2S_{\text{GS}}^{x}|~.\end{eqnarray} \tag{ 5 }$$

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\ne 1$. The 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.1.2. Exact identity.

In fact, for the JW excitation we find a surprisingly much closer identity with TB fermions: the complete time evolution of the normalised excess longitudinal magnetisation obeys

$$\begin{eqnarray}&&M(n,t)={{N}_{\text{TB}}}(n,vt)\end{eqnarray} \tag{ 6 }$$

where v  =  h is the TFI signal velocity (appendix B.2) and ${{N}_{\text{TB}}}(n,vt)$ 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}_{\text{c}}}$, i.e. in the ferromagnetic phase, but for the longitudinal magnetisation only.

The steps in ${{N}_{\text{TB}}}(n,t)$ 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 $\langle {{\hat{S}}^{z}}\rangle$ is finite in the TFI model (see appendix B) while the corresponding quantity $\langle {{c}^{\dagger}}+c\rangle$ vanishes for TB fermions. The bipartite entanglement ${{S}_{\text{ent}}}(n,t)$ 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.

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],

$$\begin{eqnarray}&&\hat{H}=-\underset{n}{\mathop \sum}(\hat{S}_{n}^{x}\hat{S}_{n+1}^{x}+\hat{S}_{n}^{y}\hat{S}_{n+1}^{y}+\Delta \hat{S}_{n}^{z}\hat{S}_{n+1}^{z}\,),\end{eqnarray} \tag{ 7 }$$

in the gapped symmetry broken phase for several $\Delta <-1$, where the ground state is also twofold degenerate. We prepare the system in the maximally symmetry broken ground state $|\Downarrow \rangle$ with staggered magnetisation $\tilde{S}_{\text{GS}}^{z}={{\left(-1\right)}^{n}}\langle \hat{S}_{n}^{z}\rangle <0$ using iDMRG and again study the evolution of a JW excitation

$$\begin{eqnarray}&&(c_{{{n}_{0}}}^{\dagger}+{{c}_{{{n}_{0}}}})|\Downarrow \rangle =\underset{n<{{n}_{0}}}{\mathop \prod}(-2\hat{S}_{n}^{x})(2\hat{S}_{{{n}_{0}}}^{z})|\Downarrow \rangle \end{eqnarray} \tag{ 8 }$$

at site n₀ inside the window (figure 4).

Zoom In Zoom Out Reset image size

Notice that due to $S_{\text{GS}}^{x}=0$ 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₀ (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 $\tau =0.01$ and maximum matrix dimensions of ${{m}_{\text{max}}}=150$ for $\Delta =-4$, ${{m}_{\text{max}}}=160$ for $\Delta =-3$ and ${{m}_{\text{max}}}=180$ for $\Delta =-2$ with discarded weights of at most $\text{O}({{10}^{-8}})$.

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

$$\begin{eqnarray}&&M(n,t):=[{{\tilde{S}}^{z}}(n,t)-\tilde{S}_{\text{GS}}^{z}]/|2\tilde{S}_{\text{GS}}^{z}|\end{eqnarray} \tag{ 9 }$$

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 $|\Delta |$ (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 $|\Delta |$. 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 $|\Delta |\approx 5$ (see inset of figure 5). It would be very interesting to understand these interaction effects between individual steps in more detail analytically.

Zoom In Zoom Out Reset image size

The asymptotic scaling function G(y) is approached differently for different $\Delta$, but the scaling exponents appear to be independent of $\Delta$ 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 $\approx 1/4$ for the bipartite entanglement entropy (figure 5).

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}^{{{\sigma}_{j}}}}$ and ${{B}^{{{\sigma}_{j}}}}$), the semi-infinite right part (matrices $R_{A}^{\sigma}$ and $R_{B}^{\sigma}$ forming this part's two-site unit cell) and the semi-infinite left part (all matrices ${{L}^{{{\sigma}_{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.

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 ${{\hat{u}}_{j,j+1}}(\tau )={{\text{e}}^{-\text{i}\tau {{{\hat{h}}}_{j,j+1}}}}$ 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 apply ${{\hat{u}}_{j,j+1}}\left(\tau \right)={{\text{e}}^{-\text{i}\tau {{{\hat{h}}}_{j,j+1}}}}$ to the bonds from $\left(\ell +1,\ell +2\right)$ until (r  −  1, r) and update matrices ${{A}^{{{s}_{j}}}}$ and ${{B}^{{{s}_{j}}}}$ contained within the CMW. The junction bonds $\left(\ell ,\ell +1\right)$ and (r, r  +  1) at the left and right boundary of the CMW are updated separately.

We first update all odd bonds $\left\{\ldots ,\left(r-3,r-2\right),\left(r-1,r\right)\right\}$ and then all even bonds $\left\{\ldots ,\left(r-4,r-3\right),\left(r-2,r-1\right)\right\}$. The junction bonds $\left(\ell ,\ell +1\right)$ 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.

Zoom In Zoom Out Reset image size

At this stage all the even and odd bonds have been updated, except for the junction bonds $\left(\ell ,\ell +1\right)$ and (r, r  +  1), i.e. the boundary matrices ${{A}^{{{s}_{\ell +1}}}}$ and ${{B}^{{{s}_{r}}}}$ 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.

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

$$\begin{eqnarray}&&\ldots {{B}^{{{s}_{r-1}}}}{{B}^{{{s}_{r}}}}R_{A}^{{{s}_{r+1}}}R_{B}^{{{s}_{r+2}}}R_{A}^{{{s}_{r+3}}}R_{B}^{{{s}_{r+4}}}\ldots .\end{eqnarray} \tag{ A.1 }$$

The evolution of the matrices $R_{A}^{{{s}_{j}}}$ and $R_{B}^{{{s}_{j}}}$ is performed by iTEBD (or variations thereof) using local operators ${{\hat{u}}^{A}}\left(\tau \right)$ and ${{\hat{u}}^{B}}\left(\tau \right)$ [47, 66], where ${{\hat{u}}^{A}}\left(\tau \right)$ acts on odd bonds and ${{\hat{u}}^{B}}\left(\tau \right)$ acts on even bonds.

In a first step, we apply an odd bond iTEBD update in the right part to get

$$\begin{eqnarray}&&R_{A\circ}^{{{s}_{A}}}R_{B\circ}^{{{s}_{B}}}=\underset{s_{A}^{\prime}s_{B}^{\prime}}{\mathop \sum}u\left(\tau \right)_{\left({{s}_{A}}{{s}_{B}}\right)\left(s_{A}^{\prime}s_{B}^{\prime}\right)}^{A}R_{A}^{s_{A}^{\prime}}R_{B}^{s_{B}^{\prime}},\end{eqnarray} \tag{ A.2 }$$

where $\circ$ 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

$$\begin{eqnarray}&&\ldots B_{\bullet}^{{{s}_{r-2}}}B_{\bullet}^{{{s}_{r-1}}}B_{\circ}^{{{s}_{r}}}R_{A\circ}^{{{s}_{r+1}}}R_{B\circ}^{{{s}_{r+2}}}\ldots ,\end{eqnarray} \tag{ A.3 }$$

where • denotes matrices having received both odd and even updates.

Special attention has to be paid to the operation of ${{\hat{u}}_{r,r+1}}\left(\tau \right)$ at the junction bond in order to update $B_{\circ}^{{{s}_{r}}}$. For this we form $\Phi_{\circ \circ}^{{{s}_{r}}{{s}_{r+1}}}:=B_{\circ}^{{{s}_{r}}}R_{A\circ}^{{{s}_{r+1}}}$ and act with ${{\hat{u}}_{r,r+1}}$ to get $\Phi_{\bullet \bullet}^{{{s}_{r}}{{s}_{r+1}}}$. In parallel we perform an even bond iTEBD update in the right part to get

$$\begin{eqnarray}&&R_{B\bullet}^{{{s}_{B}}}R_{A\bullet}^{{{s}_{A}}}=\underset{s_{A}^{\prime}s_{B}^{\prime}}{\mathop \sum}u\left(\tau \right)_{\left({{s}_{B}}{{s}_{A}}\right)\left(s_{B}^{\prime}s_{A}^{\prime}\right)}^{B}R_{B\circ}^{s_{B}^{\prime}}R_{A\circ}^{s_{A}^{\prime}}.\end{eqnarray} \tag{ A.4 }$$

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 $\Phi_{\bullet \bullet}^{{{s}_{r}}{{s}_{r+1}}}=B_{\bullet}^{{{s}_{r}}}R_{A\bullet}^{{{s}_{r+1}}}$, where both $B_{\bullet}^{{{s}_{r}}}$ and $R_{A\bullet}^{{{s}_{r+1}}}$ are right-orthogonal and $R_{A\bullet}^{{{s}_{r+1}}}$ is obtained from (A.4). We extract $B_{\bullet}^{{{s}_{r}}}$ from $\Phi_{\bullet \bullet}^{{{s}_{r}}{{s}_{r+1}}}$ by exploiting the right-orthogonality of $R_{A\bullet}^{{{s}_{r+1}}}$:

$$\begin{eqnarray}&&B_{\bullet}^{{{s}_{r}}}=\underset{{{s}_{r+1}}}{\mathop \sum}\Phi_{\bullet \bullet}^{{{s}_{r}}{{s}_{r+1}}}R_{A\bullet}^{{{s}_{r+1}}\dagger}.\end{eqnarray} \tag{ A.5 }$$

The wavefunction in MPS form is now completely updated around the right boundary and reads

$$\begin{eqnarray}&&\ldots B_{\bullet}^{{{s}_{r-2}}}B_{\bullet}^{{{s}_{r-1}}}B_{\bullet}^{{{s}_{r}}}R_{A\bullet}^{{{s}_{r+1}}}R_{B\bullet}^{{{s}_{r+2}}}\ldots .\end{eqnarray} \tag{ A.6 }$$

We note that in general, the decomposition $\Phi_{\bullet \bullet}^{{{s}_{r}}{{s}_{r+1}}}=\tilde{\mathop{B}}_{\bullet}^{{{s}_{r}}}\tilde{\mathop{R}}_{\bullet}^{{{s}_{r+1}}}$ into right-orthogonal matrices is not unique, but involves a gauge freedom $\tilde{\mathop{B}}_{\bullet}^{{{s}_{r}}}\tilde{\mathop{R}}_{\bullet}^{{{s}_{r+1}}}=B_{\bullet}^{{{s}_{r}}}{{x}^{-1}}xR_{A\bullet}^{{{s}_{r+1}}}$ with x a unitary matrix (Exploiting the right orthogonality of both $\tilde{R}_{\bullet}^{{{s}_{r+1}}}$ and $R_{A\bullet}^{{{s}_{r+1}}}$ we have $\mathbb{1}=\sum\nolimits_{s}\tilde{\mathop{R}}_{\bullet}^{{{s}_{r+1}}}\tilde{\mathop{R}}_{\bullet}^{{{s}_{r+1}}\dagger}=x\sum\nolimits_{s}R_{A\bullet}^{{{s}_{r+1}}}R_{A\bullet}^{{{s}_{r+1}}\dagger}{{x}^{\dagger}}=x{{x}^{\dagger}}$. Since x is square this also means ${{x}^{\dagger}}x=\mathbb{1}$ and thus x is unitary).

If the decomposition $\Phi_{\bullet \bullet}^{{{s}_{r}}{{s}_{r+1}}}=\tilde{\mathop{B}}_{\bullet}^{{{s}_{r}}}\tilde{\mathop{R}}_{\bullet}^{{{s}_{r+1}}}$ was carried out in the standard TEBD/tDMRG way (i.e. by means of an SVD), then a different gauge $\tilde{\mathop{R}}_{\bullet}^{{{s}_{r+1}}}\ne R_{A\bullet}^{{{s}_{r+1}}}$ and thus $\tilde{B}_{\bullet}^{{{s}_{r}}}\ne B_{\bullet}^{{{s}_{r}}}$ would result in general, since $\tilde{\mathop{R}}_{\bullet}^{{{s}_{r+1}}}$ was produced algorithmically in a different way than $R_{A\bullet}^{{{s}_{r+1}}}$. In that case, i.e. if $\tilde{\mathop{B}}_{\bullet}^{{{s}_{r}}}$ was used instead of $B_{\bullet}^{{{s}_{r}}}$, incompatible basis sets would meet at the junction bond, which would result in perturbations spreading from the boundary. By use of (A.5) we ensure that the correct gauge is chosen automatically.

This concludes one time step for the right part and right boundary. For an algorithmic summary see table A1, for a detailed graphical representation see figure A2.

Table A1. Algorithm for Method I updating the right boundary of the CMW.

(i)	Apply odd iTEBD update to get $R_{A\circ}^{{{s}_{A}}}$, $R_{B\circ}^{{{s}_{B}}}$.
(ii)	Use $R_{A\circ}^{{{s}_{A}}}$ to form $\Phi_{\circ \circ}^{{{s}_{r}}{{s}_{r+1}}}=B_{\circ}^{{{s}_{r}}}R_{A\circ}^{{{s}_{r+1}}}$
(iii)	Apply ${{\hat{u}}_{r,r+1}}$ to get $\Phi_{\bullet \bullet}^{{{s}_{r}}{{s}_{r+1}}}$.
(iV)	Apply even iTEBD update to get $R_{B\bullet}^{{{s}_{B}}}$, $R_{A\bullet}^{{{s}_{A}}}$.
(V)	Use $R_{A\bullet}^{{{s}_{A}}}$ to obtain $B_{\bullet}^{{{s}_{r}}}=\sum\limits_{{{s}_{r+1}}}\Phi_{\bullet \bullet}^{{{s}_{r}}{{s}_{r+1}}}R_{A\bullet}^{{{s}_{r+1}}\dagger}$.

Note: For a graphical representation see figure A2.

Zoom In Zoom Out Reset image size

The procedure can also be easily translated to the left boundary exploiting left orthogonality, where the translation invariance of the left-orthogonal matrices ${{L}^{{{s}_{j}}}}$ is then required.

Method I is also applicable when ${{H}_{\text{R}}}$ is time dependent, e.g. in case of a global quench.

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}^{{{s}_{j}}}}$ 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 ${{\hat{H}}_{\triangleleft,\ell +1}}:={{\hat{H}}_{\text{L}}}+{{\hat{h}}_{\ell ,\ell +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}_{\ell +1}}}}$. Note that for this method the matrices ${{L}^{{{s}_{j}}}}$ in the left part need not be translation invariant.

Since matrices ${{L}^{{{s}_{j}}}}$ are not changed during this update, we rewrite the wavefunction in MPS form after the CMW update in terms of the auxiliary basis states $|{{a}_{\ell}}\rangle =\sum\nolimits_{{{s}_{j\leqslant \ell}}}{{\left(\ldots {{L}^{{{s}_{\ell -1}}}}{{L}^{{{s}_{\ell}}}}\right)}_{{{a}_{\ell}}}}|\ldots {{s}_{\ell -1}}{{s}_{\ell}}\rangle$ connecting ${{L}^{{{s}_{\ell}}}}$ and ${{A}^{{{s}_{\ell +1}}}}$:

$$\begin{eqnarray}&&\Psi({{a}_{\ell}}|s_{\ell +1}^{\circ}s_{\ell +2}^{\bullet}\ldots ;t)={{(A_{\circ}^{{{s}_{\ell +1}}}A_{\bullet}^{{{s}_{\ell +2}}}\ldots )}_{{{a}_{\ell}}}},\end{eqnarray} \tag{ A.7 }$$

where ${{a}_{\ell}}$ is the left index of matrix $A_{\circ}^{{{s}_{\ell +1}}}$. The right-hand side of (A.7) is formally just the semi-infinite product of all matrices to the right of site $\ell$. The overall state vector after the CMW update can thus also be written

$$\begin{eqnarray}&&|\Psi(t)\rangle =\underset{{{a}_{\ell}}}{\mathop \sum}\underset{{{s}_{\ell +1}}\ldots}{\mathop \sum}\Psi\left({{a}_{\ell}}|s_{\ell +1}^{\circ}s_{\ell +2}^{\bullet}\ldots ;t\right)|{{a}_{\ell}}\rangle |{{s}_{\ell +1}}\ldots \rangle .\end{eqnarray} \tag{ A.8 }$$

Here • and ° again mark sites which have received a complete and incomplete update, respectively. Notice that in Method II the basis $\left\{|{{a}_{\ell}}\rangle \right\}$ will remain unchanged at all times. For a graphical representation of (A.7) see figure A3(a).

Zoom In Zoom Out Reset image size

We now need a renormalised representation $H_{\triangleleft,\ell +1}^{\text{eff}}$ of ${{\hat{H}}_{\triangleleft,\ell +1}}$ formulated in the block-spin basis $\left\{|{{a}_{\ell}}{{s}_{\ell +1}}\rangle \right\}$. A possible method to calculate $H_{\triangleleft,\ell +1}^{\text{eff}}$ is outlined in appendix A.5. Once we have such a renormalised expression we can determine the renormalised time evolution operator for the left part

$$\begin{eqnarray}&&U(\tau )_{\triangleleft,\ell +1}^{\text{eff}}:=\exp (-\text{i}\tau H_{\triangleleft,\ell +1}^{\text{eff}}),\end{eqnarray} \tag{ A.9 }$$

where we use the same small time step τ as for the Suzuki–Trotter updates. This time evolution operator is then used to update $\Psi\left({{a}_{\ell}}|s_{\ell +1}^{\circ}s_{\ell +2}^{\bullet}\ldots ;t\right)$. However, as it is a unitary operator defined in the block-spin basis $\left\{|{{a}_{\ell}}{{s}_{\ell +1}}\rangle \right\}$, it only updates $A_{\circ}^{{{s}_{\ell +1}}}$ and we get

$$\begin{eqnarray}&&A_{\bullet}^{{{s}_{\ell +1}}}=\underset{s_{\ell +1}^{\prime}}{\mathop \sum}U(\tau )_{\triangleleft,\ell +1}^{\text{eff},{{s}_{\ell +1}}s_{\ell +1}^{\prime}}A_{\circ}^{s_{\ell +1}^{\prime}}.\end{eqnarray} \tag{ A.10 }$$

This concludes one time step for the left boundary. Notice that the matrices in the left part are not updated as all changes in the left part are compressed into the boundary matrix ${{A}^{{{s}_{\ell +1}}}}$ with constant basis $\left\{|{{a}_{\ell}}\rangle \right\}$. In this sense the update is non-adaptive.

Notice also that $U_{\triangleleft,\ell +1}^{\text{eff}}$ breaks the structure of the even–odd 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_{\circ}^{{{s}_{\ell +1}}}$. For an algorithmic summary see table A2, for a graphical representation of this update see figure A4.

Table A2. Algorithm for Method II updating the left boundary of the CMW.

(i)	Perform only once for each CMW position:
	(a) Determine renormalised expression $H_{\triangleleft,\ell +1}^{\text{eff}}$ of ${{\hat{H}}_{\triangleleft,\ell +1}}$, formulated in the block-spin basis $\left\{|{{a}_{\ell}}{{s}_{\ell +1}}\rangle \right\}$.
	(b) Calculate $U(\tau )_{\triangleleft,\ell +1}^{\text{eff}}=\exp (-\text{i}\tau H_{\triangleleft,\ell +1}^{\text{eff}})$.
(ii)	Update $A_{\circ}^{s_{\ell +1}^{\prime}}$ using $U_{\triangleleft,\ell +1}^{\text{eff}}$ in each time step to get $A_{\bullet}^{{{s}_{\ell +1}}}=\sum\limits_{s_{\ell +1}^{\prime}}U\left(\tau \right)_{\triangleleft,\ell +1}^{\text{eff},{{s}_{\ell +1}}s_{\ell +1}^{\prime}}A_{\circ}^{s_{\ell +1}^{\prime}}$.

Note: For a graphical representation see figure A4.

Zoom In Zoom Out Reset image size

For determining $H_{\triangleleft,\ell +1}^{\text{eff}}$ used in Method II, we first assume a left part that is semi-infinite. Consider $\hat{H}$ in MPO form [68–72]

$$\begin{eqnarray}&&\hat{H}=\underset{\boldsymbol{s}}{\mathop \sum}\underset{j=-\infty}{\overset{\infty}{\mathop \prod}}W_{\left[\,j\right]}^{{{s}_{j}}s_{j}^{\prime}}|\boldsymbol{s}\rangle \langle \boldsymbol{s}^{{\prime}}|,\end{eqnarray} \tag{ A.11 }$$

where $W_{\left[\,j\right]}^{{{s}_{j}}s_{j}^{\prime}}$ are matrices of some dimension ${{d}_{W}}\times {{d}_{W}}$ containing operator elements ${{O}^{{{s}_{j}}s_{j}^{\prime}}}$. This decomposition can also be written in operator form, where we define ${{\hat{W}}_{\left[\,j\right]}}:=\sum\nolimits_{{{s}_{j}}s_{j}^{\prime}}W_{\left[\,j\right]}^{{{s}_{j}}s_{j}^{\prime}}|{{s}_{j}}\rangle \langle s_{j}^{\prime}|$ as ${{d}_{W}}\times {{d}_{W}}$ matrices containing operators. We can then simply write $\hat{H}=\prod\nolimits_{j=-\infty}^{\infty}{{\hat{W}}_{\left[\,j\right]}}.$ For finite size (or semi-infinite) operators, this product of MPO matrices is terminated by d_W-dimensional operator-valued boundary vectors ${{\hat{w}}_{\langle \left[j\right]}}$ and (or) ${{\hat{w}}_{\left[j\right]\rangle}}$.

An example for an MPO decomposition for the transverse field Ising (TFI) Hamiltonian ${{\hat{H}}_{\text{TFI}}}=-\sum\nolimits_{j}\hat{S}_{j}^{x}\hat{S}_{j+1}^{x}-h\sum\nolimits_{j}\hat{S}_{j}^{z}$ is given by

$$\begin{eqnarray}{{\hat{W}}_{\left[\,j\right]}}=\left[\begin{array}{c} {{{\hat{\mathbb{1}}}}_{j}} & 0 & 0 \\ -\hat{S}_{j}^{x} & 0 & 0 \\ -h\hat{S}_{j}^{z} & \hat{S}_{j}^{x} & {{{\hat{\mathbb{1}}}}_{j}} \end{array}\right]\end{eqnarray} \tag{ A.12 }$$

$$\begin{eqnarray}{{\hat{w}}_{\langle \left[\,j\right]}}=\left[\begin{array}{c} -h\hat{S}_{j}^{z} & \hat{S}_{j}^{x} & {{{\hat{\mathbb{1}}}}_{j}} \end{array}\right]\qquad {{\hat{w}}_{\left[\,j\right]\rangle}}={{\left[\begin{array}{c} {{{\hat{\mathbb{1}}}}_{j}} & -\hat{S}_{j}^{x} & -h\hat{S}_{j}^{z} \end{array}\right]}^{T}}.\end{eqnarray} \tag{ A.13 }$$

For the TFI Hamiltonian we thus have d_W  =  3.

We can express ${{\hat{H}}_{\triangleleft,\ell +1}}$ in terms of an MPO as

$$\begin{eqnarray}&&{{\hat{H}}_{\triangleleft,\ell +1}}=\ldots {{\hat{W}}_{\left[\ell -1\right]}}{{\hat{W}}_{\left[\ell \right]}}{{\hat{w}}_{\left[\ell +1\right]\rangle}},\end{eqnarray} \tag{ A.14 }$$

where we have terminated the semi-infinite product of MPO matrices with the boundary vector ${{\hat{w}}_{\left[\ell +1\right]\rangle}}$.

In order to get $H_{\triangleleft,\ell +1}^{\text{eff}}$ we use matrices ${{L}^{{{s}_{j}}}}$ to renormalise ${{\hat{H}}_{\triangleleft,\ell +1}}$. For this, consider the ${{d}_{W}}\times {{d}_{W}}$ dimensional MPO transfer matrix defined as

$$\begin{eqnarray}&&F_{\left[\,j\right]}^{{{b}_{j-1}}{{b}_{j}}}:=\underset{{{s}_{j}}s_{j}^{\prime}}{\mathop \sum}W_{\left[\,j\right],{{b}_{j-1}}{{b}_{j}}}^{{{s}_{j}}s_{j}^{\prime}}{{\bar{L}}^{{{s}_{j}}}}\otimes {{L}^{s_{j}^{\prime}}},\end{eqnarray} \tag{ A.15 }$$

containing ${{m}^{2}}\times {{m}^{2}}$ matrices, where m is the matrix dimension of the MPS matrices ${{L}^{{{s}_{j}}}}$ and ${{\bar{L}}^{{{s}_{j}}}}$ denotes the complex conjugate of ${{L}^{{{s}_{j}}}}$.

$H_{\triangleleft,\ell +1}^{\text{eff}}$ can then be written as

$$\begin{eqnarray}&&H_{\triangleleft,\ell +1}^{\text{eff},{{s}_{\ell +1}}s_{\ell +1}^{\prime}}=\underset{{{b}_{\ell}}=1}{\overset{{{d}_{W}}}{\mathop \sum}}{{\left(\underset{j\leqslant \ell}{\mathop \prod}{{F}_{\left[\,j\right]}}\right)}_{{{b}_{\ell}}}}w_{\left[\ell +1\right]\rangle ,{{b}_{\ell}}}^{{{s}_{\ell +1}}s_{\ell +1}^{\prime}}\end{eqnarray} \tag{ A.16 }$$

$$\begin{eqnarray}\begin{array}{*{20}{l}}&=\underset{{{b}_{\ell}}=1}{\overset{{{d}_{W}}}{\mathop \sum}}E_{\left[\ell \right]}^{{{b}_{\ell}}}w_{\left[\ell +1\right]\rangle ,{{b}_{\ell}}}^{{{s}_{\ell +1}}s_{\ell +1}^{\prime}},\end{array}\end{eqnarray} \tag{ A.17 }$$

where we have defined the semi-infinite product

$$\begin{eqnarray}&&E_{\left[\ell \right]}^{{{b}_{\ell}}}:={{\left(\underset{j\leqslant \ell}{\mathop \prod}{{F}_{\left[\,j\right]}}\right)}_{{{b}_{\ell}}}}.\end{eqnarray} \tag{ A.18 }$$

$H_{\triangleleft,\ell +1}^{\text{eff},{{s}_{\ell +1}}s_{\ell +1}^{\prime}}$ is then a set of $m\times m$ matrices labelled by ${{s}_{\ell +1}}$ and $s_{\ell +1}^{\prime}$ and ${{E}_{\left[\ell \right]}}$ can be understood as a d_W-dimensional vector containing $m\times m$ matrices. For a graphical representation of these steps see figure A4(1). Note that the vector element $E_{\left[j\right]}^{1}$ accumulates the renormalised Hamiltonian containing all sites $k\leqslant j$ (see e.g. [69]).

To determine (A.17) we need a way to calculate the semi-infinite matrix product ${{E}_{\left[\ell \right]}}$. For the moment we consider the case where both ${{\hat{H}}_{L}}$ and the matrices ${{L}^{{{s}_{j}}}}$ are translation invariant. In this case F_[ j] is also translation invariant and ${{E}_{\left[\ell \right]}}$ 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}_{\left[\ell \right]}}$, which is inspired by standard DMRG formulations. For this we relax the condition of semi-infinity for the left Hamiltonian ${{\hat{H}}_{\triangleleft,\ell +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 ${{\hat{H}}_{\triangleleft,\ell +1}}$ in MPO form is thus contracted also on the left side by ${{\hat{w}}_{\langle \left[k\right]}}$ for some $k\ll \ell$. We first compute an initial E_[k] as

$$\begin{eqnarray}&&E_{\left[k\right],{{b}_{k}}}^{\text{init}}=\underset{{{s}_{k}}s_{k}^{\prime}}{\mathop \sum}w_{\langle \left[k\right],{{b}_{k}}}^{{{s}_{k}}s_{k}^{\prime}}{{L}^{{{s}_{k}}\dagger}}{{L}^{s_{k}^{\prime}}},\end{eqnarray} \tag{ A.19 }$$

exploiting the left-orthogonality of the matrices ${{L}^{{{s}_{j}}}}$. For a graphical representation see figure A3(b).

We can now iteratively calculate ${{E}_{\left[\,j+1\right]}}={{E}_{\left[\,j\right]}}{{F}_{\left[\,j+1\right]}}$ many times until this process has converged. As a convergence criterion we can use, e.g. the ground state energy per site of the renormalised Hamiltonian which is accumulated in $E_{\left[\,j\right]}^{1}$. Using the converged result as an approximation for ${{E}_{\left[\ell \right]}}$ we can then easily determine $H_{\triangleleft,\ell +1}^{\text{eff}}$ from (A.17).

In the case where MPS matrices ${{L}^{{{s}_{j}}}}$ and/or MPOs $W_{\left[\,j\right]}^{{{s}_{j}}s_{j}^{\prime}}$ are site dependent for some sites $k\leqslant j\leqslant \ell$, we can first calculate E_[k] up to site k approximately as described above and then calculate the finite product

$$\begin{eqnarray}&&{{E}_{\left[\ell \right]}}=E_{\left[k\right]}^{\text{init}}\underset{j=k+1}{\overset{\ell}{\mathop \prod}}{{F}_{\left[\,j\right]}}.\end{eqnarray} \tag{ A.20 }$$

Notice that we can in principle even define the left part to be finite altogether, with site dependent matrices ${{L}^{{{s}_{j\leqslant \ell}}}}$ and/or site dependent MPOs $W_{\left[\,j\right]}^{{{s}_{j}}s_{j}^{\prime}}$, such that ${{E}_{\left[\ell \right]}}=\prod\nolimits_{j=1}^{\ell}{{F}_{\left[\,j\right]}}$ is also a finite product. In this case one would have to specify left boundary conditions. In our simulations we do not consider this case.

If Method II is used at the right boundary, we use the uniform matrices ${{R}^{{{s}_{j}}}}$ to construct a renormalised expression for ${{\hat{H}}_{r,\triangleright}}:={{\hat{h}}_{r,r+1}}+{{\hat{H}}_{\text{R}}}=\sum\nolimits_{j=r}^{\infty}{{\hat{h}}_{j,j+1}}$.

Generally the computational effort for calculating $H_{\triangleleft,\ell +1}^{\text{eff}}$ is dictated by the computational effort for calculating ${{E}_{\left[\ell \right]}}$. 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}^{{{s}_{j}}}}$.

In practice, it takes about 75 renormalisation steps for the TFI model at h  =  0.45 (m₀  =  30) and about 100 steps for the XXZ model at J_z  =  −2 (m₀  =  88) for convergence in energy up to an accuracy of 10⁻¹⁵, where m₀ 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.

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 no boundary update is used at the left boundary, the matrix ${{A}^{{{s}_{\ell +1}}}}$ is discarded. If Method II is used, we incorporate ${{A}^{{{s}_{\ell +1}}}}$ into the left part by using it to calculate a renormalised expression for ${{\hat{H}}_{\triangleleft,\ell +2}}:={{\hat{H}}_{\triangleleft,\ell +1}}+{{\hat{h}}_{\ell +1,\ell +2}}$. More precisely, we construct ${{F}_{\left[\ell +1\right]}}$ as defined in (A.15) using ${{A}^{{{s}_{\ell +1}}}}$

$$\begin{eqnarray}&&F_{\left[\ell +1\right]}^{{{b}_{\ell}}{{b}_{\ell +1}}}:=\underset{{{s}_{\ell +1}}s_{\ell +1}^{\prime}}{\mathop \sum}W_{\left[\ell +1\right],{{b}_{\ell}}{{b}_{\ell +1}}}^{{{s}_{\ell +1}}s_{\ell +1}^{\prime}}{{\bar{A}}^{{{s}_{\ell +1}}}}\otimes {{A}^{s_{\ell +1}^{\prime}}}.\end{eqnarray} \tag{ A.21 }$$

With ${{E}_{\left[\ell \right]}}$ from earlier calculations we can then construct ${{E}_{\left[\ell +1\right]}}={{E}_{\left[\ell \right]}}{{F}_{\left[\ell +1\right]}}$, calculate $H_{\triangleleft,\ell +2}^{\text{eff}}$ as defined in (A.17) and determine $U(\tau )_{\triangleleft,\ell +2}^{\text{eff}}=\exp (-\text{i}\tau H_{\triangleleft,\ell +2}^{\text{eff}})$.

At the right boundary we introduce $R_{A}^{{{s}_{r+1}}}$ as a new rightmost matrix ${{B}^{{{s}_{r+1}}}}$.

After the window has been moved by a single site according to these steps, we redefine $\ell \leftarrow \ell +1$, $r\leftarrow r+1$ (and we exchange labels $A\leftrightarrow B$ in the case of iTEBD).

Notice that for the left boundary the dimension of the block basis $\left\{|{{a}_{\ell}}\rangle \right\}$ 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}_{\text{max}}}$.

The TFI model

$$\begin{eqnarray}&&{{\hat{H}}_{\text{TFI}}}=-\underset{j}{\mathop \sum}\hat{S}_{j}^{x}\hat{S}_{j+1}^{x}-h\underset{j}{\mathop \sum}\hat{S}_{j}^{z}\end{eqnarray} \tag{ B.1 }$$

can be solved exactly [59, 60] by first transforming to spinless fermionic operators $c_{j}^{\dagger}$, c_j via a Jordan–Wigner (JW) transformation [61]

$$\begin{eqnarray}&&\hat{S}_{j}^{+}=\underset{n<j}{\mathop \prod}(1-2c_{n}^{\dagger}{{c}_{n}})c_{j}^{\dagger},\qquad \hat{S}_{j}^{-}=\underset{n<j}{\mathop \prod}(1-2c_{n}^{\dagger}{{c}_{n}}){{c}_{j}},\end{eqnarray} \tag{ B.2 }$$

where $\hat{S}_{j}^{+}$ and $\hat{S}_{j}^{-}$ are the spin raising and lowering operators. With $\hat{S}_{j}^{x}=\frac{1}{2}(\hat{S}_{j}^{+}+\hat{S}_{j}^{-})$ and $\hat{S}_{j}^{z}=\hat{S}_{j}^{+}\hat{S}_{j}^{-}-\frac{1}{2}$ the Hamiltonian becomes

$$\begin{eqnarray}&&{{\hat{H}}_{\text{TFI}}}=-\frac{1}{4}\underset{j}{\mathop \sum}(c_{j}^{\dagger}-{{c}_{j}})(c_{j+1}^{\dagger}+{{c}_{j+1}})-h\underset{j}{\mathop \sum}\left(c_{j}^{\dagger}{{c}_{j}}-\frac{1}{2}\right).\end{eqnarray} \tag{ B.3 }$$

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 $\text{O}(1/L)$ where L is the system size).

A subsequent Bogoliubov transformation [73] to fermionic operators ${{\eta}_{k}}$, $\eta _{k}^{\dagger}$ in momentum space

$$\begin{eqnarray}&&{{c}_{j}}=\frac{1}{\sqrt{2\pi}}\int_{-\pi}^{\pi}\text{d}k\,{{\text{e}}^{\text{i}kj}}({{a}_{k}}{{\eta}_{k}}-\text{i}{{b}_{k}}\eta _{-k}^{\dagger})\end{eqnarray} \tag{ B.4 }$$

then diagonalises the Hamiltonian. The coefficients a_k and b_k are real and satisfy

$$\begin{eqnarray}&&{{a}_{-k}}={{a}_{k}},\qquad {{b}_{-k}}=-{{b}_{k}},\qquad a_{k}^{2}+b_{k}^{2}=1\end{eqnarray} \tag{ B.5 }$$

and can be determined as

$$\begin{eqnarray}&&{{a}_{k}}=\frac{{{\varepsilon}_{k}}-\frac{1}{2}\cos (k)-h}{\sqrt{2{{\varepsilon}_{k}}\left({{\varepsilon}_{k}}-\frac{1}{2}\cos (k)-h\right)}}\end{eqnarray} \tag{ B.6 }$$

$$\begin{eqnarray}&&{{b}_{k}}=-\frac{\frac{1}{2}\sin (k)}{\sqrt{2{{\varepsilon}_{k}}\left({{\varepsilon}_{k}}-\frac{1}{2}\cos (k)-h\right)}}\end{eqnarray} \tag{ B.7 }$$

$$\begin{eqnarray}&&{{\varepsilon}_{k}}=\sqrt{\frac{1}{4}+h\cos (k)+{{h}^{2}}}.\end{eqnarray} \tag{ B.8 }$$

The Hamiltonian then reads

$$\begin{eqnarray}&&{{\hat{H}}_{\text{TFI}}}=\frac{1}{2\pi}\int_{-\pi}^{\pi}\text{d}k\,{{\varepsilon}_{k}}\left(\eta _{k}^{\dagger}{{\eta}_{k}}-\frac{1}{2}\right)\end{eqnarray} \tag{ B.9 }$$

and the ground state corresponds to the vacuum state $|0\rangle$ in terms of the fermionic operators ${{\eta}_{k}}$ and $\eta _{k}^{\dagger}$.

The propagation of a signal induced on top of the ground state $|0\rangle$ 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 ${{\varepsilon}_{k}}$ created by $\eta _{k}^{\dagger}$. In this picture, the maximum velocity ${{v}_{\text{max}}}$ of the signal can be exactly calculated as the maximum of the group velocity

$$\begin{eqnarray}&&{{v}_{k}}:=\frac{\text{d}{{\varepsilon}_{k}}}{\text{d}k}=\frac{h}{2}\frac{\sin (k)}{{{\varepsilon}_{k}}}.\end{eqnarray} \tag{ B.10 }$$

A short calculation shows that v_k takes its extrema at $\cos (k)=2h$ and $\cos (k)=\frac{1}{2h}$, which gives

$$\begin{eqnarray}{{v}_{\text{max}}}=\left\{\begin{array}{*{35}{l}} h, & h\leqslant {{h}_{\text{crit}}} \\ {{h}_{\text{crit}}}, & h\geqslant {{h}_{\text{crit}}} \end{array}\right. \end{eqnarray} \tag{ B.11 }$$

where ${{h}_{\text{crit}}}=\frac{1}{2}$.

Consider a Jordan–Wigner (JW) excitation at site $\ell$ on top of the thermodynamic limit ground state $|0\rangle$, defined as

$$\begin{eqnarray}&&|\psi {{\rangle}_{\ell}}=(c_{\ell}^{\dagger}+{{c}_{\ell}})|0\rangle \end{eqnarray} \tag{ B.12 }$$

$$\begin{eqnarray}&&=\frac{1}{\sqrt{2\pi}}\int_{-\pi}^{\pi}\text{d}k\,{{\text{e}}^{-\text{i}k\ell}}\left({{a}_{k}}+\text{i}{{b}_{k}}\right)\eta _{k}^{\dagger}|0\rangle .\end{eqnarray} \tag{ B.13 }$$

Using the results from the previous sections for the TFI model, the time evolution of the magnetisation in z after such an excitation

$$\begin{eqnarray}&&{{\langle {{\hat{S}}^{z}}(n,t)\rangle}_{\ell}}={{\langle \psi {{|}_{\ell}}c_{n}^{\dagger}(t){{c}_{n}}(t)|\psi \rangle}_{\ell}}-\frac{1}{2}.\end{eqnarray} \tag{ B.14 }$$

can be calculated analytically.

Solving the Heisenberg equation of motion for ${{\eta}_{k}}(t)$ yields ${{\eta}_{k}}(t)={{\text{e}}^{-\text{i}{{\varepsilon}_{k}}t}}{{\eta}_{k}}$. Using (B.4) then allows one to write

$$\begin{eqnarray}&&{{c}_{n}}(t)=\frac{1}{\sqrt{2\pi}}\int_{-\pi}^{\pi}\text{d}k\,{{\text{e}}^{\text{i}kn}}({{\alpha}_{k}}(t)\,{{\eta}_{k}}+{{\beta}_{k}}(t)\,\eta _{-k}^{\dagger})\end{eqnarray} \tag{ B.15 }$$

with ${{\alpha}_{k}}(t)={{a}_{k}}{{\text{e}}^{-\text{i}{{\varepsilon}_{k}}t}}$ and ${{\beta}_{k}}(t)=-\text{i}{{b}_{k}}{{\text{e}}^{\text{i}{{\varepsilon}_{k}}t}}$. Plugging (B.15) and (B.13) into (B.14) yields after some calculation

$$\begin{eqnarray}&&\begin{array}{|c|} \hline \langle{\hat{S}^{z}(n,t)}\rangle_{\ell} = S^{z}_{\rm GS} + \frac{I^{A}_{\ell}(n,t) - I^{B}_{\ell}(n,t)}{(2\pi)^2},\\ \hline \end{array}\end{eqnarray} \tag{ B.16 }$$

where

$$\begin{eqnarray}&&S_{\text{GS}}^{z}=\langle 0|\hat{S}_{j}^{z}|0\rangle =\frac{1}{2\pi}\int_{-\pi}^{\pi}\text{d}k\,|{{b}_{k}}{{|}^{2}}-\frac{1}{2}\end{eqnarray} \tag{ B.17 }$$

is the ground state magnetisation and

$$\begin{eqnarray}&&I_{\ell}^{A}(n,t)=\left|\int_{-\pi}^{\pi}\text{d}k\,{{\text{e}}^{-\text{i}\left[k\left(\ell -n\right)+{{\varepsilon}_{k}}t\right]}}\left({{a}_{k}}+\text{i}{{b}_{k}}\right){{a}_{k}}\right|^{2}\end{eqnarray} \tag{ B.18 }$$

$$\begin{eqnarray}&&I_{\ell}^{B}(n,t)=\left|\int_{-\pi}^{\pi}\text{d}k\,{{\text{e}}^{-\text{i}\left[k\left(\ell -n\right)-{{\varepsilon}_{k}}t\right]}}\left({{a}_{k}}+\text{i}{{b}_{k}}\right){{b}_{k}}\right|^{2}.\end{eqnarray} \tag{ B.19 }$$

In appendix D we use (B.16) to compare with results obtained from a CMW simulation.

The XXZ model defined by the Hamiltonian

$$\begin{eqnarray}&&{{H}_{\text{XXZ}}}=-\underset{j}{\mathop \sum}\hat{S}_{j}^{x}\hat{S}_{j+1}^{x}+\hat{S}_{j}^{y}\hat{S}_{j+1}^{y}+\Delta \hat{S}_{j}^{z}\hat{S}_{j+1}^{z}\end{eqnarray} \tag{ C.1 }$$

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 $\Delta <-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₀(x), which for $\Delta <-1$ satisfies the following integral equation

$$\begin{eqnarray}\begin{array}{*{20}{l}}{{g}_{0}}(x)&+&\frac{\sinh (2\Phi)}{2\pi}\int_{-\pi}^{\pi}\frac{{{g}_{0}}\left({{x}^{\prime}}\right)\text{d}{{x}^{\prime}}}{\cosh (2\Phi)-\cos \left(x-{{x}^{\prime}}\right)}\end{array}\end{eqnarray} \tag{ C.2 }$$

$$\begin{eqnarray}\begin{array}{*{20}{l}}&=&\frac{2\sinh (\Phi)}{\cosh (\Phi)-\cos (x)},\end{array}\end{eqnarray} \tag{ C.3 }$$

where $\cosh (\Phi)=-\Delta$.

The solution to this integral equation is given by

$$\begin{eqnarray}&&{{g}_{0}}(x)=\frac{2K\left({{m}_{0}}\right)}{\pi}\text{dn}\left(\frac{K\left({{m}_{0}}\right)}{\pi}x,{{m}_{0}}\right),\quad -\pi <x<\pi .\end{eqnarray} \tag{ C.4 }$$

Here $\text{dn}(x,m)$ is a Jacobian elliptic function [76], K(m) the complete elliptic integral of the first kind

$$\begin{eqnarray}&&K(m)=\int_{0}^{\frac{\pi}{2}}\frac{\text{d}x}{\sqrt{1-m\,{{\sin}^{2}}\,x}}\end{eqnarray} \tag{ C.5 }$$

and the parameter m₀ the solution of the equation

$$\begin{eqnarray}&&\frac{K\left({{m}_{0}}\right)}{K\left(1-{{m}_{0}}\right)}=\frac{\pi}{\Phi}.\end{eqnarray} \tag{ C.6 }$$

The root density g₀(x) can then be used to calculate various quantities such as the ground state energy and elementary excitations.

To calculate the maximum signal velocity ${{v}_{\text{max}}}$ as a function of $\Delta$ we first determine the dispersion relation ${{\varepsilon}_{k}}$ for the elementary excitations. As for the TFI model in appendix B.2 we then obtain ${{v}_{\text{max}}}$ as the maximum of the group velocity ${{v}_{k}}=\frac{\text{d}{{\varepsilon}_{k}}}{\text{d}k}$.

The dispersion of the elementary excitations is given by [75]

$$\begin{eqnarray}&&{{\varepsilon}_{k}}=\frac{1}{2}\sinh (\Phi)\left[{{g}_{0}}\left({{x}_{0}}(k)\right)-{{g}_{0}}\left(\pi \right)\right]+G\left(\Delta \right),\end{eqnarray} \tag{ C.7 }$$

where $G\left(\Delta \right)$ is the finite energy gap present in this phase and x₀(k) has to be determined by inverting

$$\begin{eqnarray}&&k\left({{x}_{0}}\right)=\frac{1}{2}\left(\pi -\int_{0}^{{{x}_{0}}}{{g}_{0}}(x)\text{d}x\right),\quad -\pi \leqslant {{x}_{0}}\leqslant \pi \end{eqnarray} \tag{ C.8 }$$

for a given momentum k.

From this we can calculate the group velocity

$$\begin{eqnarray}&&{{v}_{k}}=\frac{\text{d}{{\varepsilon}_{k}}}{\text{d}k}=\frac{\text{d}{{\varepsilon}_{k}}}{\text{d}{{x}_{0}}}{{\left(\frac{\text{d}k}{\text{d}{{x}_{0}}}\right)}^{-1}}\end{eqnarray} \tag{ C.9 }$$

where we need

$$\begin{eqnarray}&&\frac{\text{d}{{\varepsilon}_{k}}}{\text{d}{{x}_{0}}}=\frac{1}{2}\sinh (\Phi)\frac{\text{d}{{g}_{0}}}{\text{d}{{x}_{0}}},\qquad \frac{\text{d}k}{\text{d}{{x}_{0}}}=-\frac{{{g}_{0}}\left({{x}_{0}}\right)}{2}.\end{eqnarray} \tag{ C.10 }$$

Using some properties of Jacobian elliptic functions [76] and defining ${{K}_{0}}=K\left({{m}_{0}}\right)$ we get

$$\begin{eqnarray}&&\frac{\text{d}{{g}_{0}}}{\text{d}{{x}_{0}}}=-2\,{{m}_{0}}{{\left(\frac{{{K}_{0}}}{\pi}\right)}^{2}}\text{sn}\left(\frac{{{K}_{0}}{{x}_{0}}}{\pi},{{m}_{0}}\right)\text{cn}\left(\frac{{{K}_{0}}{{x}_{0}}}{\pi},{{m}_{0}}\right)\end{eqnarray} \tag{ C.11 }$$

where $\text{sn}(x,m)$ and $\text{cn}(x,m)$ are also Jacobian elliptic functions. Defining

$$\begin{eqnarray}&&S\left({{x}_{0}},{{K}_{0}},{{m}_{0}}\right):=\text{sn}\left(\frac{{{K}_{0}}{{x}_{0}}}{\pi},{{m}_{0}}\right)\text{cn}\left(\frac{{{K}_{0}}{{x}_{0}}}{\pi},{{m}_{0}}\right)\end{eqnarray} \tag{ C.12 }$$

we can then write

$$\begin{eqnarray}&&\begin{array}{|c|} \hline v_{k} = \frac{2\,m_{0}\sinh(\Phi)}{g_{0}(x_{0})}\left(\frac{K_{0}}{\pi}\right)^{2}S(x_{0},K_{0},m_{0}).\\ \hline \end{array}\end{eqnarray} \tag{ C.13 }$$

We determine the maximum of this function numerically to get ${{v}_{\text{max}}}$ as a function of $\Delta$. The values of ${{v}_{\text{max}}}$ for various interaction strengths $\Delta <-1$ are listed in table C1 and are obtained with a numerical precision of 10⁻⁸.

Table C1. Values for the maximum signal velocity ${{v}_{\text{max}}}$ in the XXZ model for various values of the interaction strength $\Delta <-1$.

$\Delta$	${{v}_{\text{max}}}$	$\Delta$	${{v}_{\text{max}}}$
−1.5	1.781 734 04	−3.5	1.959 201 77
−2.0	1.875 595 02	−4.0	1.968 758 00
−2.5	1.920 144 92	−4.5	1.975 312 55
−3.0	1.944 491 13	−5	1.980 002 06

Note: These values are obtained from numerically finding the maximum of (C.12) with a numerical precision of 10⁻⁸.

In our simulations we use a Trotter step size of $\tau =0.002$ and a maximum matrix dimension ${{m}_{\text{max}}}=120$. The unscaled 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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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 $S_{\text{ent}}^{\text{TFI}}(n,t)$ after a JW excitation at h  =  0.45 in figure E3. It is smaller than the entanglement entropy for tight binding fermions after a domain wall excitation at all times. In the latter case, $S_{\text{ent}}^{\text{TB}}(n,t)$ approaches the asymptotic scaling function H(y) without any scaling in time. In fact, $S_{\text{ent}}^{\text{TFI}}(n,t)$ decreases in time, whereas $S_{\text{ent}}^{\text{TB}}(n,t)$ approaches H(y) from below. The exact relation between the scaled excess longitudinal magnetisation $M(n,t)=({{S}^{x}}(n,t)-S_{\text{GS}}^{x})/|2S_{\text{GS}}^{x}|$ and the fermion density N_TB(n, t) described in the main text therefore does not carry over to the entanglement entropy.

Zoom In Zoom Out Reset image size

For the XXZ simulations we use a Trotter step size of $\tau =0.01$ and maximum matrix dimensions of ${{m}_{\text{max}}}=150$ for $\Delta =-4$, ${{m}_{\text{max}}}=160$ for $\Delta =-3$ and ${{m}_{\text{max}}}=180$ for $\Delta =-2$. We show a representative plot of the unscaled staggered magnetisation ${{\tilde{S}}^{z}}(n,t)={{\left(-1\right)}^{n}}{{S}^{z}}(n,t)$ at $\Delta =-3$ after a JW excitation on top of the infinite system ground state versus absolute position n at various times 800  <  t  <  960 in figure E4. Again we observe a ballistic propagation of the signal front as well as magnetisation steps near the front as in the TFI model case. For the XXZ model an additional micro step structure appears due to 'pairing' of neighbouring sites, which is due to the spinon nature of the elementary excitations created by the signal (see section 3.2).

Zoom In Zoom Out Reset image size

The scaling behaviour of the larger step structure is investigated in detail in the main text. The overall shape of the unscaled staggered magnetisation ${{\tilde{S}}^{z}}(n,t)$ 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.

We again consider the TFI model at h  =  0.45 after a JW excitation in the infinite system ground state. We use a non-moving window with N  =  50 and maximum bond dimension ${{m}_{\text{max}}}=120$, where the ground state MPS representation has bond dimension m₀  =  30. The time evolution of the bipartite entanglement entropy ${{S}_{\text{ent}}}(n,t)$ and the magnetisation S^x(n, t) can be seen in figure F1. The signal reaches the boundaries at $t\approx 40$ and reflections start to emerge at $t\approx 90$.

Zoom In Zoom Out Reset image size

We compare the magnetisation S^x(n, t) of this simulation with the magnetisation $S_{\text{ref}}^{x}(n,t)$ of the reference simulation of appendix D and show their absolute difference $\Delta M(n,t)=\,\left|{{S}^{x}}(n,t)-S_{\text{ref}}^{x}(n,t)\right|$ in figure F2, where subplot (a) shows $\Delta 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 $\Delta M(n,t)$ versus position n inside the non-moving window at various times t.

Zoom In Zoom Out Reset image size

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\approx 50$. The deviation at both boundaries then increases exponentially further until $t\approx 100$, where it becomes of the order $\text{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.

We also consider the S  =  1 bilinear, biquadratic chain at the AKLT point [63] defined by the Hamiltonian

$$\begin{eqnarray}&&\hat{H}=\underset{j}{\mathop \sum}{{\mathbf{\hat{S}}}_{j}}\cdot {{\mathbf{\hat{S}}}_{j+1}}+\frac{1}{3}{{({{\mathbf{\hat{S}}}_{j}}\cdot {{\mathbf{\hat{S}}}_{j+1}}\,)}^{2}}.\end{eqnarray} \tag{ F.1 }$$

The ground state is a valence bond state and has an exact MPS representation with bond dimension m₀  =  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 ${{m}_{\text{max}}}=200$. The time evolution of the bipartite entanglement entropy ${{S}_{\text{ent}}}(n,t)$ and the magnetisation S^z(n, t) can be seen in figure F3.

Zoom In Zoom Out Reset image size

Here the signal impacting at $t\approx 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₀  =  2, which is much too small for the matrices to absorb this information for a long time span.