Quenched topological boundary modes can persist in a trivial system

Topological boundary modes (TBM) can occur at the spatial interface between a topological and gapped trivial phase and exhibit a wavefunction that exponentially decays in the gap. Here we argue that this intuition fails for a temporal boundary between a prequench topological phase that possess TBM eigenstates and a postquench gapped trivial phase that does not possess any eigenstates in its gap. In particular, we find that characteristics of states (e.g., probability density) prepared in a topologically non-trivial system can persist long after it is quenched into a gapped trivial phase with spatial profiles that appear frozen over long times postquench. After this near-stationary window, TBM profiles decay slowly in a power-law fashion. This behavior highlights the unusual features of nonequilibrium protocols enabling quenches to extend and control topological states.

Topological boundary modes (TBM) can occur at the spatial interface between a topological and gapped trivial phase and exhibit a wavefunction that exponentially decays in the gap. Here we argue that this intuition fails for a temporal boundary between a prequench topological phase that possess TBM eigenstates and a postquench gapped trivial phase that does not possess any eigenstates in its gap. In particular, we find that characteristics of states (e.g., probability density) prepared in a topologically non-trivial system can persist long after it is quenched into a gapped trivial phase with spatial profiles that appear frozen over long times postquench. After this near-stationary window, TBM profiles decay slowly in a power-law fashion. This behavior highlights the unusual features of nonequilibrium protocols enabling quenches to extend and control topological states.
Here we argue that vestiges of spatially localized TBMs initially prepared in a gapped topologically non-trivial system persist long after it is quenched into a gapped trivial and uniform phase, see Fig. 1b. In particular, we find that even though TBMs are not an eigenstate of the postquench phase and exist inside the postquench bulk gap, characteristics such as a spatially peaked TBM probability density (PD) persist and appear frozen over a long and tunable time window postquench.
This near-stationary nature of the TBM postquench PD (Fig. 1b, Fig. 2b) is particularly surprising given the rapid evolution of the postquench wavefunction. Indeed, this contrasts with recurrent Loschmidt-echo type responses found in many quenched systems [25][26][27] wherein observables oscillate as a function of time. Further, even as translational symmetry in the postquench hamiltonian is maintained, the postquench PD remains spatially localized (Fig. 1b) retaining its prequench profile. It decays as a power law at very long times.
As we explain below, the persistence of PD arises from an interplay of TBM spatial localization and its fast dynamical evolution postquench. It is a striking example of how quantum systems pushed far out-of-equilibrium can possess very unusual properties with no analogue in equilibrium [28][29][30]. We anticipate TBM persistence can In all panels PD snapshots are plotted using Eq. (4) taken at logarithmic time spacings. We also used Φ(x) = sech(x/α) and N = (π √ 8α) −1 in Eq. (3).
Probability density and quench protocol -The persistence of TBMs can be most easily illustrated by a timedependent Dirac Hamiltonian (two-band) that undergoes a quench at time t = 0 ( Fig. 1, bottom): where v is the Dirac velocity, r = (x, y), p = (p x , p y ), σ i , i = 1, 2, 3 are Pauli matrices, and time-varying mass In the following, we shall refer to the pre-(post-) quench Hamiltonian asĤ i (Ĥ f ). For t < 0 before the quench, we have M (x) satisfying M (x < 0) < 0 and M (x > 0) > 0, such that it describes a mass domain wall along the yaxis (x = 0). Due to the jump in sgn M (x), the domain wall carries a unit topological charge, and thus supports gapless topological boundary modes (TBM) |ψ py that linearly disperse as v p y . Their spatial wavefunction is well localized at the domain wall x = 0: where N is a normalization constant, and Φ(x) decays exponentially on both sides. The decay length, α, defines the spatial extent of Φ(x).
Here we will consider a TBM state |ψ py prepared in a prequench bulk gap with |p y | < α −1 . At t = 0, the system is quenched via Eq. (2) and the mass parameter becomes uniform in space. As a result, |ψ py no longer exist as eigenstates of the postquench Hamiltonian and evolve as exp[−iĤ f t/ ]|ψ py . TBM |ψ py postquench characteristics at t > 0 are most saliently captured by its probability density (PD) ρ t,x = ψ py |e iĤ f t/ P r e −iĤ f t/ |ψ py ,P r = |r r|, (4) with explicit p y dependence suppressed hereafter for brevity. Postquench, the initial boundary-localized TBM state is no longer an eigenstate of the new spatialtranslation invariant HamiltonianĤ f ; instead its temporal evolution can be understood from the Larmor precession of its projected components in the postquench eigenbasis ofĤ f . Due to the tightly localized (in x) profile of the initial TBM, many postquench p x eigenstates are accessed. As a result, their destructive interference and multi-frequency Larmor precession can lead to dephasing [29] and decay of the initial TBM. Indeed, for mα/( v) 1, we find through direct numerical integration of Eq. (4), the PD generically decays exponentially in time, Fig. 2a.
However, when mα/( v) 1, we find that the PD no longer exhibits fast decay. Direct numerical integration of Eq. (4) yields an unusual regime wherein the initial TBM state freezes maintaining its original localized (in Log-linear plot of dynamical evolution of postquench probability density ρ t,x for small mα/( v) = 0.1. This exhibits a fast exponential decay as displayed by the linear slope quenching the PD rapidly. b. Log-Log plot of postquench probability density ρ t,x for large mα/( v) = 20. In contrast to (a), this exhibits a very slow decay and has a evolution characterized by two distinct regimes: (blue shaded) a slow frozen regime where probability density hardly evolves, and (orange shaded) a power law melting regime where the probability density decays as t −1 (indicated by dashed line). For both panels, curves from top to bottom (red to green) indicate various positions away from the domain wall x/α = 0 → 1.2 in steps of 0.15. PD was obtained from direct numerical integration of Eq. (4) using prequench parameters the same as Fig. 1.
x) profile ( Fig. 1b and Fig. 2b) at short times. As we explain below, this frozen TBM stays locked over times 0 ≤ t τ mα/ v controlled by the size of the postquench gap m; here the characteristic time τ = α/2v. At longer times, the TBM melts, displaying a much slower powerlaw decay. As we now explain, this arises from a near lock-step interference process between p x components enforced in the mα/( v) 1 regime. To more concretely understand the freezing and subsequent melting of the TBM mode, we analyze the solution of ρ t,x obtained through direct evaluation of Eq. (4): where S(x, t) = S xx + S yŷ + S zẑ , andŝ = Ψ (0) |σ|Ψ (0) The quantities C(x, t) and S(x, t) represent the even and odd temporal components of the postquench dynamical evolution of the PD and can be written as C( where j = 0, x, y, z, and the j-dependent weights as f 0 = 1 and (f x , f y , f z ) = d(p)/|d(p)|. Here ε p = 2 v 2 |p| 2 + m 2 is the (postquench) Dirac eigenenergy andΦ(p x ) is the Fourier transform of TBM profile Φ(x). In obtaining Eq. (5), we expanded exp[−iĤ f t/ ] and P r in Eq. (4) as superpositions of the postquench energy eigenstates, and applied the Pauli matrix identity σ i σ j = δ ij + i ijk σ k , see full details in Supplementary Information, SI.
Initially at t = 0, S(x, t = 0) = 0 and C(x, t = 0) reduces to the spatial profile Φ(x) of the topological boundary state yielding Eq. (5) that follows the profile of the TBM mode. For t > 0, however, both the arguments of S(x, t) and C(x, t) oscillate with time capturing the Larmor precession of the TBM projected components in the postquench eigenbasis. Indeed, given the wide range of p x eigenstates that the TBM state projects to, the temporal evolution of PD in Eq. (5) involves (momentumspace) nonlocal interference between multiple distinct p x momentum modes, and can lead to a complex spatiotemporal evolution of the PD, see SI.
When mα/( v) 1, the dynamical phase factor exp[iε p t/ ] in Eq. (6) rapidly oscillates as a function of p x . As a result, it suppresses the large momentum contributions to C(x, t) (and S(x, t)). Since the tight localization of the initial TBM depends on large p x momentum contributions, the initial state is rapidly eroded. On a physical level, this can be understood as an intrinsic dephasing phenomena where the multiple frequency oscillation get out of phase destroying the coherence of the initial state. In real space, this manifests as two wavetrains moving almost uniformly in opposite directions with a group velocity v g = v p x / ∆ 2 + 2 v 2 p 2 x ≈ ±v for small ∆ 2 = m 2 + v 2 2 p 2 y (see Fig. 1a) Frozen TBMs and Wick-rotated diffusion -However, when mα/( v) 1, the Larmor precession of the postquench state for various p x components oscillate in near lock-step, see below. Indeed, in this limit, the oscillation frequency In the latter, we focus on small p y < α −1 m. As the tight localization of the TBM state arises from large p x ∼ α −1 components, when mα/( v) 1, the differences in the frequencies are severely suppressed. This near-lock-step oscillation dramatically slows the erosion of the TBM state and as we will now discuss, at short times can even halt it. To see this explicitly, we take ε p ≈ ∆ + ( vp x ) 2 /(2∆) in Eq. (6) and directly integrate: with the postquench profile is an imaginary Gaussian kernel that results from the gapped dispersion in the large mα/( v) 1 limit, see SI. This kernel G t (x) is the mathematical embodiment of the random interference processes between different momentum sectors, with its imaginary Gaussian nature paralleling Brownian motion scattering in a diffusive medium; it is fundamentally distinct from the dynamical behavior expected when the system is quenched into a gapless hamiltonian, e.g., m = 0 above. We note that F x,y contributions are small (suppressed in the mα/( v) 1 limit) and do not contribute substantially to ρ t,x profile [40].
For short times 0 < t mα 2 /(2 v 2 ) = τ mα/( v), the imaginary gaussian G t (x) in Eqs. (7) and (8) exhibit spatial oscillations far more rapid than Φ(x + x ) varies in x. Here τ = α/2v is an intrinsic timescale. As a result, at these short times F 0,z ≈ N Φ(x)e i∆t/ , so that C(x, t) ≈ N Φ(x)cos(∆t/ ) and S z ≈ N Φ(x)sin(∆t/ ). Substituting into Eq. (5) we obtain the PD: We note PD in this regime is flat to infinite order in vt/α, Fig. 2b (blue shaded); any temporal decay in this frozen regime is non-analytic and slower than any power law. As a result, we term this a region of frozen PD. This is particularly surprising since the prequench TBM state is not an eigenstate of the postquench HamiltonianĤ f . Perhaps even more striking is the fact that the period over which the TBM's PD is frozen can be controlled by the postquench gap size m; increasing m gives a wider frozen window for TBM PD, Eq. (9). At long times t mα 2 /(2 v 2 ), the situation is dramatically different with the imaginary gaussian G t (x) in Eqs. (7) and (8) slowly varying over large ranges of x. As a result, the frozen states "melts" with its spatial profile spreading slowly out according to the convolution Ψ t (x) and an overall PD amplitude decaying as t −1 : This t −1 power-law melting decay (dashed line) conforms with that found from direct numerical integration in Fig. 2b (orange shaded), and even exhibits a data collapse for PD taken at different values of x. In this regime, the spatial extent/width of |Ψ t (x)| 2 grows as t thereby conserving the total probability density over all space. To understand this long-time behavior more intuitively, we note that G t (x) can be interpreted as a 1D heat kernel corresponding to a Wick-rotated "diffusion" process (i.e. evolution by Schrödinger's equation involving the id/dt instead of d/dt operator) with diffusion constant −i /2∆. While scattering processes in ordinary diffusion processes lead to a "Gaussian-blurred" distribution characterized by a (real) Gaussian decaying kernel, the imaginary Gaussian kernel G t (x) represents the rapid interference effects from multiple non-coherent contributions in Schrödinger evolution. A large ∆ results in slow spreading of the initial state, which by Eq. (8) is asymptotically governed by power-law decay, instead of Gaussian decay as in the more familiar real-time diffusion scattering processes. Physically, the slow melting TBM behavior can be likened to that of slow diffusion found in classical systems with large inertial masses; the large mα/( v) 1 limit corresponds to the regime where postquench modes are energetically inaccessible.
Pseudo-spin precession and creep current -The frozen-in-time PD profile described in Eq. (9) hides the fact that the TBM modes are not eigenstates of the postquench HamiltonianĤ f . Are all other observables similarly frozen in time as well? To further interrogate the TBM postquench, we consider its local (spatiallyresolved) pseudo-spin expectation with explicit p y dependence suppressed hereafter for notational brevity. σ t,x encodes the spinor-wavefunction information of the TBM. prequench, σ is aligned alonĝ s in theê 2 direction; hereê 1,2,3 denote the x, y, z directions in a Bloch sphere. However, postquench, the eigenstates ofĤ f generically possess spinor components in all three-directions. As a result, the dynamical evolution of σ t,x postquench involves a complex intertwining of precession and interference between wave components composing the spatially localized profile of the TBM state. Indeed, σ t,x possesses a dynamics that generically departs from that of the Bloch equation, see SI [40].
To see this, we directly evaluate Eq. (11) in the same fashion as Eq. (5) above producing the closed form [40] The first term tracks the persistence of the initial pseudospin directionŝ, while the second term represents a "pure" precession contribution; the other terms correspond to additional dynamical contributions from S, which arises physically from projection onto postquench eigenstates. Eq. (12) in fact describes the dynamical solution of any two-component state with an initial spatially inhomogeneous profile. In our case, Eq. (12) contains information on how the localization of them TBM interplays with precession effects. This generically yields a complex spatio-temporal and spin-dependent evolution, see explicit example in SI [40]. However, when mα/ v 1, the pseudo-spin expectation σ t,x of the TBM mode can act like a localized block spin with TBM pseudo-spin expectation at each x precessing in unison. Similar to the PD described above, this also leads to a frozen regime, where the magnitude of the | σ t,x | pseudo-spin is preserved for a sizeable period of time, Fig. 3a. Intuitively, in this limit, the d = (v p x , 0, m) precession axis points strongly towards in theê 3 direction, causing the initial pseudospinŝ to precess largely in theê 1 -ê 2 plane. Since the precession axis does not shift much across different p x and p y sectors due to large m, destructive interference is minimized and the quenched TBM states periodically return to their initial configurations: | σ t,x | is preserved. Indeed, in this regime, various positions in the TBM have σ x,y t,x that  Fig. 2. b. σx and σy oscillation at x = 0 (red) and x/α = 1.2 (brown). c Trajectory of σ t,x for x = 0 for mα/( v) = π exhibiting a decay in the pseudo-spin density over time. d The decay over each cycle leads to a creep current: a non-vanishing drift current of the TBM. Creep current estimated by averaging σx,y t,x over a period twice. oscillate in phase, Fig. 3b. At longer times t τ mα/ v, | σ t,x | melts, and decays as a power-law t −1 , Fig. 3a. See SI [40] where the power-law is derived.
Since −1 ∂Ĥ f /∂p = v(σ x , σ y ), the pseudo-spin oscillations indicate a cyclotron-type motion of the TBM mode; in the large mα/ v 1, the oscillations are largely locked to a frequency 2m/h. Interestingly, when the magnitude of the pseudo-spin decays, the value of σ x or σ y does not come back to itself after a full revolution, Fig. 3c. This decay in pseudo-spin leads to a creep current: an uncompensated drift current of the TBM (see estimate in Fig. 3d) that drifts along the domain wall.
The persistence of TBM characteristics can form the basis for a strategy to extend TBM features long after the topological system that supports it is gone: namely, by quenching to a large gap wherein postquench eigenstates are energetically inaccessible. Since the persistent TBM PD arises from the fast postquench Larmor precession enforced by large m, we expect that it is robust against inelastic scattering with energies far smaller than the gap scale as well as slowly varying disorder with typical lengths longer than the TBM width. The frozen and melting quenched PD regimes can be readily prepared and measured in ultra-cold atom optical lattices [18,19,21,31,32] with its pseudospin components of the wavepacket independently extracted [31,32], and its real-space dynamics tracked [19].
Perhaps most striking, is the contrast in behavior of TBMs across a sharp temporal boundary (as realized in the quench at t = 0) against those found at spatial boundaries. In the latter, boundary states (topological or otherwise) that exist in an energy gap typically exponentially decay over space into a gapped insulating bulk. Whereas in the former, TBM PD persists, freezing for a long window then decaying in a power-law fashion even when there are no eigenstates in trivial gap. This underscores the stark difference between sharp spatial and temporal interfaces and displays the dichotomy between equilibrium and out-of-equilibrium responses.

