#  Entanglement dynamics in confining spin chains

Stefano Scopa  SISSA and INFN, via Bonomea 265, 34136 Trieste, Italy  Pasquale
Calabrese  SISSA and INFN, via Bonomea 265, 34136 Trieste, Italy
International Centre for Theoretical Physics (ICTP), I-34151, Trieste, Italy
Alvise Bastianello  Department of Physics, Technical University of Munich,
85748 Garching, Germany  Institute for Advanced Study, 85748 Garching, Germany
Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4,
D-80799 München, Germany

###### Abstract

The confinement of elementary excitations induces distinctive features in the
non-equilibrium quench dynamics. One of the most remarkable is the suppression
of entanglement entropy which in several instances turns out to oscillate
rather than grow indefinitely. While the qualitative physical origin of this
behavior is clear, till now no quantitative understanding away from the field
theory limit was available. Here we investigate this problem in the weak
quench limit, when mesons are excited at rest, hindering entropy growth and
exhibiting persistent oscillations. We provide analytical predictions of the
entire entanglement dynamics based on a Gaussian approximation of the many-
body state, which captures numerical data with great accuracy and is further
simplified to a semiclassical quasiparticle picture in the regime of weak
confinement. Our methods are valid in general and we apply explicitly to two
prototypical models: the Ising chain in a tilted field and the experimentally
relevant long-range Ising model.

##  I Introduction

The relentless advances in cold atom experiments made available extremely
tunable platforms where synthetic phases of matter can be engineered with high
precision. Quantum simulators Feynman (2018); Georgescu  _et al._ (2014) are
on the verge of replacing classical devices in the challenge of understanding
and faithfully describing many-body and strongly interacting quantum systems,
promising a pathway to simulations that are far beyond the reach of classical
computers. A prominent example is the physics of strongly coupled gauge
theories routinely probed at large hadron colliders. Even though the promised
land of an accurate quantum simulator for such complex systems is getting
closer to a fast pace, it nowadays remains still low on the horizon. Recent
times witnessed an increasing interest in simpler, yet challenging, toy models
where the capability of quantum simulators can be put at test against the
state of the art analytical and numerical techniques. The physics of
confinement canonically belongs to high-energy physics in $3+1$ dimensions:
Quarks cannot exist in isolation and strong interaction confines them into
composite particles, such as mesons and baryons. However, confined excitations
can also be realised in simple one-dimensional condensed matter settings McCoy
and Wu (1978), which are within reach of current experimental venues Martinez
_et al._ (2016); Tan  _et al._ (2021). Admittedly, these realizations crudely
oversimplify the complicate processes taking place in hadron colliders, but
nevertheless capture many of their salient features. An ideal testbed for the
physics of confinement is the one dimensional Ising spin chain in a tilted
field Delfino  _et al._ (1996); Rutkevich (2008); Fonseca and Zamolodchikov
(2003); Kormos  _et al._ (2017). In the presence of a pure transverse field,
the Ising model is exactly solvable and equivalent to non interacting
fermions. Of particular interest is the $\mathbb{Z}_{2}-$spontaneously broken
phase, featuring two degenerate ground states of opposite magnetization. In
this regime, topological excitations have the form of domain walls (or kinks)
interpolating between the two vacua.
The nature of excitations dramatically changes in the presence of an
additional longitudinal field that explicitly breaks the aforementioned
symmetry by inducing a Zeeman-like gap in energy between the two lowest
levels, which split in a lower and higher energy states, dubbed true vacuum
and false vacuum respectively. This means that the spin chain pays energy
every time it visits the false vacuum and that domain walls experience a
potential that grows linearly in their separation. As a result, the kinks
(playing the analogous role of quarks) cannot be pulled infinitely far apart
without breaking energy conservation and get confined in composite
excitations, which are identified as mesons.

The appealing simplicity of the model is self-evident and, together with other
comparably simple spin chains, has been used to probe several aspects of
confinement physics. Already the pioneering paper Kormos  _et al._ (2017)
shows how confinement strongly suppresses the spreading of correlation after a
quantum quench. Shortly after, connections with quantum scarring Robinson  _et
al._ (2019); James  _et al._ (2019) and lattice gauge theories Kogut (1979);
Lerose  _et al._ (2020); Borla  _et al._ (2020); Chanda  _et al._ (2020); Pai
and Pretko (2020) have been unveiled. In general, confinement has been
understood to have striking consequence on transport Lerose  _et al._ (2020);
Mazza  _et al._ (2019); Bastianello  _et al._ (2021) and features anomalously
slow thermalization due to a strong suppression of meson creation/annihilation
mediated by the Schwinger mechanism Schwinger (1951); Sinha  _et al._ (2021);
Lagnese  _et al._ (2021a); Tortora  _et al._ (2020); Rigobello  _et al._
(2021); Verdel  _et al._ (2020); Myers  _et al._ (2017); Cubero and Robinson
(2019); Vovrosh and Knolle (2021). Notably, recent works closely mimic the
large hadron physics by studying scattering events of mesonic particles Surace
and Lerose (2021); Karpov  _et al._ (2020); Milsted  _et al._ (2021); Vovrosh
_et al._ (2021). The stepping stone for these investigations is the
possibility of efficiently simulate physics in one dimension, thanks to tensor
network techniques Schollwöck (2011). Numerical investigations backed up
qualitative or semi quantitative analytical results, contributing to clarify
the global picture. In contrast, the strongly interacting nature of these
models makes fully analytical treatments scarce, especially in nonequilibrium
settings.

This paper aims to analytically quantify the information spreading in confined
spin chains after a quantum quench by focusing on the dynamics of
entanglement. In the framework of confinement, analytical predictions for the
entanglement in the Ising chain close to the critical point have been obtained
in Ref. Castro-Alvaredo  _et al._ (2020) in the field theory limit, but the
far more common gapped scenario remains unexplored and motivates our
investigation. In the following, we first build on the standard quasiparticle
picture Calabrese and Cardy (2006, 2005); Calabrese (2020); Alba and Calabrese
(2017) by developing a semiclassical approach able to capture the behavior of
entanglement. Such semiclassical method is valid in the regime of large
quantum numbers, which is achieved with a small confining force. Away from
this regime, quantum effects lead to clear interference patterns in the
entanglement entropy, as already observed in Ref. Kormos  _et al._ (2017). In
order to capture this features, we develop a fully-quantum analytical
treatment that reproduces the numerical data with great accuracy.

Figure 1: Illustration of the quasiparticle picture in the confined model:
Small quenches mostly excite pairs of domain walls with zero total momentum,
which are then confined into mesons with zero velocity. In the dilute regime,
mesons behave as stationary extended particles with semiclassical trajectory
$d_{t}(k)$ in Eq. (5). The entanglement entropy measures the shared
information among the two subsystems $A\cup B$ mediated by the mesons sat at
the boundaries between the two regions.

##  II Model and quench protocol

As a prototypical model of confined spin chain, we shall focus on the Ising
chain in a tilted magnetic field with Hamiltonian

$\hat{H}=-J\sum_{j\in\mathbb{Z}}\left(\hat{\sigma}_{j}^{z}\hat{\sigma}_{j+1}^{z}+h_{\perp}\hat{\sigma}_{j}^{x}+h_{\parallel}\hat{\sigma}_{j}^{z}\right),$
(1)

where $h_{\parallel}$ (resp. $h_{\perp}$) are longitudinal (resp. transverse)
components of the magnetic field and $J$ is the overall energy scale that we
set hereafter to one. We denote with $\hat{\sigma}_{j}^{\alpha}$
($\alpha=x,y,z$) the standard Pauli operators acting on site $j$. For this
model, we consider the following quench protocol. At times $t\leq 0$, we
prepare the Ising chain in the ferromagnetic ground state by setting the
transverse field to $0\leq\bar{h}_{\perp}<1$ at zero longitudinal field
$h_{\parallel}=0$. For $t>0$, the system is then brought out of equilibrium by
the sudden variation of the transverse field $\bar{h}_{\perp}\to h_{\perp}$
and by turning on the longitudinal field $h_{\parallel}\neq 0$.
For $h_{\parallel}=0$, the model is equivalent to non interacting spinless
fermions and it allows for a well-known exact solution, see Appendix A for a
short summary.

In particular, the Hamiltonian (1) at $h_{\parallel}=0$ is diagonalized in
terms of fermionic modes
$\\{\hat{\gamma}(k),\hat{\gamma}^{\dagger}(q)\\}=\delta(k-q)$ obeying the
dispersion law $\epsilon(k)=2\sqrt{(h_{\perp}-\cos k)^{2}+\sin^{2}k}$. These
excitations are semiclassically interpreted as ballistic particles with
velocity $v(k)=\partial_{k}\epsilon(k)$. Deep in the ferromagnetic phase
$h_{\perp}\ll 1$, the model shows two degenerate ground states (${\rm
GS}_{\uparrow}$ and ${\rm GS}_{\downarrow}$, respectively) with opposite
magnetization $\langle\hat{\sigma}^{z}\rangle=\pm 1$ and the excitations can
be identified as sharp domain walls interpolating between the two ground
states. At finite transverse field $h_{\perp}<1$, the ground-state
magnetization is properly renormalized
$\langle\hat{\sigma}^{z}\rangle=\pm\bar{\sigma}$, with
$\bar{\sigma}=(1-h_{\perp}^{2})^{1/8}$ and the fermions are still identified
as dressed domain walls.
We further consider the confined model in Eq. (1) in the limit of small
quenches $|h_{\perp}-\bar{h}_{\perp}|,|h_{\parallel}|\ll 1$, where the
excitations of the gas are made of dilute mesons with zero total momentum, as
depicted in Fig. 1. As matter of convention, we choose to describe low-energy
excitations on the ground state with positive magnetization ($\ket{{\rm
GS}_{\uparrow}}\equiv\ket{0}$). Indeed, the ${\rm GS}_{\downarrow}$ sector has
the same description after replacing $h_{\parallel}\to-h_{\parallel}$. To this
end, we consider the matrix elements of $\hat{H}$ on the eigenbasis of the
pure transverse part, which can be labeled with asymptotic states
$|\\{k_{i}\\}_{i=1}^{n}\rangle$ and we focus on the form factors of the
longitudinal Pauli operator $\hat{\sigma}^{z}_{0}$, see Appendix A for their
explicit expression. These matrix elements can be divided in two classes:
those inducing interactions within a sector with a fixed number of fermions
and those violating fermion-number conservation.

