# Entanglement echo and dynamical entanglement transitions

Kim Pöyhönen  Teemu Ojanen  Computational Physics Laboratory, Physics Unit,
Faculty of Engineering and Natural Sciences, Tampere University, P.O. Box 692,
FI-33014 Tampere, Finland  Helsinki Institute of Physics P.O. Box 64,
FI-00014, Finland

###### Abstract

We formulate dynamical phase transitions in subsystems embedded in larger
quantum systems. Introducing the entanglement echo as an overlap of the
initial and instantaneous entanglement ground states, we show its analytic
structure after a quench provides natural definition of dynamical phase
transitions in the subsystem. These transitions come in two varieties, the
entanglement-type transitions and the bulk-type Loschmidt transitions. The
entanglement-type transitions arise from periodic reorganization of quantum
correlations between the subsystem and its environment, manifesting in
instantaneous entanglement ground state degeneracies. Furthermore, the
entanglement echo distinguishes the direction of the quench, resolves
spatially distinct dynamical phase transitions for non-uniform quenches and
give rise to sharply-defined transitions for mixed initial states. We propose
an experimental probe to identify entanglement-type transitions through
temporal changes in subsystem fluctuations.

_Introduction–_ The rapidly growing field of dynamical quantum phase
transitions aims to uncover general principles in nonequilibrium many-body
dynamics and explore the parallels between dynamics and critical phenomena
Heyl  _et al._ (2013); Karrasch and Schuricht (2013); Andraschko and Sirker
(2014); Vajna and Dóra (2014, 2015); Budich and Heyl (2016); Zvyagin (2016);
Halimeh and Zauner-Stauber (2017). While there is no direct relation between
far-from-equilibrium dynamics and equilibrium phases of matter, recent efforts
have revealed a wealth of connections between them Heyl (2018). Moreover, in
the modern age of quantum simulation and synthetic designer systems,
theoretical predictions are directly stimulating new experimental directions
Jurcevic  _et al._ (2017); Zhang  _et al._ (2017); Fläschner  _et al._ (2018);
Guo  _et al._ (2019); Tian  _et al._ (2019); Wang  _et al._ (2019); Tian  _et
al._ (2020); Xu  _et al._ (2020). Unifying themes across various subfields,
dynamical phase transitions have an extraordinarily wide appeal in current
research.

In the present work, we develop the theory of dynamical phase transitions of
subsystems of larger many-body systems following a sudden quench,
schematically illustrated in Fig. 1 (a). In the diagnostics of quantum
correlations in many-body systems, the entanglement spectrum has become an
invaluable tool Li and Haldane (2008); Fidkowski (2010). Recently, it has also
found applications in far-from-equilibrium systems Torlai  _et al._ (2014);
Canovi  _et al._ (2014); Jafari  _et al._ (2019); De Nicola  _et al._ (2021);
Gong and Ueda (2018); Gong  _et al._ (2019); Lu and Yu (2019); Pastori  _et
al._ (2020); Surace  _et al._ (2020). Here we consider a bipartite system and
introduce the entanglement echo
$\mathcal{E}(t)=\langle\lambda_{0}(0)|\lambda_{0}(t)\rangle$ as an overlap of
the initial and instantaneous ground states $|\lambda_{0}\rangle$ of the
entanglement Hamiltonian of a subsystem. We show that a vanishing entanglement
echo at time $t_{c}$ provides a natural definition of dynamical phase
transitions in a subsystem. The entanglement echo contains essential
information on quantum correlations that is not captured by the much-studied
Loschmidt echo Heyl  _et al._ (2013) and signals novel observable properties.

Figure 1: (a): Dynamics of subsystem $A$ embedded in a larger system (seen
here for 1d and 2d geometries) display singular features after a sudden quench
through a critical point.(b): Entanglement echo $\mathcal{E}(t)$ resolves two
types of dynamical criticality, an entanglement-type transition (top) and a
bulk-type transition (bottom). (c): Entanglement echo can spatially resolve
several dynamical transitions for a single quench configuration (solid and
dashed lines). (d): Entanglement-type transitions can be probed by monitoring
temporal change in the subsystem observables such as number of particles.

By solving dynamical entanglement transitions in 1d and 2d topological lattice
models, we demonstrate substantial conceptual advances in the theory of
dynamical phase transitions. Most importantly i) the zeros of the entanglement
echo exhibit two types of dynamical criticality as depicted in Fig. 1 (b), the
usual Loschmidt-type bulk transitions and the entanglement-type transitions
which indicate periodic redistribution of quantum correlations between the
subsystems and have no closed system analogy ii) inhomogeneous systems or
macroscopically non-uniform quenches give rise to distinct robust spatially-
resolved dynamical phase transitions as illustrated in Fig. 1 (c) iii) the
entanglement-type transition can be probed by monitoring the temporal
behaviour of subsystem fluctuations (Fig. 1 (d)) which also gives rise to
oscillating entanglement entropy. In addition, the entanglement echo provides
a natural framework to study entanglement transitions of systems in mixed
states and non-unitary evolution.

_Entanglement echo–_ To formulate a subsystem’s dynamical phase transitions,
we consider a bipartioning of a time-evolving system into two subsystems $A$
and $B$. The properties of the subsystem $A$ is encoded in the reduced density
matrix
$\rho_{A}(t)=\mathrm{Tr}_{B}\rho(t)=\sum_{i}\lambda_{i}(t)|\lambda_{i}(t)\rangle\langle\lambda_{i}(t)|$,
obtained by tracing out the degrees of freedom corresponding to the subsystem
$B$ from the full density matrix $\rho$ (representing a pure or mixed state).
The reduced density matrix can be parametrized by the entanglement Hamiltonian
$H_{E}$ defined by $\rho_{A}=\frac{e^{-H_{E}}}{Z}$, where $Z$ ensures the
normalization $\mathrm{Tr}_{A}\,\rho_{A}=1$. The dominant contribution to
$\rho_{A}$ comes from the state $|\lambda_{0}\rangle$ with the largest
eigenvalue $\lambda_{0}$, corresponding to the ground state of the
entanglement Hamiltonian. The entanglement ground state calculated for a many-
body ground state typically encodes universal information about the phase,
such as topology and low-lying excitations. The significance of the
entanglement ground state points to its potential importance also in far-from-
equilibrium systems. Thus, we define the entanglement echo by

$\displaystyle\mathcal{E}(t)=\langle\lambda_{0}(0)|\lambda_{0}(t)\rangle$ (1)

which measures the overlap between the initial and instantaneous entanglement
ground states during temporal evolution. If the entanglement ground state is
degenerate in the thermodynamic limit, the echo can be defined as the overlap
with the degenerate subspaces.

As depicted in Fig. 1 (a), we consider quench protocols where at $t=0$ the
state of the whole system is prepared to a known initial state, such as the
ground state or finite-temperature state of a pre-quench Hamiltonian. Then,
the Hamiltonian of the system is instantaneously modified to the post-quench
form. Analogous to the Loschmidt echo
$\mathcal{L}(t)=\langle\Psi(0)|\Psi(t)\rangle$, the vanishing of which defines
dynamical phase transitions for the full system, we define dynamical phase
transitions for a subsystem in terms of the entanglement echo. We regard the
subsystem $A$ as undergoing a dynamical phase transition at time $t_{c}$ if
the entanglement echo vanishes $\mathcal{E}(t_{c})=0$. It is convenient to
define the entanglement rate function
$\Gamma(t)=-\ln|\mathcal{E}(t)|^{2}/\Omega_{A}$, where $\Omega_{A}$ is the
characteristic size of the subsystem $A$. Dynamical phase transitions are
clearly seen in the non-analytic behaviour of $\Gamma(t)$.

Here we describe two methods of calculating the entanglement echo. In the case
of pure initial states and unitary evolution, the state of the system can be
expanded
$|\Psi(t)\rangle=\sum_{\mu\nu}M_{\mu\nu}(t)|\psi_{\mu}^{A}\rangle|\psi_{\nu}^{B}\rangle$,
where $|\psi_{\mu}^{A}\rangle,|\psi_{\nu}^{B}\rangle$ form a complete basis of
each subsystem. The singular-value decomposition of matrix $M_{\mu\nu}(t)$
leads to the Schmidt decomposition of a state as

$\displaystyle|\Psi(t)\rangle=\sum_{i}\lambda_{i}(t)^{1/2}|\lambda_{i}(t)\rangle|\lambda_{i}^{B}(t)\rangle,$
(2)