SUPPLEMENTARY INFORMATION FOR "QUENCHED TOPOLOGICAL BOUNDARY MODES CAN PERSIST IN A TRIVIAL SYSTEM" POSTQUENCH DYNAMICAL EVOLUTION
In this section, we review how the topological boundary mode (TBM) in Eq. (3) of the main text and its characteristics (namely probability density (PD) and pseudospom expectation) evolve postquench. For generality, we do not make any assumption on the form of the initial state (unless otherwise stated), other than it can be spatially inhomogeneous in the x-direction. We shall also assume the most general post-quench Hamiltonian that is translation-invariant in both directions. Hence our following results are applicable for a broader class of quench settings than those discussed in the main text. The time-evolved wavefunction postquench can be directly evaluated as ψ py (r, t) = r|e −iĤ f t |ψ py .
Writing exp(−iĤ f t) = p,± e ∓iεpt |ε ± p ε ± p |, and noting that |ε ± p ε ± p | = 1 2 I ±d(p) · σ |p p|, where r|p = e ip·r is the wavefunction of the momentum eigenstate, andd(p) = d(p)/|d(p)|, we obtain ψ py (r, t) = px,r e ip·(r−r ) (cos ε p t−id(p)·σ sin ε p t) r |ψ py (S-1) In obtaining Eq. (S-1) we summed across the ± bands as well as inserted the resolution of the identity r |r r |. Summing across r , p x by writing out r|ψ py = px NΦ(p x )e ipxx e ipyy |Ψ (0) we obtain the compact form where C and S are the same as the main text. We reproduce these here for the convenience of the reader In Eq. (S-2), we see that ψ(t, r) comprises two pseudospinor contributions, one proportional to cos ε p t in the direction of the original spinor |Ψ (0) , and another proportional to i sin ε p t for which |Ψ (0) has undergone â d(p) · σ spinor rotation.

