$\mathbb{I}\otimes\mathbb{I}$$\tau^{z}\otimes\tau^{z}$$\mathbb{I}\otimes\mathbb{I}$$\sigma^{z}\otimes\mathbb{I}$$\tau^{x}\otimes\mathbb{I}$$\tau^{y}\otimes\mathbb{I}$$\tau^{x}\otimes\tau^{x}$$\tau^{y}\otimes\tau^{y}$$\tau^{z}\otimes\tau^{x}$$\tau^{z}\otimes\tau^{y}$$\tau^{z}\otimes\mathbb{I}$$\mathbb{I}\otimes\tau^{z}$$\tau^{x}\otimes\tau^{y}$$\tau^{y}\otimes\tau^{x}$$\tau^{x}\otimes\tau^{z}$$\tau^{y}\otimes\tau^{z}$$\mathbb{I}\otimes\tau^{x}$$\mathbb{I}\otimes\tau^{y}$$\mathfrak{A}_{A}^{+(\sigma)}\\!\otimes\mathbb{I}$$\mathfrak{A}_{\tilde{A}}^{(\tau)}$$\tilde{\mathfrak{A}}_{\tilde{A}}^{+(\tau)}$$\mathcal{D}^{-1}_{\mathbb{Z}_{2}}(\ell)$
Figure 11: A graphical representation of the mapping (78) when the subsystem
$A$ consists of a single site. Operators $P_{1,2,3}$ are colored and
underlined. They anticommute with the operators in the correspondingly colored
circles: their effect is to remove the undesired terms in the density matrix
(80).

To sum up, we have started with enlarging the subsystem by one site
$\rho_{A}(\Psi)\rightarrow\rho_{A}(\Psi)\otimes\frac{\mathbb{I}}{2}$; this has
allowed us to apply the duality transformation (71), which alters the meaning
of spatial interval. We have then found a representation of
$\rho_{A}(\Psi)\otimes\frac{\mathbb{I}}{2}$ in terms of reduced density
matrices in the semilocal theory with the new notion of spatial subsystem. The
final object, $\bar{\rho}^{(\tau)}_{\tilde{A}}(\Psi)$ is an average of density
matrices selecting the part that is symmetric under three spin-flip
transformations. Finally, the entanglement entropies can be expressed in terms
of $\bar{\rho}^{(\tau)}_{\tilde{A}}(\Psi)$ by applying Eq. (72):

$S_{\alpha}(\ell,\Psi)=\frac{\log\mathrm{tr}[(\bar{\rho}^{(\tau)}_{\tilde{A}}(\Psi))^{\alpha}]}{1-\alpha}-\log
2\,,$ (82)

where $\log 2$ cancels the effect of the spurious double degeneracy introduced
when the subsystem was enlarged.

##### Bounds on the entropies of the dual RDMs.

We show here that, up to $\mathcal{O}(1)$ corrections, the averaged density
matrix $\bar{\rho}_{\tilde{A}}^{(\tau)}(\Psi)$ can be replaced by
$\rho_{\tilde{A}}^{(\tau)}(\Psi)$ in (82). Since no special property of the
density matrix and of the $\mathbb{Z}_{2}$ symmetry we average over is
required, we consider a generic density matrix $\rho$ (Hermitian positive
semidefinite operator of trace one) and the symmetry generated by a Hermitian
involution $P$; Eq. (81) is just a subsequent application of such an average.

We define the symmetric density matrix

$\bar{\rho}=\frac{\rho+P\rho P}{2},$ (83)

and we call $\bar{S}_{\alpha}$ its Rényi entropies ($S_{\alpha}$ indicates
instead the entropies of $\rho$).

Since $\bar{\rho}$ can be interpreted as the sum of two density matrices
weighted with the probability distribution $\\{1/2,1/2\\}$, general properties
of the von Neumann entropy (see Section $12.2$ of Ref. [102]) imply

$S_{1}\leq\bar{S}_{1}\leq S_{1}+\log 2\,$ (84)

For instance, the lower bound follows directly from the concavity of the von
Neumann entropy. We obtain an analogous result for the Rényi entropies, which
satisfy

$S_{\alpha}\leq\bar{S}_{\alpha}\leq S_{\alpha}+\frac{\alpha\log
2}{\alpha-1}\,;$ (85)

see Appendix E.1. Returning to Eq. (82), these bounds imply

$S_{\alpha}(\ell,\Psi)=S_{\alpha}^{(\tau)}(\ell+1,\Psi)+\mathcal{O}(\ell^{0})\,,$
(86)

where we have defined the entropy of the subsystem $\tilde{A}$ in the dual
representation

$S^{(\tau)}_{\alpha}(\ell+1,\Psi)=\frac{\log\mathrm{tr}[(\rho^{(\tau)}_{\tilde{A}}(\Psi))^{\alpha}]}{1-\alpha}.$
(87)

###  V.2 Excess of entropy

##### The protocol.

We consider a system invariant under $\mathcal{P}^{z}$ and compare two quench
protocols:

  1. (I)

nonequilibrium time evolution of a translationally invariant state
$\ket{\Psi(t)}=e^{-i\boldsymbol{\bf H}t}\ket{\Psi(0)}$;

  2. (II)

nonequilibrium time evolution of the same state as in I perturbed by a
local(ised) unitary operator $\ket{\Psi_{r}(t)}=e^{-i\boldsymbol{\bf
H}t}\boldsymbol{\bf\sigma}_{r}^{x}\ket{\Psi(0)}$.

The question is: do the Rényi entropies of subsystems discriminate between the
two protocols in the limit of infinite time? Inspired by Refs. [42, 43], we
define the excess of entropy as

$\Delta_{r}S_{\alpha}(\ell,t)=S_{\alpha}(\ell,\Psi_{r}(t))-S_{\alpha}(\ell,\Psi(t))\,,$
(88)

where the subsystem of length $\ell$ is associated with the sites
$\\{1,\ldots,\ell\\}$, that is to say, the reduced density matrix to be
investigated reads

$\rho_{A}(\Psi)=\frac{1}{2^{\ell}}\kern-10.00002pt\sum_{\\{\alpha\\}_{\ell}\atop\alpha_{j}=0,x,y,z}\kern-10.00002pt\braket{\Psi}{\prod_{j=1}^{\ell}\boldsymbol{\bf\sigma}_{j}^{\alpha_{j}}}{\Psi}\,\bigotimes_{j=1}^{\ell}\sigma^{\alpha_{j}}.$
(89)

Here $\Psi$ can be either $\Psi_{r}(t)$ or $\Psi(t)$.

In Section IV.3 we have discussed the limit of infinite time in the quench
protocol I, showing that, in the quasilocal theory, the state locally relaxes
to a projected $\mathbb{Z}_{2}$-semilocal generalised Gibbs ensemble. In order
to compute the excess of entropy, we should also take the limit of infinite
time in the protocol II, which goes however beyond our assumption of
translational invariance. Although a theory to describe the infinite-time
limit in inhomogeneous systems with semilocal charges is still missing, this
problem can be solved in the dual representation, where the spin flip is
mapped into a topological excitation resulting in a domain-wall initial state
[31]. Specifically, if the initial state is $\ket{\Uparrow}_{\sigma}$, it
follows from the first equation of (38) that after flipping the $r$-th spin
the state in the dual representation becomes
$\ket{\cdots\downarrow_{r-1}\downarrow_{r}\uparrow_{r+1}\uparrow_{r+2}\cdots}_{\tau}$
(or
$\ket{\cdots\uparrow_{r-1}\uparrow_{r}\downarrow_{r+1}\downarrow_{r+2}\cdots}_{\tau}$,
depending on how the symmetry is broken). The problem then moves to
identifying the stationary state describing the limit of infinite time after a
quench from a domain-wall state in a canonical theory, which is a much more
studied situation [103, 104, 105, 106, 107, 108]. The drawback is that we need
to express the observable under investigation in the dual representation; for
the Rényi entropies this corresponds to using Eq. (82).

The main effect of the presence of semilocal charges (namely, odd charges in
the dual representation) is that the stationary states emerging in the
protocols I and II are macroscopically different 111111The existence of odd
charges implies that the stationary state after the quench from
$\ket{\Uparrow}_{\tau}$ is macroscopically different from that after the
quench from $\ket{\Downarrow}_{\tau}$; by symmetry, the nonequilibrium
stationary state emerging from
$\ket{\cdots\downarrow_{r-1}\downarrow_{r}\uparrow_{r+1}\uparrow_{r+2}\cdots}_{\tau}$
has to be macroscopically different from the two of them. As far as we can
see, however, the two macrostates have generally the same entropy per unit
length, therefore the excess of entropy is expected to be subextensive. In the
examples that we consider, the excess of entropy still grows with $\ell$,
therefore we can also replace Eq. (82) by Eq. (86) without affecting the
asymptotic behaviour. This important simplification allows us to express the
excess of entropy in the dual representation as

$\lim_{t\rightarrow\infty}\Delta_{r}S_{\alpha}(\ell,t)=S^{(\tau)}_{\alpha}(\ell,\boldsymbol{\bf\rho}^{(\tau)}_{\mathrm{NESS}})-S^{(\tau)}_{\alpha}(\ell,\boldsymbol{\bf\rho}^{(\tau)}_{\infty})+O(1)\,,$
(90)

where

$\displaystyle\boldsymbol{\bf\rho}^{(\tau)}_{\mathrm{NESS}}=$
$\displaystyle\lim_{B\rightarrow\infty}\lim_{t\rightarrow\infty}\mathrm{tr}_{\overline{B}}[e^{-i\boldsymbol{\bf
H}^{(\tau)}t}(\ket{\Downarrow\Uparrow}\\!\bra{\Downarrow\Uparrow})_{\tau}e^{i\boldsymbol{\bf
H}^{(\tau)}t}]$ (91) $\displaystyle\boldsymbol{\bf\rho}^{(\tau)}_{\infty}=$
$\displaystyle\lim_{B\rightarrow\infty}\lim_{t\rightarrow\infty}\mathrm{tr}_{\overline{B}}[e^{-i\boldsymbol{\bf
H}^{(\tau)}t}(\ket{\Uparrow}\\!\bra{\Uparrow})_{\tau}e^{i\boldsymbol{\bf
H}^{(\tau)}t}]\,.$

In order to gain some intuition about the behaviour of the excess of entropy,
we follow the suggestion of Ref. [31] and consider the simplified case in
which $\boldsymbol{\bf H}^{(\tau)}$ commutes with $\boldsymbol{\bf S}^{z}$
(and hence $\ket{\Uparrow}_{\tau}$ is an eigenstate). We also assume that
$\ket{\Uparrow}_{\tau}$ could be interpreted as the vacuum of stable
excitations, whose number is linear in $\boldsymbol{\bf S}^{z}$ (this
characterises a class of integrable systems). In this situation the
nonequilibrium stationary state
$\boldsymbol{\bf\rho}^{(\tau)}_{\mathrm{NESS}}$ can be interpreted as a Fermi
sea or split Fermi seas [110], which can be described by a conformal field
theory (see also Ref. [111]). Since in a conformal field theory the Rényi
entropies have a logarithmic dependence on the subsystem’s length [112], we
can expect the excess of entropy to grow logarithmically in $\ell$. Although
not a-priori evident (but somehow consistent with the observations of Ref.
[31]), the log behaviour survives a symmetry-preserving local unitary
transformation of the initial state, which transforms this protocol into a
genuine global quench. This will be derived in the next section, in which we
will investigate systems that are dual to $\mathbb{Z}_{2}$-symmetric
generalised XY models [113].

####  V.2.1 Example: Dual generalised XY model

In this section we work out the excess of entropy in a generalisation of the
dual XY model that preserves the noninteracting structure and the existence of
a family of semilocal conservation laws. We do not focus only on the dual XY
model because, as will be clear later, the dual XY model is rather special and
one could read the outcome of the investigation as an indication of a
universality that instead is not present. In its most general form the
Hamiltonian reads

$\boldsymbol{\bf
H}=\sum_{\alpha\in\\{x,y\\}}\sum_{n=1}^{\infty}J_{\alpha}^{(n)}\sum_{j\in\mathbb{Z}}\boldsymbol{\bf
h}_{j}^{(\alpha,n)}\,,$ (92)

where the coupling constants $J_{\alpha}^{(n)}$ are assumed to approach zero
in the limit $n\rightarrow\infty$ exponentially fast in $n$. The operators
$\boldsymbol{\bf h}_{j}^{(\alpha,n)}$ are Hermitian and local, with support
including $j$ and extending over $n+1$ sites; their explicit form in the
$\sigma$ representation is not essential for the discussion but we report it
for the sake of completeness:

$\displaystyle\boldsymbol{\bf h}^{(x,1)}_{j}=$
$\displaystyle\boldsymbol{\bf\sigma}_{j-1}^{x}\boldsymbol{\bf\sigma}_{j+1}^{x}$
$\displaystyle=\boldsymbol{\bf\tau}_{j}^{x}\boldsymbol{\bf\tau}_{j+1}^{x}$
(93) $\displaystyle\boldsymbol{\bf h}^{(y,1)}_{j}=$
$\displaystyle-\boldsymbol{\bf\sigma}_{j-1}^{x}\boldsymbol{\bf\sigma}_{j}^{z}\boldsymbol{\bf\sigma}_{j+1}^{x}$
$\displaystyle=\boldsymbol{\bf\tau}_{j}^{y}\boldsymbol{\bf\tau}_{j+1}^{y}$
$\displaystyle\boldsymbol{\bf h}^{(\alpha,n)}_{j}=$
$\displaystyle\boldsymbol{\bf\sigma}_{j-1}^{x}\boldsymbol{\bf\sigma}_{j}^{\alpha}\left.\begin{cases}\prod_{l=1}^{\tfrac{n+1}{2}-\alpha}\boldsymbol{\bf\sigma}^{z}_{j+2(l-1)+\alpha}\boldsymbol{\bf\sigma}_{j+n-1}^{\alpha}&n\text{
odd}\\\
\prod_{l=1}^{\tfrac{n}{2}-1}\boldsymbol{\bf\sigma}^{z}_{j+2(l-1)+\alpha}\boldsymbol{\bf\sigma}_{j+n-1}^{3-\alpha}&n\text{
even}\end{cases}\right\\}\boldsymbol{\bf\sigma}_{j+n}^{x}$
$\displaystyle=\boldsymbol{\bf\tau}_{j}^{\alpha}\prod_{l=1}^{n-1}\boldsymbol{\bf\tau}_{l}^{z}\begin{cases}\boldsymbol{\bf\tau}_{j+n}^{\alpha}&n\text{
odd}\\\ \boldsymbol{\bf\tau}_{j+n}^{3-\alpha}&n\text{ even}\end{cases}\qquad
n>1\,,$

where $\alpha\in\\{x\equiv 1,y\equiv 2\\}$ and we abused the notation by using
the letter or the number depending on what is convenient in each case. The
dual XY model corresponds to $J_{\alpha}^{(n)}=\delta_{n1}J_{\alpha}$. More
generally, the model with Hamiltonian (92) is the dual of a generalised XY
model (introduced by Suzuki in Ref. [113]) with
$\mathcal{P}^{x}_{\tau}[\boldsymbol{\bf H}]=\boldsymbol{\bf H}$ and
$\mathcal{P}^{z}_{\tau}[\boldsymbol{\bf H}]=\boldsymbol{\bf H}$. In the
explicit examples, we will only consider the effect of having also
$J_{\alpha}^{(2)}$ and $J_{\alpha}^{(4)}$ different from zero (besides
$J_{\alpha}^{(1)}$). However, since it is not convenient to write the explicit
dependence on the coupling constants, we will still treat the dual generalised
XY model in (almost) full generality.

When written in the $\tau$ representation, the Hamiltonian has the form (22)
with $\boldsymbol{\bf
a}^{x,y}_{\ell}=(\Pi_{j<\ell}\boldsymbol{\bf\tau}^{z}_{j})\boldsymbol{\bf\tau}^{x,y}_{\ell}$.
As done in Section III.2.1, it is convenient to express the matrices
$\mathcal{H}_{\ell,n}$ characterising the Hamiltonian in terms of the symbol
$\mathcal{H}(p)$ as $\mathcal{H}_{\ell,n}=\int\tfrac{{\rm
d}k}{2\pi}e^{i(\ell-n)p}\mathcal{H}(p)$, where

$\mathcal{H}^{\dagger}(p)=\mathcal{H}(p)=-\mathcal{H}^{T}(-p)\,.$ (94)