where the sum contains at most min(dim $A$, dim $B$) terms Ekert and Knight
(1995). The entanglement ground state at time $t$ can be readily read off from
(2), allowing a direct evaluation of the entanglement echo (1).

For non-interacting fermions, the evaluation of the entanglement echo
simplifies. Pioneered by Peschel Peschel (2003); Peschel and Eisler (2009),
the entanglement spectrum for free fermions in a Gaussian state can be
obtained from the correlation matrix
$\mathcal{C}_{lm}^{\sigma\sigma^{\prime}}=\langle\hat{c}^{\dagger}_{l\sigma}\hat{c}_{m\sigma^{\prime}}\rangle$,
where fermion operators $\hat{c}_{m\sigma}$ annihilate particles with spin
$\sigma$ and $l,m$ label positions in the subsystem $A$. The eigenstates
$|\xi_{i}\rangle$ and eigenvalues $\xi_{i}\in[0,1]$ of the correlation matrix
can be regarded as the eigenstates and occupation probabilities of a single-
particle entanglement Hamiltonian. To evaluate the entanglement echo, we first
need to compute the dynamical correlation matrix
$\mathcal{C}_{lm}^{\sigma\sigma^{\prime}}(t)=\langle\hat{c}^{\dagger}_{l\sigma}(t)\hat{c}_{m\sigma^{\prime}}(t)\rangle$
and diagonalize it. The value $\xi=\frac{1}{2}$ marks the Fermi level of the
entanglement Hamiltonian, so the entanglement ground state
$|\lambda_{0}(t)\rangle$ is a Slater determinant constructed from the states
satisfying $\frac{1}{2}\leq\xi(t)\leq 1$. In the second-quantized notation, it
can be expressed as
$|\lambda_{0}(t)\rangle=\prod_{\xi_{i}(t)\geq\frac{1}{2}}\hat{c}_{\xi_{i}(t)}^{\dagger}|0\rangle$.
Then, the entanglement echo
$\mathcal{E}(t)=\langle\lambda_{0}(0)|\lambda_{0}(t)\rangle$ becomes

$\displaystyle\mathcal{E}(t)=\det{\langle\xi_{i}(0)|\xi_{j}(t)\rangle},$ (3)

where the single-particle states satisfy
$\xi_{i}(0),\xi_{j}(t)\geq\frac{1}{2}$. This formula applies to zero- as well
as to finite-temperature pre-quench states.

_Dynamical entanglement transitions in 1d and 2d–_ Now we demonstrate
dynamical entanglement phase transitions in solvable two-band Fermi systems.
Our analysis applies to arbitrary spatial dimensions but we focus of 1d and 2d
topological lattice models. We consider pre- and post quench Hamiltonians of
the form

$H^{i/f}=\sum_{l,m}\hat{\boldsymbol{c}}_{l}^{\dagger}\left[\boldsymbol{d}^{i/f}_{lm}\cdot\boldsymbol{\sigma}\right]\hat{\boldsymbol{c}}_{m},$
(4)

where $\boldsymbol{\sigma}=(I,\sigma_{x},\sigma_{y},\sigma_{z})$ is a vector
of Pauli matrices and the set of matrices
$\boldsymbol{d}^{i/f}=(d_{0}^{i/f},d_{x}^{i/f},d_{y}^{i/f},d_{z}^{i/f})$
determine the specific form of pre $(H^{i})$ and post $(H^{f})$ quench
Hamiltonians. The spinor operator
$\hat{\boldsymbol{c}}_{l}^{\dagger}=(\hat{c}^{\dagger}_{l\uparrow},\hat{c}^{\dagger}_{l\downarrow})$
creates fermions at site $l$. For translationally invariant systems and
quenches, matrices $\boldsymbol{d}^{i/f}$ become diagonal in a $n$-dimensional
quasimomentum space. With a minor modification, which extends spinors
$\hat{\boldsymbol{c}}$ to the Nambu space, model (6) also describes quenches
in topological superconductors and solvable spin chains Heyl  _et al._ (2013);
Najafi  _et al._ (2019); Zamani  _et al._ (2020) and spin liquids Schmitt and
Kehrein (2015). In Sec. I of the supplementary material (SI), we have derived
the expression for the dynamic correlation matrix for the model (6) for
spatially dependent parameters and quench protocols in zero and finite-
temperature initial states.

Figure 2: (a): Entanglement echo rate of 1d system (top) for quench $m=1.5\to
m=0.3$ and the system size $L=100$ (total) $L_{A}=30$ (subsystem). The jump
singularities, which arise from the crossing of the entanglement spectrum
(bottom), do not coincide with the cusps singularities of the Loschmidt echo
of the full system. (b): Same as (a) but for the opposite quench $m=0.3\to
m=1.5$. The two echos agree (apart from finite size effects which vanish in
the thermodynamic limit). (c): Momentum-resolved entanglement echo rate (top)
for Chern insulator quench $m=0.5\to m=-0.5$ and the system size $L=100$
(total) $L_{A}=30$ (subsystem). The jump singularities correspond to $k_{y}$
values for which the entanglement spectrum (bottom) exhibits temporal gap
closings at $t_{c}$. (d): Same as (c) but for the opposite quench $m=-0.5\to
m=0.5$.

We first consider a 1d topological insulator defined by
$\boldsymbol{d}(k)=(0,\sin k,0,m-\cos k)$. This model belongs to the Altland-
Zirnbauer class BDI Altland and Zirnbauer (1997); Schnyder  _et al._ (2008)
and exhibits a nontrivial phase for $|m|<1$ and trivial phase for $|m|>1$. In
Fig. 2 (a)-(b) we have illustrated the dynamical phase transitions when the
system is quenched through the critical point $m=1$. The entanglement echo
distinguishes whether the quench is performed from the trivial to the
topological phase or vice versa. In the former case, which we dub as an
entanglement-type transition, the entanglement echo displays periodic jump
discontinuities as seen in shown in Fig. 2 (a). The entanglement spectrum
reveals that the jumps arise from stroboscopic level crossing signalling an
instantaneous entanglement ground state degeneracy. In the latter case (Fig. 2
(b)), which we call a bulk-type transition, the entanglement echo exhibits
cusps and agrees with the Loschmidt rate function. As shown in Sec. V in the
SI, the stroboscopic entanglement ground state degeneracies persists also to
finite-temperature initial states. In contrast to the Loschmidt echo, which
does not offer a natural generalization with sharply-defined transitions at
finite temperatures Bhattacharya and Dutta (2017); Heyl and Budich (2017);
Sedlmayr  _et al._ (2018a); Abeling and Kehrein (2016), the entanglement-type
transitions remain well-defined. While the entanglement echo for the bulk-type
transitions reduce to the Loschmidt echo of the total system, as seen in Fig.
2 (a) and discussed in Sec. III of the SI, the analytic structure and the
critical times of entanglement-type transitions do not coincide with the
Loschmidt transitions. The entanglement echo quantifies a temporal
reorganization of quantum correlations between the two subsystems and captures
essential information not contained in the Loschmidt echo.

Two-dimensional systems exhibit similar bulk- and entanglement-type
transitions as 1d systems. By considering the geometry shown in Fig. 1 (a),
the subsystem $A$ can be chosen as a segment in the $x$ direction so that the
perpendicular momentum $k_{y}$ remains a good quantum number. The reduced
density matrix decouples to blocks labelled by $k_{y}$, and the entanglement
echo can be decomposed from the echoes of each block as
$\mathcal{E}(t)=\prod_{k_{y}}\mathcal{E}(k_{y},t)$. The $k_{y}$-resolved
partial echoes can be obtained by diagonalizing the momentum-resolved
correlation matrix derived in Sec. I of the SI. After obtaining the
eigenfunctions $|\xi_{i}(k_{y},t)\rangle$, the partial echos can be calculated
by applying Eq. (3). In fact, quench dynamics are conveniently analysed in
terms of partial echos $\mathcal{E}(k_{y},t)$. Here we consider Chern
insulators defined by $\boldsymbol{d}(k)=(0,\sin k_{x},\sin k_{y},m-\cos
k_{x}-\cos k_{y})$ which exhibits three distinct topological phases with Chern
numbers $C=-1$ (when $-2<m<0$), $C=1$ ($0<m<2$) and $C=0$ (when $|m|>2$). As
in 1d case, when the system is quenched through a critical point, the
entanglement echo shows non-analytic behaviour which depends on the direction
of the quench. In addition, the entanglement echo rate function may exhibit
either a cusp or jump singularity depending on $k_{y}$. This is illustrated in
Fig. 2 (c)-(d) for transitions between $C=\pm 1$ phases. At times when the
rate function shows a jump singularity for specific $k_{y}$ values, the
instantaneous entanglement spectrum exhibits temporal gap closing for the
corresponding $k_{y}$. Also, the momentum for which the gap closing takes
place changes when the direction of the quench is inverted. Thus, the
entanglement echo in both 1d and 2d systems reveal two distinct dynamical
phase transitions and, in contrast to the Loschmidt echo, makes a qualitative
distinction in which direction the critical point is crossed.

