Detecting topology through dynamics in interacting fermionic wires

We describe a protocol to read out the topological invariant of interacting 1D chiral models, based on measuring the mean chiral displacement of time-evolving bulk excitations. We present analytical calculations and numerical Matrix Product State simulations of interacting Su-Schrieffer-Heeger (SSH) chains, demonstrating how the mean chiral displacement allows to distinguish between topological insulator, trivial insulator and symmetry-broken phases. Finally, we provide an experimental blueprint for realizing a model displaying these three phases and describe how to detect those.

Since the discovery of sharp resistivity jumps in Quantum Hall devices and their explanation in terms of quantized invariants, topology played an ever increasing role in modern condensed matter [1,2]. Transport experiments are often employed in solid state systems to observe quantized responses, but in ultracold atoms, polaritonic and photonic systems other avenues are generally pursued, such as direct imaging of edge states, detection of anomalous displacement or interferometry [3,4].
Topological invariants are generally defined as integrals over momentum space. In this paper we will be focusing specifically on one-dimensional chiral models, for which the relevant invariant is the winding number of a vector in the two-dimensional space dictated by the symmetry [5]. In presence of interactions, the winding number can be defined as the chiral winding of the zero-frequency component of the Fouriertransformed single-particle imaginary-time Green's function g ≡ G(k, ω = 0), via [6,7] γ = tr dk 4πi Γg −1 ∂ k g , with g −1 the matrix inverse of g, and Γ the chiral symmetry operator, which anti-commutes with the Bloch Hamiltonian if this symmetry is present [8]. For non-interacting systems with Hamiltonian H 0 (k) the Green's function is G(k, ω) = [iω − H 0 ] −1 , so that g = H −1 0 and Eq. (1) reduces to the usual winding of the Bloch Hamiltonian.
Accessing momentum space to compute such invariants in experiments is possible [9], but still not always evident. Obvious examples are systems with broken translational invariance. There is therefore a high demand for efficient detection methods of topological phases of matter working directly in real space. One such method, based on spectral projectors, was introduced by Bianco and Resta for clean and regular 2D Chern insulators [10], and later employed to describe quasicrystals [11]. A real-space analysis of non-Hermitian models was instead performed in Refs. [12,13].
Another method which proved very successful to detect the topology of single-particle Bloch bands in 1D chiral models is the read-out of the mean chiral displacement (MCD), which can be defined in a many-body context as In the above, ΓX = x,τ,τ ,s,s c † x,τ,s (xΓ τ s,τ s )c x,τ ,s is defined via c x,τ,s -the annihilation operator of a fermionic particle with sublattice index τ and spin s acting on the unit-cell at position x (lattice spacing is set to unity). In this Letter, we only consider stationary Hamiltonians H such that time-evolution is readily given by U (t) = e −iHt/ . The mean chiral displacement thus measures the (chiral) propagation of an excitation created at the central unit-cell x = 0 after an evolution time t. In Refs. [14][15][16] it was shown theoretically and confirmed experimentally that the MCD of a localized single-particle excitation above the vacuum state converges in the long-time limit to the chiral invariant of one dimensional static, periodically-driven, and even disordered systems. Since then, the MCD has found multiple applications in very different systems [17][18][19][20][21][22][23].
In this work, we extend the study of the MCD to the case of many-particle systems. First, we show that the single-particle (disconnected) part of the MCD, measured after sudden creation of local excitations, coincides precisely with the many-body winding defined in Eq. (1). This implies that measuring the MCD on top of a halffilled Fermi sea provides a direct and exact detection of the winding. Then, we present nonperturbative matrix product states (MPS) simulations of a versatile 1D chiral model. In a weakly-correlated configuration, the connected part of the MCD remains negligible, and the MCD keeps converging reliably to the winding number. Switching instead to strong many-body correlations, we show that the MCD can also signal emergent symmetrybreaking long-range ordered phases. We conclude by discussing a blueprint scheme for the experimental detection of the MCD on a strongly-correlated model.
Mean chiral displacement as a topological marker.-Let us start by clarifying the relation between the many-body topological invariant γ defined via single-particle Green's functions and the MCD. For this purpose, let us consider a system in its ground state for times t < 0, and excite it by adding one particle at (x = 0, t = 0). The time dependence of the MCD evaluated on the perturbed ground state at half filling can be expanded (for details see Supplemental Material [24]) as with local density n x,τ,s = c † x,τ,s c x,τ,s , averages . . . evaluated on the unperturbed ground state, ξ(t) the connected (many-body) part of C(t) and G the one-body Green's function matrix, with components iG τ s,τ s (x, t) = c x,τ,s (t)c † 0,τ ,s (0) (t > 0). The second sum is constant over time, so that the time-dependence of the MCD is given by with ∆ξ(t) = ξ(t)−ξ(0 − ) the many-body contribution to the time-dependent MCD, where 0 − indicates the time before creating the excitation. For later convenience, we define the one-body contribution ∆C d (t) = ∆C − ∆ξ. In a non-interacting many-body system (e.g., a Fermi sea), the many-body part ∆ξ of the correlator is zero and the only remainder is the one-body part. When a single particle is added to a half-filled Fermi sea, the MCD is remarkably quantized, and equal to the topological invariant (see Fig. 1(c), and a detailed proof in [24]). When instead excitations are generated within an interacting system, or within a Fermi sea away from half-filling, scattering processes between conduction and valence band give rise to damped oscillations, or more complex behaviors. To filter away oscillatory contributions, it is useful to define the time-averaged MCD ∆C ≡ 1 T T 0 dt ∆C(t). In the following we show explicitly by nonperturbative MPS simulations, presented in Figs. 2 and 3, that the time-average of the one-body (disconnected) term provides an approximation of the winding number γ defined in real space and real time. Moreover, we show that the MCD is also a useful marker of spontaneous symmetry-breaking in trivial insulators, as the latter gives rise to a divergent many-body term ξ(t).
Model: Interacting SSH chains.-To illustrate our findings we focus on a specific 1D model, the interacting Su-Schrieffer-Heeger (SSH) fermionic chain, which is shown schematically in Fig. 1(a). The Hamiltonian of the model reads H = H 0 + H ⊥ + H . The kinetic part H 0 is a tight binding model with two sublattices (τ ∈ {A, B}) and internal spin-1/2 (s ∈ {↑, ↓}) degree of freedom denoted by The intra-cell and inter-cell tunneling amplitudes are denoted by J and J , respectively. In the following, we will set = 1, and time will be measured in units of 1/|J |. On top of the kinetic terms we consider Hubbard on-site interactions between ↑ and ↓ spins, and interactions which act between two identical spins located on the A and B sites of the same unit-cell, Note that H 0 is fully decoupled in the spin-sector, and, when considering only H 0 +H , it is redundant to account for both spin orbitals.
The non-interacting case.-Let us first review the topological properties of H 0 . This 1D Hamiltonian has chiral symmetry, which places it in the AIII class of topological insulators [25]. If moreover J and J are real, H 0 features particle-hole and time-reversal symmetries, and the model belongs to the more restrictive class BDI. For each of the two disconnected SSH chains, the Bloch Hamiltonian winds once or not at all around the origin, so that γ = 2 for |J/J | < 1 and γ = 0 for |J/J | > 1. This is particularly easy to see in two extreme limits. When J = 0, unit-cells are disconnected and G As,Bs (x, t) is zero unless x = 0, thus g As,Bs (k) is momentum-independent and with zero winding. In the other extreme case J = 0, we have g As,Bs = g 0 δ s,s e ik , which (summing over the spin degrees of freedom) yields a winding of 2 [7]. In Fig. 1(b) we evaluate the dynamics of excitations on top of a halffilled Fermi sea by exact diagonalization of the quadratic Hamiltonian (see [24] for details), and we show that the corresponding MCD captures correctly the winding.
The short-ranged case.-Let us now consider the Hamiltonian H sr = H 0 + H ⊥ , also known as Peierls-Hubbard model, whose topological properties were first studied in Ref. [7]. When |J/J | < 1 the system remains topological for every value of U ⊥ because the interaction terms in H ⊥ do not couple edge with bulk modes and therefore zero-energy edge excitations remain unaffected [7]. At |J| = |J |, the system is a 1D Hubbard model: while a charge gap is immediately opened by U ⊥ > 0, the spin sector remains gapless and thus provides the phase boundary between topological and trivial insulator. Its dimerized-and bond-limits (J = 0 and J = 0, respectively) are short-range correlated and thus we expect that ∆C ≈ ∆C d will oscillate around the winding number γ. In Fig. 2 we confirm this by computing the MCD and its disconnected part with MPS simulations. Our findings agree with the results of Ref. [7], and here we also provide the missing link to an observable which is easily accessible in ongoing experiments based on fermionic quantum gas microscopes [26][27][28][29] or optical tweezers [30].
Effective interacting spin-model.-To proceed further, we introduce an alternative representation of the Hamiltonian. This will allow us to illustrate the crucial differences between H and H ⊥ interactions on the topological properties of the underlying system. in which σ i denotes the Pauli matrix i ∈ {x, y, z}. These are particularly useful in the limit J = 0, because the individual terms of H become expressions in S i x,s that all clearly commute with each link-density n x+1,A,s + n x,B,s separately. For a translationally invariant system at halffilling, this implies n x+1,A,s + n x,B,s = 1 ∀x, s, with effective Hamiltonian As a consequence, H 0 + H ⊥ results in a product state with anti-aligned pseudospins at each unit-cell, which implies a many-body state with short-range correlations. Strong correlations are instead induced by the interplay between H 0 and H , resulting in two copies of a transverse-field Ising model with quasi long-range ferromagnetic ordering if U > 4|J | [31,32]. A more detailed derivation of this model may be found in the Supplemental Material [24].
The long-ranged case.-Let us now consider H lr = H 0 + H . The two spin states are effectively decoupled, so that we focus only on a single chain, for which the winding γ may be either 0 or 1. As shown in the Supplemental material, an appropriate rotation maps H lr to the Creutz model considered in Ref. [31,32]. Its phase diagram features an interaction-induced transition from topological (TOI) to trivial insulator (TRI), which is readily detectable by the disconnected part of the MCD -we present MPS simulations thereof in Fig. 3(a). Correspondingly, the time-averaged MCD shown in Fig. 3(b) goes smoothly to zero in the trivial region, while it saturates to a finite value in the topological region. Interestingly, strong U interactions can also induce an exotic interaction-induced phase transition to a ferromagnetic phase, where the Z 2 symmetry is spontaneously broken (SB). Indeed the excitation generates the appearance of an anti-aligned domain in the local A, B density, whose size grows linearly in time. In turn, this implies that the corresponding MCD grows quadratically in time -as it is clearly visible in our simulations in the "SB" region [see gray and black traces in Fig. 3 Experimental blueprint for H lr .-The intriguing phase diagram shown in Fig. 3 is readily accessible through experiments with ultracold atoms trapped in optical lattices. Indeed, in Fig. 4 we present a simple proposal yielding the long-ranged Hamiltonian H lr = H 0 + H , based on two-component fermions in a tilted lattice, exposed to RF and two-photon Raman transitions. In this implementation, spin (↓ / ↑) and sublattice (A/B) degrees of freedom are effectively identified. To perform the meaurement, the system should be prepared in the FIG. 4. Proposal to engineer the long-range Hamiltonian H lr = H0 + H with two-component ultracold fermions in a tilted optical lattice. A Zeeman term splits the degeneracy between the ↑ and ↓ states. A linear tilt is used to hinder particles from naturally hopping between lattice sites. The desired hopping processes are reintroduced by i) RF transitions between the ↑ / ↓ states (green), and ii) two-photon Raman transitions (yellow). To create a hole excitation, a blast beam tightly-focused on x = 0 transfers one particle to an internal state which is uncoupled to the Hamiltonian. model's ground state, i.e., at half-filling. Subsequently a localized excitation |e is created at position x = 0 and arbitrary sublattice. The easiest strategy here would be to create a quasi-hole, by "blasting" away one particle by means of a focused excitation. Successive readouts of the spin-resolved average positions . . e denotes an average over the excited state created by |e ) are used to compute the MCD, which in this case very simply coincides with the difference C(t) = ±(X ↓ − X ↑ ). Here ± depends on whether |e adds or removes a particle [24] Conclusions.-In this work, we showed that reading out the mean chiral displacement of a localized excitation probes the topological invariant of chiral 1D many-body systems and signals the presence of symmetry-broken phases. Furthermore, we have identified a simple technique to realize and detect an interaction-driven transition to symmetry-protected chiral order by means of a dynamical experiment with ultracold atoms. Intriguing perspectives of this work are a generalization of the MCD to higher dimensions or to other symmetry-protected topological phases, and an investigation of the effects of disorder, losses and temperature on top of interactions. the project PASQUANS, too. M.R. acknowledges support from the Alexander von Humboldt foundation and the kind hospitality of the BEC Center in Trento, where part of the manuscript writing was conducted. P.M. acknowledges supports by the "Ramón y Cajal" program, and by the Spanish MINECO (FIS2017-84114-C2-1-P). The MPS simulations were run on the Mogon cluster of the Johannes Gutenberg-Universität (made available by the CSM and AHRP), with a code based on a flexible Abelian Symmetric Tensor Networks Library, developed in collaboration with the group of S. Montangero at the University of Ulm (now moved to Padua). * pietro.massignan@upc.edu