The symbol of $\boldsymbol{\bf H}$ — (92) — has the additional property of
being a smooth matrix-valued function of $p$ satisfying

$\mathcal{H}(p+\pi)=\sigma^{z}\mathcal{H}(p)\sigma^{z}\,.$ (95)

Under close scrutiny, we realise that this condition represents evenness under
$\mathcal{P}^{x}_{\tau}$, which is required in order for the Hamiltonian
density to remain (strongly) quasilocal also in the $\sigma$ representation.
For the sake of simplicity, we are going to restrict ourselves to reflection
symmetric systems ($J_{x}^{(2n)}=J_{y}^{(2n)}$), whose symbol has the
additional property

$(J_{x}^{(2n)}=J_{y}^{(2n)}\quad\forall
n)\quad\Rightarrow\quad\mathrm{tr}[\mathcal{H}(p)]=0\,.$ (96)

Incidentally, the dispersion relation of the quasiparticle excitations, which
can be identified with the positive eigenvalues of the symbol, can be written
as

$\varepsilon(p)=\sqrt{\frac{\mathrm{tr}[\mathcal{H}(p)^{2}]}{2}}$ (97)

As mentioned in Section III.2.1, the symbol of the Hamiltonian generates the
time evolution of the symbol of the correlation matrix (24) and of the
operators (in the Heisenberg representation) with a quadratic fermionic
representation.

####  V.2.2 Free-fermion techniques for the entropies

Since $\boldsymbol{\bf H}$ is noninteracting in the dual representation, if
the initial state $\ket{\Psi(0)}^{(\tau)}$ is Gaussian, Wick’s theorem implies
that every expectation value can be written in terms of the correlation matrix
$\Gamma_{t=0}$ capturing the expectation values of quadratic operators. As
long as the subsystem is connected, this statement holds true even if both the
operators and the correlation matrix are restricted to the subsystem. As a
result, the density matrices $\rho_{\tilde{A}}^{(\tau)}(\Psi(t))$ and
$\rho_{\tilde{A}}^{(\tau)}(\Psi_{r}(t))$ are Gaussian at any time, that is to
say

$\rho_{\tilde{A}}^{(\tau)}(\Psi)\propto\exp\Bigl{(}\frac{1}{4}\sum_{m,n=1}^{\ell+1}\vec{a}_{m}^{T}\mathcal{W}_{m,n}(\Psi)\vec{a}_{n}\Bigr{)}$
(98)

with $\vec{a}_{m}=\begin{pmatrix}a^{x}_{m}&a^{y}_{m}\end{pmatrix}^{T}$,
$a^{x,y}_{m}$ being the Majorana (Jordan-Wigner) fermions defined in
$\tilde{A}=\\{1,\ldots,\ell+1\\}$:

$a^{x,y}_{m}=(\tau^{z})^{\otimes(m-1)}\otimes\tau^{x,y}\otimes\mathbb{I}^{(\ell+1-m)}\,.$
(99)

The Gaussian structure is very useful also to compute the Rényi entropies,
which can be expressed in terms of the correlation matrix as [114, 115]

$S_{\alpha}^{(\tau)}(\ell+1;\Psi)=\frac{\log\det\left[\left(\frac{I+\Gamma_{\tilde{A}}}{2}\right)^{\alpha}+\left(\frac{I-\Gamma_{\tilde{A}}}{2}\right)^{\alpha}\right]}{2(1-\alpha)}\,,$
(100)

where $\Gamma_{\tilde{A}}$ is the correlation matrix of $\ket{\Psi}$
restricted to $\tilde{A}$. Therefore, using the approximation (90), as
$\ell\to\infty$, the unbounded part of the excess of entropy can be expressed
in terms of the correlation matrices of
$\boldsymbol{\bf\rho}_{\mathrm{NESS}}^{(\tau)}$ and
$\boldsymbol{\bf\rho}_{\infty}^{(\tau)}$. They are reported in the following,
together with the correlation matrix of the initial state, $\Gamma_{0}(p)$, to
which they are related.

##### Correlation matrix of $\ket{\Psi(0)}^{(\tau)}$.

In the numerical investigations we have focussed on the initial state
$\ket{\Uparrow}_{\tau}$, which is one of the two representations of
$\ket{\Uparrow}_{\sigma}$ with the hidden symmetry broken — see Section IV.2.
Its correlations matrix has the symbol $\Gamma_{0}(p)=\sigma^{y}$.

##### Correlation matrix of $\boldsymbol{\bf\rho}_{\infty}^{(\tau)}$.

We have already reviewed how to obtain the correlation matrix of
$\boldsymbol{\bf\rho}_{\infty}^{(\tau)}$ after a global quench — see Section
III.2.1. In our specific case with reflection symmetry (and one-site shift
invariance), the symbol of the correlation matrix of
$\boldsymbol{\bf\rho}_{\infty}^{(\tau)}$ can be written as

$\Gamma_{\infty}(p)=\frac{\mathrm{tr}[\Gamma_{0}(p)\mathcal{H}(p)]\mathcal{H}(p)}{2\varepsilon^{2}(p)}\,,$
(101)

where we used that the only matrix commuting with $\mathcal{H}(p)$ different
from the identity (which is excluded because the initial state is reflection
symmetric and hence $\mathrm{tr}[\Gamma_{0}(p)]=0$ — see also (96)) is
proportional to $\mathcal{H}(p)$. From $\Gamma_{\infty}(p)$ we can then obtain
$\Gamma_{\tilde{A};\infty}$ by Fourier transforming $\Gamma_{\infty}(p)$ and
restricting the indices to the subsystem $\tilde{A}$.

##### Correlation matrix of $\boldsymbol{\bf\rho}_{\mathrm{NESS}}^{(\tau)}$.

We can not apply the same technique for the correlation matrix of
$\boldsymbol{\bf\rho}_{\mathrm{NESS}}^{(\tau)}$: the time evolved domain wall
$\ket{\Psi_{r}(t)}^{(\tau)}$ is not translationally invariant, so the symbol
is not even defined. This problem can be exactly solved by redefining the
symbol so as to incorporate the inhomogeneities of the Gaussian state [103].
As generally expected in bipartitioning protocols in integrable systems,
however, $\ket{\Psi_{r}(t)}^{(\tau)}$ locally relaxes along rays at fixed
$\zeta=(x-r)/t$ [116, 103]. It is then easier (and asymptotically correct in
the limit $t\rightarrow\infty$) to associate a translationally invariant
macrostate with each ray. Being translationally invariant, the macrostate is
characterised by a symbol as defined before. Finally, a generalised
hydrodynamic equation provides the connection between the macrostates at
different rays [117, 118], and one can easily extract the symbol associated
with the ray $\zeta=0$, which corresponds to the infinite time limit that we
indicated with $\boldsymbol{\bf\rho}_{\mathrm{NESS}}^{(\tau)}$. Specifically,
Appendix E.2 shows

$\Gamma_{\rm
NESS}(p)=-\mathrm{sgn}[v(p)]\frac{\mathrm{tr}[\Gamma_{0}(p)\mathcal{H}(p)]\mathbb{I}}{2\varepsilon(p)}\,,$
(102)

where $v(p)=\varepsilon^{\prime}(p)$ is the velocity of the quasiparticle
excitations.

We now have all the ingredients to work out Eq. (90) and to obtain analytical
asymptotic results for large $\ell$. In the next paragraph we show how, for
$\alpha=2,3,\ldots$, one can also obtain exact numerical results.

##### Exact numerical evaluation of the entropies.

In order to evaluate the Rényi entropies exactly, we use Eq. (82), which is
written in terms of the averaged density matrix
$\bar{\rho}^{(\tau)}_{\tilde{A}}(\Psi)$ — (81). The method is explained in
details in Appendix E.4 and is based on the realisation that each term of the
sum in Eq. (81) is Gaussian. The Rényi entropies are expressed as traces of
powers of $\bar{\rho}^{(\tau)}_{\tilde{A}}(\Psi)$, which are linear
combinations of Gaussians. The trace of a product of Gaussians can be
expressed as the square root of a determinant of a $2(\ell+1)\times 2(\ell+1)$
matrix [119, 120] and can therefore be computed efficiently. There is only a
subtlety related to the sign ambiguity of the square root. For a product of
two Gaussians (which appears in the expression for the second Rényi entropy)
the sign problem is resolved by using the property that the trace of a product
of two positive semidefinite Hermitian operators is non-negative [121]. For
higher-order Rényi entropies the sign is fixed in a less straightforward way,
as discussed in Ref. [119].

####  V.2.3 Results

Figure 12: Time evolution of the excess of the second Rényi entropy, generated
in the XY model (19) and in the dual XY model (27). For large times the excess
in the XY model reaches a zero value, while in the dual XY model the its
asymptotic value is nonzero. It is in agreement with the result obtained via
the method used also in Fig. 13 (here the latter results are presented as
dashed lines). The spin flip occurs in the middle of the subsystem $A$
($r=5,|A|=10$). Figure 13: Exact results for the excess of the second Rényi
entropy after the quench, for different values of model’s parameters. The
excess grows logarithmically with the subsystem size (the $x$-axis is in the
log scale). The analytical results for the coefficients $b_{\alpha}$ (i.e.,
the slopes), given by Eq. (105), are in agreement with the numerical data and
indicated in the plot.

Both with and without the local spin flip in the initial state, the quench is
global and the Rényi entropies grow linearly in time until they saturate to an
extensive value. In noninteracting models this was established in Refs. [63,
122, 103]; it was later generalised to interacting integrable models in Refs.
[123, 124, 125]. Less is known in regards to the subleading behaviour.

For the quench from the initial state $\ket{\Uparrow}_{\tau}$ we find

$\lim_{t\rightarrow\infty}S_{\alpha}^{(\tau)}(\ell+1;\Psi(t))=a_{\alpha}\ell+\mathcal{O}(1)\,,$
(103)

with positive constants $a_{\alpha}$, whose values can be determined
analytically using the integral representation for the entropies [114] and the
Szegő-Widom theorem [126, 127]. Their values are reported in Appendix E.3.

For the quench from the domain-wall state (the dual representation of the
state $\boldsymbol{\bf\sigma}_{r}^{x}\ket{\Uparrow}_{\sigma}$) we instead find
a subleading term that is logarithmic in $\ell$:

$\lim_{t\rightarrow\infty}S_{\alpha}^{(\tau)}(\ell;\Psi_{r}(t))=a_{\alpha}\ell+b_{\alpha}\log\ell+\mathcal{O}(1),$
(104)

with positive constants $b_{\alpha}$. This result is obtained using the
asymptotic formula for the determinant of a block Toeplitz matrix with
piecewise-continuous symbols, conjectured and checked numerically in Refs.
[128, 129]. The values of $b_{\alpha}$ are related to the discontinuities of
the symbol (102), which coincide with the zeroes of the velocity $v(p)$ of the
quasiparticle excitations. For instance, for $\alpha=2$ we have

$b_{2}=\frac{2}{\pi^{2}}\\!\sum_{p\in(-\pi,\pi]:\,v(p)=0}\mathrm{Arg}^{2}\left(1+i\frac{\mathrm{tr}[\Gamma(p;0)\mathcal{H}(p)]}{2\varepsilon(p)}\right)\,,$
(105)

where $\mathrm{Arg}$ denotes (the principal value of) the argument, and the
sum is over the zeroes of $v(p)$. For general $\alpha$ the expression can be
found in Appendix E.3.

In the symmetric generalised XY model there are always zeroes of $v(p)$ at the
momenta $p=-\pi/2,0,\pi/2,\pi$, and additional ones can appear depending on
the coupling constants $J_{\alpha}^{(n)}$. In the special case of the XY model
these four zeros are the only ones and we find

$b^{\rm(XY)}_{\alpha}=\frac{1}{6}\left(1+\frac{1}{\alpha}\right)\,.$ (106)

For more general coupling constants $J_{\alpha}^{(n)}$ this expression is only
a lower bound, i.e., $b_{\alpha}\geq b^{\rm(XY)}_{\alpha}$.

From Eqs. (104) and (105) it follows that the excess (88) of Rényi entropies
grows logarithmically with $\ell$:

$\lim_{t\rightarrow\infty}\Delta_{r}S_{\alpha}(\ell,t)=b_{\alpha}\log\ell+\mathcal{O}(1)\,.$
(107)

This prediction is compared against exact numerical data for $\alpha=2$: Fig.
12 shows the time evolution of $\Delta_{r}S_{\alpha}(\ell,t)$, while in Fig.
13 its dependence on $\ell$ is presented.

Finally, the von Neumann entropy is obtained by taking the limit $\alpha\to
1^{+}$; in the dual XY model it is given by $b_{1}^{\rm(XY)}=1/3$.

In the next section we conclude our preliminary study of the signatures of
semilocal order after global quenches by showing that the order melts down
over large time scales even when the initial state is not symmetric.

###  V.3 Melting of the order

Just as a continuous phase transition is characterised by the behaviour of the
system when the critical point is approached, so is the symmetry protected
topological order when the symmetry is only approximate. We are familiar with
an example of this: in Fig. 5 we showed that, after a quench from the tilted
state $\ket{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}$, the time
scale on which local observables relax towards the predictions of the
(generalised) Gibbs ensemble grows as $\theta$ approaches $0$. In the limit
$\theta\to 0$ the relaxation time diverges: the canonical form of a quasilocal
ensemble, valid for any finite $\theta$, yields a wrong prediction when
$\theta=0$. At the symmetric point, instead, the correct prediction is given
by the $\mathbbm{Z}_{2}$-semilocal ensemble, and the corresponding relaxation
time is again finite. This implies that the limits $t\to\infty$ and $\theta\to
0$ do not commute (note that the existence of noncommuting limits is one of
the hallmarks of spontaneous symmetry breaking [57] as well as of
prethermalisation/prerelaxation behaviours after global quenches [45, 46, 47,
48]).

Signatures of semilocal order thus exist also in states that weakly break the
symmetry with which the order is associated. In particular, we envisage the
existence of a scaling limit with $t\to\infty$ and $\theta\to 0$, in which
local observables exhibit some form of prerelaxation. This scaling limit is
discussed in Section V.3.1.

Symmetric initial states exhibiting semilocal order can be regarded as ground
states of prequench Hamiltonians $\boldsymbol{\bf H}_{0}$ (see Section IV.2).
Since string order does not survive nonzero temperatures [12], the
prerelaxation behaviour mentioned above is expected also for sufficiently low,
but nonzero temperatures. We will see in Section V.3.2 that the melting of the
order arising at small nonzero temperatures has essentially the same origin as
the prerelaxation occurring for small symmetry-breaking unitary perturbations
of the initial state.

####  V.3.1 Slow relaxation after weak symmetry breaking

Being interested in the time scales on which semilocal order disappears after
a weak symmetry breaking, we consider an initial state prepared close to the
symmetric point, i.e.,
$\ket{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}$ for small
$\theta>0$.

Since
$\braket{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}{\boldsymbol{\bf\sigma}_{\ell}^{y}}{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}=0$
and
$\braket{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}{\boldsymbol{\bf\sigma}_{\ell}^{x}}{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}=\sin\theta\to
0$ as $\theta\to 0$, the local operators to consider in order to explore the
melting of the order consist solely of products of
$\boldsymbol{\bf\sigma}^{z}$. This includes, in particular, the strings
$\prod_{\ell=0}^{n-1}\boldsymbol{\bf\sigma}_{j+\ell}^{z}$, which, when
evaluated in a particular state, can be regarded as truncations of a semilocal
order parameter, corresponding to the limit $n\rightarrow\infty$. They satisfy

$\braket{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}{\prod_{\ell=0}^{n-1}\boldsymbol{\bf\sigma}_{j+\ell}^{z}}{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}=\cos^{n}\theta\,.$
(108)

The inverse duality transformation (38) maps operators composed solely of
$\boldsymbol{\bf\sigma}^{z}_{j}$ into operators consisting of an even number
of $\boldsymbol{\bf\tau}_{j}^{z}$. Specifically, the strings appearing in Eq.
(108) are mapped into
$\boldsymbol{\bf\tau}^{z}_{j}\boldsymbol{\bf\tau}^{z}_{j+n}$, which are the
4-fermion operators $-\boldsymbol{\bf a}^{x}_{j}\boldsymbol{\bf
a}^{y}_{j}\boldsymbol{\bf a}^{x}_{j+n}\boldsymbol{\bf a}^{y}_{j+n}$. To
describe the expectation values of such operators we introduce an effective
density matrix $\boldsymbol{\bf\rho}_{\rm eff}$ that is only asked to capture
the expectation value of 4-fermion operators $\boldsymbol{\bf O}$:

$\braket{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}{\boldsymbol{\bf
O}^{(\sigma)}(t)}{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}=\mathrm{tr}[\boldsymbol{\bf\rho}_{\rm
eff}(t)\boldsymbol{\bf O}]\,,$ (109)

where $\boldsymbol{\bf O}^{(\sigma)}(t)=e^{i\boldsymbol{\bf
H}^{(\sigma)}t}\boldsymbol{\bf O}e^{-i\boldsymbol{\bf H}^{(\sigma)}t}$ is the
time evolution of $\boldsymbol{\bf O}$ in the Heisenberg picture, and
$\boldsymbol{\bf H}^{(\sigma)}$ is the dual XY model’s Hamiltonian, given in
Eq. (27).

For each $t$ we can expand $\boldsymbol{\bf O}^{(\sigma)}(t)$ in the operator
basis $\\{\boldsymbol{\bf e}_{\alpha}\\}$ as $\boldsymbol{\bf
O}^{(\sigma)}(t)=\sum_{\alpha}({\rm tr}[\boldsymbol{\bf O}(t)\boldsymbol{\bf
e}_{\alpha}]/{\rm tr}\boldsymbol{\bf I})\,\boldsymbol{\bf e}_{\alpha}$, where
$\boldsymbol{\bf e}_{\alpha}$ can be chosen, for example, to be strings of
Pauli matrices of various lengths. The time evolution preserves the number of
fermions, so the expansion of $\boldsymbol{\bf O}^{(\sigma)}(t)$ retains only
the basis elements $\boldsymbol{\bf e}_{\alpha}$ corresponding to $4$-fermion
operators. Moreover, since in Eq. (109) we project onto the tilted initial
state, we can keep only the basis elements for which
$\braket{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}{\boldsymbol{\bf
e}_{\alpha}}{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}\neq 0$. The
latter again correspond to strings of $\boldsymbol{\bf\sigma}^{z}_{j}$, whence
we obtain

$\displaystyle\begin{aligned}
\braket{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}{\boldsymbol{\bf
O}^{(\sigma)}(t)}{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}=&\sum_{j\in\mathbbm{Z}}\sum_{n\geq
1}\braket{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}{\prod_{\ell=0}^{n-1}\boldsymbol{\bf\sigma}_{j+\ell}^{z}}{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}\\\
&\times\frac{\mathrm{tr}[\boldsymbol{\bf
O}^{(\tau)}(t)\boldsymbol{\bf\tau}_{j}^{z}\boldsymbol{\bf\tau}_{j+n}^{z}]}{\mathrm{tr}\boldsymbol{\bf
I}}\,.\end{aligned}$ (110)

The effective density matrix for small $\theta$ is thus 121212In case the
reader is not comfortable with replacing a density matrix with a traceless
operator, we point out that the result does not change if a term proportional
to the identity is added to $\boldsymbol{\bf\rho}_{\rm eff}$.

$\displaystyle\boldsymbol{\bf\rho}_{\rm eff}(t)\sim e^{-i\boldsymbol{\bf
H}^{(\tau)}t}\Big{[}\frac{1}{{\rm tr}\boldsymbol{\bf
I}}\sum_{j\in\mathbbm{Z}}\sum_{n\geq
1}e^{-n\theta^{2}/2}\boldsymbol{\bf\tau}_{j}^{z}\boldsymbol{\bf\tau}_{j+n}^{z}\Big{]}e^{i\boldsymbol{\bf
H}^{(\tau)}t}\,,$ (111)

where we used Eq. (108) with its right-hand side rewritten, for small
$\theta$, as $\cos^{n}\theta=e^{-n\theta^{2}/2}[1+\mathcal{O}(\theta^{4})]$.

The right-hand side of Eq. (111) is in $\tau$ representation: $\boldsymbol{\bf
H}^{(\tau)}$ corresponds to an XY Hamiltonian (19), so the standard free-
fermionic techniques apply for calculation of the expectation values — see
Appendix F. For example, the truncated string order parameter
$Z_{n}=\braket{\prod_{\ell=0}^{n-1}\boldsymbol{\bf\sigma}_{j+\ell}^{z}}$,
which, in the fermionic language, corresponds to $-\braket{\boldsymbol{\bf
a}^{x}_{1}\boldsymbol{\bf a}^{y}_{1}\boldsymbol{\bf
a}^{x}_{n+1}\boldsymbol{\bf a}^{y}_{n+1}}$, can be efficiently computed by
invoking Wick’s theorem: the full expression derived in the scaling limit
$t\to\infty$, $\theta\to 0$, at fixed finite $0<\theta^{2}t<\infty$ for any
$n\geq 1$ is reported in Appendix F. It can be simplified in the limit
$n\to\infty$, where it becomes

$\displaystyle Z_{\infty}\\!=\\!\frac{1}{2}\int\\!\frac{{\rm
d}^{2}k}{(2\pi)^{2}}(e^{-|v(k_{1})-v(k_{2})|\frac{\theta^{2}t}{2}}\\!+\\!e^{-|v(k_{1})+v(k_{2})|\frac{\theta^{2}t}{2}})\frac{(J_{x}\\!+\\!J_{y})^{4}\cos^{2}k_{1}\cos^{2}k_{2}}{[J_{x}^{2}\\!+\\!J_{y}^{2}\\!+\\!2J_{x}J_{y}\cos(2k_{1})][J_{x}^{2}\\!+\\!J_{y}^{2}\\!+\\!2J_{x}J_{y}\cos(2k_{2})]}\,.$
(112)

From this expression one can read off the time scale $\xi$ on which the order
melts: $\xi\sim\theta^{-2}$. In fact this time scale is inherent to the
behaviour of the effective density matrix $\boldsymbol{\bf\rho}_{\rm eff}$
(see Appendix F) and is therefore exhibited also by the time evolution of
$Z_{n}$ for finite $n$. Fig. 14 corroborates it numerically for $n=1$, i.e.,
in the case of $Z_{1}=\braket{\boldsymbol{\bf\sigma}_{\ell}^{z}}$.

We warn the reader that setting $\theta^{2}t=0$ in Eq. (112) does not
correspond to the expectation value of the string order parameter at the
initial time, it is rather the starting value from which the prerelaxation
behaviour ensues. This is because the derivation is based on asymptotics for
large $t>0$ and small $\theta>0$. Despite this, setting $\theta^{2}t=0$ in Eq.
(112) correctly reproduces the large time value of the string order parameter
for $\theta=0$.

Figure 14: Numerically computed (iTEBD) expectation value
$\braket{\boldsymbol{\bf\sigma}^{z}_{\ell}}$ plotted versus $\theta^{2}t$ is
consistent with a scaling prediction derived in Appendix F. In the limit
$\theta\to 0$ the scaling prediction yields
$\braket{\boldsymbol{\bf\sigma}^{z}_{\ell}}\approx 0.5626$, which is also the
GGE prediction for the expectation value of
$\braket{\boldsymbol{\bf\sigma}_{\ell}^{z}\boldsymbol{\bf\sigma}_{\ell+1}^{z}}$
in the XY model prepared in the initial state $\ket{\Uparrow}$ — see Fig. 7.
The prediction coming from a GGE with the standard pseudolocal charges of the
dual XY model is instead zero. The inset is in logarithmic scale.

####  V.3.2 Melting at low finite temperatures

The simplest example of a finite-temperature initial state reads

$\boldsymbol{\bf\rho}(0)=\frac{e^{\beta\sum_{\ell}\boldsymbol{\bf\sigma}_{\ell}^{z}}}{{\rm
tr}[e^{\beta\sum_{\ell}\boldsymbol{\bf\sigma}_{\ell}^{z}}]}\,.$ (113)

In the limit $\beta\to\infty$ this thermal state reproduces the symmetric
state with semilocal order, as considered in Section IV.2. For any finite
$\beta$ $\boldsymbol{\bf\rho}(0)$ contains only terms composed of products of
$\boldsymbol{\bf\sigma}^{z}_{j}$. Hence, focusing again only on 4-fermion
operators, the only relevant operator with a nonzero expectation value is
again $\boldsymbol{\bf a}_{\ell}^{x}\boldsymbol{\bf
a}_{\ell}^{y}\boldsymbol{\bf a}_{\ell+n}^{x}\boldsymbol{\bf a}_{\ell+n}^{y}$,
i.e., a fermionic representation of a string of Pauli matrices
$\boldsymbol{\bf\sigma}^{z}_{j}$. This time its expectation value reads

$-\braket{\boldsymbol{\bf a}_{\ell}^{x}\boldsymbol{\bf
a}_{\ell}^{y}\boldsymbol{\bf a}_{\ell+n}^{x}\boldsymbol{\bf
a}_{\ell+n}^{y}}=\tanh^{n}\beta\,.$ (114)

For large $\beta$ and small $\theta$, this is equivalent to setting the
rotation angle in the tilted initial state
$\ket{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}$ to be

$\theta_{\beta}\approx 2e^{-\beta}\,.$ (115)

Having established a connection with the tilted initial state, we can retrace
the steps in the previous section and observe that the nontrivial scaling
limit is now $t\to\infty$, $\beta\to\infty$, at fixed finite $4e^{-2\beta}t$.
The time scale $\xi$ on which the order melts at a low, but finite
temperature, after the system has been prepared in the state
$\boldsymbol{\bf\rho}(0)$, is thus
$\xi\sim\theta_{\beta}^{-2}=(e^{\beta}/2)^{2}$.

Finally, we observe a similarity with the findings of Ref. [47], in which an
analogous scaling behaviour was observed at low temperature in a phase in
which a symmetry is spontaneously broken at zero temperature. The reader can
understand this connection as a manifestation of the duality between the
quasilocal theory and the even $\mathbb{Z}_{2}$-semilocal theory, the latter
being characterised by a spontaneous symmetry breaking at zero temperature —
cf. (38).

##  VI Discussion

We have shown that symmetry-protected topological phases of matter can emerge
also at infinite times after quenches of global Hamiltonian parameters in
isolated many-body quantum systems. We have traced this phenomenon back to the
existence of conserved operators (semilocal charges) that do not belong to the
natural theory in which local observables are represented by local operators
and their quasilocal completion. From that perspective, the topological
character of the system is manifest: different theories of local observables
are associated with different concepts of locality for the operators; being
outside the quasilocal theory, semilocal charges can keep memory of
information that is not contained in the bulk of the system.

We aimed at reconciling the observation that semilocal charges affect the
stationary values of local operators with the established belief that the
those values are solely determined by the local conserved quantities. This has
led us to introduce statistical ensembles sensitive to string (semilocal)
order, which we called semilocal (generalised) Gibbs ensembles.

Finally, we have reported signatures of semilocal order in the entropies of
subsystems and in the time evolution of systems that are not exactly
symmetric. In particular, we have found that a perturbation with a finite
support affects the entanglement entropies of subsystems with an exceptional
leading correction scaling as the logarithm of the subsystem’s length.

### Open problems

  * •

The description we have proposed relies on translational invariance. If we
break it, the state can be globally non-symmetric still exhibiting the local
symmetries characterising semilocal order, in the sense that the two-point
functions of semilocal operators could approach nonzero values at intermediate
distances. In addition, since the expectation values of semilocal charges are
also affected by local inhomogeneities, the latter could look relevant despite
the Lieb-Robinson bounds ruling out their importance at space-like distances.
A theory that could capture time evolution of inhomogeneous states goes beyond
the purposes of this work and will be addressed in separate investigations.

  * •

We define the $G$-semilocal (generalised) Gibbs ensemble in an extended space
that includes operators that do not represent local observables. Having
realised the existence of an ambiguity in what should be called “local
observable” in the extended space, we propose to project back the ensemble
onto a theory of local observables. Taking into account that similar
situations seem to occur in quantum field theories, where fields that are
semilocal with respect to each other are part of the theory, we wonder whether
such a projection is really required. Perhaps a more rigorous treatment along
the lines of the algebraic formulation of a local quantum field theory [131]
could clarify this issue. Incidentally, we wonder whether the finding of Ref.
[132] could be connected with the picture presented here.

  * •

We show that, in the presence of semilocal order, the excess of entropy
triggered by a local perturbation of the initial state develops a logarithmic
dependence on the subsystem’s length. We expect that even the entanglement
properties of the symmetric translationally invariant system (without any
perturbation) should exhibit exceptional properties, but, at the moment, this
is still an open question.

  * •

Despite our description applying to more general situations, the explicit
examples that we consider are not interacting. This choice has a twofold
motivation: on the one hand, we aim at providing the cleanest examples
exhibiting the exotic phenomenology of semilocal order after global quenches;
on the other hand, interactions could introduce complications that need clues
from simpler models in order to be solved. Whether the interplay between
interactions and semilocality could enable additional phenomena is an open
problem.

  * •

We discuss semilocal order both in generic and in integrable systems.
Extending the theory to nonabelian integrable structures remains an intriguing
development, in which semilocal dynamical symmetries could play an important
role [133].

###### Acknowledgements.

We thank Kemal Bidzhiev, Saverio Bocini, Mario Collura, Fabian Essler, Ivan
Kukuljan, Leonardo Mazza, Spyros Sotiriadis, and Gabor Takacs for useful
discussions. This work was supported by the European Research Council under
the Starting Grant No. 805252 LoCoMacro.

##  Appendix A Convergence of the iTEBD data

Here we estimate the times $t$ up to which the results of the iTEBD time
evolution in Fig. 5 are virtually indistinguishable from those of the exact
calculation. By “virtually indistinguishable” we mean, that the difference
between the numerically computed and the exact expectation value
$\braket{\boldsymbol{\bf\sigma}_{\ell}^{z}}$ at any particular time should be
smaller than some accuracy, which we choose to be $0.01$ in Fig. 5.

We compare the data from calculations with maximum bond dimensions that are
increasing exponentially up to the largest one, $M_{\rm max}=1000$ (bond
dimension is expected to increase exponentially with $t$). The results of the
calculation with $M_{\rm max}=1000$ are assumed to match the exact time
evolution longer than those at smaller $M_{\rm max}$, which is why we take
them as a reference. Fig. 15 shows the differences
$\braket{\boldsymbol{\bf\sigma}_{\ell}^{z}}_{1000}-\braket{\boldsymbol{\bf\sigma}_{\ell}^{z}}_{M_{\rm
max}}$, where the subscript denotes the maximum bond dimension in the
calculation. From the times at which these differences become virtually
discernible (i.e., larger than $0.01$), we estimate that the calculations with
$M_{\rm max}=1000$ can be trusted up to times $t\sim 4.5$.

Figure 15: Differences between the results of the calculation with $M_{\rm
max}=1000$ and smaller maximum bond dimensions. Red circles denote when the
differences become noticeable: the corresponding times grow linearly.

##  Appendix B $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ duality

We report here what could be considered as a double application of the rotated
Kramers-Wannier duality transformation. We denote it by
$\mathcal{D}_{\mathbb{Z}_{2}\times\mathbb{Z}_{2}}$ and it is an algebra
homomorphism mapping the algebra of observables generated by
$\\{\boldsymbol{\bf\tau},\boldsymbol{\bf T}^{x}_{\tau;-2,0},\boldsymbol{\bf
T}^{x}_{\tau;-2,1}\\}$ into the one generated by
$\\{\boldsymbol{\bf\sigma},\boldsymbol{\bf T}^{z}_{\sigma;2,0},\boldsymbol{\bf
T}^{z}_{\sigma;2,1}\\}$. Here again $\boldsymbol{\bf\sigma}_{j}^{\alpha}$ and
$\boldsymbol{\bf\tau}_{j}^{\alpha}$ act like Pauli matrices on site $j$ and
like the identity elsewhere. The additional operators $\boldsymbol{\bf
T}_{\sigma;2s,y}^{\alpha}$ are instead defined by the following conditions

$\displaystyle[\boldsymbol{\bf T}_{\sigma;2s,y}^{\alpha}]^{2}=\boldsymbol{\bf
I}\,,\quad\boldsymbol{\bf T}_{\sigma;2s,1}^{\alpha}\boldsymbol{\bf
T}_{\sigma;2s,0}^{\alpha}=\boldsymbol{\bf
T}_{\sigma;2s,0}^{\alpha}\boldsymbol{\bf
T}_{\sigma;2s,1}^{\alpha}=\boldsymbol{\bf T}_{\sigma;s}^{\alpha}\,,$ (116)
$\displaystyle\boldsymbol{\bf T}_{\sigma;2s,y}^{\alpha}\boldsymbol{\bf
O}\boldsymbol{\bf
T}_{\sigma;2s,y}^{\alpha}=\lim_{n\to\infty}\Big{[}\prod_{\ell=0}^{n}\boldsymbol{\bf\sigma}_{s(2\ell+y)}^{\alpha}\Big{]}\boldsymbol{\bf
O}\Big{[}\prod_{\ell=0}^{n}\boldsymbol{\bf\sigma}_{s(2\ell+y)}^{\alpha}\Big{]}$