_Spatially-varying quenches–_ The Loschmidt echo characterizes dynamics of the
system as a whole and, as such, is incapable of providing spatially-resolved
information. However, the entanglement echo reveals novel dynamical
criticality in macroscopically inhomogeneous systems or spatially varying
quenches. In fact, a single quench can give rise to several spatially-resolved
dynamical phase transitions characterized by different time scales. In Fig. 3
we have illustrated a quench in 1d system, where the pre- and post-quench
configurations vary in space. The subsystem $A$ experiences a quench from a
trivial to topological phase while subsystem $B$ experiences the opposite
quench. The entanglement echo reveals that, indeed, the different parts of the
system exhibit distinct sharply-defined dynamical phase transitions. Not only
are their critical times different, but the non-analytic structure shows that
the transition in $A$ is of entanglement type and the transition in $B$ is of
bulk-type. Since it takes finite time for information to propagate through the
system Gong  _et al._ (2019), the short-time behaviour giving rise to early
dynamical phase transitions is sensitive only to local quench properties.
Thus, a system which exhibits several distinct equilibrium critical points can
display multiple spatially-resolved dynamical phase transitions in a single
quench.

Figure 3:  Top: Spatially-varying pre-quench (solid) and post-quench (dashed)
configurations for the 1d model with periodic boundary conditions. The initial
configuration consists of regions with $\mu_{1}=1.5$ and $\mu_{2}=0.3$, while
the post-quench parameters have $\mu_{3}=0.5$ and $\mu_{4}=1.7$. Bottom:
Resulting distinct transitions displayed in subsystems $A$ (red) and $B$
(blue).

_Observable consequences–_ The entanglement-type transitions arise from
instantaneous degeneracies of the entanglement ground state which persist to
finite-temperature initial states. It is natural to wonder what the observable
consequences of this are, especially in contrast to the Loschmidt-type
criticality. Far from equilibrium, the states in the entanglement spectrum are
not in simple correspondence with the physical edge modes, thus preventing the
most direct experimental probes. Here we devise a method to probe and
distinguish the entanglement-type transitions by monitoring the temporal
changes in subsystem fluctuations. The stroboscopic degeneracy of two
entanglement ground states is expected to lead to enhanced fluctuations for
observables which have different expectation values in the two states. Indeed,
we demonstrate this by considering the number of particles in the subsystem
$A$ in a setup depicted in Fig. 1 (d). The particle number operator is
$\hat{N}_{A}=\sum_{i\in
A,\alpha=\uparrow,\downarrow}\hat{c}_{i\alpha}^{\dagger}\hat{c}_{i\alpha}$,
where the summation is over the lattice sites in $A$ and spin. As shown in
Sec. IV of the SI, the time-dependent variance of the particle number is given
by

$\displaystyle\mathrm{Var}\,N_{A}(t)=\sum_{i}\left[\xi_{i}(t)-\xi_{i}^{2}(t)\right],$
(5)

where $\xi_{i}(t)$ are the eigenvalues of the correlation matrix. In a
translation-invariant system, the first term is a constant fixed by the
average density, however, the second term should reflect the pronounced
oscillations of midgap states characterizing the entanglement-type transitions
shown in Fig. 2 (a). As seen in Fig. 4 (a), the particle number variance
indeed oscillates with periodicity of the critical times. In addition to
oscillations, it shows a linear trend due to mixing of the two subsystems. The
onset time of the linear growth depends on the depth of a quench while the
oscillation period reflects the periodicity of critical times. The pronounced
oscillations, which are visible even for small subsystems down to $\sim 10$
sites, persist to finite temperature initial states and provide an
experimental signal that distinguish entanglement-type transitions from
Loschmidt transitions seen in Fig. 4 (b) and trivial quenches shown in Sec. V
in the SI.

Figure 4: Subsystem particle number variance and entanglement entropy (inset)
after a quench. (a) corresponds to the entanglement-type transition in Fig. 1
(a) and (b) to the bulk-type transition in Fig. 1 (b).

Moreover, as illustrated in Fig. 4, particle number fluctuations essentially
reflect the behaviour of the von Neumann entropy
$S(t)=-\sum_{i}\left[\xi_{i}\log\xi_{i}+(1-\xi_{i})\log(1-\xi_{i})\right]$
Zhang  _et al._ (2014). The difference in entropy oscillations Sedlmayr  _et
al._ (2018b) depending on the direction of the quench is naturally explained
by the existence of the two types of subsystem transitions discussed in our
work. Since the above discussed mechanism of the subsystem fluctuations follow
from the oscillating entanglement ground state degeneracy, it applies to
generic observables and entanglement-type dynamical phase transitions.

_Conclusion and outlook–_ In this work we formulated dynamical phase
transitions for a subsystem of a many-body system by introducing the
entanglement echo. The entanglement echo provides an appropriate
generalization of the Loschmidt echo, giving rise to several conceptual
advances as well as new observable predictions discussed in our work. In the
present work, we studied dynamical criticality resulting from a zero- and
finite-temperature initial state undergoing unitary time evolution after a
quench. Since the entanglement echo is formulated in terms of the reduced
density matrix, it can be straightforwardly employed to study mixed states and
non-unitary time evolution. In the future, it will be interesting to study
subsystem dynamics in systems subjected to measurements, the effects of
measurements on quench dynamics Kuo  _et al._ (2021) and possible measurement-
induced dynamical entanglement phase transitions Li  _et al._ (2018); Chan
_et al._ (2019); Skinner  _et al._ (2019).

_Acknowledgements–_ The authors acknowledge the Academy of Finland project
331094 for support.