Supplemental Material on
Detecting topology through dynamics in interacting fermionic wires

Andreas Haller, Pietro Massignan, and Matteo Rizzi
The additional material is organized as follows. In Section A we present a detailed derivation of the mean chiral displacement (MCD) in interacting systems. In particular, we prove that the disconnected (one-body) part is actually quantized to the many-body winding number γ for a half-filled (noninteracting) Fermi sea, while at arbitrary fillings it presents time-damped oscillations around the same invariant, so that it converges to γ in the long time limit. Incidentally, we show that this provides an exact real-space reformulation of the many-body winding number. In Section B we provide details on the derivation of the effective spin model presented in the main text. In Section C we elucidate on how the long-ranged Hamiltonian H lr presented in the text is directly connected to a particular limit of the Creutz model. We conclude by discussing the exact simulations of quadratic Hamiltonians in Section D, and the MPS simulations of interacting systems in Section E.
Let us begin by presenting some definitions, which we use throughout all sections. In this manuscript we focus specifically on chiral symmetric Hamiltonians, i.e., those which commute with an anti-unitary operator Γ, acting locally on the fermionic operators c x,τ,s as Γc x,τ,s Γ −1 = Γ τ s,τ s c † x,τ ,s . Here Γ is an even-dimensional unitary matrix squaring to the identity, Γ 2 = 1. As discussed in detail in Refs. [6, 7], their topological properties are captured by the many-body invariant where g(k) = G(k, ω = 0) is the zero-frequency component of the imaginary-time Green's function. Such an object is not straightforward to measure in common experiments, highlighting the importance of alternative formulations of the same object.
Our aim is to study the connection between the invariant γ in Eq. (S.1) and the mean chiral displacement (MCD), The MCD is therefore defined as the time-dependent expectation value of the chiral displacement operator ΓX, ΓX ≡ Moreover, chiral symmetry ensures that around half-filling it is possible to extract the same information (up to a global sign) by removing particles at time t = 0, instead of adding them. For example, we may write where the sum runs over localized hole excitations (in the specific case of H sr , these excitations may be generated for example by {c 0,A,↑ , c 0,A,↓ , c 0,B,↑ , c 0,B,↓ }).
The most general form of the many-body MCD is therefore given explicitly by with ± depending on whether the operators ψ † i create particles (+), or removes them (-), and the sum runs over a complete orthonormal set of creation/destruction operators acting on a given unit cell sufficiently far away from the boundaries of the chain (to avoid scattering from the edges for times t < T ).