A naive analysis of the Hamiltonian (1) would put the two terms on equal
footing, but effects spoiling the number conservation have been shown to be
exponentially suppressed in the weak longitudinal field Lerose  _et al._
(2020). Therefore, once the fermionic excitations are generated, they behave
as stable particles. Nevertheless, these modes experience a non trivial
interaction for $h_{\parallel}\neq 0$ due to the fermion-conserving part of
the perturbation. Such interactions are mainly of the form of a kink-antikink
pair that interact through a linear potential $2h_{\parallel}\bar{\sigma}d$,
where $d$ is their relative distance, plus negligible short range terms, see
Appendix A. Following Ref. Kormos  _et al._ (2017), one can perform a quantum
quench in the transverse field to create an initial distribution of fermionic
excitations Calabrese  _et al._ (2011) and later evolve the latter in the
presence of confinement. In this respect, the longitudinal field might appear
as a mere spectator during the preparation of the intial state as its effect
enters only at a second stage. However, its sudden activation gives non
trivial contributions to the early stage dynamics.

Let $|\psi_{0}\rangle$ be the initial state, which is let to evolve with
Hamiltonian (1) as $|\psi_{t}\rangle=e^{-\text{i}t\hat{H}}|\psi_{0}\rangle$.
We split $\hat{H}$ into two parts
$\hat{H}=\hat{H}_{\text{diag}}+\hat{H}_{\text{off-diag}}$, where
$\hat{H}_{\text{diag}}$ is the projection of $\hat{H}$ onto the number
conserving sector. Despite the fact that $\hat{H}_{\text{off-diag}}$ can be
neglected in the late time dynamics, it creates excitations within a finite
time scale after the quench. For this reason, it is useful to write
$|\psi_{t}\rangle=e^{-\text{i}t\hat{H}_{\text{diag}}}{\rm
T}\exp\left[-\text{i}\int_{0}^{t}{\rm d}\tau\,\hat{H}^{(I)}_{\text{off-
diag}}(\tau)\right]|\psi_{0}\rangle$ where $\hat{H}^{(I)}_{\text{off-diag}}$
is the off-diagonal term of the Hamiltonian in the interaction picture with
respect to the diagonal part. Under the assumption of weak longitudinal field
$h_{\parallel}\ll 1$, one can argue on a separation of time scales and
approximate the state as $|\psi_{t}\rangle\simeq
e^{-\text{i}t\hat{H}_{\text{diag}}}|\tilde{\psi}_{0}\rangle$ by defining the
renormalized initial state

$|\tilde{\psi}_{0}\rangle\simeq{\rm
T}\exp\left[-\text{i}\int_{0}^{t_{h_{\perp}}}{\rm
d}\tau\,\hat{H}^{(I)}_{\text{off-diag}}(\tau)\right]|\psi_{0}\rangle\,,$ (2)

with $t_{h_{\perp}}$ a large time scale compared with the transverse Ising
energy scale $t_{h_{\perp}}\epsilon(k)\gg 1$, but still smaller when compared
with the longitudinal field $t_{h_{\perp}}h_{\parallel}\ll 1$. In Appendix B
we substantiate on this idea and compute the renormalized initial state
$\ket{\tilde{\psi}_{0}}$ in the sought limit of small quenches using form
factor perturbation theory Delfino  _et al._ (1996); Delfino (2014); Delfino
and Viti (2017) and obtaining a simple squeezed state

$|\tilde{\psi}_{0}\rangle\propto\exp\left[\int_{0}^{\pi}{\rm
d}k\,\mathcal{K}(k)\hat{\gamma}^{\dagger}(k)\hat{\gamma}^{\dagger}(-k)\right]|0\rangle\,.$
(3)

Above, mode operators act on the post-quench vacuum
$\hat{\gamma}(k)|0\rangle=0$. The wave function $\mathcal{K}$ can be written
as

$\mathcal{K}(k)=-\text{i}\tan(\theta_{k}^{h_{\perp}}-\theta_{k}^{\bar{h}_{\perp}})-\text{i}h_{\parallel}\bar{\sigma}v(k)/\epsilon^{2}(k)\,,$
(4)

where $\theta_{k}^{h_{\perp}}$ is the Bogoliubov angle that diagonalizes the
transverse part of the Hamiltonian, see Appendix B for the detailed
derivation. Here, the energy $\epsilon(k)$ and velocity $v(k)$ are computed
with respect to the post-quench transverse field $h_{\perp}$. The squeezed
form represents an incoherent superposition of pairs of particle with opposite
momenta and density $n_{k}=|\mathcal{K}(k)|^{2}(1+|\mathcal{K}(k)|^{2})^{-1}$.
The limit of small quenches is then attained for $|\mathcal{K}(k)|^{2}\ll 1$.

##  III Entanglement dynamics

Entanglement quantifies the amount of shared quantum information in a
bipartition $A\cup B$. The most common entanglement measure is the von Neumann
entropy $S_{A}=-\text{Tr}\hat{\rho}_{A}\log\hat{\rho}$ of the reduced density
matrix $\hat{\rho}_{A}=\text{Tr}_{B}|\psi_{t}\rangle\langle\psi_{t}|$, but the
Rényi entropies $S^{(n)}_{A}=-\frac{1}{1-n}\log(\text{Tr}\hat{\rho}_{A}^{n})$
recently became experimentally accessible Islam  _et al._ (2015); M.  _et al._
(2016); Elben  _et al._ (2018); Brydges  _et al._ (2019); Elben  _et al._
(2020); Neven  _et al._ (2021); Vitale  _et al._ (2021).

The quasiparticle picture (QP) Calabrese and Cardy (2006, 2005); Alba and
Calabrese (2017) judiciously supplements few quantum inputs with semiclassical
arguments, providing the leading order extensive part of the entanglement
entropy. Within this framework, local pairs of entangled excitations with
opposite momenta are seen as an initial source for the entanglement and,
during the dynamics, the propagating pairs that are shared by the two
subsystems become responsible for the entanglement growth. Notably, the QP
picture can be used beyond the pair structure of the post-quench state Bertini
_et al._ (2018); Bastianello and Calabrese (2018); Bastianello and Collura
(2020). In the context of confinement, these ideas give a quick qualitative
grasp on the expected behavior of entanglement. In particular, sooner or
later, the confining force bends the short time ballistic trajectories
$d_{t}(k)\simeq 2v(k)t$ and consequently mitigates the entanglement growth,
see Fig. 1 for an illustration. Focusing on a single pair created at momenta
$(-k,k)$ and overlapping position, the classical distance is Rutkevich (2008)

$d_{t}(k)=(h_{\parallel}\bar{\sigma})^{-1}[\epsilon(k)-\epsilon(k-2h_{\parallel}\bar{\sigma}(t\text{
mod }k/(h_{\parallel}\bar{\sigma})))]\,.$ (5)

In the limit of dilute mesons, each fermionic pair remains well-separated and
thus independent from others during the time evolution. Precisely, this regime
is reached if $d_{\text{max}}\int_{0}^{\pi}\frac{{\rm d}k}{2\pi}\,n_{k}\ll 1$,
where $d_{\text{max}}=4h_{\perp}/(|h_{\parallel}|\bar{\sigma})$ is the maximum
extension of a meson. In the case where only the longitudinal field is
quenched and $h_{\perp}$ remains constant, Eq. (4) leads to the simple
expression $\int_{0}^{\pi}\frac{{\rm
d}k}{2\pi}\,n_{k}=(h_{\parallel}h_{\perp}\bar{\sigma})^{2}/[16(1-h_{\perp}^{2})^{3}]$.
Hence, the bond is made explicit and it is linearly improved in the weak
longitudinal field. Conversely, at higher densities, the single-meson physics
breaks down due to collisions among mesons. We shall focus on the dilute
regime, where the standard quasiparticle prediction for the entanglement
entropy can be readily generalized to the confined model by replacing the term
$2tv(k)$ with $d_{t}(k)$ given in Eq. (5). In the case of an infinite system
divided in two semi-infinite halves $A\cup B=(-\infty,0]\cup(0,+\infty)$, the
von Neumann entropy is given by the formula

$S_{A}(t)=\int_{0}^{\pi}\frac{{\rm d}k}{2\pi}d_{t}(k)[-n_{k}\log
n_{k}-(1-n_{k})\log(1-n_{k})]\,.$ (6)

Figure 2: von Neumann entropy dynamics of the confined Ising model (1). We
show different quenches of the transverse field on different columns and
different values of $h_{\parallel}$ on different rows. In each panel, we
compare the semiclassical prediction (6) (dashed line), the quantum prediction
in Gaussian approximation (full line) and numerical data obtained with time-
dependent DMRG simulations (symbols). Notice the validity of the approach also
for quenches starting from the false vacuum (g)-(h).

In Fig. 2 we test the semiclassical prediction in Eq. (6) against tensor
network simulations for different quench parameters, finding a very good
agreement in the regime of weak longitudinal field $h_{\parallel}\ll 1$.
Notice that the applicability of our analysis holds also for negative
longitudinal fields, namely when mesons are true-vacuum excitations in a false
vacuum sea, see Fig. 2 (g)-(h). In this case, the linear potential changes
sign and becomes repulsive, but the lattice induces Bloch oscillations
preventing the domain walls to be repelled infinitely far apart (see also Ref.
Pomponio  _et al._ (2021)). By increasing the value of the longitudinal field,
deviations from Eq. (6) eventually undermine the approximation of isolated
mesons, see Appendix E for further details. However, we observe a failure of
the semiclassical approximation far before this limiting case. Indeed, in Eq.
(5) $d_{t}(k)$ is treated as a continuum variable and, a consequence, Eq. (6)
is expected to hold only in the limit $d_{\text{max}}\gg 1$. Needless to say,
a simple theory outclassing this semiclassical treatment and including quantum
effects is highly desirable.

###  III.1 Quantum effects beyond the quasiparticle picture