## References

  * Heyl  _et al._ (2013) M. Heyl, A. Polkovnikov, and S. Kehrein, Phys. Rev. Lett. 110, 135704 (2013).
  * Karrasch and Schuricht (2013) C. Karrasch and D. Schuricht, Phys. Rev. B 87, 195104 (2013).
  * Andraschko and Sirker (2014) F. Andraschko and J. Sirker, Phys. Rev. B 89, 125120 (2014).
  * Vajna and Dóra (2014) S. Vajna and B. Dóra, Phys. Rev. B 89, 161105 (2014).
  * Vajna and Dóra (2015) S. Vajna and B. Dóra, Phys. Rev. B 91, 155127 (2015).
  * Budich and Heyl (2016) J. C. Budich and M. Heyl, Phys. Rev. B 93, 085416 (2016).
  * Zvyagin (2016) A. A. Zvyagin, J. Low Temp. Phys. 42, 971 (2016).
  * Halimeh and Zauner-Stauber (2017) J. C. Halimeh and V. Zauner-Stauber, Phys. Rev. B 96, 134427 (2017).
  * Heyl (2018) M. Heyl, Rep. Prog. Phys. 81, 054001 (2018).
  * Jurcevic  _et al._ (2017) P. Jurcevic, H. Shen, P. Hauke, C. Maier, T. Brydges, C. Hempel, B. P. Lanyon, M. Heyl, R. Blatt, and C. F. Roos, Phys. Rev. Lett. 119, 080501 (2017).
  * Zhang  _et al._ (2017) J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong, and C. Monroe, Nature 551, 601 (2017).
  * Fläschner  _et al._ (2018) N. Fläschner, D. Vogel, M. Tarnowski, B. S. Rem, D.-S. Lühmann, M. Heyl, J. C. Budich, L. Mathey, K. Sengstock, and C. Weitenberg, Nat. Phys. 14, 265 (2018).
  * Guo  _et al._ (2019) X.-Y. Guo, C. Yang, Y. Zeng, Y. Peng, H.-K. Li, H. Deng, Y.-R. Jin, S. Chen, D. Zheng, and H. Fan, Phys. Rev. Applied 11, 044080 (2019).
  * Tian  _et al._ (2019) T. Tian, Y. Ke, L. Zhang, S. Lin, Z. Shi, P. Huang, C. Lee, and J. Du, Phys. Rev. B 100, 024310 (2019).
  * Wang  _et al._ (2019) K. Wang, X. Qiu, L. Xiao, X. Zhan, Z. Bian, W. Yi, and P. Xue, Phys. Rev. Lett. 122, 020501 (2019).
  * Tian  _et al._ (2020) T. Tian, H.-X. Yang, L.-Y. Qiu, H.-Y. Liang, Y.-B. Yang, Y. Xu, and L.-M. Duan, Phys. Rev. Lett. 124, 043001 (2020).
  * Xu  _et al._ (2020) K. Xu, Z.-H. Sun, W. Liu, Y.-R. Zhang, H. Li, H. Dong, W. Ren, P. Zhang, F. Nori, D. Zheng, H. Fan, and H. Wang, Science Adv. 6, eaba4935 (2020).
  * Li and Haldane (2008) H. Li and F. D. M. Haldane, Phys. Rev. Lett. 101, 010504 (2008).
  * Fidkowski (2010) L. Fidkowski, Phys. Rev. Lett. 104, 130502 (2010).
  * Torlai  _et al._ (2014) G. Torlai, L. Tagliacozzo, and G. D. Chiara, Journal of Statistical Mechanics: Theory and Experiment 2014, P06001 (2014).
  * Canovi  _et al._ (2014) E. Canovi, E. Ercolessi, P. Naldesi, L. Taddia, and D. Vodola, Phys. Rev. B 89, 104303 (2014).
  * Jafari  _et al._ (2019) R. Jafari, H. Johannesson, A. Langari, and M. A. Martin-Delgado, Phys. Rev. B 99, 054302 (2019).
  * De Nicola  _et al._ (2021) S. De Nicola, A. A. Michailidis, and M. Serbyn, Phys. Rev. Lett. 126, 040602 (2021).
  * Gong and Ueda (2018) Z. Gong and M. Ueda, Phys. Rev. Lett. 121, 250601 (2018).
  * Gong  _et al._ (2019) Z. Gong, N. Kura, M. Sato, and M. Ueda, (2019), arXiv:1904.12464 [quant-ph] .
  * Lu and Yu (2019) S. Lu and J. Yu, Phys. Rev. A 99, 033621 (2019).
  * Pastori  _et al._ (2020) L. Pastori, S. Barbarino, and J. C. Budich, Phys. Rev. Research 2, 033259 (2020).
  * Surace  _et al._ (2020) J. Surace, L. Tagliacozzo, and E. Tonni, Phys. Rev. B 101, 241107 (2020).
  * Ekert and Knight (1995) A. Ekert and P. L. Knight, American Journal of Physics 63, 415 (1995), https://doi.org/10.1119/1.17904 .
  * Peschel (2003) I. Peschel, Journal of Physics A: Mathematical and General 36, L205 (2003).
  * Peschel and Eisler (2009) I. Peschel and V. Eisler, Journal of Physics A: Mathematical and Theoretical 42, 504003 (2009).
  * Najafi  _et al._ (2019) K. Najafi, M. A. Rajabpour, and J. Viti, J. Stat. Mech. 2019, 083102 (2019).
  * Zamani  _et al._ (2020) S. Zamani, R. Jafari, and A. Langari, Phys. Rev. B 102, 144306 (2020).
  * Schmitt and Kehrein (2015) M. Schmitt and S. Kehrein, Phys. Rev. B 92, 075114 (2015).
  * Altland and Zirnbauer (1997) A. Altland and M. R. Zirnbauer, Phys. Rev. B 55, 1142 (1997).
  * Schnyder  _et al._ (2008) A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Phys. Rev. B 78, 195125 (2008).
  * Bhattacharya and Dutta (2017) U. Bhattacharya and A. Dutta, Phys. Rev. B 96, 014302 (2017).
  * Heyl and Budich (2017) M. Heyl and J. C. Budich, Phys. Rev. B 96, 180304 (2017).
  * Sedlmayr  _et al._ (2018a) N. Sedlmayr, M. Fleischhauer, and J. Sirker, Phys. Rev. B 97, 045147 (2018a).
  * Abeling and Kehrein (2016) N. O. Abeling and S. Kehrein, Phys. Rev. B 93, 104302 (2016).
  * Zhang  _et al._ (2014) Y. F. Zhang, L. Sheng, R. Shen, R. Wang, and D. Y. Xing, Journal of Physics: Condensed Matter 26, 105502 (2014).
  * Sedlmayr  _et al._ (2018b) N. Sedlmayr, P. Jaeger, M. Maiti, and J. Sirker, Phys. Rev. B 97, 064304 (2018b).
  * Kuo  _et al._ (2021) W.-T. Kuo, D. Arovas, S. Vishveshwara, and Y.-Z. You, (2021), arXiv:2103.08068 [quant-ph] .
  * Li  _et al._ (2018) Y. Li, X. Chen, and M. P. A. Fisher, Phys. Rev. B 98, 205136 (2018).
  * Chan  _et al._ (2019) A. Chan, R. M. Nandkishore, M. Pretko, and G. Smith, Phys. Rev. B 99, 224307 (2019).
  * Skinner  _et al._ (2019) B. Skinner, J. Ruhman, and A. Nahum, Phys. Rev. X 9, 031009 (2019).

SUPPLEMENTAL INFORMATION to “Entanglement Echo and dynamical entanglement
transitions”

##  I Derivation of the dynamical correlation matrix

The free fermion entanglement properties can be obtained from the correlation
matrix of a subsystem. Here we provide the detailed calculation of the
dynamical correlation matrix for two-band Fermi systems employed in the main
text, with number of additional results for completeness. As a starting point,
we assume the system is described by the two-band Hamiltonian of Eq. (4) of
the main text,

$H^{i/f}=\sum_{l,m}\hat{\boldsymbol{c}}_{l}^{\dagger}\left[\vec{d}^{i/f}_{lm}\cdot\boldsymbol{\sigma}\right]\hat{\boldsymbol{c}}_{m},$
(6)

where $\vec{d}_{lm}^{i/f}$ is a position-dependent four-component set of
matrices parametrizing the Hamiltonian,
$\boldsymbol{\sigma}=\begin{pmatrix}\sigma_{0},\sigma_{x},\sigma_{y},\sigma_{z}\end{pmatrix}$
is a vector of Pauli matrices with $\sigma_{0}$ the $2\times 2$ unit matrix,
$\hat{c}_{l}=(c_{l,\uparrow},\ c_{l,\downarrow})$ is an annihilation spinor at
position $l$, and the superscript $i/f$ refers to pre-/post-quench form of the
Hamiltonian.

###  I.1 General case

In the following, we will use $\hat{c},\psi$ and $E$ to refer to annihilation
operators, wavefunction and energies of the post-quench Hamiltonian, and
equivalently $\hat{a}$, $\phi$ and $\epsilon$ for the pre-quench Hamiltonian.
The time-dependent correlation matrix defined as

$\mathcal{C}_{lm}^{\sigma\sigma^{\prime}}(t)=\langle\hat{c}^{\dagger}_{l\sigma}(t)\hat{c}_{m\sigma^{\prime}}(t)\rangle,$
(7)

where the expectation value is taken in the ground state or a finite
temperature state of the pre-quench Hamiltonian. Rewriting the operators in
the eigenbasis of the post-quench Hamiltonian
$c_{m\sigma}=\sum_{E}\psi_{m\sigma}^{E}c_{E}$ where their time-dependence
becomes simple yields

$\mathcal{C}_{lm}^{\sigma\sigma^{\prime}}(t)=\sum_{E,E^{\prime}}(e^{-iEt}\psi_{l\sigma}^{E})^{*}e^{-iE^{\prime}t}\psi_{m\sigma^{\prime}}^{E^{\prime}}\langle\hat{c}^{\dagger}_{E}\hat{c}_{E^{\prime}}\rangle$
(8)

In order to calculate the expectation value, we express the post-quench
operators in the pre-quench eigenbasis
$\hat{c}_{E}=\sum_{\epsilon}\langle\psi^{E}||\phi^{\epsilon}\rangle\hat{a}_{\epsilon}$
to obtain