Probability density dynamics
The evolution of the probability density postquench can be obtained from direct evaluation of Eq. (4) of the main text. This amounts to finding the square amplitude of the wavefunction above, |ψ py (r, t)| 2 . This can be evaluated as where we have suppressed the explicit x, t dependence of C, S for brevity. This expression can be readily simplified by recalling the vector identity: This yields the compact expression for PD as we obtain Eq. (5) in the main text. We note, parenthetically, when d(p) describes a Dirac model, further simplifications can be made to Eq. (S-7). For example, for Dirac models ε p is even in p x . This forces C to be real, and Re[S] (Im[S]) to arise solely from the even(odd) components ofd(p x , p y ). Hence Im[S] ê 1 and Re[S] ⊥ê 1 , and the time-evolved probability density reads as with the last term nonzero only for p y = 0.

Pseudospin dynamics
The evolution of the pseudospin expectation (density) postquench can be evaluated in the same fashion as above. Using the C, S notation above, the pseudospin expectation density in Eq. (11) of the main text can be re-written compactly as This expression can be readily simplified by repeated application of the identiy σ a σ b = δ ab I + i abc σ c and recalling the identity abc ade = δ bd δ ce − δ be δ cd . Here abc is the Levita-Cevita symbol, and a, b, c, d, e indices run over x, y, z. Focussing on the a-th component of the pseudospin expectation density, we obtain Re-writing in terms of real and imaginary parts of S and C, we obtain Eq. (12) in the main text. In obtaining Eq. (12) we have noted the identity (u + iv) × (u − iv) = −2iu × v, where u and v are vectors with real components. Similar to that discussed above, simplifications to Eq. (S-10) arise when ε p is even in p x . This forces C to be real. In that case, Re[S] (Im[S]) arises from the integral of the even(odd) components ofd(p x , p y ). Ifd(p) is furthermore purely even/odd, Eq. (S-10) simplifies to σ even t,x =ŝ(|C| 2 − |S| 2 ) + 2Re[S * C] ×ŝ + 2Re[(S * ·ŝ)S] (S-11) and (S-12) In either case, the last term of Eq. (S-10) always disappears, since it requires S to have both real and imaginary parts.
Recalling that the TBM modes discussed in the main text haveŝ =ê 2 , one also obtains the squared expected pseudospin magnitude (not to be confused with ρ t,r ) From the rotational invariance of Eq. (S-10), one can also easily deduce the expressions valid for general initial states withŝ =ê 2 .