Computing the entanglement of many-body systems is an extremely hard task with
very few notable exceptions where analytical derivations are possible,
including e.g. critical models Calabrese and Cardy (2004, 2009) and non
interacting systems Casini and Huerta (2009). While in Ref. Castro-Alvaredo
_et al._ (2020), conformal invariance has been exploited for similar purposes,
we rather borrow methods usually applied to non interacting models. Let us
consider a system of spinless fermionic operators $\hat{c}_{j}$,
$\hat{c}^{\dagger}_{j}$, which will eventually be chosen as the Jordan-Wigner
fermions of the Ising chain (1), see Appendix C. The distinctive feature of a
non interacting system is the applicability of Wick’s theorem, which holds for
any Gaussian state. Gaussian states are fully characterized by their quadratic
correlations

$C=\begin{pmatrix}\langle\hat{c}_{j}\hat{c}^{\dagger}_{j^{\prime}}\rangle&&\langle\hat{c}_{j}\hat{c}_{j^{\prime}}\rangle\\\
\langle\hat{c}^{\dagger}_{j}\hat{c}^{\dagger}_{j^{\prime}}\rangle&&\langle\hat{c}^{\dagger}_{j}\hat{c}_{j^{\prime}}\rangle\end{pmatrix}$
(7)

From the knowledge of $C$, one can easily construct the correlation matrix
$C_{A}$ restricted to the subsystem $A$ and evaluate the entanglement entropy
as $S_{A}=-\text{Tr}[C_{A}\log C_{A}]$ Casini and Huerta (2009). Notice that
the trace is over the lattice sites of $A$, thus the dimension of the vector
space scales linearly with the size of $A$, rather than exponentially as the
dimension of the Hilbert space describing the subsystem. Therefore, the
entanglement entropy is obtained in a matter of seconds on a laptop.
Notably, the initial state in Eq. (3) is Gaussian. We now argue that the
dilute regime is compatible with a Gaussian-preserving dynamics. Indeed, since
Eq. (3) describes an incoherent superposition of paired fermions that evolve
independently, it is expected to capture dilute mesons provided the single-
meson dynamics is properly described by a non trivial time-dependent wave
function $\mathcal{K}_{t}(k)$. The wave function $\mathcal{K}_{t}(k)$ is
obtained from the Schrödinger equation
$\text{i}\partial_{t}|\psi_{t}\rangle=\hat{H}_{\text{off-
diag}}|\psi_{t}\rangle$ by projecting onto the two-particle subspace. The
details of this derivation are lengthy but straightforward and can be found in
Appendix A. The resulting equation is conveniently written in real space as

$\text{i}\partial_{t}W_{t}(j)=2h_{\parallel}\bar{\sigma}|j|W_{t}(j)+\sum_{\ell\in\mathbb{Z}}[T_{j-\ell}-T_{-j-\ell}]W_{t}(\ell)\,,$
(8)

where $W_{t}(j)$ is an antisymmetric wave function related to ${\cal
K}_{t}(k)$ as $W_{t}(j)=\int\frac{{\rm
d}k}{2\pi}e^{-\text{i}kj}\mathcal{K}_{t}(k)$ and $T_{j}=\int\frac{{\rm
d}k}{2\pi}e^{-\text{i}kj}2\epsilon(k)$ is the kinetic energy term. Notice that
in the regime of weak transverse field $\epsilon(k)\simeq 2-2h_{\perp}\cos k$,
$T_{j}$ reduces to a first neighbor hopping and the Schrödinger equation (8)
can be exactly diagonalized Rutkevich (2008). For finite transverse fields,
the one-body problem can be easily tackled numerically. For completeness, it
should be mentioned that additional small short-range interactions are present
in Eq. (8), see Appendix A for details. In Fig. 2 we compare tensor network
simulations against the quantum prediction obtained in Gaussian approximation,
finding an excellent agreement. Technical details on the numerical
calculations and the comparison with tensor network simulations are reported
in Appendix C.

Figure 3: (Left) Rényi-2 and (right) von Neumann entropy dynamics for the
long-range Ising model (9) with different exponent $\alpha$ on different rows.
The system is initially prepared in a ferromagnetic state with all spins up.
The quantum prediction in Gaussian approximation (full line) is compared with
numerical data (dashed line) obtained with global subspace expansion time-
dependant variational principle Yang and White (2020).

##  IV Confined spin chains in experiments

Recent experimental platforms of trapped ions show clear signature of confined
physics Tan  _et al._ (2021). In this setups, the confinement is induced by
the long-range interactions of the Ising Hamiltonian Liu  _et al._ (2019)

$\hat{H}=-\sum_{i<j}\frac{1}{|i-j|^{\alpha}}\hat{\sigma}^{z}_{i}\hat{\sigma}_{j}^{z}-h_{\perp}\sum_{i}\hat{\sigma}_{i}^{x}\,.$
(9)

For small transverse field $h_{\perp}$, the fundamental excitations of the
model (9) are domain walls, analogously to the short-range Ising chain in Eq.
(1) we previously discussed. In view of Eq. (9), a pair of domain walls
experiences an attractive potential that for $1<\alpha\leq 2$ indefinitely
grows as the distance between the particles increases. In this way, a
confining mechanisms is induced. For $\alpha>2$, the intra-kink interaction is
bounded, but it still supports deep bound states with meson-like dynamics for
$\alpha$ not too large. The methods discussed previously for the Hamiltonian
(1) can be easily generalized to the long-range Ising model (9). In Fig. 3, we
initialize the system in a product state of all spins up and we let it evolve
with the Hamiltonian (9) for small $h_{\perp}$. Due to its experimental
relevance Islam  _et al._ (2015); M.  _et al._ (2016); Elben  _et al._ (2018);
Brydges  _et al._ (2019); Elben  _et al._ (2020); Neven  _et al._ (2021);
Vitale  _et al._ (2021), we also show the results for the Rényi $S^{(2)}_{A}$
which, analogously to the Von Neumann entropy, is accessible from correlation
functions (7) (see Appendix D). In contrast with the longitudinal case, the
long-range Hamiltonian (9) induces meson-meson and boundary-meson interactions
beyond short range terms, harming the approximation of dilute, thus non
interacting, mesons. As a consequence of interactions, mesons become mobile
causing an additional linear drift of the entanglement over the single-meson
physics result. As expected, the drift is more evident as $\alpha$ is smaller,
see Fig. 3. Further details on the state preparation and on the breakdown our
approximations are left to Appendices D and E.

##  V Conclusions

We provide an analytical study of the entanglement dynamics in confined spin
chains based on a minimal set of ingredients. Precisely: (i) mesons are
incoherently excited and non interacting, hence squeezed states remain a good
approximation also at late times and (ii) Gaussian correlators are evolved
according to the dynamics projected in the two-fermion sector. Although our
results have been explicitly worked out only for Ising spin chains with short
and long range interactions, it is clear that our methods apply in general, as
e.g. in confining spin ladders Lagnese  _et al._ (2021b, 2020), Potts models
Lencsés and Takács (2015); Delfino and Grinza (2008) and tricritical Ising
Lencsés  _et al._ (2021). Natural future directions address the breakdown of
the single-meson approximation and the effect of moving _—_ and hence
scattering _—_ mesons on the entanglement dynamics. Particularly appealing is
the semiclassical limit, where entangled particles move along classical
trajectories in the spirit of Bertini  _et al._ (2018); Bastianello and
Calabrese (2018); Bastianello and Collura (2020); Moca  _et al._ (2017).
Quantifying how the mesons scatter would eventually allow us to understand
whether these confining models thermalize.

## Acknowledgements

AB acknowledges support from the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy – EXC–2111–390814868.
PC and SS acknowledge support from ERC under Consolidator Grant No. 771536
(NEMO). The authors acknowledge Mario Collura for useful discussions and
insights on the numerical simulations. SS is grateful to Lorenzo Piroli for
discussions and to Emanuele Tirrito for help with TDVP simulations. AB is
grateful to Alessio Lerose, Michael Knap, Stefan Birnkammer, Johannes Knolle,
Joseph Vovrosh and Hongzheng Zhao for several discussions and their insight on
related work. Numerical simulations have been performed with the open-source
library iTensor Fishman  _et al._ (2020).

##  Appendix A The Ising chain in a weak longitudinal field: the form factor
approach

We consider the pure transverse field Ising chain

$\hat{H}_{\text{trans}}=-\sum_{j\in\mathbb{Z}}\left(\hat{\sigma}^{z}_{j+1}\hat{\sigma}_{j}+h_{\perp}\hat{\sigma}^{x}_{j}\right)\,,$
(10)

which can be exactly diagonalized in terms of non interacting spinless
fermions as we now briefly review. For simplicity, we shall describe the Ising
model on a lattice of $N$ sites, assuming periodic boundary conditions of the
chain and eventually we take the thermodynamic limit where $N\to\infty$.
First, we introduce fermionic creation and annihilation operators
$\hat{c}^{\dagger}_{j}$, $\hat{c}_{j}$ by means of a Jordan-Wigner
transformation of the ladder Pauli operators as

$\frac{\hat{\sigma}_{j}^{x}+\text{i}\hat{\sigma}_{j}^{y}}{2}=\exp\big{(}\text{i}\pi\sum_{i<j}\hat{c}^{\dagger}_{i}\hat{c}_{i}\big{)}\hat{c}^{\dagger}_{j}\,.$
(11)

With this transformation, it is possible to divide the Ising Hamiltonian (10)
in two terms $\hat{H}_{\text{trans}}^{\pm}$ according with the parity operator
$\hat{P}=\prod_{j}\hat{\sigma}_{j}^{x}$ (notice that $\hat{P}$ has eigenvalues
$P=\pm 1$ and it commutes with the transverse Hamiltonian)

$\hat{H}_{\text{trans}}=\frac{1}{2}(1+\hat{P})\hat{H}_{\text{trans}}^{+}+\frac{1}{2}(1-\hat{P})\hat{H}_{\text{trans}}^{-}\,.$
(12)

In this form, the Ising Hamiltonian becomes quadratic in the fermionic
operators and translationally invariant in each of the two parity sectors.
Notice that the terms $\hat{H}_{\text{trans}}^{\pm}$ differ only for the
boundary terms: The parity sector $P=+1$ ($P=-1$) is characterized by
antiperiodic (periodic) boundary conditions for the fermionic operators and it
induces a quantization of momenta in half-integer (integer) multiples of
$2\pi/N$. Nevertheless, these boundary effects become negligible in the
thermodynamic limit and will be thus omitted in what follows. Hence, in both
parity sectors the Hamiltonian can be diagonalized in the momentum space as