for all local operators $\boldsymbol{\bf O}$, extended then by linearity as in
the case of $\mathcal{D}_{\mathbbm{Z}_{2}}$, in such a way that
$\boldsymbol{\bf O}_{1}+\boldsymbol{\bf O}_{2}\boldsymbol{\bf
T}_{\sigma;2s,0}^{\alpha}+\boldsymbol{\bf O}_{3}\boldsymbol{\bf
T}_{\sigma;2s,1}^{\alpha}+\boldsymbol{\bf O}_{4}\boldsymbol{\bf
T}_{\sigma;s}^{\alpha}$ is transformed by the adjoint action of
$\boldsymbol{\bf T}_{\sigma;2s,y}^{\alpha}$ into an analogous operator in
which the local operators $\boldsymbol{\bf O}_{j}$ are replaced by
$\boldsymbol{\bf T}_{\sigma;2s,y}^{\alpha}\boldsymbol{\bf
O}_{j}\boldsymbol{\bf T}_{\sigma;2s,y}^{\alpha}$. Explicitly, the duality
transformation $\mathcal{D}_{\mathbb{Z}_{2}\times\mathbb{Z}_{2}}$ reads

$\displaystyle\boldsymbol{\bf\Pi}^{x}_{\tau,-2}(j)=$
$\displaystyle\boldsymbol{\bf\sigma}_{j}^{x}\,,$ (117)
$\displaystyle\boldsymbol{\bf\tau}_{j}^{y}=$
$\displaystyle\boldsymbol{\bf\sigma}_{j-2}^{x}\boldsymbol{\bf\sigma}_{j}^{y}\boldsymbol{\bf\Pi}_{\sigma,+2}^{z}(j+2)\,,$
$\displaystyle\boldsymbol{\bf\tau}_{j}^{z}=$
$\displaystyle\boldsymbol{\bf\Pi}_{\sigma,+2}^{z}(j)\,,$

where $\boldsymbol{\bf\Pi}^{x}_{\tau,-2}(j)$ and
$\boldsymbol{\bf\Pi}_{\sigma,+2}^{z}(j)$ are defined as
$\boldsymbol{\bf\Pi}_{\sigma,s}^{\alpha}(j)$ in the sublattice of sites with
the same parity of $j$. This is obtained using $\boldsymbol{\bf
T}_{\sigma;2s,j\mod 2}^{\alpha}$ instead of $\boldsymbol{\bf
T}_{\sigma;s}^{\alpha}$. In the bulk the transformation can be expressed in
the more standard form

$\boldsymbol{\bf\tau}_{j}^{x}=\boldsymbol{\bf\sigma}_{j-2}^{x}\boldsymbol{\bf\sigma}_{j}^{x}\qquad\boldsymbol{\bf\tau}_{j}^{z}\boldsymbol{\bf\tau}_{j+2}^{z}=\boldsymbol{\bf\sigma}_{j}^{z}\,,$
(118)

which corresponds to two independent rotated Kramers-Wannier dualities on the
sublattices consisting of the sites labelled by numbers with the same parity.

Finally, the inverse duality
$\mathcal{D}^{-1}_{\mathbb{Z}_{2}\times\mathbb{Z}_{2}}$ reads

$\displaystyle\boldsymbol{\bf\sigma}_{j}^{x}=$
$\displaystyle\boldsymbol{\bf\Pi}^{x}_{\tau,-2}(j)\,,$ (119)
$\displaystyle\boldsymbol{\bf\sigma}_{j}^{y}=$
$\displaystyle\boldsymbol{\bf\Pi}^{x}_{\tau,-2}(j-2)\boldsymbol{\bf\tau}_{j}^{y}\boldsymbol{\bf\tau}_{j+2}^{z}\,,$
$\displaystyle\boldsymbol{\bf\Pi}_{\sigma,+2}^{z}(j)=$
$\displaystyle\boldsymbol{\bf\tau}_{j}^{z}\,.$

##  Appendix C Zero-temperature phases of the dual XY model

Since the (rotated) Kramers-Wannier transformation maps the Hamiltonian (27)
of the dual XY model into the Hamiltonian (19) of the XY model, we can
immediately infer that the model is noncritical for $|J_{x}|\neq|J_{y}|$ and
critical for $|J_{x}|=|J_{y}|$.

For $|J_{y}|<|J_{x}|$ it should thus be possible to connect smoothly the
ground state(s) of the model to the one(s) at the classical point $J_{y}=0$,
above which there is a finite energy gap from the first excited state(s) in
the thermodynamic limit. It is easy to see that at the classical point the
ground states exhibit magnetic order with translational invariance under a
shift by one or several sites in the bulk, the simplest being the
ferromagnetic ground state $\otimes_{j}\ket{\leftarrow}_{j}$ for $J_{x}<0$
(and the one related by the spin flip symmetry). A symmetry breaking field
[57] can then select a ground state with a particular pattern of the magnetic
order and we are in the standard Landau phase. Note that, besides symmetry
breaking, there could be other ambiguities in the zero temperature limit. For
example, for $J_{y}>0$, imposing periodic boundary conditions on a chain with
an odd number of sites frustrates the system, with consequences on the
properties of the state that survive the thermodynamic limit [134, 135, 136].

It is simple to show that the Hamiltonian of the dual XY model with
$|J_{y}|<|J_{x}|$ can be mapped through a symmetry-preserving local unitary
transformation, into an operator, with strongly quasilocal densities, that
commutes with

$\boldsymbol{\bf
H}(J_{x},0)=J_{x}\sum_{\ell}\boldsymbol{\bf\sigma}_{\ell-1}^{x}\boldsymbol{\bf\sigma}_{\ell+1}^{x}$
(120)

and shares the same ground states of $\boldsymbol{\bf H}(J_{x},0)$. The latter
states completely break the $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ symmetry of
the Hamiltonian (see Section IV.3.4), belonging therefore to a standard
ordered phase à la Landau. Specifically and just for the sake of example, the
symmetry breaking states of $\boldsymbol{\bf H}(J_{x},0)$ with clustering
properties are mapped into the analogues with $J_{y}\neq 0$ (and
$|J_{y}|<|J_{x}|$) by the transformation $e^{-i\boldsymbol{\bf W}}$, where
$\boldsymbol{\bf W}$ is a dual generalised XY model — Section V.2.1 and
Appendix E.2 — characterised by the symbol

$W(p)=-\frac{1}{2}\arctan\frac{2J_{y}\tan
p}{J_{x}+J_{y}+(J_{x}-J_{y})\tan^{2}p}\sigma^{z}\,.$ (121)

Analogously, there is a transformation mapping the Hamiltonian of the dual XY
model with $|J_{y}|<|J_{x}|$ into an operator, with strongly quasilocal
densities, that commutes with

$\boldsymbol{\bf
H}(0,J_{y})=-J_{y}\sum_{\ell}\boldsymbol{\bf\sigma}_{\ell-1}^{x}\boldsymbol{\bf\sigma}_{\ell}^{z}\boldsymbol{\bf\sigma}_{\ell+1}^{x}$
(122)

and shares the same ground state of $\boldsymbol{\bf H}(0,J_{y})$. As
discussed in Ref. [86], the latter Hamiltonian has the ground state in a
nontrivial $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ protected topological phase,
which therefore describes the entire parameter space $|J_{y}|>|J_{x}|$.

##  Appendix D Tilted initial state

###  D.1 Tilted state
$\ket{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}_{\sigma}$

The mapping of the tilted initial state (20) is perhaps most conveniently
performed on a finite lattice of size $L$. To this end we report the finite-
system inverse duality transformation

$\displaystyle\begin{aligned}
\boldsymbol{\bf\sigma}_{j}^{x}=&\boldsymbol{\bf\Pi}_{\tau,-}^{x}(j)\,,\\\
\boldsymbol{\bf\sigma}_{j}^{y}=&\begin{cases}\boldsymbol{\bf\tau}_{1}^{x}\boldsymbol{\bf\tau}_{2}^{z}&j=1\\\
\boldsymbol{\bf\Pi}^{x}_{\tau,-}(j-1)\boldsymbol{\bf\tau}_{j}^{y}\boldsymbol{\bf\tau}_{j+1}^{z}&1<j<L\\\
-\boldsymbol{\bf\tau}_{L}^{z}&j=L\,,\end{cases}\\\
\boldsymbol{\bf\sigma}_{j}^{z}=&\begin{cases}\boldsymbol{\bf\tau}_{j}^{z}\boldsymbol{\bf\tau}_{j+1}^{z}&1\leq
j<L\\\
\boldsymbol{\bf\Pi}^{x}_{\tau}\boldsymbol{\bf\tau}_{1}^{z}\boldsymbol{\bf\tau}_{L}^{z}&j=L\,,\end{cases}\end{aligned}$
(123)

where

$\displaystyle\boldsymbol{\bf\Pi}^{x}_{\tau,-}(j)=-\boldsymbol{\bf\tau}_{1}^{y}\prod_{\ell=2}^{j}\boldsymbol{\bf\tau}_{\ell}^{x},\qquad\boldsymbol{\bf\Pi}^{x}_{\tau}=\prod_{\ell=1}^{L}\boldsymbol{\bf\tau}_{\ell}^{x}\,.$
(124)

The state

$\displaystyle\ket{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}_{\sigma}\equiv\bigotimes_{j=1}^{L}\begin{pmatrix}\cos\frac{\theta}{2}\\\
\sin\frac{\theta}{2}\end{pmatrix}\,,$ (125)

is an eigenstate of all
$\boldsymbol{\bf\eta}^{(\sigma)}_{j}(\theta)=\cos\theta\boldsymbol{\bf\sigma}_{j}^{z}+\sin\theta\boldsymbol{\bf\sigma}_{j}^{x}$
for $j=1,\ldots L$, with eigenvalue $1$. Applying the inverse duality
transformation on a finite system (123), we have

$\displaystyle\boldsymbol{\bf\eta}^{(\tau)}_{j}\\!(\theta)\\!=\\!\begin{cases}\cos\theta\boldsymbol{\bf\tau}_{j}^{z}\boldsymbol{\bf\tau}_{j+1}^{z}\\!+\\!\sin\theta\boldsymbol{\bf\Pi}_{\tau,-}^{x}(j)&1\leq
j<L\\\
\cos\theta\boldsymbol{\bf\Pi}_{\tau}^{x}\boldsymbol{\bf\tau}_{1}^{z}\boldsymbol{\bf\tau}_{L}^{z}\\!-\\!i\sin\theta\boldsymbol{\bf\Pi}_{\tau}^{x}\boldsymbol{\bf\tau}_{1}^{z}&j=L\,.\end{cases}$
(126)

In the basis of all $\boldsymbol{\bf\tau}_{j}^{z}$ (their eigenvalues being
denoted by $s_{j}\in\\{-1,1\\}$) the tilted state reads

$\displaystyle\ket{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}_{\sigma}=\sum_{\underline{s}\in\\{-1,1\\}^{\times
L}}c(\underline{s})\ket{\underline{s}}_{\tau}\,.$ (127)

Applying the operator (126) we now obtain relations

$\displaystyle\begin{aligned} \frac{c(-s_{1},\ldots-s_{j},s_{j+1},\ldots
s_{L})}{c(\underline{s})}=\\\
=\begin{cases}\frac{1-s_{j}s_{j+1}\cos\theta}{is_{1}\sin\theta}&1\leq j<L\\\
s_{1}s_{L}e^{-i\theta s_{L}}&j=L\,,\end{cases}\end{aligned}$ (128)

which are solved by

$\displaystyle
c(\underline{s})=\frac{1}{\sqrt{Z}}e^{i\frac{\pi}{4}(s_{1}-s_{L})+i\frac{\theta}{2}s_{L}}e^{-\frac{\beta}{2}E(\underline{s})}\,.$
(129)

Here $E(\underline{s})=-\sum_{j=1}^{L-1}s_{j}s_{j+1}$,
$\beta=-\log\tan(\theta/2)$, and

$\displaystyle Z=\sum_{\underline{s}\in\\{-1,1\\}^{\times L}}e^{-\beta
E(\underline{s})}\,.$ (130)

We recognize that the squared absolute value $|c(\underline{s})|^{2}$ yields
the canonical ensemble probability distribution for the one-dimensional
classical Ising model with free boundary conditions: $\beta$ plays the role of
an effective inverse temperature, while $Z$ is the corresponding partition
function. In this regard a nonzero angle $\theta$ corresponds to a nonzero
temperature.

Finally, the partition function (130) can be evaluated using the transfer
matrix formalism (see, e.g., Ref. [137]). We define transfer matrix

$\displaystyle T=\frac{1+\sigma^{x}}{2}+\frac{1-\sigma^{x}}{2}\cos\theta$
(131)

with elements $T_{s,s^{\prime}}=(1/2)e^{\beta ss^{\prime}}\sin\theta$,
$s,s^{\prime}\in\\{-1,1\\}$, and use it to rewrite the partition function:

$\displaystyle
Z=\left(\frac{2}{\sin\theta}\right)^{L-1}\sum_{s_{1},s_{L}\in\\{-1,1\\}}\left[T^{L-1}\right]_{s_{1},s_{L}}.$
(132)

The powers of the transfer matrix are readily obtained, i.e.,
$T^{n}=(1+\sigma^{x})/2+[(1-\sigma^{x})/2]\cos^{n}\theta$, whence

$Z=2\left(\frac{2}{\sin\theta}\right)^{L-1}.$ (133)

###  D.2 Decay of correlations in the state
$\ket{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}_{\sigma}$

The transfer matrix (131) enables an efficient method for calculation of local
operators that act in the bulk of the spin chain, i.e., far from the
boundaries, which will eventually be sent to infinity as the thermodynamic
limit is taken. The expectation value of a local operator $\boldsymbol{\bf O}$
acting on the sites $j,\ldots j+r-1$, where $r\geq 1$ is its range, reads

$\displaystyle\braket{\boldsymbol{\bf
O}}=\frac{1}{Z}\sum_{\underline{s},\underline{s}^{\prime}}\bra{\underline{s}}\boldsymbol{\bf
O}\ket{\underline{s^{\prime}}}e^{\frac{\beta}{2}[E(\underline{s})+E(\underline{s}^{\prime})]}.$
(134)

We have used $s_{1}=s^{\prime}_{1}$ and $s_{L}=s^{\prime}_{L}$ to cancel the
phases coming from the coefficients (129). This is possible since spins in
positions $1$ and $L$ are not affected by the action of $\boldsymbol{\bf O}$:

$\displaystyle\begin{aligned} \bra{\underline{s}}\boldsymbol{\bf
O}\ket{\underline{s^{\prime}}}=&\prod_{\ell=1}^{j-1}\delta_{s_{\ell},s^{\prime}_{\ell}}\prod_{m=j+r}^{L}\delta_{s_{m},s^{\prime}_{m}}\\\
&\times\bra{s_{j},\ldots s_{j+r-1}}\boldsymbol{\bf
O}\ket{s^{\prime}_{j},\ldots s^{\prime}_{j+r-1}}\,.\end{aligned}$ (135)

In fact, this formula allows us to simplify Eq. (134) even further. Defining a
$2\times 2$ matrix with elements

$\displaystyle\begin{aligned}
O_{s_{j-1},s_{j+r}}=&\left(\frac{\sin\theta}{2}\right)^{r+1}\sum_{s_{j},\ldots
s_{j+r-1},\atop s^{\prime}_{j},\ldots s^{\prime}_{j+r-1}}\bra{s_{j},\ldots
s_{j+r-1}}\boldsymbol{\bf O}\ket{s^{\prime}_{j},\ldots s^{\prime}_{j+r-1}}\\\
&\times\exp\left\\{\beta\left[s_{j-1}(s_{j}+s^{\prime}_{j})+s_{j+r}(s_{j+r-1}+s^{\prime}_{j+r-1})+\sum_{\ell=j}^{j+r-2}(s_{\ell}s_{\ell+1}+s^{\prime}_{\ell}s^{\prime}_{\ell+1})\right]\right\\},\end{aligned}$
(136)

we can rewrite the expectation value of $\boldsymbol{\bf O}$ as