$\displaystyle\mathcal{C}_{lm}^{\sigma\sigma^{\prime}}(t)$
$\displaystyle=\sum_{E,E^{\prime}}(e^{-iEt}\psi_{l\sigma}^{E})^{*}e^{-iE^{\prime}t}\psi_{m\sigma^{\prime}}^{E^{\prime}}\sum_{\epsilon,\epsilon^{\prime}}\langle\phi^{\epsilon}|\psi^{E}\rangle\langle\psi^{E^{\prime}}|\phi^{\epsilon^{\prime}}\rangle\langle\hat{a}^{\dagger}_{\epsilon}\hat{a}_{\epsilon^{\prime}}\rangle=\langle\Psi_{m\sigma^{\prime}}(t)|M|\Psi_{l\sigma}(t)\rangle$
(9)

where we have defined

$\displaystyle|\Psi_{n\sigma}(t)\rangle$
$\displaystyle=\sum_{E}(e^{-iEt}\psi_{n\sigma}^{E})^{*}\ket{\psi^{E}}$ (10)
$\displaystyle M$
$\displaystyle=\sum_{\epsilon}n_{F}(\epsilon)\ket{\phi^{\epsilon}}\bra{\phi^{\epsilon}}.$
(11)

In the above, $n_{F}(\epsilon)=\left[e^{\epsilon/(k_{B}T)}+1\right]^{-1}$ is
the Fermi-Dirac distribution.

###  I.2 Translation invariant case

In the translation invariant case, we can write the dynamical correlation
matrix completely analytically. To streamline notation, we present the
derivation for the 1d case and state the result for higher dimensions in the
end. If translation invariance is preserved throughout the quench, it is
convenient to begin by moving to $k$ space:

$\mathcal{C}_{lm}^{\sigma\sigma^{\prime}}(t)=\langle\hat{c}^{\dagger}_{l\sigma}(t)\hat{c}_{m\sigma^{\prime}}(t)\rangle=\frac{1}{L}\sum_{k}e^{-i(l-m)k}\langle\hat{c}^{\dagger}_{k\sigma}(t)\hat{c}_{k\sigma^{\prime}}(t)\rangle,$
(12)

where, as before, the expectation value is evaluated in the ground state or a
finite temperature state of the pre-quench Hamiltonian. The spatial indices
$l,m$ label positions in the reduced subsystem $A$ while the whole system is
assumed to have periodic boundary conditions.

We transform the operators to the basis of the post-quench Hamiltonian energy
eigenstates as per
$\hat{c}_{k\sigma}=u_{k\sigma+}\hat{c}_{\omega+}+u_{k\sigma-}\hat{c}_{\omega-}$
where $u_{k\sigma\pm}$ are the eigenfunctions of the post-quench Hamiltonian.
Again, the post-quench operators have simple exponential time dependence, so
we obtain

$\displaystyle\mathcal{C}_{lm}^{\sigma\sigma^{\prime}}(t)$
$\displaystyle=\frac{1}{L}\sum_{k}e^{-i(l-m)k}\langle(u_{k\sigma+}^{*}\hat{c}^{\dagger}_{\omega+}(t)+u_{k\sigma-}^{*}\hat{c}^{\dagger}_{\omega-}(t))(u_{k\sigma^{\prime}+}\hat{c}_{\omega+}(t)+u_{k\sigma^{\prime}-}\hat{c}_{\omega-}(t))\rangle$
$\displaystyle=\frac{1}{L}\sum_{k}e^{-ilk+imk}\left[u_{k\sigma+}^{*}u_{k\sigma^{\prime}+}\langle\hat{c}^{\dagger}_{\omega+}\hat{c}_{\omega+}\rangle+u_{k\sigma-}^{*}u_{k\sigma^{\prime}-}\langle\hat{c}^{\dagger}_{\omega-}\hat{c}_{\omega-}\rangle\right.$
$\displaystyle\left.\qquad+u_{k\sigma+}^{*}u_{k\sigma^{\prime}-}e^{i(\omega_{k}^{+}-\omega_{k}^{-})t}\langle\hat{c}^{\dagger}_{\omega+}\hat{c}_{\omega-}\rangle+u_{k\sigma-}^{*}u_{k\sigma^{\prime}+}e^{i(\omega_{k}^{-}-\omega_{k}^{+})t}\langle\hat{c}^{\dagger}_{\omega-}\hat{c}_{\omega+}\rangle\right],$
(13)

where $\omega^{\pm}_{k}$ are the eigenvalues of the two bands of the post-
quench Hamiltonian. The expectation values are to be taken over the pre-quench
state, while the operators here are in the post-quench state format. To
evaluate the expectation values we hence convert the operators to the
eigenbasis of the pre-quench Hamiltonian

$\hat{c}_{\omega\pm}=\langle u_{k\pm}|v_{k+}\rangle\hat{a}_{\epsilon+}+\langle
u_{k\pm}|v_{k-}\rangle\hat{a}_{\epsilon-}(0),$ (14)

where we label the pre-quench states with $v$ and $\epsilon$ analogously to
the $u$ and $\omega$ of the post-quench states. We then have

$\langle\hat{c}^{\dagger}_{\omega\pm}\hat{c}_{\omega\pm^{\prime}}\rangle=\langle
v_{k,-}|u_{k,\pm}\rangle\langle u_{k\pm^{\prime}}|v_{k,-}\rangle
n_{F}(\epsilon_{k}^{-})+\langle v_{k,+}|u_{k\pm}\rangle\langle
u_{k\pm^{\prime}}|v_{k,+}\rangle n_{F}(\omega_{k}^{+}),$ (15)

Hence, we have

$\displaystyle\mathcal{C}_{lm}^{\sigma\sigma^{\prime}}(t)=\frac{1}{L}\sum_{k}e^{-i(l-m)k}$
$\displaystyle\left[u_{k\sigma+}^{*}u_{k\sigma^{\prime}+}\left(|\langle
u_{k,+}|v_{k,-}\rangle|^{2}n_{F}(\epsilon_{k}^{-})+|\langle
u_{k,+}|v_{k,+}\rangle|^{2}n_{F}(\epsilon_{k}^{+})\right)\right.$
$\displaystyle+u_{k\sigma-}^{*}u_{k\sigma^{\prime}-}\left(|\langle
u_{k,-}|v_{k,-}\rangle|^{2}n_{F}(\epsilon_{k}^{-})+|\langle
u_{k,-}|v_{k,+}\rangle|^{2}n_{F}(\epsilon_{k}^{+})\right)$ $\displaystyle
u_{k\sigma+}^{*}u_{k\sigma^{\prime}-}e^{i\Delta\omega_{k}t}\left(\langle
v_{k,-}|u_{k,+}\rangle\langle u_{k,-}|v_{k,-}\rangle
n_{F}(\epsilon_{k}^{-})+\langle v_{k,+}|u_{k,+}\rangle\langle
u_{k,-}|v_{k,+}\rangle n_{F}(\epsilon_{k}^{+})\right)$
$\displaystyle\left.+u_{k\sigma-}^{*}u_{k\sigma^{\prime}+}e^{-i\Delta\omega_{k}t}\left(\langle
v_{k,-}|u_{k,-}\rangle\langle u_{k,+}|v_{k,-}\rangle
n_{F}(\epsilon_{k}^{-})+\langle v_{k,+}|u_{k,-}\rangle\langle
u_{k,+}|v_{k,+}\rangle n_{F}(\epsilon_{k}^{+})\right)\right].$

where $\Delta\omega\equiv\omega_{k}^{+}-\omega_{k}^{-}$. Now, we notice that
we can regroup things in terms of the vectors

$\displaystyle x_{k\sigma}(t)$ $\displaystyle\equiv
u_{k\sigma+}e^{-i\omega^{+}_{k}t}\langle
u_{k,+}|v_{k,-}\rangle+e^{-i\omega_{k}^{-}t}v_{k\sigma-}\langle
u_{k,-}|v_{k,-}\rangle$ (16) $\displaystyle y_{k\sigma}(t)$
$\displaystyle\equiv u_{k\sigma+}e^{-i\omega^{+}_{k}t}\langle
u_{k,+}|v_{k,+}\rangle+e^{-i\omega_{k}^{-}t}u_{k\sigma-}\langle
u_{k,-}|v_{k,+}\rangle,$ (17)

and finally obtain