$\hat{H}_{\text{trans}}^{\pm}=\int{\rm d}k\
\epsilon(k)\hat{\gamma}_{\pm}^{\dagger}(k)\hat{\gamma}_{\pm}(k)+\
\text{const.}$ (13)

with single-particle energy

$\epsilon(k)=2\sqrt{(h_{\perp}-\cos(k))^{2}+\sin^{2}k}.$ (14)

The new set of operators $\hat{\gamma}_{\pm}$ is defined by the Bogoliubov
rotation

$\begin{pmatrix}\hat{c}_{j}\\\\[4.0pt]
\hat{c}^{\dagger}_{j}\end{pmatrix}=\int\frac{{\rm
d}k}{\sqrt{2\pi}}\;U_{\theta_{k}}\begin{pmatrix}\hat{\gamma}_{\pm}(k)\\\\[4.0pt]
\hat{\gamma}_{\pm}^{\dagger}(-k)\end{pmatrix}$ (15)

where

$U_{\theta_{k}}=\begin{pmatrix}\cos\theta_{k}&&\text{i}\sin\theta_{k}\\\
\text{i}\sin\theta_{k}&&\cos\theta_{k}\end{pmatrix}$ (16)

and $\theta_{k}$ is the Bogoliubov angle

$\theta_{k}=-\frac{1}{2\text{i}}\log\left(\frac{h_{\perp}-e^{\text{i}k}}{(\cos
k-h_{\perp})^{2}+\sin^{2}k}\right)\,.$ (17)

With our conventions for the Fourier transform, the new modes obey standard
anticommutation relations
$\\{\hat{\gamma}_{\pm}(k),\hat{\gamma}^{\dagger}_{\pm}(q)\\}=\delta(k-q)$ for
any value of $\theta_{k}$. In each of the two parity sectors, the Hilbert
space is understood as a Fock space built on the vacua $|0_{\pm}\rangle$ that
are annihilated by the respective mode operators
$\hat{\gamma}_{\pm}(k)|0_{\pm}\rangle=0$. For $h_{\perp}<1$, the model
undergoes spontaneous symmetry breaking in the thermodynamic limit. This can
be easily understood in the classical limit $h_{\perp}\to 0$, where there are
two degenerate ground states $|{\rm
GS}_{\uparrow}\rangle=|...\uparrow\uparrow\uparrow...\rangle$ and $|{\rm
GS}_{\downarrow}\rangle=|...\downarrow\downarrow\downarrow...\rangle$. Through
linear combinations of these two vacua, one forms $P$-symmetric and
antisymmetric pairs, which are then identified with the vacua
$|0_{\pm}\rangle=(|{\rm GS}_{\uparrow}\rangle\pm|{\rm
GS}_{\downarrow}\rangle)/\sqrt{2}$ with opposite value of magnetization
$\bar{\sigma}=\pm 1$. At finite $h_{\perp}$, the degenerate ground states are
not simple product states any longer and the longitudinal magnetization gets
renormalized to the value $\bar{\sigma}\equiv\langle{\rm
GS}_{\uparrow}|\sum_{j=1}^{N}\hat{\sigma}^{z}_{j}/N|{\rm
GS}_{\uparrow}\rangle=-\langle{\rm
GS}_{\downarrow}|\sum_{j=1}^{N}\hat{\sigma}^{z}_{j}/N|{\rm
GS}_{\downarrow}\rangle=(1-h_{\perp}^{2})^{1/8}$ Rutkevich (2008).

###  A.1 The form factor approach to the longitudinal field

We now wish to consider the effect of the longitudinal field $h_{\parallel}$
in Eq. (1) on the low-energy excitations above one of the degenerate vacua.
For the sake of concreteness, we choose the ground state with positive
magnetization $|{\rm GS}_{\uparrow}\rangle$. As a first step, one builds the
proper excitation basis on the Fock space of the transverse part, by
considering the symmetric combination of the Fock spaces with opposite parity.
More precisely, let $|(\\{k_{i}\\}_{i}^{N})_{\pm}\rangle$ be a multiparticle
state in the $P=\pm 1$ sectors, then one defines
$|\\{k_{i}\\}_{i=1}^{N}\rangle=\frac{1}{\sqrt{2}}(|(\\{k_{i}\\}_{i}^{N})_{+}\rangle+|(\\{k_{i}\\}_{i}^{N})_{-}\rangle)$.
In this notation, the state with no excitations is the desired ground state,
i.e., $|0\rangle\equiv|{\rm GS}_{\uparrow}\rangle$. Notice that the analysis
at finite system size $N$ is more involved due to the different quantization
conditions of the even and odd parity sectors. The action of the confined
Ising Hamiltonian in Eq. (1) on this basis is

$\hat{H}|\\{k_{i}\\}_{i=1}^{N}\rangle=\left(\sum_{i=1}^{N}\epsilon(k_{i})\right)|\\{k_{i}\\}_{i=1}^{N}\rangle+h_{\parallel}\sum_{M=1}^{\infty}\frac{1}{M!}\int_{-\pi}^{\pi}\frac{{\rm
d}^{M}q}{(2\pi)^{M}}\,2\pi\delta\left(\sum_{j=1}^{M}q_{j}-\sum_{i=1}^{N}k_{i}\right)|\\{q_{j}\\}_{j=1}^{M}\rangle\langle\\{q_{j}\\}_{j=1}^{M}|\hat{\sigma}^{z}_{0}|\\{k_{i}\\}_{i=1}^{N}\rangle\,.$
(18)

In this expression, the form factors
$\langle\\{q_{j}\\}_{j=1}^{M}|\hat{\sigma}^{z}_{0}|\\{k_{i}\\}_{i=1}^{N}\rangle$
that appear were previously determined in Ref. Rutkevich (2008). By a repeated
use of Wick theorem, one can be express the latter in terms of the following
two particle form factors:

$\langle
k_{1}k_{2}|\hat{\sigma}^{z}_{0}|0\rangle=\bar{\sigma}F(k_{1},k_{2}|);$ (19a)
$\langle
0|\hat{\sigma}^{z}_{0}|k_{1},k_{2}\rangle=\bar{\sigma}F(|k_{1},k_{2});$ (19b)
$\langle k_{1}|\sigma_{0}^{z}|k_{2}\rangle=\bar{\sigma}F(k_{1}|k_{2}),$ (19c)

where

$\begin{split}&F(k_{1},k_{2}|)=-F^{*}(|k_{1},k_{2})\\\\[5.0pt]
&\quad=\frac{1}{1-\exp[\text{i}(k_{1}+k_{2})]}\frac{\epsilon(k_{1})-\epsilon(k_{2})}{\sqrt{\epsilon(k_{1})\epsilon(k_{2})}}\,\end{split}$
(20)

and

$F(k_{1}|k_{2})=\frac{1}{1-\exp[\text{i}(k_{1}-k_{2})]}\frac{\epsilon(k_{1})+\epsilon(k_{2})}{\sqrt{\epsilon(k_{1})\epsilon(k_{2})}}\,.$
(21)

Notice that $F(|k_{1},k_{2})$ is always finite while $F(k_{1}|k_{2})$ develops
a kinematic singularity for $k_{1}=k_{2}$ that has to be properly regularized
by a symmetric shift in the imaginary axis

$\begin{split}&\frac{2}{1-\exp[\text{i}(k_{1}-k_{2})]}\to\left[\frac{2}{1-\exp[\text{i}(k_{1}-k_{2})]}\right]_{\text{reg}}\\\\[5.0pt]
&=\frac{1}{1-\exp[\text{i}(k_{1}-k_{2})]+\text{i}0^{+}}+\frac{1}{1-\exp[\text{i}(k_{1}-k_{2})]-\text{i}0^{+}}\,.\end{split}$
(22)

###  A.2 The dynamics in the two particle sector

The singular part of $F(k_{1}|k_{2})$ is responsible for the long-range
confinement, as it is clarified by the dynamics in the two-fermion sector that
we now analyze. Hence, we now focus on the matrix elements of the form
$\langle q_{1},q_{2}|\hat{H}|k_{1},k_{2}\rangle$. In particular, the matrix
element of the longitudinal Pauli operators reads

$\begin{split}&\bar{\sigma}^{-1}\langle
q_{1}q_{2}|\hat{\sigma}_{0}^{z}|k_{1},k_{2}\rangle=F(q_{1},q_{2}|)F(|k_{1},k_{2})\\\\[3.0pt]
&\qquad+F(q_{1}|k_{2})F(q_{2}|k_{1})-F(q_{1}|k_{1})F(q_{2}|k_{2})\,.\end{split}$
(23)

We further specialize to the zero-momentum sector where one sets
$(k_{1},k_{2})=(k,-k)$ and similarly $(q_{1},q_{2})=(q,-q)$, due to momentum
conservation. In this case, we obtain $F(|k,-k)=\text{i}v(k)/\epsilon(k)$,
while the singularity in Eq. (22) requires a more careful treatment. The
calculation are lengthy but straightforward, hence we report below the final
result only (see Ref. Rutkevich (2008) for further details). We define the
two-body interaction in the momentum space as
$V(q,k)\equiv\bar{\sigma}^{-1}\langle q,-q|\hat{\sigma}_{0}^{z}|k,-k\rangle$
and we divide the regular part from the singular part as

$\begin{split}&V(q,k)=V^{\text{reg}}(q,k)+\\\\[4.0pt]
&\quad\left[\frac{2}{1-\exp[\text{i}(q+k)]}\right]_{\text{reg}}\left[\frac{2}{1-\exp[\text{i}(-q-k)]}\right]_{\text{reg}}\\\\[4.0pt]
&\quad-\left[\frac{2}{1-\exp[\text{i}(q-k)]}\right]_{\text{reg}}\left[\frac{2}{1-\exp[\text{i}(-q+k)]}\right]_{\text{reg}}\,.\end{split}$
(24)

The double pole singularity is responsible for the linear confinement, while
$V^{\text{reg}}(q,k)$ is not singular and gives rise to a short-range
interaction. Next, we introduce the wave function ${\cal K}(k)$ and we write
the two-body state as