$\displaystyle\braket{\boldsymbol{\bf
O}}=\frac{1}{2}\sum_{s_{1},s_{L}}\left[T^{j-2}OT^{L-j-r}\right]_{s_{1},s_{L}}\,.$
(137)

Assuming $\theta\notin\\{0,\pi\\}$, the thermodynamic limit suppresses the
powers of $\cos\theta$, originating in the second term of Eq. (131). We
finally obtain

$\displaystyle\braket{\boldsymbol{\bf
O}}\\!=\\!\frac{1}{2}\\!\sum_{s_{1},s_{L}}\\!\left[\frac{1+\sigma^{x}}{2}O\frac{1+\sigma^{x}}{2}\right]_{s_{1},s_{L}}\\!=\\!\frac{1}{2}\\!\begin{pmatrix}1&1\end{pmatrix}\\!O\\!\begin{pmatrix}1\\\
1\end{pmatrix}.$ (138)

Consider now two operators $\boldsymbol{\bf O}_{1}$ and $\boldsymbol{\bf
O}_{2}$, acting on regions that are separated for $d$ sites. Using the same
argument as above, the thermodynamic limit now yields

$\displaystyle\begin{aligned} \braket{\boldsymbol{\bf O}_{1}\boldsymbol{\bf
O}_{2}}=&\left[\frac{1}{2}\\!\begin{pmatrix}1&1\end{pmatrix}\\!O_{1}\\!\begin{pmatrix}1\\\
1\end{pmatrix}\right]\left[\frac{1}{2}\\!\begin{pmatrix}1&1\end{pmatrix}\\!O_{2}\\!\begin{pmatrix}1\\\
1\end{pmatrix}\right]\\\ &+\mathcal{O}(\cos^{d}\theta)\,,\end{aligned}$ (139)

implying

$\braket{\boldsymbol{\bf O}_{1}\boldsymbol{\bf O}_{2}}-\braket{\boldsymbol{\bf
O}_{1}}\braket{\boldsymbol{\bf O}_{1}}=\mathcal{O}(\cos^{d}\theta)$ (140)

for the connected 2-point correlation function. The exponential decay of the
correlations between any distant local operators in the $\tau$-representation
signifies clustering properties of the tilted state
$\ket{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}_{\sigma}$ for even
local and even semilocal operators in the $\sigma$-representation.

Note that, for diagonal local operators $\boldsymbol{\bf O}$ (which are
written solely in terms of $\\{\boldsymbol{\bf\tau}_{j}^{z}\\}$), Eq. (134)
becomes

$\displaystyle\braket{\boldsymbol{\bf
O}}=\frac{\sum_{\underline{s}}\braket{\underline{s}}{\boldsymbol{\bf
O}}{\underline{s}}e^{\beta E(\underline{s})}}{\sum_{\underline{s}}e^{\beta
E(\underline{s})}}\,,$ (141)

i.e., a thermal expectation value. In the infinite system such operators
include, for example, $\boldsymbol{\bf\Pi}^{z}_{\sigma,+}(j)$ (the latter is
mapped into $\boldsymbol{\bf\tau}_{j}^{z}$ by the inverse duality
transformation (119)). Eq. (140) can then simply be seen as a consequence of
the absence of finite-temperature phase transitions in the one-dimensional
classical Ising model.

##  Appendix E Excess of entropy

###  E.1 Bounds on Renyi entropies

Here we prove Eq. (85) by showing

$2^{-\alpha}\mathrm{tr}(\rho^{\alpha})\leq\mathrm{tr}(\bar{\rho}^{\alpha})\leq\mathrm{tr}(\rho^{\alpha}).$
(142)

Taking the logarithm and multiplying the inequalities by $1/(1-\alpha)$ yields
Eq. (85).

We prove both inequalities in Eq. (142) using the inequality of Ault [138,
139], which says that, given arbitrary $n\times n$ complex matrices
$A_{1},A_{2},\ldots,A_{m}$, one has

$\mathrm{tr}\Big{(}\frac{B+B^{\dagger}}{2}\Big{)}\leq\frac{1}{m}\mathrm{tr}\Big{(}\sum_{j=1}^{m}(A_{j}A_{j}^{\dagger})^{\frac{m}{2}}\Big{)}$
(143)

for $B=A_{1}A_{2}\cdots A_{m}$.

To show the upper bound in (142) we write

$\mathrm{tr}{\bar{\rho}^{\alpha}}=\mathrm{tr}\big{(}\frac{\rho\bar{\rho}^{\alpha-1}+\mathrm{h.c.}}{2}\big{)},$
(144)

which follows from the definition (83) of $\bar{\rho}$, using the cyclic
property of the trace and the invariance of $\bar{\rho}$ under the Hermitian
involution $P$. The inequality of Ault for $A_{1}=\rho,\
A_{2}=A_{3}=\ldots=A_{\alpha}=\bar{\rho}$ then yields

$\mathrm{tr}\bar{\rho}^{\alpha}\leq\frac{1}{\alpha}\left[\mathrm{tr}(\rho^{\alpha})+(\alpha-1)\mathrm{tr}(\bar{\rho}^{\alpha})\right].$
(145)

Rearranging the terms proves the second inequality in Eq. (142).

To prove the lower bound in Eq. (142) we introduce projectors $P_{\pm}=(I\pm
P)/2$ and write

$\bar{\rho}=P_{+}\rho P_{+}+P_{-}\rho P_{-}.$ (146)

Note that the operators $P_{+}\rho P_{+}$ and $P_{-}\rho P_{-}$ are Hermitian
and positive semidefinite. Since $P_{+}P_{-}=0$ it follows

$\bar{\rho}^{\alpha}=\left(P_{+}\rho P_{+}\right)^{\alpha}+\left(P_{-}\rho
P_{-}\right)^{\alpha}.$ (147)

Since $\mathrm{tr}\left[\left(\rho
P_{s}\right)^{\alpha}\right]=\mathrm{tr}\left[\left(P_{s}\rho
P_{s}\right)^{\alpha}\right]$, we then see

$0\leq\mathrm{tr}\left[\left(\rho
P_{s}\right)^{\alpha}\right]\leq\mathrm{tr}\bar{\rho}^{\alpha},\quad s=\pm
1\,.$ (148)

On the other hand using $P_{+}+P_{-}=I$ we have

$\mathrm{tr}\rho^{\alpha}=\sum_{s_{1,\ldots\alpha=\pm
1}}\mathrm{tr}\left(P_{s_{1}}\rho P_{s_{2}}\rho\ldots
P_{s_{\alpha}}\rho\right)\,,$ (149)

which can be rewritten, using the cyclic property of the trace, as

$\displaystyle\begin{aligned}
\mathrm{tr}\rho^{\alpha}\\!=\\!\sum_{s_{1,\ldots\alpha}=\pm
1}\\!\mathrm{tr}\left[\\!\frac{\left(\sqrt{\rho}P_{s_{1}}\\!\sqrt{\rho}\right)\cdots\left(\sqrt{\rho}P_{s_{\alpha}}\\!\sqrt{\rho}\right)\\!+\\!\mathrm{h.c.}}{2}\\!\right]\,.\end{aligned}$
(150)

Here, $\sqrt{\rho}$ denotes the square root of the positive semidefinite
operator $\rho$. Applying the inequality of Ault with
$A_{j}=\sqrt{\rho}P_{s_{j}}\sqrt{\rho}$ for $j=1,\ldots\alpha$ now yields

$\mathrm{tr}\rho^{\alpha}\leq\sum_{s_{1,\ldots\alpha}=\pm
1}\frac{1}{\alpha}\sum_{j=1}^{\alpha}\mathrm{tr}\left[\left(\rho
P_{s_{j}}\right)^{\alpha}\right].$ (151)

Finally, using (148) we obtain

$\mathrm{tr}\rho^{\alpha}\leq 2^{\alpha}\mathrm{tr}\bar{\rho^{\alpha}},$ (152)

which proves the lower bound in Eq. (142).

###  E.2 Correlation matrix after the quench

At time $t$ the correlation matrix symbol for the initial states
$\ket{\Uparrow}_{\tau}$ and $\ket{\Downarrow}_{\tau}$ reads, respectively,

$\Gamma(p)=\pm e^{-i\mathcal{H}(p)t}\sigma^{y}e^{i\mathcal{H}(p)t},$ (153)

where
$\mathcal{H}(p)=f_{0}(p)I+f_{1}(p)\sigma^{x}+f_{2}(p)\sigma^{y}+f_{3}(p)\sigma^{z}$
is the symbol of the Hamiltonian. For the generalized XY model we have
explicitly

$\displaystyle
f_{0}(p)=2\sum_{k=1}^{\infty}(J_{x}^{(2k)}-J_{y}^{(2k)})\sin(2kp),$ (154)
$\displaystyle
f_{1}(p)=2\sum_{k=0}^{\infty}(J_{x}^{(2k+1)}-J_{y}^{(2k+1)})\sin\left[(2k+1)p\right],$
(155) $\displaystyle
f_{2}(p)=-2\sum_{k=0}^{\infty}(J_{x}^{(2k+1)}+J_{y}^{(2k+1)})\cos\left[(2k+1)p\right],$
(156) $\displaystyle
f_{3}(p)=-2\sum_{k=1}^{\infty}(J_{x}^{(2k)}+J_{y}^{(2k)})\sin(2kp)\,,$ (157)

and focusing on the reflection symmetric system ($J_{x}^{(2k)}=J_{y}^{(2k)}$
$\forall k$) we have the simplification $f_{0}(p)=0$. For $t\to\infty$ due to
dephasing the symbol can be effectively replaced by its time average, equal to

$\Gamma(p)=\pm\frac{f_{2}(p)}{\varepsilon^{2}(p)}\mathcal{H}(p),$ (158)

which gives Eq. (101) in the main text.

Let us determine the root density $\rho(p)$ starting from this expression,
since we are going to need it to derive eq. (102). It can be shown [103] that
for translationally invariant and reflection symmetric free fermionic models
the root density $\rho(p)$ is related to the correlation matrix symbol
$\Gamma(p)$ through the relation

$\begin{split}&V(p)\Gamma(p)V^{\dagger}(p)=\\\
&4\pi\rho_{o}(p)I+(4\pi\rho_{e}(p)-1)\sigma^{y}+4\pi\psi_{R}(p)\sigma^{z}-4\pi\psi_{I}(p)\sigma^{x}.\end{split}$
(159)

Here $\rho_{o}(p)$ and $\rho_{e}(p)$ are the odd and the even part of
$\rho(p)$ respectively, $\psi(p)=\psi_{R}(p)+i\psi_{I}(p)$ is a field related
to the off-diagonal elements of the density matrix, that can be neglected in
the limit $t\to\infty$, while $V(p)$ is a unitary transformation related to
Bogoliubov rotation. The action of the latter on the symbol of the Hamiltonian
reads

$V(p)\mathcal{H}(p)V(p)^{\dagger}=\varepsilon(p)\sigma^{y}.$ (160)

In the limit $t\to\infty$ we have thus simply

$\Gamma(p)=4\pi\rho_{o}(p)I+(4\pi\rho_{e}(p)-1)\frac{\mathcal{H}(p)}{\varepsilon(p)},$
(161)

so, comparing with (158), we find

$4\pi\rho(p)-1=\pm\frac{f_{2}(p)}{\varepsilon(p)},$ (162)

where we have $\rho(p)=\rho_{e}(p)$ because the state is reflection symmetric.

When the initial state is the domain wall state
$\ket{\ldots\downarrow\downarrow\downarrow_{r}\uparrow\uparrow\uparrow\ldots}_{\tau}$,
as in bipartitioning protocols in general, the system reaches stationarity
along rays $\zeta=x/t$, where $x$ is the distance from the junction [116,
103]. Along a given ray one can still describe the observables by the root
density $\rho_{\zeta}(p)$ or the correlation matrix symbol
$\Gamma_{\zeta}(p)$, through (161), which now holds along the ray. The late-
time dynamics is described by the hydrodynamic equation [117, 118, 103]

$\partial_{t}\rho_{\zeta}(p)+v(p)\partial_{x}\rho_{\zeta}(p)=0,$ (163)

where $v(p)=\varepsilon^{\prime}(p)$ is the velocity of excitations. The
equation implies that $\rho_{\zeta}(p)$ depends only on the difference
$x-v(p)t=t[\zeta-v(p)]$.

The subsystem $A$ corresponds to $\zeta=0$, since the domain wall is placed
inside. For $t\to\infty$ and $\zeta=0$ we have
$x-v(p)t\to-\mathrm{sgn}[v(p)]\infty$. However, the same value of $x-v(p)t$ is
achieved for $\zeta\to-\mathrm{sgn}[v(p)]\infty$ and $t\to\infty$. The latter
corresponds simply to the evolution starting from the state
$\ket{\Uparrow}_{\tau}$ or $\ket{\Downarrow}_{\tau}$, depending on the sign of
the velocity. From (162) we conclude that the zero ray root density is given
by

$4\pi\rho_{\zeta=0}(p)-1=-\mathrm{sgn}[v(p)]\frac{f_{2}(p)}{\varepsilon(p)}.$
(164)

Finally, from (161) we find the correlation matrix symbol for the subsystem
$A$,

$\Gamma_{\zeta=0}(p)=-\mathrm{sgn}[v(p)]\frac{f_{2}(p)}{\varepsilon(p)}I.$
(165)

###  E.3 Asymptotic results for Renyi entropies

The Renyi entropyi can be expressed as the integral [114]

$S_{\alpha}^{(\tau)}(\ell,\Psi)=\lim_{\varepsilon\to 0^{+}}\frac{1}{2\pi
i}\oint_{C_{\varepsilon}}F_{\alpha}(z)\frac{d}{dz}\log\det\left[zI-\Gamma\right]dz,$
(166)

where

$F_{\alpha}(z)=\frac{1}{2(1-\alpha)}\log\left[\left(\frac{1+z}{2}\right)^{\alpha}+\left(\frac{1-z}{2}\right)^{\alpha}\right]$
(167)

and the integration curve $C_{\varepsilon}$ goes around the interval $[-1,1]$
on the real axis, at distance $\varepsilon$ from it. The information on the
state $\ket{\Psi(t)}$ is, of course, encoded in the correlation matrix
$\Gamma$.

A block Toeplitz matrix $T_{\ell}(\mathcal{M})$ with a $2\times 2$ matrix
symbol $\mathcal{M}(p)$ is defined by the elements

$\left[T_{\ell}(\mathcal{M})\right]_{2j+m,2l+n}=\int_{-\pi}^{\pi}\mathcal{M}(p)e^{ip(j-l)}\frac{dp}{2\pi}$
(168)

for $0\leq j,l\leq\ell-1$, $0\leq m,n\leq 1$. For a piecewise continuous
symbol, such that $\det\mathcal{M}(p)$ is non-zero and has zero winding
number, and jump discontinuities at points $p_{r}$ for $r=1,2,\ldots,R$, we
use the asymptotic formula for large $\ell$,

$\begin{split}&\log\det T_{\ell}(\mathcal{M})\\\
&=\int_{-\pi}^{\pi}\log\det\mathcal{M}(p)\frac{dp}{2\pi}+\log(\ell)\sum_{r=1}^{R}\mathcal{B}_{r}+O(1).\end{split}$
(169)

Without discontinuities, the second term is not present. The leading term
grows linearly with $\ell$ and its value comes from the Szegő-Widom theorem
[127, 126]. It has been conjectured and checked numerically in [128, 129] that
with discontinuities there is a subdominant term that grows logarithmically
with $\ell$. Denoting by $\mu_{r,j}^{\pm}$ for $j=1,2$ the eigenvalues of the
limits $\lim_{p\to p_{r}^{\pm}}\mathcal{M}(p)$, the conjecture says that the
logarithmic coefficients $\mathcal{B}_{r}$ are given by

$\mathcal{B}_{r}=\frac{1}{4\pi^{2}}\sum_{j=1}^{2}\left(\log\frac{\mu_{r,j}^{-}}{\mu_{r,j}^{+}}\right)^{2}.$
(170)

Note that the formula generalizes straightforwardly to a symbol of arbitrary
size.

The asymptotic results for the Renyi entropies are obtained by using the
asymptotic formula (169) to compute the determinant of the Toeplitz matrix
$zI-\Gamma$ in (166). The latter is associated to the symbol
$\mathcal{M}(p)=zI-\Gamma(p)$, where $\Gamma(p)$ is the symbol of the
correlation matrix, given by (101) for a quench from the state of all spin up
and by (102) for a quench from a domain wall state. In this way we find that
the leading term both in (103) and in (104) is given by