FROZEN AND MELTING REGIMES
In this section, we provide a detailed description of how both frozen and melting regimes arise in the limit of large mα/ v 1. Before we proceed, we note C and S can be expressed as C(x, t) = [(F 0 (x, t) + F 0 (x, −t)]/2 and S x,y,z = [F x,y,z (x, t) − F x,y,z (x, −t)]/2i with where j = 0, x, y, z, and the j-dependent weights as f 0 = 1 and (f x , f y , f z ) = d(p)/|d(p)|. Since the dynamics of C(x, t) and S(x, t) are controlled by F j (x, t), we will analyze the temporal dynamics of F j (x, t) in the limit of mα/ v 1. In so doing, we first note that in the large postquench mass limit we can write where ∆ = (m 2 + 2 v 2 p 2 y ) 1/2 . We note parenthetically, when the postquench mass is large and by focussing on small p y < α −1 m/ v, we have ∆ ≈ m. The tight localization of the TBM state arises from large p x ∼ α −1 components. When mα/( v) 1, the differences in their precession frequency are very much suppressed and can be well approximated by Eq. (S-15). Applying Eq. (S-15) into the precession frequency of Eq. (S-14) we obtain where we have taken the principal value, and integrating out p x by completing the square in the argument of the exponential function, we obtain (S-19) where we have changed dummy variables x−x → x . We note that in obtaining the PD profile in both frozen and melting regimes, F x,y contributions are small (suppressed in the large mα/ v 1 limit) and are negligible in ρ x,t . This is because while the integrand in F z is proportional to m/∆ → 1, the integrand in F x , F y are proportional to p x , p y respectively, which are much smaller than the gap size m.

Frozen Regime
The imaginary Gaussian kernal in Eq. (S-19) plays a crucial role in the dynamical evolution of the probability density. At short times, G t (x) exhibit spatial oscillations far faster than the spatial variations of Φ(x). Indeed, Φ(x) varies significantly for x varying on order of the TBM width, α. In contrast, G t (x) changes rapidly on length scales of order 2t v 2 /∆, and for short times, this length scale can be far shorter than α. As a result, for short time windows (controllable by postquench m) we have ∆α 2 /(2t v 2 ) 1. In this short time window Eq. (S-20), G t → e −iπ/4 v √ 2π tδ(∆x), and Φ(x + x ) is effectively replaced by a constant Φ(x) times a factor containing t 1/2 , and we can approximate yielding F 0,z (x, t) = N e i∆t/ . As a result, at these short times F 0,z (x, t) ≈ N Φ(x)e i∆t/ , so that C(x, t) ≈ N Φ(x)cos(∆t/ ) and S z ≈ N Φ(x)sin(∆t/ ), where we have noted that F x,y are suppressed for large mα/ v. Substituting into Eq. (5) in the main text we obtain the Frozen profile This frozen (i.e. near-stationary in time) profile closely conforms to that found from a direct numerical evaluation of the probability density found in Fig.2b of the main text.

Melting regime
In the opposite regime to Eq. (S-20) (maintaining mα/ v 1), t ∆α 2 2 v 2 = mα v τ , the imaginary Gaussian kernal in Eq. (S-19) varies slowly in x. This will tend to spread F 0,z (x, t) and hence the density of states out in space. The amplitude of F 0,z (x, t) decays as t −1/2 as in Eq. (S-18). We similarly obtain C(x, t) ≈ [∆/(hv 2 t)] 1/2 N cos(∆t/ + π/4)Ψ t (x) and S z ≈ [∆/(hv 2 t)] 1/2 N sin(∆t/ + π/4)Ψ t (x). Note that S x , S y S z because the values of p x , p y that contribute to the integral giving F x , F y are much smaller than ∆ in the large mα/ v 1 limit. Plugging this into Eq. (5) in the main text we obtain, the power-law melting profile found in Eq. (10) of the main text.
Interestingly, in this regime the imaginary Gaussian kernal in Eq. (S-19) varies slowly in x far more slowly than Φ(x + x ) varies in position. As a result, for x α close to the domain wall Ψ t ≈ Φ(x + x ) = N −1 . In this limit, we find that the probability density at different x α collapse onto each other and decay as displaying a universal x-independent power-law decay. Indeed, this collapse of PD at different x positions in the long time limit is seen in a direct numerical evaluation in Fig. 2b of the main text. Indeed, we expect that at sufficiently long times, when Ψ t (x) has spread out, larger x > α within the width of Ψ t (x) will also similarly exhibit the universal x-independent power-law decay.

Illustration: exactly solvable example
While we have obtained approximate expressions for Ψ t (x) [and hence, F j (x, t)] above, it is instructive to illustrate the "frozen" and "melting" regime behavior through an exactly solvable model. To do so we consider the case where the initial TBM state Φ(x) takes the form of a Gaussian: Φ(x) = e −x 2 /2α 2 , where α is the standard deviation. We can then obtain by direct integrating Eq. (S-19) Substituting in Eq. (S-18), we obtain F j (x, t) ≈ g j e i∆t N exp(−x 2 /2α 2 ) ∝ Φ(x)e i∆t that leads to the frozen (near-stationary in time) PD profile as discussed above.
In the opposite limit t Recalling that √ 2πα 2 = N −1 is the inverse normalization constant for a Gaussian profile Φ(x) = e −x 2 /2α 2 , and substituting into Eq. (S-18), we similarly obtain the power-law decay in the PD described in Eq. (S-23) (in the limit x α). Note that the Gaussian profile above is special in that it does not spread due to its covariance under convolution by a Gaussian Kernel. For more general profiles, we will expect spread effects like what shown in Fig. 1 of the main text, as well as that shown in Figs. S-1 and S-2 of the Supplementary Information.

PRECESSION AND DYNAMICAL EQUATIONS OF MOTION
Although we have already provided the full solution to the PD and pseudospin expectation, it is still instructive to understand how the dynamical equations governing these quantities are related to the well-known spin precession equation. Due to the lack of translation invariance in the initial state, the evolution is not diagonal in momentum space, and the equations of motion will involve projectors. With |ψ py (t, r) = |r r|e −iĤ f t |Ψ (0) , the rate of change of the PD can be directly evaluated from Eq. (4) of the main text and takes the form whereP r = |r r| is projection operator onto the realspace state |r at position r,d = p d(p)|p p| is the projection operator onto the spinor eigenstate corresponding to d(p) and ... t denotes an expectation evaluated with the evolved initial state at time t (after the quench).
The projectionP enters because the expectation was evaluated at the particular position r, not across the entire wavefunction ψ py (t, r). Importantly, the simultaneous presence of operatorsP andĤ f , which are respectively diagonal in real-space and momentum-space, leads to a nontrivial imaginary part of the expectation Im P rĤf t , even though Im Ĥ f t = 0 due to the Hermiticity ofĤ f . Indeed the rate of change d ρ t,r dt at the specific position r can be construed as a measure of the "non-Hermiticity" of P rĤf , which is manifested through the imaginary part of its eigenvalues. This exhibits how the probability density at position r changes over time evolving to other positions.
Analogously, the rate of change of the pseudospin expectation is given by This is reminiscent of the usual precession equation with RHS given by 2d × σ, except that theP r projector now appear prominently. Comparing Eqs. S-27 and S-28, we see that the PD and pseudospin evolutions are both governed byP rd . Explicitly, the effect ofP r is to render this operator product non-local in momentum space: (S-29) Due to the spatial inhomogeneity of the real-space profile of the initial state, both p and p are not connected by a delta function, resulting in a more complicated expression involving both p and p . This can be interpreted as an interference mechanism between different momentum components, which gives rise to time dependence in the PD as well as the pseudospin expectation. As we discuss in the main text, this interference leads to a "Wickrotated" diffusion process, by which slow decay emerges in the melting regime of both ρ x,t and σ x,t arise.

Relation to Bloch dynamics for single spins
One can relate the above results to the much simpler well-known precession equation by defining a magnetization i.e. spectral pseudospin denstiy functionS(p, t) via x dx is the total pseudospin expectation for the entire system. With the spatial inhomogeneity integrated over, we indeed obtain an equation of motion that is local in momentum.

ILLUSTRATION: DYNAMICAL EVOLUTION OF PSEUDOSPIN EXPECTATION DENSITY
In this section we show the dynamical evolution of the full pseudospin expectation density which, as discussed in the main text, involves a complex intertwining of precession and interference between wave components composing the spatially localized profile of the TBM state. Concentrating on quenches into a Dirac-type Hamiltonian as discussed in the main text, we plot the dynamical evolution in Figs. S-1 and S-2 for special but important cases p y = 0, m = 0 and m = 0, p y = 0 respectively. When m = p y = 0 (exactly solved in the final Supplementary section), a state prepared as a TBM prequench with σ ∝ê 2 starts at x = 0 splits into two oppositely traveling wavepackets with σ ∝ ±ê 1 . This light-cone like spreading without attenuation is the result of the absence of any time scale/energy scale imposed by nonzero m or p y .
When p y = 0, the pseudospin expectation never acquires anyê 3 component. As the postquench gap m increases, attenuation occurs and the spins takes longer to relax to ±ê 1 , although the integrated PD over all space of course remains conserved at unity. As shown in Fig. S-1(b) and (c), the attenuation is asymmetric, and the initial decay on one side can be even slower than the case of m = 0. However, when mα/ v increases above unity, the decay becomes rapid while the group velocity and spin relaxation becomes very slow. Finally, at large mα/ v 1, the wavepacket becomes very spatially localized and behaves almost like a single Bloch spinor with negligible spatial inhomogeneity [ Fig. S-1(e)]. This is also the regime where it exhibits frozen and then melting behavior about the initial domain wall x = 0, as discussed extensively in the main text, as well as in the following section. This is akin to an imaginary damped precessing pseudospin.
When m = 0 instead, the attenuation occurs symmetrically and the evolved pseudospin cants out of the plane with nonzeroê 3 component [ Fig. S-2]. This occurs because the precession axis no longer all lie within the same plane across different momentum sectors.
Interestingly, when p y α increases above unity [ Fig. S-2], the decay becomes very rapid like in the m = 0 case. At large p y α, the pseudospin evolution proceeds similarly as in the large m, p y = 0 case, exhibiting a frozen and then melting behavior. This is because such behavior results fundamentally from large α∆/ v = α m 2 / 2 v 2 + p 2 y , which can result from either large mα/ v or large p y α (or both).
FIG. S-1: Time evolution of σ t,x for mα/ v = 0, py = 0 from an initial x = 0 spike, with profiles at vt/α = 1, 2, 3 being increasingly spread out or attenuated. The pseudospin direction is color coded according to the color wheel in (f). For the py = 0 case at hand, the evolution is restricted to theê1-ê2 plane. While it initially spreads out unattenuated for small mα/ v, for larger mα/ v attenuation occurs and the spins takes longer to relax to ±ê1. Finally for large mα/ v 1, the system exhibits frozen followed by melting behavior, with diminishing group velocity and attenuation.
FIG. S-2: Time evolution of σ t,x for pyα = 0, mα/ v = 0 from an initial x = 0 spike, with profiles at vt/α = 1, 2, 3 being increasingly spread out or attenuated. The main differences with the mα/ v = 0, py = 0 case (Fig. S-1 above) is that the attenuation occurs symmetrically and the spin evolution cants out of the plane with nonzeroê3 component, indicated by darker shades (see color legend in Fig. S-1(f)). At large pyα, the system also exhibits frozen followed by melting behavior.

PSEUDOSPIN EXPECTATION DENSITY IN THE MELTING REGIME FOR ARBITRARY py
In the following, we analyze the pseudospin expectation density in the melting regime when mα/ v 1. While in earlier sections we chose p y small, in this section we relax the assumption that |p y | < α −1 , in order to extend our scope to arbitrary electronic fillings. As we shall see, power-law behavior shall still emerge asymptotically even when p y is comparable to ∆. We shall temporarily work in units of α, , v = 1 for notational brevity.

Pseudospin Expectation σ t,x
We now explicitly derive the pseudospin expectation behavior in the regime ∆ = m 2 + p 2 y α −1 , specifically to leading order in ∆ −1 . Firstly, note thatd(p x ) can be approximated aŝ where ip x has been replaced by a derivative operator on the profile Φ(x). For simplicity of notation in the following derivations, we shall introduce -32) and the corresponding convolution with Φ (x):