A: DECOMPOSITION OF THE MEAN CHIRAL DISPLACEMENT THROUGH CONTRACTIONS
We use the following contractions to decompose the correlator of interest, i.e., the MCD in which ξ(t) denotes the many-body (connected) part of the MCD. The one-body (disconnected) part of the MCD is therefore defined by C d = C − ξ. Note that, for a noninteracting system, the contractions correspond to the application of Wick's theorem and, obviously, ξ = 0. In the above decomposition, we made use of the particle number conservation to neglect expectation values with unequal number of creation and annihilation operators. Furthermore, since t > 0, we can readily identify the single-particle (retarded) Green's functions iG τ s,τ s (x, t) ≡ c x,τ,s (t)c † 0,τ ,s (0) and c 0,τ ,s (0)c † x,τ,s (t) = c x,τ,s (t)c † 0,τ ,s (0) * = −iG * τ s,τ s (x, t). Using the definition of the chiral displacement and the matrix definitions of both G and Γ (their indices being spin and sublattice ones), we obtain: which is Eq.
(3) of the main text. It is important to notice that the first sum does not depend on time, because the average ΓX is intended over the state before the perturbation, which is assumed to be an eigenstate of H (either the true vacuum or the half-filled system, in our work here).
Let us now focus on the time-dependent one-body part of the correlator C, i.e. ∆C d (t) ≡ x x · tr G † ΓG (x, t). Using the Fourier transform G(x, t) = (2π) −1 dk e ikx G(k, t) = (2π) −2 dk dω e i(kx−ωt) G(k, ω) and performing an integration by parts (to substitute x → i∂ k ) we arrive at The identity ∂ k GG −1 = 0 allows one to write ∂ k G(k,ω) = −G(k,ω)∂ k G −1 (k,ω) G(k,ω). Upon commuting ΓG(k,ω) = −G(k, −ω)Γ one finds To proceed further, we consider now a non-interacting Fermi sea in a lattice with unit-cells containing two sites only (like the famous SSH model, for example). The lattice filling is controlled by the Fermi momentum k F , and the Green's function may be written as which allows for the simplification ∂ k G −1 (k,ω) = ∂ k H 0 (k) ≡ H 0 , that we use to good advantage in the following two subsections, where we explicitly perform the integrations in ω andω. Before proceeding with the calculation of ∆C d (t), let us recall that in an appropriate basis, the Hamiltonian may be cast in a completely off-diagonal form to anti-commute with the chiral operator Γ = σ z : with ε k = |h x (k)| 2 + |h y (k)| 2 , and n is a vector of unit norm. In this case g(k) = H −1 0 (k) = H 0 (k)/ε 2 k and thus the integral leading to the invariant in Eq. (S.1) simplifies to since |w k | = 1, so that ∂ k (w k w * k ) = 0.
A.1: Adding one particle on top of the vacuum When a particle is added to a completely empty lattice, as considered originally in Refs. [14][15][16], the chemical potential coincides with the bottom of the band, so that the Green's function reads simply G 0 (k, ω) = 1/(ω − H(k) + i0 + ), and all its poles are in the lower half of the complex plane. Making use of Cauchy's residue theorem and closing the integration contour in the lower (upper) semi-plane for G (G † ), one finds In the last step, we used [Γ, 1/H 0 (k)] = 0, then H 0 = 1/g, and finally g∂ k (g −1 ) = ∂ k ln(g −1 ) = −g −1 ∂ k g. Due to the presence of the time t inside the phase factor e i2H0(k)t , the integrand of the last term on the r.h.s. changes sign very rapidly when t → ∞, so that its integral is increasingly damped, and as time grows the whole expression converges smoothly to the chiral invariant γ introduced in Eq. (S.1). In terms of the matrix element w k an the energy eigenvalue ε k introduced in Eq. (S.12), the matrix expression Eq. (S.14) simplifies notably to the scalar result which shows more explicitly how the result is the chiral winding γ, plus a damped oscillatory term (guaranteed to be purely real). The latter result, derived here by means of a many-body formalism, coincides exactly with the one found in Refs. [14]. A lengthy but straightforward calculation (in analogy to the one presented in Ref. [14]) allows one to show that for a system with unit-cells containing two sites only C(t) = 2C 1/2 (t), where C 1/2 (t) ≡ ΓX ψ is the MCD computed over a single excitation, instead of over a collection of them as in Eq. (S.2) (provided of course the system has global spin-rotation invariance; else, a sum on spins should also be performed). This convenient feature was used repeatedly in the main text.
A.2: Adding one particle to a half-filled Fermi sea As the lattice is increasingly filled, the poles of G (G † ) which are below the Fermi momentum move to the upper (lower) plane, so that they exit the integration contours defined above and do not contribute to the final result. In particular, at half filling one obtains the very simple and appealing result ∆C d (t) = 1 4πi dk tr Γg −1 (k)∂ k [g(k)] = γ , (S.16) directly quantized to the chiral invariant (i.e., without any oscillatory term). The latter constitutes an analytical proof that the disconnected time-dependent part ∆C d = x x tr(G † ΓG)(x, t) of the MCD, evaluated on a half-filled non-interacting system, is the real-space and real-time equivalent of the many-body winding defined in Eq. (S.1). This feature was also explicitly demonstrated numerically in Fig. 1(c) of the main text. Furthermore, a lengthy but otherwise straightforward calculation shows that the oscillatory term in Eqs. (S.14)-(S.15) in the previous subsection is caused by inter-band scattering processes between conduction and valence bands which are Fermi-blocked at half-filling. Interactions re-enable such scattering processes, introducing again oscillatory contributions on top of the quantized winding number. with the Pauli matrices being σ i , i ∈ {x, y, z}, and densities read The bare kinetic term transforms to If we define the local pseudospin operators as it becomes clear that the term dictated by J commutes with all local densities S 0 x,s , whereas J-contributions exchange particles between adjacent links and thus do not commute with S 0 x,s . Interactions result (up to irrelevant chemical potentials and boundary terms) in We now assume a fully translational-invariant system at half filling, which implies setting the link density S 0 x,s = 1/2 at every link, and arrive thus at the effective model presented in Eq. (9) of the main text.