$\mathcal{C}_{lm}^{\sigma\sigma^{\prime}}(t)=\frac{1}{L}\sum_{k}e^{-i(l-m)k}\left[x_{k\sigma}^{*}(t)x_{k\sigma^{\prime}}(t)n_{F}(\epsilon_{k}^{-})+y_{k\sigma}^{*}(t)y_{k\sigma^{\prime}}(t)n_{F}(\epsilon_{k}^{+})\right].$
(18)

The result can also be brought to a form more similar to the general case
solution above by defining

$|\Psi_{n\sigma,k}(t)\rangle=\sum_{\eta=\pm
1}e^{-ikn}\left(e^{-i\omega_{k}^{\eta}t}u_{k,\sigma,\eta}\right)^{*}\ket{u_{k,\eta}}\qquad\qquad
M_{k,k^{\prime}}=\delta_{k,k^{\prime}}\sum_{\eta=\pm
1}n_{F}(\epsilon_{k}^{\eta})\ket{v_{k,\eta}}\bra{v_{k^{\prime},\eta}},$ (19)

in terms of which we have

$\mathcal{C}_{lm}^{\sigma\sigma^{\prime}}(t)=\sum_{k,k^{\prime}}\langle\Psi_{m\sigma^{\prime},k}|M_{k,k^{\prime}}|\Psi_{l\sigma,k^{\prime}}(t)\rangle.$
(20)

The final result for 1d case can be generalized to arbitrary dimension
straightforward along the lines presented above. The result for $n$ spatial
dimensions is

$\mathcal{C}_{\mathbf{l}\mathbf{m}}^{\sigma\sigma^{\prime}}(t)=\frac{1}{\Omega}\sum_{\mathbf{k}}e^{-i(\mathbf{l}-\mathbf{m})\cdot\mathbf{k}}\left[x_{\mathbf{k}\sigma}^{*}x_{\mathbf{k}\sigma^{\prime}}n_{F}(\epsilon_{\mathbf{k}}^{-})+y_{\mathbf{k}\sigma}^{*}y_{\mathbf{k}\sigma^{\prime}}n_{F}(\epsilon_{\mathbf{k}}^{+})\right],$
(21)

where $\Omega$ is the volume of the full system and summation runs over
$n$-dimensional Brillouin zone.

The main advantage with this approach over Eq. (9) is the ability to obtain
explicit analytical expressions for eigenstates and eigenvalues in $k$-space.
While the correlation matrix in itself is gauge invariant, care must be taken
to ensure the gauge used in the eigenvectors is consistent. We will assume the
Hamiltonians (pre- and post-quench) in $k$-space are of the form

$H^{i/f}_{k}=d^{i/f}_{0}\sigma_{0}+\vec{d}^{i/f}\cdot\boldsymbol{\sigma}.$
(22)