(S-33)
This gives us where we have defined γ(x, t) = ∆t + η(x) and γ (x, t) = ∆t + η (x), and Θ t (x) = |Θ t (x)|e iη(x) and Θ t (x) = |Θ t (x)|e iη (x) , i.e. tan η(x) = Im[Θ t (x)]/Re[Θ t (x)]. Then upon subtitution in Eq. S-10, and noting that Im[S] = 0, we get whereŝ =ê 2 is the original pseudospin polarization of the state prepared in the TBM. In final expression, we see that the component in the direction ofŝ originates not just from the convolution Θ t (x) involving the initial state profile, but also Θ t (x) which involves the convolution with the spatial derivative of the initial state profile. This is the full solution of the dynamical evolution in the ∆ α −1 limit, with the first term of the final line indicating the extent of pseudospin expectation persistence.
To make further headway, we shall define which contains information about the nontrivialty of the imaginary Gaussian kernel. Firstly, in the frozen regime of small t, the kernel tends to a delta function, and Θ t (x) becomes genuinely the derivative of Θ t (x). Hence for small t, m t (x)| t≈0 → − d dx log e − x M (x )dx = M (x), the initial mass gap distribution. Hence m t (x) can be regarded as a "time-relaxed" mass gap distribution which equilibrates to zero after some time. But for large t, the definition Eq. S-36 permits no further simplification due to imaginary Gaussian Kernel. Nevertheless, the large t limit can still be reduced as follows: If we recall the assumption |p y | < α −1 i.e. that the occupied states are in the pre-quench gap, ∆ α −1 also implies that p y ∆, m ≈ ∆ and m t (x) < α −1 ∆. Hence Eq. S-37 further simplifies to which is the precession behavior in the melting regime as discussed in the main text. To summarize the above derivations, we first identified the large ∆ condition for slow power-law decay of the key quantities C and S. Next, the full expression for σ t,x is derived in this limit (Eq. S-35). The large t limit is then taken, resulting in a simplification to Eq. S-37. Finally, the assumption |p y | < α −1 on the initial state occupancy is invoked to result in the even simpler expression Eq. S-38, as borne out in the melting regime of Fig. 3 of the main text. In this final section, we provide pedagogical examples of the analytically tractable limit of m = p y = 0. It was shown in Fig. S-1(a) that in this special limit, the initial wavepacket does not decay at all, but instead "splits" into two and depart in equal and opposite velocities. Here we shall show how this behavior can be derived.
The post-quench Hamiltonian p x σ x is gapless with eigenenergy ε px = |p x |, and all its bulk states are energetically accessible from any initial TBM state. Its exact linear dispersion also admits an exact analytic solution for the special choice of initial TBM state given by our ansatzΦ(p x ) = πα sech παpx 2 . The post-quench state evolution can be computed from Eq. S-1 (still working with units v = 1): py=0 (x, t) = N px,py,± δ py,0 e ip·x e ∓i p tΦ (p x ) I ±d(p) · σ 2ŝ On the first line, we have projected the initial spinor into the upper and lower bands, where time evolution is given by phase rotations. This is then simplified in terms of a sum of p x , which will be converted into an integral below: where X = 2 πα x, T = 2 πα t are the spatial and temporal displacements in units of π 2 α. The integral is, with this fortuitous choice of profileΦ(p x ), analytically tractable with contour integration.
At long t, we observe manifest exponential decay behavior ψ m=0 py=0 (x, t) ∼ e −t/α with a timescale of α. This scale arises from the spatial decay of the initial mass profile m t (x), which should govern the post-quench dynamics because there is no other scale introduced. In the massive (m = 0) post-quench case, m will also determine the post-quench dynamics. Unlike in the special large m case, the initial state generically decays exponentially according to timescales set by the energy scales involved.
One can directly compute the pseudospin expectation of this m = p y = 0 case: takes a nice symmetric form containing the "light cone" rays t ± x, which govern the speed of propagation of the evolved wavefunction. At long time t, however, we observe that ϑ is still zero at the original location x = 0 of the TBM state, although the bulk states should have γ = ± π 2 . One can also directly obtain the pseudospin expectation amplitude | σ t,x | = sech 2 t − x α + sech 2 t + x α (S-43) which is manifestly a superposition of contributions from two oppositely traveling pulse trains at x ± t. This result (Eq. S-43) can also be obtained be obtained through Eq. S-13 withd(p x ) = (p x , 0, 0)/|p x |. Alternatively, one may also notice thatd(p x ) is odd, and proceed via Eq. S-12. Since S is hence purely imaginary, σ t,x lies in thê e 1 -ê 2 plane. Withŝ ∝ê 2 ⊥ S, only the first two terms of Eq. S-12 survives, and its evaluation is simple.