$|\psi_{2-\text{fermions}}\rangle=\frac{1}{2}\int_{-\pi}^{\pi}{\rm
d}k\,\mathcal{K}(k)|k,-k\rangle\,,$ (25)

where the notation is chosen to match the expansion of the squeezed state in
Eq. (3). Furthermore, due to the long-range nature of confinement, it is
convenient to move to coordinate space by defining the real-space wave
function $W(j)$ as

$W(j)=\int\frac{{\rm d}k}{2\pi}e^{-\text{i}kj}\mathcal{K}(k)\,.$ (26)

Hence, projecting the Schrödinger equation to the two-particles state and by
evaluating the action of the Hamiltonian in Eq. (1) on the wave function
$W(j)$, one obtains

$\text{i}\partial_{t}W(j)=2h_{\parallel}\bar{\sigma}|j|W(j)+\sum_{\ell>0}\Big{[}T_{j-\ell}-T_{-j-\ell}-T_{j+\ell}+T_{-j+\ell}+\frac{h_{\parallel}\bar{\sigma}}{2}V^{\text{reg}}_{j,\ell}\Big{]}W(\ell)\,,$
(27)

with $T_{j}=\int_{-\pi}^{\pi}\frac{{\rm d}k}{2\pi}e^{-\text{i}kj}2\epsilon(k)$
and

$V^{\text{reg}}_{j,\ell}=\int\frac{{\rm d}k{\rm
d}q}{(2\pi)^{2}}e^{-\text{i}kj}e^{\text{i}q\ell}V^{\text{reg}}(k,q)\,.$ (28)

Notice that $V^{\text{reg}}_{j,\ell}$ is a short range interaction that
vanishes whenever $j$ and $\ell$ are large, regardless their relative
difference. While in principle both the linear potential and the short term
interactions are of order ${\cal O}(h_{\parallel})$, the long-ranged nature of
the first makes $V^{\text{reg}}_{j,\ell}$ negligible in most of practical
cases.

##  Appendix B State preparation thorough a quantum quench

In this section, we characterize the state obtained after a quantum quench in
the transverse and the longitudinal magnetic field. Let $|\psi_{0}\rangle$ be
the initial state, set as the ground state of the model for a certain value of
transverse field $\bar{h}_{\perp}$ at zero longitudinal coupling. To begin
with, we express $\ket{\psi_{0}}$ in the basis of the post-quech fermions,
which amounts to solve the quench in the transverse field Calabrese  _et al._
(2011). From Eq. (15), one can relate the pre- and post-quench modes
$\hat{\bar{\gamma}}(k)$ and $\hat{\gamma}(k)$ respectively by a simple
rotation
$\hat{\bar{\gamma}}(k)=\cos(\theta_{k}^{h_{\perp}}-\theta_{k}^{\bar{h}_{\perp}})\hat{\gamma}(k)+\text{i}\sin(\theta_{k}^{h_{\perp}}-\theta_{k}^{\bar{h}_{\perp}})\hat{\gamma}^{\dagger}(-k)$.
Imposing that $\hat{\bar{\gamma}}(k)|\psi_{0}\rangle=0$, one obtains a simple
equation for the post-quench modes, which has the following solution

$|\psi_{0}\rangle=\frac{1}{\sqrt{\mathcal{N}}}\exp\left[\int_{0}^{\pi}{\rm
d}k\,K(k)\hat{\gamma}^{\dagger}(k)\hat{\gamma}^{\dagger}(-k)\right]|0\rangle\,,$
(29)

where
$K(k)=-\text{i}\tan(\theta_{k}^{h_{\perp}}-\theta_{k}^{\bar{h}_{\perp}})$ and
$\mathcal{N}$ is a normalization factor. However, as anticipated in Sec. II,
the activation of a longitudinal magnetic field has non trivial consequences
on the state that we now address by computing the renormalized initial state

$|\tilde{\psi}_{0}\rangle\simeq{\rm
T}\exp\left[-\text{i}\int_{0}^{t_{h_{\perp}}}{\rm
d}\tau\,\hat{H}^{(I)}_{\text{off-diag}}(\tau)\right]|\psi_{0}\rangle.$ (30)

In the limit of small quench where $K(k)$ is small, the initial state (29) can
be expanded as

$|\psi_{0}\rangle\propto|0\rangle+\int_{0}^{\pi}{\rm
d}k\,K(k)|k,-k\rangle+\dots$ (31)

and, at leading order in $h_{\parallel}\ll 1$, one obtains from (30)

$\begin{split}|\tilde{\psi}_{0}\rangle&=\frac{1}{\sqrt{\mathcal{N}}}\Big{(}|0\rangle+\int_{0}^{\pi}{\rm
d}k\,K(k)|k,-k\rangle\\\ &\qquad-\text{i}\int_{0}^{t_{h_{\perp}}}{\rm d}\tau\
\hat{H}^{(I)}_{\text{off-diag}}(\tau)|0\rangle+\dots\Big{)}\,.\end{split}$
(32)

In principle, $\hat{H}^{(I)}_{\text{off-diag}}(\tau)|0\rangle$ can couple to a
sector with an arbitrary even number of fermions, however in the small quench
limit $h_{\perp}\ll 1$ the most important contribution comes from the two-
particle sector. This can be argued on the basis of two facts: (i) as we will
see, the cross ratio of a particles pair creation is associated with the form
factors $F(k_{1},k_{2}|)$ in Eq. (20), which for small transverse field scale
as $\sim h_{\perp}$, hence couplings to the sector with $2n$ fermions are of
order $\sim h_{\perp}^{n}$; (ii) any process is multiplied by the inverse of
the energy of the final state, which means that couplings to sectors with more
fermions are further suppressed. Motivated by this reasoning, we now wish to
compute the quantity $\int_{0}^{t_{h_{\perp}}}{\rm d}\tau\langle
k_{1},k_{2}|\hat{H}^{(I)}_{\text{off-diag}}(\tau)|0\rangle$. As a further
approximation, it is expected that for $h_{\parallel}\ll 1$ fermions are
locally created on a much shorter time scale compared to the one over which
the confining force shows appreciable effects. Therefore, we write the
Hamiltonian in the interaction picture as $\hat{H}^{(I)}_{\text{off-
diag}}(\tau)\simeq e^{\text{i}\tau H_{\text{trans}}}\hat{H}_{\text{off-
diag}}e^{-\text{i}\tau\hat{H}_{\text{trans}}}$, neglecting the term of
confinement. Under this program, we finally obtain the expression

$\begin{split}&-\text{i}\int_{0}^{t_{h_{\perp}}}{\rm d}\tau\langle
k_{1},k_{2}|\hat{H}^{(I)}_{\text{off-diag}}(\tau)|0\rangle\simeq\\\
&-\text{i}\int_{0}^{t_{h_{\perp}}}{\rm d}\tau
h_{\parallel}\bar{\sigma}e^{\text{i}\tau(\epsilon(k_{1})+\epsilon(k_{2}))}\left[\sum_{j}e^{\text{i}(k_{1}+k_{2})j}\right]F(k_{1},k_{2}|)\\\\[4.0pt]
&\stackrel{{\scriptstyle\epsilon(k)t_{h_{\perp}}\gg
1}}{{=}}-\text{i}2\pi\delta(k_{1}+k_{2})h_{\parallel}\bar{\sigma}\frac{v(k)}{\epsilon^{2}(k)}.\end{split}$
(33)

In the last passage, the large time limit is performed in the distribution
sense. Plugging this result in Eq. (32), we obtain
$|\tilde{\psi}_{0}\rangle\simeq\frac{1}{\sqrt{\mathcal{N}}}\left(|0\rangle+\int_{0}^{\pi}{\rm
d}k\,\mathcal{K}(k)|k,-k\rangle+\dots\right)$ with a renormalized wave
function
$\mathcal{K}(k)=-\text{i}\tan(\theta_{k}^{h_{\perp}}-\theta_{k}^{\bar{h}_{\perp}})-\text{i}h_{\parallel}\bar{\sigma}v(k)/\epsilon^{2}(k)$
(cf with Eq. (4)).

Finally, we can ri-exponentiate the state $\ket{\tilde{\psi}_{0}}$ to obtain
the squeezed form

$|\tilde{\psi}_{0}\rangle\simeq\frac{1}{\sqrt{\mathcal{N}}}\exp\left[\int_{0}^{\pi}{\rm
d}k\,\mathcal{K}(k)\hat{\gamma}^{\dagger}(k)\hat{\gamma}^{\dagger}(-k)\right]|0\rangle\,.$
(34)

Notice that Eq. (34) reduces to the linearized result above after a series
expansion to the first order. Higher order terms in the squeezed form (34) are
justified from the physical assumption that pairs of fermions are incoherently
created across the system. In Fig. 4, we put to test our quantum prediction in
Gaussian approximation obtained with the renormalized state (34) and that for
an initial state obtained from the wave function
$K(k)=-\text{i}\tan(\theta^{\bar{h}_{\perp}}_{k}-\theta^{h_{\perp}}_{k})$ of a
transverse field quench. The comparison of the curves with tensor network
simulations undoubtedly shows the relevance of the longitudinal field
corrections on the initial state for a correct description of the entanglement
dynamics.

Figure 4: Comparison of the tensor network simulations (symbols) with the
single-meson prediction obtained for two different initial wave function: (i)
the wave function
$K(k)=-\text{i}\tan(\theta^{\bar{h}_{\perp}}_{k}-\theta^{h_{\perp}}_{k})$ for
the transverse field quench with no additional corrections (dash-dot line);
(ii) the wave function ${\cal
K}(k)=-\text{i}\tan(\theta^{\bar{h}_{\perp}}_{k}-\theta^{h_{\perp}}_{k})-\text{i}h_{\parallel}\bar{\sigma}v(k)/\epsilon^{2}(k)$
in Eq. (4), obtained with a renormalized initial state (full line). The figure
clearly shows that only a properly renormalized initial state
$\ket{\tilde{\psi}_{0}}$ gives a good description of the entanglement
dynamics.

##  Appendix C Computation of entanglement from Gaussian correlations