where $i/f$ labels pre-quench (initial) and post-quench (final) parameters. If
it can be guaranteed that $d^{i/f}_{3}\neq-d^{i/f}\equiv-|\vec{d}^{i/f}|$ for
$k\in]0,2\pi[$, a suitable gauge choice is

$\ket{u_{k+}}=\frac{1}{\sqrt{2d^{f}(d^{f}+d^{f}_{3})}}\begin{pmatrix}d^{f}_{3}+d^{f}\\\
d^{f}_{1}+id^{f}_{2}\end{pmatrix}\qquad\ket{u_{k-}}=\frac{1}{\sqrt{2d^{f}(d^{f}+d^{f}_{3})}}\begin{pmatrix}d^{f}_{1}-id^{f}_{2}\\\
-(d^{f}_{3}+d^{f})\end{pmatrix},$ (23)

where, as previously, $\ket{u_{k\pm}}$ is a post-quench eigenstate; the same
holds for the pre-quench eigenstates $\ket{v_{k\pm}}$, but with $d^{i}$
instead. This gauge works e.g. for the 1d topological insulator
$\boldsymbol{d}(k)=(0,\sin k,0,m-\cos k)$ when $m>-1$. If, on the other hand,
it is known that $m<1$, we can instead use the gauge

$\ket{u_{k+}}=\frac{1}{\sqrt{2d^{f}(d^{f}-d^{f}_{3})}}\begin{pmatrix}d^{f}_{1}-id^{f}_{2}\\\
d^{f}-d^{f}_{3}\end{pmatrix}\qquad\ket{u_{k-}}=\frac{1}{\sqrt{2d^{f}(d^{f}-d^{f}_{3})}}\begin{pmatrix}d^{f}_{3}-d^{f}\\\
d^{f}_{1}+id^{f}_{2}\end{pmatrix}.$ (24)

For the Hamiltonians used in this work, these expressions are ill defined at
$k=0$, but an unambiguous limit exists and is be used instead.

###  I.3 Partial translation invariance

If the system is translationally invariant either post- or pre-quench but not
both, a mixed $k$-space approach can be used. Eq. (20) can be applied
directly, but by replacing the part that corresponds to the system which is
not translationally invariant with the Fourier transform of its realspace
solution. In other words, if the pre-quench system is not translation
invariant, we substitute

$M_{k,k^{\prime}}=\sum_{\epsilon}n_{F}(\epsilon)\ket{\phi^{\epsilon}_{k}}\bra{\phi^{\epsilon}_{k^{\prime}}},$
(25)

while if the post-quench system breaks the invariance, we substitute

$\Psi_{n,\sigma,k}(t)=\sum_{E}\left[e^{-iEt}\psi_{n\sigma}^{E}\right]^{*}\ket{\psi^{E}_{k}}.$
(26)

In this way, the analytical expressions derived in the previous subsection can
be used for the Hamiltonian that does preserve translation invariance.
Explicitly, using the notation of the previous two cases:

  1. 1.

If the pre-quench system breaks translational symmetry, but the post-quench
system does not, we have

$\mathcal{C}_{lm}^{\sigma\sigma^{\prime}}(t)=\frac{1}{L}\sum_{k,k^{\prime},\epsilon}e^{-ikl}x_{k,\sigma,\epsilon}^{*}(t)n_{F}(\epsilon)x_{k^{\prime},\sigma^{\prime},\epsilon}(t)e^{ik^{\prime}m},$
(27)

where

$x_{k,\sigma,\epsilon}(t)=e^{-i\omega^{+}_{k}t}u_{k,\sigma,+}\langle\boldsymbol{u}_{k,+}|\phi_{\epsilon,k}\rangle+e^{-i\omega^{-}_{k}t}u_{k,\sigma,-}\langle\boldsymbol{u}_{k,-}|\phi_{\epsilon,k}\rangle.$
(28)

  2. 2.

If the post-quench system breaks translational symmetry, while the pre-quench
system does not, we have

$\mathcal{C}_{lm}^{\sigma\sigma^{\prime}}(t)=\sum_{k}(\Psi^{k}_{m,\sigma^{\prime}}(t))^{\dagger}M_{k}\Psi^{k}_{l,\sigma}(t)$
(29)

with
$\Psi_{n,\sigma,k}(t)=\sum_{E}\left[e^{-iEt}\psi_{n\sigma}^{E}\right]^{*}\ket{\psi^{E}_{k}}$
and
$M_{k}=n_{F}(\epsilon_{k}^{-})\ket{v_{k-}}\bra{v_{k-}}+n_{F}(\epsilon_{k}^{+})\ket{v_{k+}}\bra{v_{k+}}$

###  I.4 Momentum-resolved dynamical correlation matrix

The entanglement spectrum and the entanglement echo for translation invariant
systems beyond 1d are conveniently analyzed as a function of momenta which are
conserved by the entanglement cut defining the subsystem. For example, on a
torus we can perform the subsystem partitioning in $x$ direction and maintain
translation symmetry in perpendicular $y$ direction. The corresponding
entanglement spectrum can be calculated from a correlation matrix when $k_{y}$
is a good quantum number. This $k_{y}$-resolved matrix is defined as
$\mathcal{C}_{lm}^{\sigma\sigma^{\prime}}(k_{y},t)=\langle\hat{c}^{\dagger}_{l,k_{y}\sigma}(t)\hat{c}_{m,k_{y}\sigma^{\prime}}(t)\rangle$
where $l,m$ label the reduced system coordinates in $x$ direction. Here the
partially transformed fermion operators $\hat{c}_{m,k_{y}\sigma}$ annihilate
particles at position $m$ with transverse momentum $k_{y}$ and spin $\sigma$.
Following the derivation of the translation invariant two-band systems above,
we obtain

$\displaystyle\mathcal{C}_{lm}^{\sigma\sigma^{\prime}}(k_{y},t)=\frac{1}{L_{x}}\sum_{k_{x}}e^{-i(l-m)k_{x}}\left[x_{k_{x}\sigma}^{*}(k_{y},t)x_{k_{x}\sigma^{\prime}}(k_{y},t)n_{F}(\epsilon_{k_{x}}^{-}(k_{y}))+y_{k_{x}\sigma}^{*}(k_{y},t)y_{k_{x}\sigma^{\prime}}(k_{y},t)n_{F}(\epsilon_{k_{x}}^{+}(k_{y}))\right]$
(30)

Here $L_{x}$ is the length of the full system in $x$ direction and the
summation runs over all quasimomenta $k_{x}$.

Diagonalizing the matrix (30) for a fixed $k_{y}$ provides eigenfunctions of
the momentum-resolved single-particle entanglement Hamiltonian. The momentum-
resolved entanglement echo is then obtained by the determinant formula Eq. (3)
in the main text. Repeating the calculation for each $k_{y}$ provides the full
single-particle entanglement spectrum and all the partial entanglement echos.

##  II Loschmidt echo

For the completeness, here we derive the formula for the Loschmidt echo for
general position-dependent quench. The Loschmidt echo is defined as
$M(t)=|\mathcal{L}(t)|^{2}$, where $\mathcal{L}(t)$ is the Loschmidt amplitude

$\displaystyle\mathcal{L}(t)=\langle\Psi_{0}|U(t)|\Psi_{0}\rangle=\langle
0|\hat{a}_{\epsilon_{1}}\ldots\hat{a}_{\epsilon_{N-1}}\hat{a}_{\epsilon_{N}}U(t)\tilde{\hat{a}}_{\epsilon_{N}}^{\dagger}\hat{a}_{\epsilon_{N-1}}^{\dagger}\ldots\hat{a}_{\epsilon_{1}}^{\dagger}|0\rangle,$
(31)

where $\tilde{\hat{c}}_{\epsilon_{i}}^{\dagger}$ creates a particle in an
eigenstate of the pre-quench Hamiltonian. We convert the operators to the
post-quench basis using
$\hat{a}_{\epsilon_{i}}(t)=\sum_{E_{i}}\langle\epsilon_{i}|E_{i}\rangle\hat{c}_{E_{i}}$,
and get

$\displaystyle\mathcal{L}(t)=\sum_{E_{i},E_{i}^{\prime}}$
$\displaystyle\langle\epsilon_{N}|E_{N}^{\prime}\rangle\langle
E_{N}|\epsilon_{N}\rangle
e^{-iE_{N}t}\times\ldots\times\langle\epsilon_{1}|E_{1}^{\prime}\rangle\langle
E_{1}|\epsilon_{1}\rangle e^{-iE_{1}t}$ $\displaystyle\times\langle
0|\hat{c}_{E^{\prime}_{1}}\ldots\hat{c}_{E^{\prime}_{N-1}}U(t)\hat{c}_{E^{\prime}_{N}}\hat{c}_{E_{N}}^{\dagger}\hat{c}_{E_{N-1}}^{\dagger}\ldots\hat{c}_{E_{1}}^{\dagger}|0\rangle.$
(32)

By inserting $1=U^{\dagger}(t)U(t)$ between each set of creation operators, we
obtain a product of terms of the form
$U(t)c^{\dagger}_{E_{i}}U^{\dagger}(t)=c^{\dagger}_{E_{i}}(-t)$, where the
time dependence can then be easily obtained to get

$\displaystyle\mathcal{L}(t)=$
$\displaystyle\sum_{E_{i},E_{i}^{\prime}}\langle\epsilon_{N}|E_{N}^{\prime}\rangle\langle
E_{N}|\epsilon_{N}\rangle
e^{-iE_{N}t}\times\ldots\times\langle\epsilon_{1}|E_{1}^{\prime}\rangle\langle
E_{1}|\epsilon_{1}\rangle e^{-iE_{1}t}$ $\displaystyle\times\langle
0|\hat{c}_{E^{\prime}_{1}}\ldots\hat{c}_{E^{\prime}_{N-1}}\hat{c}_{E^{\prime}_{N}}\hat{c}_{E_{N}}^{\dagger}\hat{c}_{E_{N-1}}^{\dagger}\ldots\hat{c}_{E_{1}}^{\dagger}|0\rangle.$
(33)

The operator expectation value vanishes unless the set $E_{n}^{\prime}$ is a
permutation of the set $E_{n}$ in which case it gives the sign of the
permutation. Thus, we can write the amplitude as a sum over the set $E_{i}$
and permutations $P_{i}$ as

$\displaystyle\mathcal{L}(t)=$
$\displaystyle\sum_{E_{i},P_{i}}\mathrm{sign}P_{i}\langle\epsilon_{N}|E_{P_{i}N}\rangle\langle
E_{N}|\epsilon_{N}\rangle
e^{-iE_{N}t}\times\ldots\times\langle\epsilon_{1}|E_{P_{i}1}\rangle\langle
E_{1}|\epsilon_{1}\rangle e^{-iE_{1}t},$ (34)

which remains unchanged if we switch the permutation in the ket vectors into
the bra vectors

$\displaystyle\mathcal{L}(t)=$
$\displaystyle\sum_{E_{i},P_{i}}\mathrm{sign}P_{i}\langle\epsilon_{P_{i}N}|E_{N}\rangle\langle
E_{N}|\epsilon_{N}\rangle
e^{-iE_{N}t}\times\ldots\times\langle\epsilon_{P_{i}1}|E_{1}\rangle\langle
E_{1}|\epsilon_{1}\rangle e^{-iE_{1}t}.$ (35)

Finally, this can be written as

$\displaystyle\mathcal{L}(t)=\det{\mathcal{A}},$ (36)

where

$\displaystyle\left[\mathcal{A}\right]_{ij}=\sum_{E}\langle\epsilon_{i}|E\rangle\langle
E|\epsilon_{j}\rangle e^{-iEt}.$ (37)

In the above, the label $E$ refers to the post-quench Hamiltonian while
$\epsilon$ refers to the pre-quench $H$. Typically rather than the echo
$M(t)=|\mathcal{L}(t)|^{2}$, one studies the Loschmidt rate function given by

$\lambda(t)=-\log{|\mathcal{L}(t)|^{2}}/\Omega,$ (38)

where $\Omega$ is the volume of the system. Hence, zeros of the Loschmidt echo
correspond to divergences in the rate function.

The general formula above simplifies for two-band systems. When the pre-quench
Hamiltonian has a filled lower-band $|\phi_{\epsilon_{i}}^{-}\rangle$ and the
post-quench Hamiltonian is translation invariant, we have

$\displaystyle\left[\mathcal{A}\right]_{ij}=\sum_{k}\left[\langle{\phi_{\epsilon_{i}}^{-}}|u_{k}^{+}\rangle\langle{u_{k}^{+}}|\phi_{\epsilon_{j}}^{-}\rangle
e^{i\omega_{k}^{-}t}+\langle{\phi_{\epsilon_{i}}^{-}}|u_{k}^{-}\rangle\langle{u_{k}^{-}}|\phi_{\epsilon_{j}}^{-}\rangle
e^{i\omega_{k}^{+}t}\right],$ (39)

where the summation is over the momentum eigenstates, and, again,
$|u_{k}^{\pm}\rangle$ are eigenstates of the post-quench Hamiltonian.
Furthermore, if also the initial state is translation invariant, the general
formula reduces to the well-known form

$\displaystyle\mathcal{L}(t)=\prod_{k}\left[\cos(d^{f}_{k}t)+i\sin(d^{f}_{k}t)\cos\theta_{k}\right],$
(40)

where $\theta_{k}$ is the polar angle between $\vec{d}^{i}$ and $\vec{d}^{f}$
on the Bloch sphere and the product is over $n$-dimensional Brillouin zone.

##  III Critical times of entanglement transitions

As noted in the main text, while the critical times of bulk-type transitions
coincide with those obtained from the Loschmidt amplitude, this is not the
case for entanglement-type transition. In fact, this deviation from the
Loschmidt critical times is a parameter-dependent phenomenon, as seen from
Fig. 5. For $m>0$ the entanglement transition occurs with shorter critical
time, with the opposite for negative $m$; at the special point $m=0$, the
entanglement critical time coincides with the Loschmidt value. While we have
kept the pre-quench $m$ constant in the figure, this occurs even if it is
varied e.g. so that each quench corresponds to a constant shift in $m$
instead. Generically, then, an entanglement-type topological transition will
have a different critical time than would the analogous Loschmidt transition.
However, the precise shift depend on the details of the quench; in fact, as
seen in Fig. 5(d), when disorder is included, the critical times of
entanglement transitions split due to the lifting of the translation symmetry
protected four-fold degeneracy of the midgap crossings to two different
crossings. Regardless, the discontinuities in the entanglement echo and the
crossings in the entanglement spectrum remain as telltale signs of an
entanglement-type transitions.

Figure 5: (a): Entanglement echo rate (top) and spectrum (bottom) of 1d system
for the quench $m=1.5\to m=0.3$ and system size $L=100$ (total) $L_{A}=30$
(subsystem). The dashed vertical lines indicate Loschmidt critical times. (b)
Same, but for the quench $m=1.5\to m=0.0$. (c) Same, but for the quench
$m=1.5\to m=-0.3$. (d) Same, but where $m$ on each site is picked from a
uniform distribution $m=1.5\pm 0.3$ pre-quench and $m=\pm 0.3$ post-quench.

##  IV particle number fluctuations

Here we will derive the formula for the time-dependent fluctuations of the
number of particles in the subsystem $A$ in terms of the single-particle
entanglement spectrum evaluated in a gaussian pre-quench state. The starting
point is the particle number operator of the subsystem $A$

$\displaystyle\hat{N}_{A}=\sum_{i\in
A,\sigma=\uparrow,\downarrow}\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma},$
(41)