$\begin{split}&a_{\alpha}=\\\
&\frac{1}{1-\alpha}\int_{-\pi}^{\pi}\log\left[\left(\frac{1+\frac{f_{2}(p)}{\varepsilon(p)}}{2}\right)^{\alpha}+\left(\frac{1-\frac{f_{2}(p)}{\varepsilon(p)}}{2}\right)^{\alpha}\right]\frac{dp}{2\pi},\end{split}$
(171)

where $f_{2}(p)=\mathrm{tr}[\sigma^{y}\mathcal{H}(p)]$ for the generalized XY
model is given by eq. (156).

The discontinuities $p_{1},p_{2},\ldots,p_{R}\in(-\pi,\pi]$ of the symbol
(102) are the zeroes of the velocity, i.e. points for which $v(p_{r})=0$. From
the asymptotic formula (169) the logarithmic coefficient in (104) is given by

$b_{\alpha}=\lim_{\varepsilon\to 0^{+}}\sum_{r=1}^{R}\frac{1}{2\pi
i}\oint_{C_{\varepsilon}}F_{\alpha}(z)\frac{d}{dz}\log\mathcal{B}_{r}(z)dz,$
(172)

where

$\mathcal{B}_{r}(z)=\frac{1}{2\pi^{2}}\left(\log\frac{z-\frac{f_{2}(p_{r})}{\varepsilon(p_{r})}}{z+\frac{f_{2}(p_{r})}{\varepsilon(p_{r})}}\right)^{2}.$
(173)

To evaluate the integral we use the partial integration,

$b_{\alpha}=-\lim_{\varepsilon\to 0^{+}}\sum_{r=1}^{R}\frac{1}{2\pi
i}\oint_{C_{\varepsilon}}\frac{d}{dz}\left[F_{\alpha}(z)\right]\mathcal{B}_{r}(z)dz,$
(174)

and deform the contour around the poles of

$\frac{d}{dz}F_{\alpha}(z)=\frac{\alpha}{2(1-\alpha)}\frac{(1+z)^{\alpha-1}-(1-z)^{\alpha-1}}{(1+z)^{\alpha}+(1-z)^{\alpha}},$
(175)

which are given by

$z_{j}=i\tan\left[\frac{(2j-1)\pi}{2\alpha}\right],\quad
j=1,2,\ldots,\alpha,\quad j\neq\frac{\alpha+1}{2}.$ (176)

As $z\to z_{j}$ we have

$\frac{d}{dz}F_{\alpha}(z)=\frac{1}{2(1-\alpha)(z-z_{j})}\big{(}1+O(z-z_{j})\big{)}$
(177)

so we obtain

$\begin{split}&b_{\alpha}=-\frac{1}{4\pi^{2}(\alpha-1)}\\\
&\times\sum_{r=1}^{R}\sum_{\begin{subarray}{c}j=1,2,\ldots,\alpha\\\
j\neq\frac{\alpha+1}{2}\end{subarray}}\left(\log\frac{i\tan\left[\frac{(2j-1)\pi}{2\alpha}\right]-\frac{f_{2}(p)}{\varepsilon(p)}}{i\tan\left[\frac{(2j-1)\pi}{2\alpha}\right]+\frac{f_{2}(p)}{\varepsilon(p)}}\right)^{2}.\end{split}$
(178)

Since the term inside the logarithm is just a phase factor, the squared
logarithm is real and negative. Because of the prefactor every term in the sum
brings thus a positive contribution to $b_{\alpha}$. For instance, for
$\alpha=2$ this expression simplifies to eq. (105) in the main text.

The velocity of excitations is given by

$v(p)=\frac{\sum_{j=1}^{3}f_{j}(p)f_{j}^{\prime}(p)}{\varepsilon(p)}.$ (179)

For the XY chain the velocity has four zeroes, given by $p_{1}=-\pi/2,\
p_{2}=0,\ p_{3}=\pi/2,\ p_{4}=\pi$. Since $f_{2}(p)=0$ for $p=p_{1},p_{3}$,
these two points do not contribute in (178). As a matter of fact, since
$f_{2}(p)=0$ there is no discontinuity in the symbol at these points. For
$p=p_{2},p_{4}$ we have $f_{2}(p)/\varepsilon(p)\in\\{-1,1\\}$ so the
contribution of these two points is the same. From (178) we find the
logarithmic coefficient (106). A general model under consideration contains
the aforementioned four zeroes of the velocity but it can contain also
additional ones. Therefore, in general the logarithmic coefficient for the XY
chain is the lower bound. Namely, $b_{\alpha}\geq b^{(XY)}_{\alpha}$.

###  E.4 Exact Procedure

In this section we show how to compute the second Renyi entropy exactly for
systems mappable to free fermions, starting from eq. (81). Let us denote
shortly $\rho_{\tau}\equiv\rho^{(\tau)}_{\tilde{A}}(\Psi)$ and
$\bar{\rho}\equiv\bar{\rho}^{(\tau)}_{\tilde{A}}(\Psi)$.

The density matrix $\rho_{\tau}$ can be expressed as a Gaussian in Majorana
fermions (99), following the standard techniques [115, 140, 119] based on Wick
theorem [141]. Namely, defining the $2(\ell+1)\times 2(\ell+1)$ correlation
matrix with the elements

$\begin{split}\Gamma_{2j+\alpha,2l+\beta}=\delta_{j,l}\delta_{\alpha,\beta}-\bra{\Psi(t)}\boldsymbol{\bf
a}_{j}^{\alpha}\boldsymbol{\bf a}_{l}^{\beta}\ket{\Psi(t)},\\\ 0\leq
j,l\leq\ell,\;1\leq\alpha,\beta\leq 2,\end{split}$ (180)

(with the usual identification $x\equiv 1,y\equiv 1$) and the matrix function

$W(\Gamma)=\log[(I+\Gamma)(I-\Gamma)^{-1}]$ (181)

we have

$\rho_{\tau}=\frac{1}{\mathcal{Z}(\Gamma)}\exp\left(\frac{\vec{a}^{\mathrm{T}}W(\Gamma)\vec{a}}{4}\right),$
(182)

where
$\vec{a}=(a_{1}^{x},a_{1}^{y},a_{2}^{x},a_{2}^{y},\ldots,a_{\ell+1}^{y})$. The
normalization is given by the function

$\mathcal{Z}(\Gamma)=\pm\left[\det\left(\frac{I-\Gamma}{2}\right)\right]^{-\frac{1}{2}},$
(183)

which can be computed from the formula (188) that we introduce later. There is
a sign ambiguity in the formula (183), which can be avoided in computing the
second Renyi entropy.

The effects of the transformations $P_{i}$ on the Gaussian in eq. (81) are
completely determined by their effects on the two-point products of Majorana
fermions, i.e. by $P_{i}a_{j}^{x,y}a_{l}^{x,y}P_{i}$. We thus find

$P_{i}\rho_{\tau}P_{i}=\frac{1}{\mathcal{Z}(\mathcal{F}_{i}(\Gamma))}\exp\left(\frac{\vec{a}^{\mathrm{T}}W(\mathcal{F}_{i}(\Gamma))\vec{a}}{4}\right),$
(184)

where we define the functions $\mathcal{F}_{i}$ by their action on the matrix
elements of the correlation matrix:

$\displaystyle\Gamma_{j,l}\stackrel{{\scriptstyle\mathcal{F}_{1}}}{{\longrightarrow}}(-1)^{\delta_{j,1}+\delta_{j,2}+\delta_{l,1}+\delta_{l,2}}\Gamma_{j,l}$
(185)
$\displaystyle\Gamma_{j,l}\stackrel{{\scriptstyle\mathcal{F}_{2}}}{{\longrightarrow}}(-1)^{\delta_{j,2\ell+1}+\delta_{j,2\ell+2}+\delta_{l,2\ell+1}+\delta_{l,2\ell+2}}\Gamma_{j,l}$
(186)
$\displaystyle\Gamma_{j,l}\stackrel{{\scriptstyle\mathcal{F}_{3}}}{{\longrightarrow}}(-1)^{\lfloor\frac{j}{2}\rfloor+\lfloor\frac{l}{2}\rfloor}\Gamma_{j,l}.$
(187)

Of course, for products of different transformations, such as in
$P_{i}P_{j}\rho_{\tau}P_{j}P_{i}$, we just have to consider the compositions
such as $\mathcal{F}_{i}\left(\mathcal{F}_{j}(\cdot)\right)$.

We have thus expressed the density matrix as a sum of Gaussian fermionic
operators, similarly to the problem of the entanglement entropy of two
disjoint blocks in the XY chain [119]. The Renyi entropies can thus be
computed by exploring the formula for the trace of a product of Gaussians,
given by [120]

$\begin{split}&\mathrm{tr}\left(e^{\frac{\vec{a}^{T}W_{1}\vec{a}}{4}}e^{\frac{\vec{a}^{T}W_{2}\vec{a}}{4}}\ldots
e^{\frac{\vec{a}^{T}W_{n}\vec{a}}{4}}\right)\\\
&=\pm\sqrt{\det\left(I+e^{W_{1}}e^{W_{2}}\ldots e^{W_{n}}\right)},\end{split}$
(188)

for antisymmetric matrices $W_{1},W_{2},\ldots,W_{n}$. The formula has a sign
ambiguity, well known in such problems, that has to be resolved in some
independent way. The normalization (183) follows from formula (188) with
$n=1$.

In computing the second Renyi entropy it is easy to avoid the sign ambiguity,
based on fact that the trace of a product of two positive semidefinite
Hermitian operators is non-negative [121]. Now, from the invariance
$P_{j}\bar{\rho}P_{j}=\bar{\rho}$, for $j=1,2,3$, that follows from (81), we
get the simplification

$\mathrm{tr}\big{[}\bar{\rho}^{2}\big{]}=\mathrm{tr}\big{[}\rho_{\tau}\bar{\rho}\big{]}.$
(189)

Using (81) again, we can write this as a sum of eight non-negative terms

$\begin{split}&\mathrm{tr}\big{[}\bar{\rho}^{2}\big{]}=\\\
&2^{-3}\sum_{j_{1},j_{2},j_{3}=0,1}\mathrm{tr}\left[\rho_{\tau}\left(P_{3}^{j_{3}}P_{2}^{j_{2}}P_{1}^{j_{1}}\rho_{\tau}P_{1}^{j_{1}}P_{2}^{j_{2}}P_{3}^{j_{3}}\right)\right].\end{split}$
(190)

Each of the eight terms can be evaluated using (184) and the formula (188) for
the trace of a product of Gaussians. Using the identity

$\frac{\det\left(I+e^{W(\Gamma_{1})}e^{W(\Gamma_{2})}\right)}{\det\left(I+e^{W(\Gamma_{1})}\right)\det\left(I+e^{W(\Gamma_{2})}\right)}=\det\left(\frac{I+\Gamma_{1}\Gamma_{2}}{2}\right),$
(191)

and eq. (82) we get then the second Renyi entropy,

$\begin{split}&S_{2}(\ell,\Psi)=2\log
2-\log\bigg{[}\mathcal{G}\big{(}\Gamma,\Gamma\big{)}\\\
&+\sum_{i=1}^{3}\mathcal{G}\big{(}\Gamma,\mathcal{F}_{i}\left(\Gamma\right)\big{)}+\sum_{1\leq
i<j\leq
3}\mathcal{G}\big{(}\Gamma,\mathcal{F}_{i}\left(\mathcal{F}_{j}(\Gamma)\right)\big{)}\\\
&+\mathcal{G}\big{(}\Gamma,\mathcal{F}_{1}\left(\mathcal{F}_{2}\left(\mathcal{F}_{3}\left(\Gamma\right)\right)\right)\big{)}\bigg{]},\end{split}$
(192)

where we have defined the function

$\mathcal{G}(\Gamma_{1},\Gamma_{2})=\sqrt{\det\left(\frac{I+\Gamma_{1}\Gamma_{2}}{2}\right)}.$
(193)

##  Appendix F Melting of the order and scaling limits

This section provides details on the calculation of expectation values of
local observables (e.g., strings of $\boldsymbol{\bf\sigma}^{z}_{j}$) in a
time-evolved weakly tilted state
$\ket{\rotatebox[origin={c}]{45.0}{$\Rightarrow$}_{\theta}}$. The first part
of the section discusses the effective density matrix given in Eq. (111), from
which expectation values of semilocal order parameters can be calculated; the
latter are discussed in the second part.

###  F.1 Effective density matrix

In the fermionic representation we have

$\displaystyle\begin{aligned} e^{-i\boldsymbol{\bf
H}^{(\tau)}t}\boldsymbol{\bf\tau}_{j}^{z}e^{i\boldsymbol{\bf
H}^{(\tau)}t}\\!=\\!&\frac{1}{2}\\!\sum_{j_{1},j_{2}}\\!\int_{-\pi}^{\pi}\\!\frac{{\rm
d}^{2}k}{(2\pi)^{2}}e^{i[k_{1}(j_{1}-j)+k_{2}(j-j_{2})]}\\\
&\times\vec{\boldsymbol{\bf
a}}^{T}_{j_{1}}\\!\cdot\\!\left[e^{-i\mathcal{H}(k_{1})t}\sigma^{y}e^{i\mathcal{H}(k_{2})t}\right]\\!\cdot\\!\vec{\boldsymbol{\bf
a}}_{j_{2}}\,,\end{aligned}$ (194)

where $\vec{\boldsymbol{\bf a}}_{j}=\begin{pmatrix}\boldsymbol{\bf
a}^{x}_{j}&\boldsymbol{\bf a}^{y}_{j}\end{pmatrix}^{T}$ and $\mathcal{H}(k)$
is the symbol of the Hamiltonian, defined in Eq. (23). Using this
representation in the expression (111) for the effective density matrix,
performing the sums over $j$ and $n$, and changing the integration variables
as $k_{1}\to k_{2}-k_{1}$, $k_{2}\to k_{2}$, $k_{3}\to k_{3}$, $k_{4}\to
k_{3}-k_{4}$, we obtain

$\displaystyle\begin{aligned} \boldsymbol{\bf\rho}_{\rm
eff}(t)\\!\sim\\!&\frac{1}{4\mathrm{tr}\boldsymbol{\bf
I}}\\!\sum_{j_{1,2,3,4}}\\!\int\\!\frac{{\rm
d}^{3}\\!k}{(2\pi)^{3}}\\!\frac{e^{i[(j_{4}\\!-\\!j_{1})k_{1}\\!+\\!(j_{1}\\!-\\!j_{2})k_{2}\\!+\\!(j_{3}\\!-\\!j_{4})k_{3}]}}{e^{\frac{\theta^{2}}{2}+ik_{1}}\\!-\\!1}\\\
&\times\vec{\boldsymbol{\bf
a}}^{T}_{j_{1}}\\!\cdot\\![e^{-i\mathcal{H}(k_{2}-k_{1})t}\sigma^{y}e^{i\mathcal{H}(k_{2})t}]\\!\cdot\\!\vec{\boldsymbol{\bf
a}}_{j_{2}}\\\ &\times\vec{\boldsymbol{\bf
a}}^{T}_{j_{3}}\cdot\\![e^{-i\mathcal{H}(k_{3})t}\sigma^{y}e^{i\mathcal{H}(k_{3}-k_{1})t}]\\!\cdot\\!\vec{\boldsymbol{\bf
a}}_{j_{4}}\,,\end{aligned}$ (195)

for the effective density matrix.

Being interested in the expectation value of local operators at large times
after a weak symmetry breaking, we will assume that $t$ and $\theta$ are the
only large and small parameter, respectively, in Eq. (195). We can assume this
because local operator $\boldsymbol{\bf O}$ in
$\mathrm{tr}[\boldsymbol{\bf\rho}_{\rm eff}(t)\boldsymbol{\bf O}]$ is
characterised by indices $j_{1}$, $j_{2}$, $j_{3}$, and $j_{4}$, confined to
some finite interval. In order to obtain the asymptotic behaviour of the
integral in Eq. (195) one now has to consider the pole in $k_{1}$. The first
step consists of deforming the integration path in $k_{1}$ into a piecewise
linear curve with lines parallel to the axes, in such a way that the
horizontal lines give exponentially decaying contributions.

(a)$k_{1}^{\rm R}$$k_{1}^{\rm I}$$-\pi$$\pi$$i\frac{\theta^{2}}{2}$

(b)$k_{1}^{\rm R}$$k_{1}^{\rm I}$$-\pi$$\pi$$i\frac{\theta^{2}}{2}$