By approximating the time-evolved state $|\psi_{t}\rangle$ with a squeezed
state, it is possible compute the entanglement entropy with the help of well
known formulas for free Fermi gases Casini and Huerta (2009). In order to do
this, we first diagonalize the two-particle problem in the zero-momentum
sector and subsequently we solve the Schrödinger equation (27). Since we
explicitly used the translational invariance of the two-particle state, we
found more convenient to work with a finite-size system having $N$ sites and
periodic boundaries. Hence, we compute the correlation matrix in Eq. (7)

$C=\int\frac{{\rm
d}k}{2\pi}\frac{e^{\text{i}k(j-j^{\prime})}}{1+|\mathcal{K}_{t}(k)|^{2}}U_{\theta_{k}}\begin{pmatrix}1&&\mathcal{K}_{t}(k)\\\
\mathcal{K}^{*}_{t}(k)&&|\mathcal{K}_{t}(k)|^{2}\end{pmatrix}U^{\dagger}_{\theta_{k}}\,,$
(35)

where we used
$\langle\hat{\gamma}^{\dagger}(q)\hat{\gamma}(k)\rangle=\delta(k-q)|\mathcal{K}(k)|^{2}/(1+|\mathcal{K}(k)|^{2})$,
$\langle\hat{\gamma}(q)\hat{\gamma}(-k)\rangle=\delta(k-q)\mathcal{K}(k)/(1+|\mathcal{K}(k)|^{2})$
and we replaced the integration over momenta with a Fourier series.
It is convenient to organize the correlations in the matrix form

$C=\begin{pmatrix}\mathbbm{1}-D&O\\\
O^{\dagger}&D\end{pmatrix}\,,\,\,\,D_{i,j}=\langle\hat{c}^{\dagger}_{i}\hat{c}_{j}\rangle\,,\,\,\,O_{ij}=\langle\hat{c}_{i}\hat{c}_{j}\rangle$
(36)

from which the entanglement entropy can be conveniently computed. Following
Ref. Casini and Huerta (2009), one notices that if the density matrix
$\hat{\rho}$ is Gaussian, then the reduced density matrix $\hat{\rho}_{A}$ is
also Gaussian and it is uniquely identified by the reduced correlation matrix
$C_{A}$. For the sake of concreteness, let us assume the subsystem $A$ is an
interval $j\in\\{1,...,L\\}$, but of course the method is more general. As a
next step, one considers a Bogoliubov rotation from the fermionic fields
$\hat{c}_{j}$, $\hat{c}_{j}^{\dagger}$ to modes $\hat{f}_{j}$,
$\hat{f}_{j}^{\dagger}$ that diagonalize the reduced density matrix

$\begin{pmatrix}\hat{c}_{1}\\\ \hat{c}_{2}\\\ ...\\\ \hat{c}_{L}\\\
\hat{c}^{\dagger}_{1}\\\ \hat{c}^{\dagger}_{2}\\\ ...\\\
\hat{c}^{\dagger}_{L}\end{pmatrix}=U\begin{pmatrix}\hat{f}_{1}\\\
\hat{f}_{2}\\\ ...\\\ \hat{f}_{L}\\\ \hat{f}^{\dagger}_{1}\\\
\hat{f}^{\dagger}_{2}\\\ ...\\\ \hat{f}^{\dagger}_{L}\end{pmatrix}$ (37)

where $U$ is a $2L\times 2L$ matrix. In order to preserve the canonical
commutation relations of the new fermionic modes, $U$ must be unitary. We fix
$U$ by requiring that the correlation matrix $C_{A}$ is brought in the
diagonal form

$U^{\dagger}C_{A}U=\begin{pmatrix}\mathbbm{1}-\Lambda&O\\\
O^{\dagger}&\Lambda\end{pmatrix}$ (38)

with $\Lambda$ a $L\times L$ diagonal matrix. Hence, $\hat{\rho}_{A}$ is also
diagonal in the new fermionic basis

$\hat{\rho}_{A}\propto\exp\left(\sum_{i=1}^{L}\log\left(\frac{\Lambda_{ii}}{1-\Lambda_{ii}}\right)\hat{f}^{\dagger}_{i}\hat{f}_{i}\right).$
(39)

From this result, the von Neumann and Rényi entropies of density matrices of
the form (39) are easily obtained as

$\begin{split}&S_{A}=-\text{Tr}_{A}(\hat{\rho}_{A}\log\hat{\rho}_{A})\\\
&\quad=-\sum_{i}\Big{\\{}\Lambda_{ii}\log\Lambda_{ii}+(1-\Lambda_{ii})\log(1-\Lambda_{ii})\Big{\\}}\\\
&\quad=-\text{Tr}_{A}[C_{A}\log C_{A}]\end{split}$ (40)

and

$S_{A}^{(n)}=\frac{1}{1-n}\log(\text{Tr}_{A}\
\hat{\rho}_{A}^{n})=\frac{1}{1-n}\sum_{i}\log[\Lambda_{ii}^{n}+(1-\Lambda_{ii})^{n}].$
(41)

Crucially, we assume periodic boundary conditions on the fermion basis when
computing (35), which imply a quantization of momenta as $k=2\pi n/N$. This
introduces some subtleties in the computation of the Von Neumann and Rényi
entropies, as we now discuss. As already mentioned in Appendix A, choosing
periodic boundary conditions amounts to the choice of the parity sector with
$P=-1$. Namely, the ground state correlation functions (obtained setting
$\mathcal{K}=0$ in Eq. (35)) are not those of $|{\rm GS}_{\uparrow}\rangle$,
but rather those of the antisymmetric combination
$|0^{-}\rangle=\frac{1}{\sqrt{2}}(|{\rm GS}_{\uparrow}\rangle-|{\rm
GS}_{\downarrow}\rangle)$. Before considering excited states and dynamics, it
is instructive to compare the entanglement entropy obtained from the states
$|{\rm GS}_{\uparrow}\rangle$ and $|0^{-}\rangle$. More precisely, let us
consider a bipartition $A\cup B$ where $A$ is an interval of length $L<N$ and
let $S^{-}_{L}$ ( resp. $S^{\uparrow}_{L}$) be the entanglement entropy (which
can be either von Neumann or Rényi, the forthcoming argument holds in both
cases) of such a bipartition on the state $|0^{-}\rangle$ (resp. $|{\rm
GS}_{\uparrow}\rangle$). With simple symmetry arguments between up and down
spins, it is then easy to show that

$S^{-}_{L}=S^{\uparrow}_{L}+\log 2\,.$ (42)

As a simple check, one can set $h_{\perp}=0$ and obtain Eq. (42) from the two
degenerated ground states, which now take the form of product states with all
spins up or down. As a next step, we consider the excitations and their
dynamics. It is important to stress that, depending on the initial choice of
the ground state ($|{\rm GS}_{\uparrow}\rangle$ or $|{\rm
GS}_{\downarrow}\rangle$), one obtains very different excitations and,
consequently, a different dynamics. However, by noticing that the Schrödinger
equation (27) describes the dynamics of excitations over the state $|{\rm
GS}_{\uparrow}\rangle$, it is easy to see that using $\mathcal{K}_{t}$ in Eq.
(35), we are not evolving $|0^{-}\rangle$, but we are rather separately
evolving $|{\rm GS}_{\uparrow}\rangle$ with the Hamiltonian (1) and $|{\rm
GS}_{\downarrow}\rangle$ with a flipped sign in the longitudinal field
$h_{\parallel}\to-h_{\parallel}$. Therefore, the identity Eq. (42) remains
valid also at finite times. We then use the methods explained in Sec. III.1 to
compute $S^{-}_{A}$ at all times from the correlation matrix in Eq. (35), and
we obtain $S^{\uparrow}_{A}$ by subtracting the constant offset $\log 2$. At
this point, we have the entanglement entropy of a finite interval $L$ embedded
in a system with periodic boundary conditions. On the other hand, we performed
tensor network simulations to compute the half-size entanglement entropy of a
finite system with open boundary conditions. To compare the two, we notice
that the entanglement entropy obeys the area law, from which it follows that

$S_{\text{half-sys; OBC}}=\frac{1}{2}\lim_{L\to\infty}S_{L}^{\uparrow}\,.$
(43)

In practice, saturation is obtained as soon as $L$ exceeds the maximum meson
size $d_{\text{max}}=4h_{\perp}/(|h_{\parallel}|\bar{\sigma})$.

##  Appendix D Confinement in the long-range Ising chain

The dynamics of entanglement in the long-range Ising chain in Eq. (9) can be
obtained with the same strategy employed for the tilted Ising chain in Eq.
(1). Following Ref. Liu  _et al._ (2019), we focus on the regime of weak
transverse field $h_{\perp}\ll 1$, where the excitations of the model (9) are
identified by domain walls. Restricting to the two domain walls subspace, the
effective Hamiltonian is given by

$\begin{split}\hat{H}_{\text{eff}}|j_{1},j_{2}\rangle&=-h_{\perp}\Big{[}|j_{1}+1,j_{2}\rangle+|j_{1}-1,j_{2}\rangle+|j_{1},j_{2}+1\rangle\\\\[4.0pt]
&+|j_{1},j_{2}-1\rangle\Big{]}+U(|j_{1}-j_{2}|)|j_{1},j_{2}\rangle\,,\end{split}$
(44)

where
$|j_{1},j_{2}\rangle\equiv|\uparrow...\uparrow_{j_{1}}\downarrow..\downarrow_{j_{2}-1}\uparrow_{j_{2}}...\rangle$
and, from hereafter, we consider an ordered set of coordinates with hard core
constraint $j_{1}<j_{2}$ for the wave function. The two-kink interaction is

$U(j)=4j\zeta(\alpha)-4\sum_{1\leq\ell<n}\sum_{1\leq
r\leq\ell}\frac{1}{r^{\alpha}}\,,$ (45)

with $\zeta(\alpha)$ the Riemann-zeta function. Notice that the potential
$U(j)$ diverges at large distances whenever $1<\alpha<2$. We are now
interested to the dynamics of a ferromagnetic state with Hamiltonian (9) in
the presence of a small, but finite, transverse field $h_{\perp}$.
Equivalently, we consider a small homogeneous quench in the transverse field.
Therefore, one first considers the dynamics in the zero-momentum sector by
defining $|j\rangle=\sum_{j^{\prime}}|j,j+j^{\prime}\rangle$. Similarly to the
case of the Ising chain in a tilded field, we introduce the wave function
$W_{t}(j)$ of two domain walls in the zero-momentum sector, satisfying the
Schrödinger equation