with expectation value $\langle\hat{N}_{A}\rangle(t)=\sum_{i\in
A,\sigma=\uparrow,\downarrow}\text{Tr}_{A}[\rho_{A}(t)\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}]=\sum_{i\in
A,\sigma=\uparrow,\downarrow}\langle\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}\rangle=\mathrm{tr}\,{\mathcal{C}}(t)=\sum_{i}\xi_{i}(t)$.
This quantity is fixed by the total number of particles in a homogeneous
system. However, the subsystem particle number fluctuations display a more
interesting behaviour as they are sensitive to non-local correlations over the
whole subsystem. The variance of the particle number is obtained through

$\displaystyle\mathrm{Var}\,N_{A}(t)=\langle\hat{N}_{A}^{2}\rangle-\langle\hat{N}_{A}\rangle^{2}=\sum_{\begin{subarray}{c}i,j\in
A\\\
\sigma,\sigma^{\prime}=\uparrow,\downarrow\end{subarray}}\langle\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}\hat{c}_{j\sigma^{\prime}}^{\dagger}\hat{c}_{j\sigma^{\prime}}\rangle-\langle\hat{N}_{A}\rangle^{2}$
$\displaystyle=\sum_{\begin{subarray}{c}i,j\in A\\\
\sigma,\sigma^{\prime}=\uparrow,\downarrow\end{subarray}}\langle\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}\rangle\langle\hat{c}_{j\sigma^{\prime}}^{\dagger}\hat{c}_{j\sigma^{\prime}}\rangle+\langle\hat{c}_{i\sigma}^{\dagger}\hat{c}_{j\sigma^{\prime}}\rangle\langle\hat{c}_{i\sigma}\hat{c}_{j\sigma^{\prime}}^{\dagger}\rangle-\langle\hat{N}_{A}\rangle^{2}=\sum_{\begin{subarray}{c}i,j\in
A\\\
\alpha,\sigma^{\prime}=\uparrow,\downarrow\end{subarray}}\langle\hat{c}_{i\sigma}^{\dagger}\hat{c}_{j\sigma^{\prime}}\rangle\langle\hat{c}_{i\sigma}\hat{c}_{j\sigma^{\prime}}^{\dagger}\rangle,$
(42)

where Wick’s theorem was employed to factor the four-operator expectation
value. Thus,

$\displaystyle\mathrm{Var}\,N_{A}(t)$
$\displaystyle=\sum_{\begin{subarray}{c}i,j\in A\\\
\sigma,\sigma^{\prime}=\uparrow,\downarrow\end{subarray}}\langle\hat{c}_{i\sigma}^{\dagger}\hat{c}_{j\sigma^{\prime}}\rangle\langle\hat{c}_{i\sigma}\hat{c}_{j\sigma^{\prime}}^{\dagger}\rangle=\sum_{\begin{subarray}{c}i,j\in
A\\\
\sigma,\sigma^{\prime}=\uparrow,\downarrow\end{subarray}}\langle\hat{c}_{i\sigma}^{\dagger}\hat{c}_{j\sigma^{\prime}}\rangle\left(\delta_{ij}^{\sigma\sigma^{\prime}}-\langle\hat{c}_{j\sigma^{\prime}}^{\dagger}\hat{c}_{i\sigma}\rangle\right)$
$\displaystyle=\mathrm{tr}\,{\mathcal{C}}(t)-\mathrm{tr}\,\mathcal{C}^{2}(t)=\sum_{i}\left[\xi_{i}(t)-\xi_{i}^{2}(t)\right]$
(43)

The final expression (IV) is valid for zero- as well as finite-temperature
initial states and provides a direct connection between the dynamical
entanglement spectrum and experimentally accessible quantities.

##  V Entanglement-type dynamical phase transitions at finite temperature

In this section we will illustrate the behaviour of the entanglement echo for
finite-temperature initial states. In particular, it is shown that the
entanglement echo still predicts a sharp transition which has a clear physical
interpretation. In Fig. 6 (a) we have plotted the entanglement echo for a
quench displaying entanglement-type dynamical transitions and in Figs. 6
(b)-(d) the corresponding single-particle entanglement spectra. Notably, the
echo rate is essentially temperature independent, displaying jump
singularities coinciding with the instantaneous entanglement ground state
degeneracies arising from the level crossings at $\xi=1/2$. Increasing
temperature primarily serves to compress the edges of the entanglement
spectrum towards $\xi=1/2$. Hence, the crossing are robust and do not shift
when the temperature of the initial state is varied, preserving the
discontinuities in the echo rate.

Figure 6: (a): Entanglement echo rate of 1d system (top) for quench $m=1.5\to
m=0.3$ and the system size $L=100$ (total) $L_{A}=30$ (subsystem). (b)-(d):
Corresponding single particle entanglement spectra in the order of increasing
temperature with $k_{B}T=$ 0.25, 0.5 and 0.75 in units of $t$.

The fact that the non-analytic character of entanglement-type transitions
persists to finite temperatures has observable physical consequences for the
subsystem fluctuations, as discussed in the main text. The time-dependent
variance of the particle number of a subsystem, illustrated in Fig. 7, shows
how the oscillating entanglement ground state degeneracy translates to
pronounced oscillations in entanglement-type quenches shown in (a) and (b).
The frequency of the oscillations reflects the periodicity of the critical
times, while the onset of the linear trend depends on the particular post-
quench parameters. In the special case where the post-quench Hamiltonian is
tuned close to $m=0$, the oscillations persist a remarkably long time before
the linear increase sets in. The entanglement spectrum crossings correspond to
local maxima in the particle number variance, which provides an experimental
handle to extract critical times. Due to the four-fold degeneracy of crossings
in the model, the maximum variance at $T=0$ is quantized to 1, as can be
observed at the critical times before the onset of the approximately linear
growth. An increased initial temperature serves to shift the overall variance
upwards, which obscures the quantization; however, the overall oscillatory
behaviour, while somewhat lower in amplitude, can still clearly be observed.
Remarkably, this behaviour of the particle number fluctuations is essentially
the same as that of the von Neumann entropy, which is not easily accessed
experimentally.

Figure 7: Particle number fluctuations (top) and the entanglement entropy
(bottom). (a): Entanglement-type transition $m=2.0\to m=0.05$. The dashed line
indicates first two critical times. (b): Entanglement transition $m=2.6\to
m=0.4$ (c): bulk-type transition $m=-0.3\to m=1.6$ (d): trivial quench
$m=2.4\to m=1.2$.

In comparison to the entanglement type transitions, the particle number and
entropy oscillations are suppressed for bulk-type transitions and trivial
quenches, as illustrated in Fig. 7 (c) and (d). This difference could be
employed as an experimental probe to distinguish entanglement-type transitions
from bulk transition and trivial quenches.