C: RELATION TO THE CREUTZ MODEL
The Creutz model was first introduced in Ref. [33], and more recently studied in detail in Ref. [31,32]. Its kinetic part is defined as operator in this basis takes the form Γ = σ y . It is more convenient to work in a basis in which Γ = σ z , which can be achieved by rotating the former coordinate system by an angle π/2 around the σ x axis. Then, σ y → σ z , σ x → σ x and σ z → −σ y , which leads us to a rotated Creutz model The interacting rotated Creutz model at χ = π. We represent hopping processes by arrows indicating the direction of movement with corresponding colored amplitude. Interactions are represented by colored boxes.
The new lattice for χ = π displayed in Fig. S2 is slightly simpler, in the sense that the horizontal hoppings have been replaced by complex spin-mixing terms. At the fine-tuned point m = 0 and t = g = χ/π it is possible to find a connection to the spinless SSH chain: By tilting (and stretching) the lattice of Fig. S2 at g = t, it becomes apparent that we are dealing with a spinless SSH chain in which ↑ and ↓ play the role of A and B sublattices, subject to a density-density interaction which acts within the unit-cells (see Fig. S3). Indeed, this fine-tuned point results in the nontrivial phase diagram studied in Ref. [31], which was later studied in the context of Wilson-Hubbard matter in high energy physics [32]. We thus arrive at the Hamiltonian H lr studied in the main text.