$\text{i}\partial_{t}W_{t}(j)=-2h_{\perp}\big{[}W_{t}(j+1)+W_{t}(j-1)\big{]}+U(j)W_{t}(j)\,.$
(46)

In the last expression, one should restrict to $j>0$ with the condition
$W_{t}(j=0)=0$, due to the aforementioned choice of coordinates. In
particular, a convenient way to automatically encode the hard core constraint
is to consider an antisymmetric extension of the wave function to negative
coordinates, i.e. imposing that $W_{t}(-j)=-W_{t}(j)$. When performing this
operation, the potential is then symmetrically extended as $U(-j)=U(j)$. In
other words, we are building a fermionic extension of the original wave
function that is very useful for our scopes, as we now discuss. Let
$\Phi_{n}(j)$ be the normalized eigenfunctions in the antisymmetric sector of
energy $E_{n}$, thus satisfying the time-independent Schrödinger equation

$E_{n}\Phi_{n}(t)=-2h_{\perp}\big{[}\Phi_{n}(j+1)+\Phi_{n}(j-1)\big{]}+U(j)\Phi_{n}(j)\,.$
(47)

For fixed initial conditions, the wave function evolves as
$W_{t}(j)=\sum_{n}e^{-\text{i}E_{n}t}c_{n}\Phi_{n}(j)$ with time independent
coefficients $c_{n}$. However, in the case at hand, such coefficients become
time dependent $c_{n}\to c_{n}(t)$ since domain walls are _dynamically_
created from the fully-polarised state due to the transverse field. To compute
$c_{n}(t)$, we proceed similarly to the short range Ising case and split the
Hamiltonian (9) as $\hat{H}=\hat{H}_{\text{diag}}+\hat{H}_{\text{off-diag}}$,
where $\hat{H}_{\text{off-diag}}$ is obtained by subtracting from
$-h_{\perp}\sum_{i}\hat{\sigma}_{i}^{x}$ the projected component acting within
the sector with a conserved number of domain walls. Then, we move to the
interaction picture and consider first order perturbation theory by computing
the matrix element of $-\text{i}\int_{0}^{t}{\rm
d}\tau\,\hat{H}^{(I)}_{\text{off-diag}}(\tau)$ between the polarised state
$|...\uparrow\uparrow\uparrow...\rangle$ and the eigenfunctions $\Phi_{n}$,
reaching the simple expression

$W_{t}(j)=\sum_{n}e^{-\text{i}E_{n}t}c_{n}(t)\Phi_{n}(j)$ (48)

with coefficients

$c_{n}(t)=h_{\perp}\frac{1-e^{\text{i}tE_{n}}}{E_{n}}\sqrt{2}\Phi_{n}(1)\,.$
(49)

The energies $E_{n}$ and eigenfunctions $\Phi_{n}$ are obtained numerically by
solving Eq. (47). As a last step, we now wish to build the many-body wave
function. In particular, we look for an effective description in terms of
fermionic Gaussian states such that the formulas of Appendix C for the
entanglement remain still valid. However, in the long range setting, this
requires to consider fictitious fermionic degrees of freedom instead of the
standard fermionic excitations of the model (i.e., those obtainable through a
Jordan-Wigner transformation). More precisely, for small transverse field
$|h_{\perp}|\ll 1$, we notice that the fundamental excitations of the long
range model are domain walls, similarly to those arising in the short range
Ising model discussed previously. Therefore, one can introduce fictitious
fermionic degrees of freedom for the long range problem as the domain wall
excitations in the short range Ising model, and use these fictitious fermionic
degrees of freedom to apply the strategy of Appendix C for the calculation of
the entanglement. Following this program, we write the two-body state in
momentum space as a squeezed state of mode operators $\hat{\gamma}(k)$

$|\tilde{\psi}_{0}\rangle\propto\exp\left[\int_{0}^{\pi}{\rm
d}k\,\mathcal{K}(k)\hat{\gamma}^{\dagger}(k)\hat{\gamma}^{\dagger}(-k)\right]|0\rangle\,,$
(50)

where $\mathcal{K}(k)$ is defined as the Fourier transform of the real space
wave function $W_{t}(j)$. Finally, the real space correlation functions of the
fictitious fermions must be computed. To this end, one can simply use the
expression in Eq. (35) for the Ising model with the Bogoliubov angle
$\theta_{k}=-\frac{1}{2\text{i}}\log(-e^{\text{i}k})$, which is nothing else
than the limit $h_{\perp}\to 0$ of Eq. (17). From this result, one can obtain
the entanglement entropy with the techniques of Appendix C.

##  Appendix E Breakdown of the single-meson description

Figure 5: (a) Large quench in the transverse field: $\bar{h}_{\perp}=0\to
h_{\perp}=0.5$, at $h_{\parallel}=0.1$. Inset: exact density of quasiparticles
$n_{k}$ (full line) significantly deviates from its approximation (dashed-dot
line) in the dilute regime. In this regime, both analytical predictions for
the Von Neumann entropy (quasiparticle picture – dashed line; quantum Gaussian
treatment – full line) fail to predict the observed numerical behavior
(symbols) at finite times. (b) Example of Schwinger mechanism:
$\bar{h}_{\perp}=0.5\to h_{\perp}=0.25$ and strong confinement
$h_{\parallel}=0.4$. The analytical predictions (quasiparticle picture –
dashed line; quantum Gaussian treatment – thick full line) significantly
deviates from the numerical data (thin full line). The difference $\delta
S_{A}$ between the numerical data and the quantum prediction clearly shows a
linear growth of entanglement due to the spontaneous creation/annihilation of
new particles. Figure 6: (a) Single spin flip energy profile of a finite-size
chain with $N=150$ sites for different values of $\alpha$. For $\alpha\lesssim
2$, the profile significantly deviates from a flat one and leads to a non-
negligible boundary-meson interaction term. In each of the curves, we subtract
the max value of the energy for a better comparison. (b) Rényi-2 entropy
dynamics for the long-range Ising model (see Eq. (9)) with exponent
$\alpha\leq 2$ obtained with variational methods (dashed line). The emergence
of meson-meson and meson-boundary interactions produces a linear drift which
is not captured by the quantum prediction in single-meson approximation (full
line). In each of the right panels, we set the transverse field $h_{\perp}$
such that $h_{\perp}/U(1)\ll 1$.

In Sec. III, we derived an analytical prediction for the entanglement dynamics
of confined spin chains under the assumption of small quenches, which
physically corresponds to mesons that are well-separated and thus
noninteracting. This approximation is found in excellent agreement with the
numerical data for a wide set of quench parameters, see Fig. 2. Nevertheless,
corrections appear when the typical distance among mesons becomes comparable
with their size or with the interaction range, as we now discuss. To begin
with, we consider such deviations in the short-range Ising model with tilted
magnetic field (see Eq. (1)), where the assumption of well-separated mesons is
easier to fulfill than in the long range case. For this model, inter-mesons
interactions are possible only through contact, therefore corrections are
present only when the typical size of the meson, here estimated with the
semiclassical maximum size
$d_{\text{max}}=4h_{\perp}/(h_{\parallel}\bar{\sigma})$, becomes comparable
with the average spacing, the latter being the inverse of the meson density
$n_{\text{tot}}=\int_{0}^{\pi}\frac{{\rm d}k}{2\pi}\,n_{k}$. In Fig. 5 (a) we
consider a large quench with $d_{\text{max}}n_{\text{tot}}\simeq 0.3$.
Semiclassically, pairs of fermions are originated at the same point in space,
hence the size of mesons starts from zero and increases with time.
Consequently, we observe an initial agreement with great accuracy between the
tensor network data and the single-meson prediction, but the two drift away
for $t>t^{*}$, with $t^{*}\simeq 4$ with this choice of parameters. The regime
of dilute meson can be better attained by either diminishing the density
$n_{\text{tot}}$ or either reducing the size of mesons $d_{\text{max}}$ by
increasing the longitudinal field. However, a too strong longitudinal field
leads as well to a breakdown the single-meson approximation due to the
Schwinger effect, through which the initial mesons decay in couples of lighter
mesons at a constant rate. In Fig. 5 (b), we probed such quasiparticles
production by considering a quench with strong confinement. We observe a
linear growth of the deviation from the stable meson approximation.
We now refer to Fig. 6 for the long range case, where finite-size corrections
to the single-meson predictions are particularly visible for small exponents
$\alpha\lesssim 2$. As an example, in Fig. 6 (a), we consider the energy of an
isolated spin flip in a system of finite length $L=150$. The finite size of
the system induces a non trivial potential landscape and the meson (or, more
precisely, the deep boundstate for $\alpha>2$) is initially created at rest
and later accelerated due to the effective inhomogeneous potential. This
causes a scattering of the excitations which results in a linear growth of
entanglement on intermediate time scales, see Fig. 6 (b). Such finite size
effects are suppressed at larger $\alpha$, see Fig. 3.
Although, strictly speaking, $\alpha>2$ does not induce confining dynamics, we
experienced that setting $2.5\leq\alpha\leq 3$ is a good compromise between
having deep bound states and relatively small boundary effects on the sizes we
analyzed. Furthermore, we notice that, even in the ideal scenario of an
infinite system, similar long-tailed interactions will show up among mesons
(or bound states). Hence, we conclude that the dilute regime is attained when
the distance among the excitations is large compared to the $\alpha-$dependent
interaction range. In Fig. 6 (b) we consider the entanglement growth after
transverse field quenches for values of $\alpha$ smaller than those considered
in Sec. IV. In particular, corrections with respect to our predictions are
more relevant as $\alpha$ is reduced, in agreement with the phenomenology we
depicted.