Figure 16: The integration path over the interval $[-\pi,\pi]$ on the real
axis is deformed depending on the sign of $v(k_{3}-k_{1}^{\rm
R})s_{2}-v(k_{2}-k_{1}^{\rm R})s_{1}$, to make the exponent in (199) negative
on the horizontal lines. Panel (a): A contour deformed without crossing the
pole. Panel (b): Deforming the contour across the pole, its contribution has
to be included.

We start from the relation

$e^{-i\mathcal{H}(k_{2}-k_{1})t}=\sum_{s_{1}=\pm
1}P_{s_{1}}(k_{2}-k_{1})e^{-i\varepsilon(k_{2}-k_{1})s_{1}t},$ (196)

where $P_{s}(k)$ are the projectors on the eigenstates of $\mathcal{H}(k)$
with eigenvalues $\pm\varepsilon(k)$:

$P_{s}(k)=\frac{1}{2}\Big{(}I+s\frac{\mathcal{H}(k)}{\varepsilon(k)}\Big{)}\,.$
(197)

Denoting the real and imaginary part of $k_{1}$ by $k_{1}^{\rm R},k_{1}^{\rm
I}$ respectively, assuming the imaginary part $k_{1}^{\rm I}$ to be small, and
expanding in it, we obtain

$\displaystyle\begin{aligned}
e^{-i\mathcal{H}\\!(k_{2}\\!-\\!k_{1})t}\\!&\sim\\!\sum_{s_{1}}\\!P_{s_{1}}\\!(k_{2}\\!-\\!k_{1})e^{-i\varepsilon(k_{2}\\!-\\!k_{1}^{\rm
R})s_{1}t\\!-\\!k_{1}^{\rm I}v(k_{2}\\!-\\!k_{1}^{\rm R})s_{1}t},\\\
e^{i\mathcal{H}\\!(k_{3}\\!-\\!k_{1})t}\\!&\sim\\!\sum_{s_{2}}\\!P_{s_{2}}\\!(k_{3}\\!-\\!k_{1})e^{i\varepsilon(k_{3}\\!-\\!k_{1}^{\rm
R})s_{2}t+k_{1}^{\rm I}v(k_{3}\\!-\\!k_{1}^{\rm R})s_{2}t}.\end{aligned}$
(198)

where $v(k):=\varepsilon^{\prime}(k)$. The expressions in each line of Eq.
(198) are multiplied in Eq. (195). In particular, the real parts of each
exponential in Eq. (198) together yield a factor

$e^{k_{1}^{\rm I}\left[v(k_{3}-k_{1}^{\rm R})s_{2}-v(k_{2}-k_{1}^{\rm
R})s_{1}\right]t}.$ (199)

Assuming $\theta^{2}t$ to be fixed, we deform the integration path in $k_{1}$
so that $\theta^{2}\ll|k_{1}^{\rm I}|\ll 1$ and we choose the sign of
$k_{1}^{\rm I}$ so that the factor multiplying $t$ in Eq. (199) is always
negative and the horizontal lines in the contour are exponentially suppressed
(the choice of the integration contour thus depends also on $s_{1},s_{2}$).
The contours that we work with are therefore determined by the sign of
$v(k_{3}-k_{1}^{\rm R})s_{2}-v(k_{2}-k_{1}^{\rm R})s_{1}$, and are presented
in Fig. 16.

The pole is at $k_{1}^{\rm R}=0$ and there are two possible placements of the
obtained contour with respect to the pole, depending on the exponent in eq.
(199):

  1. 1.

For $v(k_{3})s_{2}-v(k_{2})s_{1}>0$ we have $k_{1}^{\rm I}<0$ on the deformed
contour when $k_{1}^{\rm R}=0$. In this case we do not have a pole
contribution — see Fig. 16(a).

  2. 2.

For $v(k_{3})s_{2}-v(k_{2})s_{1}<0$ we have $k_{1}^{\rm I}>0$ on the deformed
contour when $k_{1}^{\rm R}=0$. In such a case the contour has to be extended
across the pole — see Fig. 16(b). In order for the integral along the extended
contour to reproduce the correct result, the pole contribution has to be
subtracted.

In the first case the integral is exponentially suppressed along the entire
contour and it yields zero. In the second case only the subtracted pole
contribution remains of the integral. This can be described simply by
including the factor $\theta_{\rm H}(v(k_{2})s_{1}-v(k_{3})s_{2})$ in the
integral, where $\theta_{\rm H}$ is the Heaviside step function. Assuming
$\theta^{2}(\theta^{2}t)\to 0$, we then include into Eq. (195) the pole
contribution

$\displaystyle\begin{aligned}
e^{-i\mathcal{H}(k_{2}-i\frac{\theta^{2}}{2})t}\\!=\\!&\sum_{s_{1}=\pm
1}\\!P_{s_{1}}\\!(k_{2})e^{-i\mathcal{H}(k_{2})s_{1}t\\!-\\!v(k_{2})\frac{\theta^{2}}{2}s_{1}t}\\\
&\times(1+\mathcal{O}(\theta^{2}))\,,\end{aligned}$ (200)

and an analogous one for $\exp[i\mathcal{H}(k_{3}-i\theta^{2}/2)t]$, to get

$\begin{split}\boldsymbol{\bf\rho}_{\rm
eff}(t)\\!\sim\\!&\frac{1}{4\,\mathrm{tr}\boldsymbol{\bf
I}}\\!\sum_{j_{1},j_{2},\atop j_{3},j_{4}}\\!\int\\!\frac{{\rm d}k_{2}{\rm
d}k_{3}}{(2\pi)^{3}}\\!\sum_{s_{1},s_{2}=\pm
1}e^{i[(j_{1}\\!-\\!j_{2})k_{2}\\!+\\!(j_{3}\\!-\\!j_{4})k_{3}]}\\\ &\times
e^{-\frac{\theta^{2}}{2}t(v(k_{2})s_{1}\\!-\\!v(k_{3})s_{2})}\theta_{H}(v(k_{2})s_{1}-v(k_{3})s_{2})\\\
&\times\vec{\boldsymbol{\bf
a}}^{T}_{j_{1}}\cdot[P_{s_{1}}(k_{2})\Gamma_{t}(k_{2})]\cdot\vec{\boldsymbol{\bf
a}}_{j_{2}}\\\ &\times\vec{\boldsymbol{\bf
a}}^{T}_{j_{3}}\cdot[P_{s_{2}}(k_{3})\Gamma_{t}(k_{3})]\cdot\vec{\boldsymbol{\bf
a}}_{j_{4}}\,.\end{split}$ (201)

Here we denoted
$\Gamma_{t}(k)=e^{-i\mathcal{H}(k)t}\sigma^{y}e^{i\mathcal{H}(k)t}$. In the
limit $t\to\infty$ we can replace the latter by its time average

$\bar{\Gamma}(k)=-\frac{2(J_{x}+J_{y})\cos
k}{\varepsilon^{2}(k)}\mathcal{H}(k)\,,$ (202)

ending up with

$\displaystyle\begin{aligned} \boldsymbol{\bf\rho}_{\rm
eff}(t)\\!\sim\\!&\frac{1}{\mathrm{tr}\boldsymbol{\bf I}}\\!\int\\!\frac{{\rm
d}^{2}k}{(2\pi)^{2}}\frac{(J_{x}\\!+\\!J_{y})^{2}\cos k_{1}\cos
k_{2}}{\varepsilon(k_{1})\varepsilon(k_{2})}\\\ &\times\sum_{s=\pm
1}se^{-|v(k_{1})\\!-\\!sv(k_{2})|\frac{\theta^{2}t}{2}}\\\
&\times\sum_{j_{1},j_{2}}e^{i(j_{1}\\!-\\!j_{2})k_{1}}\vec{\boldsymbol{\bf
a}}^{T}_{j_{1}}\cdot
P_{\mathrm{sgn}[v(k_{1})\\!-\\!sv(k_{2})]}(k_{1})\cdot\vec{\boldsymbol{\bf
a}}_{j_{2}}\\\
&\times\sum_{j_{3},j_{4}}e^{i(j_{3}\\!-\\!j_{4})k_{2}}\vec{\boldsymbol{\bf
a}}^{T}_{j_{3}}\cdot
P_{\mathrm{sgn}[sv(k_{1})\\!-\\!v(k_{2})]}(k_{2})\cdot\vec{\boldsymbol{\bf
a}}_{j_{4}}\,.\end{aligned}$ (203)

For sufficiently large $\theta^{2}t$ this generically decays as
$1/(\theta^{2}t)$.

###  F.2 Semilocal order parameter

Having obtained the effective reduced density matrix, we are now ready to
investigate the expectation values of the string
$\prod_{\ell=0}^{n-1}\boldsymbol{\bf\sigma}_{j+\ell}^{z}$, which corresponds
to a product of four Majorana fermions. For $n=1$ this is the simplest local
operator sensitive to the semilocal order [31], while in the limit
$n\rightarrow\infty$ it can be regarded as a string-order parameter.

Let us then consider
$Z_{n}=\braket{\prod_{\ell=0}^{n-1}\boldsymbol{\bf\sigma}_{j+\ell}^{z}}=-\braket{\boldsymbol{\bf
a}_{1}^{x}\boldsymbol{\bf a}_{1}^{y}\boldsymbol{\bf
a}^{x}_{n+1}\boldsymbol{\bf a}^{y}_{n+1}}$. This can be computed from Eq.
(203) using Wick’s theorem:

$\displaystyle\begin{aligned} Z_{n}\\!=\\!&\int\\!\frac{{\rm
d}^{2}k}{(2\pi)^{2}}\frac{2(J_{x}\\!+\\!J_{y})^{2}\cos k_{1}\cos
k_{2}}{\varepsilon(k_{1})\varepsilon(k_{2})}\sum_{s,s^{\prime}=\pm
1}se^{-|v(k_{1})-sv(k_{2})|\frac{\theta^{2}t}{2}}\sin^{2}\Bigl{[}n\frac{k_{1}+s^{\prime}k_{2}}{2}\Bigr{]}\\\
&\times\Bigl{(}-s^{\prime}+4s\frac{(J_{x}+J_{y})^{2}\cos k_{1}\cos
k_{2}+s^{\prime}(J_{x}-J_{y})^{2}\sin k_{1}\sin
k_{2}}{\varepsilon(k_{1})\varepsilon(k_{2})}\Bigr{)}.\end{aligned}$ (204)

This expression simplifies in the limit $n\rightarrow\infty$ (bear in mind
that this limit has been taken after the limit of large time and small angle
at fixed $\theta^{2}t$), in which we obtain Eq. (112).