D: TIME-EVOLUTION OF QUADRATIC HAMILTONIANS
Here we focus on generic non-interacting systems, describe in detail the modelling of single-particle excitations and their time-evolution by exact diagonalization of quadratic Hamiltonians. Thus, we assume a generic Hamiltonian which can be cast into the form in which h ij is a L × L matrix and c j are fermionic annihilation operators at lattice position j. We choose to rotate into the eigenbasis of H 0 , i.e.
is the diagonalized Hamiltonian of eigenmodes d k = U kj c j with eigenenergies ε k , which we choose to sort by magnitude ε k ≤ ε k+1 . In this basis, the ground state is of the form (S.27) in which N denotes the total number of fermions. Applying the definition of the ground state and rotating the basis, we find the single-particle Green's function We want to study the dynamics of a generic excitation of holes (or particles). In particular, where we introduced an (unnormalized) vector Q ∈ R L which is modeling the annihilation (creation) of local modes c ( †) x . For example, a single excitation at position L/2 would imply a Q with one nonzero entry Q L/2 = 1. The normalization constant Ω can be explicitly calculated from ψ ± | ψ ± and evaluates to with the excitation vector P = U Q and the subset P + = {p, p > N } for the creation of particles, or P − = {p, p ≤ N } for the creation of holes. We are interested in the time-evolution of |ψ ± , which is given after application of the time-evolution operator If we take a momentum eigenstate of the Hamiltonian d † k |0 and apply A(t), we find the time-evolved state and we note that we may interpret time-evolution as a rotation of the diagonal basis according to t † k = A k d † k . The resulting set of operators t † k are creators of time-evolved eigenmodes at momentum k. Note, that the time-evolved excitations read with the time-evolved ground state Now, we have all ingredients to calculate the Green's function of the time-evolved excited state which can be used to compute any non-interacting expectation value (for example, its diagonal entries give the values of the local density). We note the static part G 0 ij , which yields a potential constant (in time) in the MCD. This constant can be amended by the quantity ∆C = C(t) − C(t = 0 − ), with C(t = 0 − ) the (static) MCD of the ground state right before the creation of the local excitation.
To test finite-size scaling of the MCD, in Fig. S4 we compare the evolution of excitations in finite chains with variable lengths to the fully translational invariant and infinite sized limits derived in Section A. As one may have naively expected, the local densities (and therefore the MCD results) are completely identical for times in which the excitation is spreading only through the bulk of the system.