## References

  * Feynman (2018) R. P. Feynman, in  _Feynman and Computation_ (CRC Press, Boca Raton, 2018) pp. 133–153.
  * Georgescu  _et al._ (2014) I. M. Georgescu, S. Ashhab, and F. Nori, Rev. Mod. Phys. 86, 153 (2014).
  * McCoy and Wu (1978) B. M. McCoy and T. T. Wu, Phys. Rev. D 18, 1259 (1978).
  * Martinez  _et al._ (2016) E. A. Martinez, C. A. Muschik, P. Schindler, D. Nigg, A. Erhard, M. Heyl, P. Hauke, M. Dalmonte, T. Monz, P. Zoller, and R. Blatt, Nature 534, 516 (2016).
  * Tan  _et al._ (2021) W. L. Tan, P. Becker, F. Liu, G. Pagano, K. S. Collins, A. De, L. Feng, H. B. Kaplan, A. Kyprianidis, R. Lundgren, W. Morong, S. Whitsitt, A. V. Gorshkov, and C. Monroe, Nature Phys. 17, 742 (2021).
  * Delfino  _et al._ (1996) G. Delfino, G. Mussardo, and P. Simonetti, Nucl. Phys. B 473, 469 (1996).
  * Rutkevich (2008) S. B. Rutkevich, J. Stat. Phys. 131, 917 (2008).
  * Fonseca and Zamolodchikov (2003) P. Fonseca and A. Zamolodchikov, J. Stat. Phys. 110, 527 (2003).
  * Kormos  _et al._ (2017) M. Kormos, M. Collura, G. Takács, and P. Calabrese, Nature Phys. 13, 246 (2017).
  * Robinson  _et al._ (2019) N. J. Robinson, A. J. A. James, and R. M. Konik, Phys. Rev. B 99, 195108 (2019).
  * James  _et al._ (2019) A. J. A. James, R. M. Konik, and N. J. Robinson, Phys. Rev. Lett. 122, 130603 (2019).
  * Kogut (1979) J. B. Kogut, Rev. Mod. Phys. 51, 659 (1979).
  * Lerose  _et al._ (2020) A. Lerose, F. M. Surace, P. P. Mazza, G. Perfetto, M. Collura, and A. Gambassi, Phys. Rev. B 102, 041118 (2020).
  * Borla  _et al._ (2020) U. Borla, R. Verresen, F. Grusdt, and S. Moroz, Phys. Rev. Lett. 124, 120503 (2020).
  * Chanda  _et al._ (2020) T. Chanda, J. Zakrzewski, M. Lewenstein, and L. Tagliacozzo, Phys. Rev. Lett. 124, 180602 (2020).
  * Pai and Pretko (2020) S. Pai and M. Pretko, Phys. Rev. Research 2, 013094 (2020).
  * Mazza  _et al._ (2019) P. P. Mazza, G. Perfetto, A. Lerose, M. Collura, and A. Gambassi, Phys. Rev. B 99, 180302 (2019).
  * Bastianello  _et al._ (2021) A. Bastianello, U. Borla, and S. Moroz, “Fragmentation and emergent integrable transport in the weakly tilted ising chain,” (2021), arXiv:2108.04845 .
  * Schwinger (1951) J. Schwinger, Phys. Rev. 82, 664 (1951).
  * Sinha  _et al._ (2021) A. Sinha, T. Chanda, and J. Dziarmaga, Phys. Rev. B 103, L220302 (2021).
  * Lagnese  _et al._ (2021a) G. Lagnese, F. M. Surace, M. Kormos, and P. Calabrese, Phys. Rev. B 104, L201106 (2021a).
  * Tortora  _et al._ (2020) R. J. V. Tortora, P. Calabrese, and M. Collura, EPL 132, 50001 (2020).
  * Rigobello  _et al._ (2021) M. Rigobello, S. Notarnicola, G. Magnifico, and S. Montangero, Phys. Rev. D 104, 114501 (2021).
  * Verdel  _et al._ (2020) R. Verdel, F. Liu, S. Whitsitt, A. V. Gorshkov, and M. Heyl, Phys. Rev. B 102, 014308 (2020).
  * Myers  _et al._ (2017) R. C. Myers, M. Rozali, and B. Way, J. Phys. A: Math. Theor. 50, 494002 (2017).
  * Cubero and Robinson (2019) A. C. Cubero and N. J. Robinson, J. Stat. Mech.: Theory and Experiment 2019, 123101 (2019).
  * Vovrosh and Knolle (2021) J. Vovrosh and J. Knolle, Sci. Rep. 2021 11:1 11, 1 (2021).
  * Surace and Lerose (2021) F. M. Surace and A. Lerose, New J. Phys. 23, 062001 (2021).
  * Karpov  _et al._ (2020) P. I. Karpov, G. Y. Zhu, M. P. Heller, and M. Heyl, “Spatiotemporal dynamics of particle collisions in quantum spin chains,” (2020), arXiv:2011.11624 .
  * Milsted  _et al._ (2021) A. Milsted, J. Liu, J. Preskill, and G. Vidal, “Collisions of false-vacuum bubble walls in a quantum spin chain,” (2021), arXiv:2012.07243 .
  * Vovrosh  _et al._ (2021) J. Vovrosh, H. Zhao, J. Knolle, and A. Bastianello, “Confinement induced impurity states in spin chains,” (2021), arXiv:2108.03976 .
  * Schollwöck (2011) U. Schollwöck, Annals of Physics 326, 96 (2011).
  * Castro-Alvaredo  _et al._ (2020) O. A. Castro-Alvaredo, M. Lencsés, I. M. Szécsényi, and J. Viti, Phys. Rev. Lett. 124, 230601 (2020).
  * Calabrese and Cardy (2006) P. Calabrese and J. Cardy, Phys. Rev. Lett. 96, 136801 (2006).
  * Calabrese and Cardy (2005) P. Calabrese and J. Cardy, J. Stat. Mech.: Theory and Experiment 2005, P04010 (2005).
  * Calabrese (2020) P. Calabrese, SciPost Phys. Lect. Notes , 20 (2020).
  * Alba and Calabrese (2017) V. Alba and P. Calabrese, PNAS 114, 7947 (2017).
  * Calabrese  _et al._ (2011) P. Calabrese, F. H. L. Essler, and M. Fagotti, Phys. Rev. Lett. 106, 227203 (2011).
  * Delfino (2014) G. Delfino, J. Phys. A: Math. Theor. 47, 402001 (2014).
  * Delfino and Viti (2017) G. Delfino and J. Viti, J. Phys. A: Math. Theor. 50, 084004 (2017).
  * Islam  _et al._ (2015) R. Islam, R. Ma, P. M. Preiss, M. Eric Tai, A. Lukin, M. Rispoli, and M. Greiner, Nature 528, 77 (2015).
  * M.  _et al._ (2016) K. A. M., T. M. Eric, L. Alexander, R. Matthew, S. Robert, P. P. M., and G. Markus, Science 353, 794 (2016).
  * Elben  _et al._ (2018) A. Elben, B. Vermersch, M. Dalmonte, J. I. Cirac, and P. Zoller, Phys. Rev. Lett. 120, 050406 (2018).
  * Brydges  _et al._ (2019) T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F. Roos, Science 364, 260 (2019).
  * Elben  _et al._ (2020) A. Elben, R. Kueng, H.-Y. R. Huang, R. van Bijnen, C. Kokail, M. Dalmonte, P. Calabrese, B. Kraus, J. Preskill, P. Zoller, and B. Vermersch, Phys. Rev. Lett. 125, 200501 (2020).
  * Neven  _et al._ (2021) A. Neven, J. Carrasco, V. Vitale, C. Kokail, A. Elben, M. Dalmonte, P. Calabrese, P. Zoller, B. Vermersch, R. Kueng, and B. Kraus, npj Quantum Information 7, 152 (2021).
  * Vitale  _et al._ (2021) V. Vitale, A. Elben, R. Kueng, A. Neven, J. Carrasco, B. Kraus, P. Zoller, P. Calabrese, B. Vermersch, and M. Dalmonte, “Symmetry-resolved dynamical purification in synthetic quantum matter,” (2021), arXiv:2101.07814 .
  * Bertini  _et al._ (2018) B. Bertini, E. Tartaglia, and P. Calabrese, J. Stat. Mech.: Theory and Experiment 2018, 063104 (2018).
  * Bastianello and Calabrese (2018) A. Bastianello and P. Calabrese, SciPost Phys. 5, 33 (2018).
  * Bastianello and Collura (2020) A. Bastianello and M. Collura, SciPost Phys. 8, 45 (2020).
  * Pomponio  _et al._ (2021) O. Pomponio, M. A. Werner, G. Zarand, and G. Takacs, “Bloch oscillations and the lack of the decay of the false vacuum in a one-dimensional quantum spin chain,” (2021), arXiv:2105.00014 .
  * Calabrese and Cardy (2004) P. Calabrese and J. Cardy, J. Stat. Mech.: Theory and Experiment 2004, P06002 (2004).
  * Calabrese and Cardy (2009) P. Calabrese and J. Cardy, J. Phys. A: Math. Theor. 42, 504005 (2009).
  * Casini and Huerta (2009) H. Casini and M. Huerta, J. Phys. A: Math. Theor. 42, 504007 (2009).
  * Yang and White (2020) M. Yang and S. R. White, Phys. Rev. B 102, 094315 (2020).
  * Liu  _et al._ (2019) F. Liu, R. Lundgren, P. Titum, G. Pagano, J. Zhang, C. Monroe, and A. V. Gorshkov, Phys. Rev. Lett. 122, 150601 (2019).
  * Lagnese  _et al._ (2021b) G. Lagnese, F. M. Surace, M. Kormos, and P. Calabrese, “Quenches and confinement in a heisenberg-ising spin ladder,” (2021b), arXiv:2110.09472 .
  * Lagnese  _et al._ (2020) G. Lagnese, F. M. Surace, M. Kormos, and P. Calabrese, J. Stat. Mech.: Theory and Experiment 2020, 093106 (2020).
  * Lencsés and Takács (2015) M. Lencsés and G. Takács, J. High Energ. Phys. 2015, 1 (2015).
  * Delfino and Grinza (2008) G. Delfino and P. Grinza, Nucl. Phys. B 791, 265 (2008).
  * Lencsés  _et al._ (2021) M. Lencsés, G. Mussardo, and G. Takács, “Confinement in the tricritical ising model,” (2021), arXiv:2111.05360 .
  * Moca  _et al._ (2017) C. P. Moca, M. Kormos, and G. Zaránd, Phys. Rev. Lett. 119, 100603 (2017).
  * Fishman  _et al._ (2020) M. Fishman, S. R. White, and E. M. Stoudenmire, “The ITensor software library for tensor network calculations,” (2020), arXiv:2007.14822 .