## References

  * Chen  _et al._ [2011a] X. Chen, Z.-C. Gu, and X.-G. Wen, Phys. Rev. B 83, 035107 (2011a).
  * Wen [1990] X. G. Wen, Phys. Rev. Lett. 64, 2206 (1990).
  * Wen [1991a] X. G. Wen, Phys. Rev. B 43, 11025 (1991a).
  * Wen [1991b] X. G. Wen, Phys. Rev. Lett. 66, 802 (1991b).
  * Chen  _et al._ [2011b] X. Chen, Z.-X. Liu, and X.-G. Wen, Phys. Rev. B 84, 235141 (2011b).
  * Fendley [2016] P. Fendley, Journal of Physics A: Mathematical and Theoretical 49, 30LT01 (2016).
  * Scaffidi  _et al._ [2017] T. Scaffidi, D. E. Parker, and R. Vasseur, Phys. Rev. X 7, 041048 (2017).
  * Pollmann  _et al._ [2010] F. Pollmann, A. M. Turner, E. Berg, and M. Oshikawa, Phys. Rev. B 81, 064439 (2010).
  * den Nijs and Rommelse [1989] M. den Nijs and K. Rommelse, Phys. Rev. B 40, 4709 (1989).
  * Pollmann and Turner [2012] F. Pollmann and A. M. Turner, Phys. Rev. B 86, 125441 (2012).
  * de Groot  _et al._ [2021] C. de Groot, A. Turzillo, and N. Schuch, Symmetry protected topological order in open quantum systems (2021).
  * Roberts  _et al._ [2017] S. Roberts, B. Yoshida, A. Kubica, and S. D. Bartlett, Phys. Rev. A 96, 022306 (2017).
  * Tsomokos  _et al._ [2009] D. I. Tsomokos, A. Hamma, W. Zhang, S. Haas, and R. Fazio, Phys. Rev. A 80, 060302 (2009).
  * Rahmani and Chamon [2010] A. Rahmani and C. Chamon, Phys. Rev. B 82, 134303 (2010).
  * Halász and Hamma [2013] G. B. Halász and A. Hamma, Phys. Rev. Lett. 110, 170605 (2013).
  * Mazza  _et al._ [2014] L. Mazza, D. Rossini, M. Endres, and R. Fazio, Phys. Rev. B 90, 020301 (2014).
  * Wang  _et al._ [2015] P. Wang, W. Yi, and G. Xianlong, New Journal of Physics 17, 013029 (2015).
  * Calvanese Strinati  _et al._ [2016] M. Calvanese Strinati, L. Mazza, M. Endres, D. Rossini, and R. Fazio, Phys. Rev. B 94, 024302 (2016).
  * Zeng  _et al._ [2016] Y. Zeng, A. Hamma, and H. Fan, Phys. Rev. B 94, 125104 (2016).
  * Rademaker [2019] L. Rademaker, SciPost Phys. 7, 71 (2019).
  * Polkovnikov  _et al._ [2011] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
  * Eisert  _et al._ [2015] J. Eisert, M. Friesdorf, and C. Gogolin, Nature Phys 11, 124 (2015).
  * Gogolin and Eisert [2016] C. Gogolin and J. Eisert, Rep. Prog. Phys. 79, 056001 (2016).
  * Calabrese and Cardy [2016] P. Calabrese and J. Cardy, J. Stat. Mech. 2016, 064003 (2016).
  * Essler and Fagotti [2016] F. H. L. Essler and M. Fagotti, J. Stat. Mech. 2016, 064002 (2016).
  * Cazalilla and Chung [2016] M. A. Cazalilla and M.-C. Chung, Journal of Statistical Mechanics: Theory and Experiment 2016, 064004 (2016).
  * Caux [2016] J.-S. Caux, Journal of Statistical Mechanics: Theory and Experiment 2016, 064006 (2016).
  * Vidmar and Rigol [2016] L. Vidmar and M. Rigol, J. Stat. Mech. 2016, 064007 (2016).
  * Ilievski  _et al._ [2015a] E. Ilievski, J. De Nardis, B. Wouters, J.-S. Caux, F. H. L. Essler, and T. Prosen, Phys. Rev. Lett. 115, 157201 (2015a).
  * Ilievski  _et al._ [2016] E. Ilievski, M. Medenjak, T. Prosen, and L. Zadnik, J. Stat. Mech. 2016, 064008 (2016).
  * Fagotti [2022] M. Fagotti, Phys. Rev. Lett. 128, 110602 (2022).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008).
  * Deutsch [1991] J. M. Deutsch, Phys. Rev. A 43, 2046 (1991).
  * Srednicki [1994] M. Srednicki, Phys. Rev. E 50, 888 (1994).
  * Deutsch [2018] J. M. Deutsch, Reports on Progress in Physics 81, 082001 (2018).
  * Rigol  _et al._ [2007] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Phys. Rev. Lett. 98, 050405 (2007).
  * Zadnik and Fagotti [2021] L. Zadnik and M. Fagotti, SciPost Phys. Core 4, 10 (2021).
  * Eisler and Maislinger [2018] V. Eisler and F. Maislinger, Phys. Rev. B 98, 161117 (2018).
  * Eisler and Maislinger [2020] V. Eisler and F. Maislinger, SciPost Phys. 8, 37 (2020).
  * Bidzhiev  _et al._ [2021] K. Bidzhiev, M. Fagotti, and L. Zadnik, arXiv:2108.03230 (2021).
  * Zadnik  _et al._ [2021] L. Zadnik, S. Bocini, K. Bidzhiev, and M. Fagotti, arXiv:2111.06325 (2021).
  * Gruber and Eisler [2021] M. Gruber and V. Eisler, SciPost Phys. 10, 5 (2021).
  * Eisler [2021] V. Eisler, J. Phys. A: Math. Theor. 54, 424002 (2021).
  * Kitaev and Preskill [2006] A. Kitaev and J. Preskill, Phys. Rev. Lett. 96, 110404 (2006).
  * Bertini  _et al._ [2015] B. Bertini, F. H. L. Essler, S. Groha, and N. J. Robinson, Phys. Rev. Lett. 115, 180601 (2015).
  * Bertini and Fagotti [2015] B. Bertini and M. Fagotti, Journal of Statistical Mechanics: Theory and Experiment 2015, P07012 (2015).
  * Alba and Fagotti [2017] V. Alba and M. Fagotti, Phys. Rev. Lett. 119, 010601 (2017).
  * Durnin  _et al._ [2021] J. Durnin, M. J. Bhaseen, and B. Doyon, Phys. Rev. Lett. 127, 130601 (2021).
  * Bratteli and Robinson [1997] O. Bratteli and D. W. Robinson,  _Operator Algebras and Quantum Statistical Mechanics 1 & 2_, 2nd ed. (Berlin: Springer-Verlag, 1997).
  * Ilievski and Prosen [2013] E. Ilievski and T. Prosen, Communications in Mathematical Physics 318, 809 (2013).
  * Prosen [2014] T. Prosen, Nuclear Physics B 886, 1177 (2014).
  * Doyon [2017] B. Doyon, Commun. Math. Phys. 351, 155 (2017).
  * Lieb and Robinson [1972] E. H. Lieb and D. W. Robinson, Commun. Math. Phys. 28, 251 (1972).
  * Bravyi  _et al._ [2006] S. Bravyi, M. B. Hastings, and F. Verstraete, Phys. Rev. Lett. 97, 050401 (2006).
  * Else  _et al._ [2013] D. V. Else, S. D. Bartlett, and A. C. Doherty, Phys. Rev. B 88, 085114 (2013).
  * Duivenvoorden and Quella [2013] K. Duivenvoorden and T. Quella, Phys. Rev. B 88, 125115 (2013).
  * Beekman  _et al._ [2019] A. J. Beekman, L. Rademaker, and J. van Wezel, SciPost Phys. Lect. Notes , 11 (2019).
  * Zeng  _et al._ [2019] B. Zeng, X. Chen, D.-L. Zhou, and X.-G. Wen,  _Quantum Information Meets Quantum Matter_ (Springer, New York, 2019).
  * Jaynes [1957] E. T. Jaynes, Phys. Rev. 106, 620 (1957).
  * Calabrese and Cardy [2004] P. Calabrese and J. Cardy, Journal of Statistical Mechanics: Theory and Experiment 2004, P06002 (2004).
  * Note [1] Although a discrete symmetry could in principle be broken, we are only aware of results pointing at its impossibility in specific classes of models.
  * Note [2] There is no standard definition of “relaxation time”, but our considerations are quite independent of the details of the definition.
  * Calabrese and Cardy [2005] P. Calabrese and J. Cardy, J. Stat. Mech. 2005, P04010 (2005).
  * Bocchieri and Loinger [1957] P. Bocchieri and A. Loinger, Phys. Rev. 107, 337 (1957).
  * Cassidy  _et al._ [2011] A. C. Cassidy, C. W. Clark, and M. Rigol, Phys. Rev. Lett. 106, 140405 (2011).
  * Caux and Essler [2013] J.-S. Caux and F. H. L. Essler, Phys. Rev. Lett. 110, 257203 (2013).
  * Fagotti and Essler [2013a] M. Fagotti and F. H. L. Essler, Phys. Rev. B 87, 245107 (2013a).
  * Wouters  _et al._ [2014] B. Wouters, J. De Nardis, M. Brockmann, D. Fioretto, M. Rigol, and J.-S. Caux, Phys. Rev. Lett. 113, 117202 (2014).
  * Pozsgay  _et al._ [2014] B. Pozsgay, M. Mestyán, M. A. Werner, M. Kormos, G. Zaránd, and G. Takács, Phys. Rev. Lett. 113, 117203 (2014).
  * Prosen and Ilievski [2013] T. Prosen and E. Ilievski, Phys. Rev. Lett. 111, 057203 (2013).
  * Pereira  _et al._ [2014] R. G. Pereira, V. Pasquier, J. Sirker, and I. Affleck, J. Stat. Mech. 2014, P09037 (2014).
  * Ilievski  _et al._ [2015b] E. Ilievski, M. Medenjak, and T. Prosen, Phys. Rev. Lett. 115, 120601 (2015b).
  * Ilievski  _et al._ [2017] E. Ilievski, E. Quinn, and J.-S. Caux, Phys. Rev. B 95, 115128 (2017).
  * Pozsgay  _et al._ [2017] B. Pozsgay, E. Vernier, and M. A. Werner, Journal of Statistical Mechanics: Theory and Experiment 2017, 093103 (2017).
  * Ilievski and Quinn [2019] E. Ilievski and E. Quinn, SciPost Phys. 7, 33 (2019).
  * Rossini  _et al._ [2009] D. Rossini, A. Silva, G. Mussardo, and G. E. Santoro, Phys. Rev. Lett. 102, 127204 (2009).
  * Note [3] Suppose instead that $\mathinner{\langle{\boldsymbol{\bf K},\boldsymbol{\bf O}}\rangle}_{\lambda}=0$ for all local operators $\bf O$ (including $\boldsymbol{\bf O}=\boldsymbol{\bf h}_{j}$). Then condition $\mathinner{\langle{-\beta^{\prime}(\lambda)\boldsymbol{\bf H}+\boldsymbol{\bf K},\boldsymbol{\bf h}_{j}}\rangle}_{\lambda}=0$ would imply either $\mathinner{\langle{\boldsymbol{\bf H},\boldsymbol{\bf h}_{j}}\rangle}=0$, which is clearly false, or $\beta^{\prime}(\lambda)=0$, the latter leading to $\delta\mathinner{\langle{\boldsymbol{\bf O}}\rangle}_{\lambda}=0$ in Eq. (14).
  * Note [4] We consider the entropy per site rather than the entropy because the latter is infinite in the thermodynamic limit.
  * Note [5] If $\mathinner{\langle{\boldsymbol{\bf K},\boldsymbol{\bf k}_{j}}\rangle}$ was infinite $\beta^{\prime}(\lambda)$ would have to compensate it in order for $\delta\mathfrak{s}[\boldsymbol{\bf\rho}_{\infty}(\lambda)]$ to remain finite. The entropy density would then not be a smooth function of $\lambda$.
  * Lieb  _et al._ [1961] E. Lieb, T. Schultz, and D. Mattis, Annals of Physics 16, 407 (1961).
  * Note [6] We remind the reader that the half-infinite string of Pauli matrices is not well defined in an infinite chain; this is however irrelevant as long as only operators quadratic in Majorana fermions are considered. For a rigorous definition of Jordan-Wigner transformation see instead Ref. [91]. We will use the latter definition in Section IV.1 where we discuss duality transformation and semilocal operators.
  * Fagotti [2014] M. Fagotti, J. Stat. Mech. 2014, P03016 (2014).
  * Calabrese  _et al._ [2012] P. Calabrese, F. H. L. Essler, and M. Fagotti, Journal of Statistical Mechanics: Theory and Experiment 2012, P07022 (2012).
  * Murthy and Srednicki [2019] C. Murthy and M. Srednicki, Phys. Rev. E 100, 012146 (2019).
  * Gluza  _et al._ [2019] M. Gluza, J. Eisert, and T. Farrelly, SciPost Phys. 7, 38 (2019).
  * Smacchia  _et al._ [2011] P. Smacchia, L. Amico, P. Facchi, R. Fazio, G. Florio, S. Pascazio, and V. Vedral, Phys. Rev. A 84, 022304 (2011).
  * Gombor and Pozsgay [2021] T. Gombor and B. Pozsgay, Phys. Rev. E 104, 054123 (2021).
  * Fagotti and Essler [2013b] M. Fagotti and F. H. L. Essler, J. Stat. Mech. 2013, P07012 (2013b).
  * Pozsgay [2013] B. Pozsgay, J. Stat. Mech. 2013, P07003 (2013).
  * Note [7] In fact, it will later become clear that the time evolution of $\mathinner{\langle{\boldsymbol{\bf\sigma}^{z}_{\ell}}\rangle}$ in the dual XY model is exactly the same as that of $\mathinner{\langle{\boldsymbol{\bf\sigma}^{z}_{\ell}\boldsymbol{\bf\sigma}^{z}_{\ell+1}}\rangle}$ in the XY model shown in Fig. 7.
  * Araki and Matsui [1985] H. Araki and T. Matsui, Communications in Mathematical Physics 101, 213 (1985).
  * Weinberg [1996] S. Weinberg,  _The Quantum Theory of Fields_, Vol. 2 (Cambridge University Press, 1996).
  * Kennedy and Tasaki [1992a] T. Kennedy and H. Tasaki, Phys. Rev. B 45, 304 (1992a).
  * Kennedy and Tasaki [1992b] T. Kennedy and H. Tasaki, Communications in Mathematical Physics 147, 431 (1992b).
  * Note [8] As far as we can see, semilocal terms in $\boldsymbol{\bf W}(\tau)$, can always be replaced by local terms.
  * Grady [1982] M. Grady, Phys. Rev. D 25, 1103 (1982).
  * Note [9] Note that the projected ensemble (64) gives correct predictions for the stationary expectation values of quasilocal operators, as can be seen by mapping it into the $\tau$ representation.
  * Pozsgay [2014] B. Pozsgay, Journal of Statistical Mechanics: Theory and Experiment 2014, P09026 (2014).
  * Müller and Shrock [1985] G. Müller and R. E. Shrock, Phys. Rev. B 32, 5845 (1985).
  * Note [10] A $\mathbbm{Z}_{2}\times\mathbbm{Z}_{2}$-semilocal charge would be mapped by $\mathcal{D}_{\mathbbm{Z}_{2}}^{-1}$ into a $\mathbbm{Z}_{2}$-semilocal charge in the $\tau$ representation, but there are no such charges in the latter — see Section IV.3.1.
  * Shohat  _et al._ [1943] J. Shohat, J. Tamarkin, and A. Society,  _The Problem of Moments_, Mathematical Surveys and Monographs (American Mathematical Society, 1943).
  * Bengtsson and Zyczkowski [2006] I. Bengtsson and K. Zyczkowski,  _Geometry of Quantum States: An Introduction to Quantum Entanglement_ (Cambridge University Press, 2006).
  * Alba  _et al._ [2021] V. Alba, B. Bertini, M. Fagotti, L. Piroli, and P. Ruggiero, Journal of Statistical Mechanics: Theory and Experiment 2021, 114004 (2021).
  * Nardis  _et al._ [2022] J. D. Nardis, B. Doyon, M. Medenjak, and M. Panfil, Journal of Statistical Mechanics: Theory and Experiment 2022, 014002 (2022).
  * Bulchandani  _et al._ [2021] V. B. Bulchandani, S. Gopalakrishnan, and E. Ilievski, Journal of Statistical Mechanics: Theory and Experiment 2021, 084001 (2021).
  * Bouchoule and Dubail [2022] I. Bouchoule and J. Dubail, Journal of Statistical Mechanics: Theory and Experiment 2022, 014003 (2022).
  * Bastianello  _et al._ [2021] A. Bastianello, A. D. Luca, and R. Vasseur, Journal of Statistical Mechanics: Theory and Experiment 2021, 114003 (2021).
  * Borsi  _et al._ [2021] M. Borsi, B. Pozsgay, and L. Pristyák, Journal of Statistical Mechanics: Theory and Experiment 2021, 094001 (2021).
  * Note [11] The existence of odd charges implies that the stationary state after the quench from $\mathinner{|{\Uparrow}\rangle}_{\tau}$ is macroscopically different from that after the quench from $\mathinner{|{\Downarrow}\rangle}_{\tau}$; by symmetry, the nonequilibrium stationary state emerging from $\mathinner{|{\cdots\downarrow_{r-1}\downarrow_{r}\uparrow_{r+1}\uparrow_{r+2}\cdots}\rangle}_{\tau}$ has to be macroscopically different from the two of them.
  * Doyon  _et al._ [2017] B. Doyon, J. Dubail, R. Konik, and T. Yoshimura, Phys. Rev. Lett. 119, 195301 (2017).
  * Dubail  _et al._ [2017] J. Dubail, J.-M. Stéphan, J. Viti, and P. Calabrese, SciPost Phys. 2, 002 (2017).
  * Calabrese and Cardy [2009] P. Calabrese and J. Cardy, Journal of Physics A: Mathematical and Theoretical 42, 504005 (2009).
  * Suzuki [1971] M. Suzuki, Physics Letters A 34, 338 (1971).
  * Jin and Korepin [2004] B.-Q. Jin and V. E. Korepin, Journal of Statistical Physics 116, 79 (2004).
  * Vidal  _et al._ [2003] G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Phys. Rev. Lett. 90, 227902 (2003).
  * Bertini and Fagotti [2016] B. Bertini and M. Fagotti, Phys. Rev. Lett. 117, 130402 (2016).
  * Bertini  _et al._ [2016] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, Phys. Rev. Lett. 117, 207201 (2016).
  * Castro-Alvaredo  _et al._ [2016] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, Phys. Rev. X 6, 041065 (2016).
  * Fagotti and Calabrese [2010] M. Fagotti and P. Calabrese, Journal of Statistical Mechanics: Theory and Experiment 2010, P04016 (2010).
  * Klich [2014] I. Klich, Journal of Statistical Mechanics: Theory and Experiment 2014, P11006 (2014).
  * Bhatia [1997] R. Bhatia, Variational principles for eigenvalues, in  _Matrix Analysis_ (Springer New York, New York, NY, 1997) pp. 57–83.
  * Fagotti and Calabrese [2008] M. Fagotti and P. Calabrese, Phys. Rev. A 78, 010306 (2008).
  * Alba and Calabrese [2017] V. Alba and P. Calabrese, Proceedings of the National Academy of Sciences 114, 7947 (2017), https://www.pnas.org/doi/pdf/10.1073/pnas.1703516114 .
  * Alba  _et al._ [2019] V. Alba, B. Bertini, and M. Fagotti, SciPost Phys. 7, 5 (2019).
  * Bertini  _et al._ [2022] B. Bertini, K. Klobas, V. Alba, G. Lagnese, and P. Calabrese, arXiv:2203.17264 (2022).
  * Widom [1976] H. Widom, Advances in Mathematics 21, 1 (1976).
  * Widom [1975] H. Widom, Proceedings of the American Mathematical Society 50, 167 (1975).
  * Ares  _et al._ [2015] F. Ares, J. G. Esteve, F. Falceto, and A. R. de Queiroz, Phys. Rev. A 92, 042334 (2015).
  * Ares  _et al._ [2018] F. Ares, J. G. Esteve, F. Falceto, and A. R. de Queiroz, Phys. Rev. A 97, 062301 (2018).
  * Note [12] In case the reader is not comfortable with replacing a density matrix with a traceless operator, we point out that the result does not change if a term proportional to the identity is added to $\boldsymbol{\bf\rho}_{\rm eff}$.
  * Haag [1996] R. Haag,  _Local Quantum Physics: Fields, Particles, Algebras_ , 2nd ed. (Springer Publishing Company, Incorporated, 1996).
  * Kukuljan  _et al._ [2020] I. Kukuljan, S. Sotiriadis, and G. Takács, Journal of High Energy Physics 2020, 224 (2020).
  * Alaeian and Buča [2022] H. Alaeian and B. Buča, Exact bistability and time pseudo-crystallization of driven-dissipative fermionic lattices (2022).
  * Marić  _et al._ [2020a] V. Marić, S. M. Giampaolo, D. Kuić, and F. Franchini, New Journal of Physics 22, 083024 (2020a).
  * Marić  _et al._ [2022] V. Marić, S. M. Giampaolo, and F. Franchini, Phys. Rev. B 105, 064408 (2022).
  * Marić  _et al._ [2020b] V. Marić, S. M. Giampaolo, and F. Franchini, Communications Physics 3, 220 (2020b).
  * Mussardo [2020] G. Mussardo,  _Statistical Field Theory: An Introduction to Exactly Solved Models in Statistical Physics_, 2nd ed., Oxford Graduate Texts (Oxford University Press, Oxford, 2020).
  * Ault [1967] J. C. Ault, Journal of the London Mathematical Society s1-42, 497 (1967).
  * Marshall  _et al._ [2011] A. W. Marshall, I. Olkin, and B. C. Arnold, Matrix theory, in  _Inequalities: Theory of Majorization and Its Applications_ (Springer New York, New York, NY, 2011) pp. 297–365.
  * Peschel and Eisler [2009] I. Peschel and V. Eisler, Journal of Physics A: Mathematical and Theoretical 42, 504003 (2009).
  * Gaudin [1960] M. Gaudin, Nuclear Physics 15, 89 (1960).