E: DETAILS ON THE MPS SIMULATIONS
All simulations of interacting systems presented in the main text are performed on a system with L = 32 unit-cells, and for simplicity we took the coupling constants J, J to be real-valued, i.e. we considered BDI Hamiltonians. In case of H lr , the system is fully decoupled along the spin such that the unit-cell can be reduced to that of a spinless SSH-chain (e.g., we consider only the ↑-wire). We consider up to M = 256 Schmidt-values and simulate the timeevolution via the time-dependent variational principle (TDVP) [34], using a timestep of at most ∆t = 0.1J . To avoid boundary effects discussed here and in Section E, we restrict the time-window to t < 10J , such that the propagating wavepacket is never touching the boundaries.
In the following, we will estimate the entanglement growth in the underlying system. It is well-known that a critical system can be represented as a conformal field theory, and methods familiar from renormalization group theory provide an estimate for the asymptotic behavior of such a generic system. For instance, a global quench of a given (critical) Hamiltonian results in a continuum of quasiparticle excitations spreading through the system (with maximum velocity v), and the entanglement S A (t) between a region A( ) of the chain of length and the rest of the chain B(L − ) will scale as S A (t) ∝ t (for times smaller than the characteristic length t < t ≈ /2v) [35]. Local perturbations, instead, produce a bipartite entanglement which grows only as S A (t) ∝ log(t) for critical systems, which implies a saturation of S A (t) ≈ const. if the underlying system is gapped [36]. Moreover, if we denote by i the distance between the location of the excitation and the A|B bipartition boundary, we expect that the entanglement entropy will suddenly increase after a critical time proportional to i. In our experiment, we thus expect a rapid saturation of the entanglement entropy in time, making it a perfect candidate to target with tensor network schemes. We numerically confirm the above statements by the MPS simulations presented in Fig. S5. Away from the quantum phase transition lines, the model H lr is gapped such that MPS Ansätze are highly efficient: the figure of merit -the maximum truncated probability in the reduced density matrix of a bipartition -can be kept at machine precision for finite systems. This allows us to aim for much larger evolution times and system sizes at individual points in parameter space, which we present in Fig. S6. In particular, the MCD presents a very irregular behavior close to the phase boundaries, as was visible in Fig. 3 of the main text. Here however we show that the MCD becomes more and more regular when one considers larger and larger systems sizes and evolution times.
Although the phase diagram of H sr is simpler than H lr , the two-wire setup is computationally more expensivee.g., the ground state at half filling in the topological phase is four-fold degenerate, and the complexity of the gapless region is increased as well. We compare the numerical complexity to approximate the ground state in Fig. S7(a-b).
Combining the precise approximation of the ground state with the accurate time-evolution of a local excitation as (a/c) The bipartition is chosen such that = L/2 and the excitation is created at the interface between A|B. Panels (b/d) demonstrate that increasing the distance between the boundary and the initial location of the excitation (as indicated by the numbers 1, . . . , 5 within the panels) causes an approximately linear delay of the onset of the rapid growth of the entanglement.  Fig. 3(b-c)). The injection of the excitation induces a spontaneous symmetry breaking, which leads to the formation of the domain-wall visible in panel (c). As a consequence, the MCD displays the characteristic quadratic divergence in time shown in panel (d).
confirmed by Fig. S5, we conclude that an MPS approach yields negligible numerical inaccuracies for H lr , whereas errors need to be discussed further in case of H sr by means of a bond dimension scaling analysis. We present the corresponding MCD in Fig. S7(c): time-traces visibly differ at the tails of the graph. Nonetheless, we show in Fig. S7(d) that the time-average of the MCD scales very quickly to the desired precision. In conclusion, ∆C is quite resistant against numerical errors. All in all, this provides good confidence that the scheme and results presented in this work are resistant against usual unavoidable inaccuracies of the TN simulations [37].