# Constrained dynamics and confinement in the two-dimensional quantum Ising
model

Luka Pavešić<EMAIL_ADDRESS>Dipartimento di Fisica e Astronomia “G.
Galilei”, via Marzolo 8, I-35131 Padova, Italy  Istituto Nazionale di Fisica
Nucleare (INFN), Sezione di Padova, I-35131 Padova, Italy  Daniel Jaschke
Dipartimento di Fisica e Astronomia “G. Galilei”, via Marzolo 8, I-35131
Padova, Italy  Istituto Nazionale di Fisica Nucleare (INFN), Sezione di
Padova, I-35131 Padova, Italy  Institute for Complex Quantum Systems, Ulm
University, Albert-Einstein-Allee 11, 89069 Ulm, Germany  Simone Montangero
Dipartimento di Fisica e Astronomia “G. Galilei”, via Marzolo 8, I-35131
Padova, Italy  Istituto Nazionale di Fisica Nucleare (INFN), Sezione di
Padova, I-35131 Padova, Italy

###### Abstract

We investigate the dynamics of the quantum Ising model on two-dimensional
square lattices up to $16\times 16$ spins. In the ordered phase, the model is
predicted to exhibit dynamically constrained dynamics, leading to confinement
of elementary excitations and slow thermalization. After demonstrating the
signatures of confinement, we probe the dynamics of interfaces in the
constrained regime through sudden quenches of product states with domains of
opposite magnetization. We find that the nature of excitations can be captured
by perturbation theory throughout the confining regime, and identify the
crossover to the deconfining regime. We systematically explore the effect of
the transverse field on the modes propagating along flat interfaces and
investigate the crossover from resonant to diffusive melting of a square of
flipped spins embedded in a larger lattice.

While interacting quantum systems typically relax to equilibrium relatively
quickly, there are striking examples where this is not the case [1, 2].
Investigating the origins of the anomalously long prethermal regime of such
systems has applications in understanding exotic phases of condensed matter
[3, 4], phenomena like false vacuum decay [5] and quark confinement [6].

One possible reason for slow thermalization is the emergence of dynamical
constraints which approximately restrict the time evolution to certain
subspaces of the full Hilbert space. A prominent example are systems with
confinement [7], where pairs of elementary excitations experience attraction
which increases with increasing separation, and thus form long-lived bound
states. This is, for example, the reason why quarks at low temperatures do not
exist as isolated particles, but are bound in mesons or hadrons.

An analogous mechanism is also found in Ising spin chains [8]. At small
transverse fields, the low energy excitations are domain walls, and pairs of
them feel a confining force provided by the longitudinal field [9, 10, 11, 12,
13] or interactions beyond nearest neighbour [14]. The direct correspondence
between confined domain walls in spin systems and pairs of charges confined by
a gauge field in lattice gauge theories makes the Ising chain a useful toy
model for studying thermalization of bound states via meson production [15,
16, 17, 18]. Domains of spins anti-aligned with the longitudinal field are
bubbles of false vacuum, and probing their dynamics is crucial for
understanding false vacuum decay [19, 20].

In two dimensions (2D), the confining potential comes from the lattice itself.
An excitation has to flip neighbouring spins to propagate, thus adding domain
walls. At small transverse fields such processes are very energetically
costly, and thus strongly suppressed. Confinement has initially been studied
in ladders [21, 22], by orthogonally coupling many chains described by the
Ising field theory in the scaling limit [23, 7] and recently in a 2D lattice
with infinite tensor network methods [24]. These approaches describe the
system in the thermodynamic limit and hence cannot probe dynamics of
inhomogeneous initial states and domain walls.

In this paper, we use tree tensor networks (TTNs) to explore the constrained
dynamics of the quantum Ising model on 2D square lattices of up to $16\times
16$ spins, through quenches of product states to finite transverse fields.
Starting with the completely polarized state, we demonstrate signatures of
confinement and investigate the nature of low-energy excitations that drive
the dynamics. Then, we turn our attention to two representative types of
domain walls, a flat interface between two spin species and a square domain
embedded in a larger lattice, to explore the melting of interfaces in the
confining regime. We characterize the surface modes that are known to dominate
the dynamics in the limit of small transverse field [25, 26], and investigate
the crossover towards diffusive melting. We find that the crossover
corresponds to the dynamical phase transition from the ordered to disordered
phase and we characterize the nature of domain melting in both regimes.
Finally, we discuss the feasibility of implementing our numerical calculations
on quantum simulators and suggest an outlook of possible extensions of the
presented work.

Model and methods \- The Hamiltonian of the two-dimensional quantum Ising
model is

$H=-J\sum_{\langle
i,j\rangle}\sigma^{z}_{i}\sigma^{z}_{j}-g\sum_{i}\sigma^{x}_{i},$ (1)

with $x$\- and $z$-Pauli matrices $\sigma^{x}$ and $\sigma^{z}$, $J$ the
interaction, and $g$ the transverse field strength. Sum indices run over sites
of a $N\times N$ square lattice, with $\langle i,j\rangle$ indicating a sum
over nearest neighbours. We choose to work with periodic boundary conditions
throughout, but the physics observed here does not depend on this choice.

Throughout the paper, we investigate instantaneous quenches of product states
with spins aligned in the $z$-direction (eigenstates at $g=0$) to finite
values of $g$. We encode them into tree tensor networks [27, 28, 29] and
propagated in real time using the time-dependent variational principle (TDVP)
[30, 31], with a timestep of $0.005J$. Our results are generally sufficiently
converged for the maximal bond dimension $\chi=100$. See Supplemental Material
(SM) for technical details and convergence analysis.

Broadly speaking, the accuracy of tensor network methods is bound by
entanglement. After instantaneous quenches of product states it typically
grows, at worst linearly with time [32]. This means that the resources
required for accurate tensor network calculations grow exponentially. However,
as shown in this work and elsewhere [9, 7, 24], the physical regime
investigated here is noteworthy precisely because entanglement grows slowly in
time. This makes TTNs a particularly apt tool to attack the present problem.

We characterize the physics through various observables. Aside from the local
magnetization $\left\langle\sigma^{z}_{i}\right\rangle$, we use the spin-spin
correlation function
$C_{ij}=\langle\sigma^{z}_{i}\sigma^{z}_{j}\rangle-\langle\sigma^{z}_{i}\rangle\langle\sigma^{z}_{j}\rangle$
where $i$ and $j$ are lattice sites of the two spins. When quenching from a
product state where initially $C_{ij}=0$ for all $i\neq j$, the time evolution
of $C_{ij}$ quantifies the spread of an excitation in time.

The entanglement of two parts of the system is quantified through the von
Neumann entropy $S=-\sum_{i}\alpha_{i}^{2}\ln\alpha_{i}^{2}$, where
$\alpha_{i}$ are the singular values of the Schmidt decomposition of the
chosen bipartition. $S$ is also a measure of the spread of excitations across
the bipartition, as this process entangles the two parts of the system.

Figure 1:  Confinement in the 2D quantum Ising model. (a) Short time
dependence of magnetization $\left\langle\sigma^{z}(t)\right\rangle$ for a set
of $g/J$. (b) $\left\langle\sigma^{z}(t)\right\rangle$ for longer times and
small $g/J$. (c) Time dependence of entanglement entropy $S$ for a bipartition
of two neighbouring spins and the rest of the system. $S_{\mathrm{max}}$ is
the maximal entanglement entropy of the bipartition (in this case
$S_{\mathrm{max}}=2\ln 2$). (d) Time dependence of a horizontal cut of the
connected correlation function $C_{ij}$ where $i=4$ for a set of $g/J$.

Confinement \- Slow thermalization of a many-body system is reflected in
persistent oscillations of local observables, like the magnetization of one
spin. Additionally, a striking signature of confined dynamics is the
suppression of the light-cone spread of correlations, and slow and non-
monotonous growth of entanglement with time [9]. We start by demonstrating
these properties after a quench of a completely polarized product state
$\left|\downarrow\dots\downarrow\right\rangle$ to finite $g$. The time
dependence of the magnetization $\left\langle\sigma_{z}\right\rangle$ is
presented in Figs. 1(a,b). For small $g/J\lesssim 1.5$, we find persistent
oscillations with no signs of suppression up to large times. However, when
increasing $g$ towards $g/J\approx 2$, $\left\langle\sigma_{z}\right\rangle$
equilibrates on a timescale of $1/J$. This indicates the onset of a regime of
quick thermalization, in rough agreement with the critical value of the
transverse field for a dynamical phase transition in the same model,
$g_{c}=2.0J$ [33].

The time evolution of the entanglement entropy $S$ is shown in Fig. 1(c). We
choose a bipartition between two neighbouring spins and the rest of the
system, but all other bipartitions produce similar results. At small $g/J$, we
observe an oscillatory pattern with the system periodically returning very
close to a product state with a frequency of approximately $\frac{\pi}{2}J$.
Interestingly, the frequency is half the frequency of the oscillations of
$\left\langle\sigma_{z}\right\rangle$. The oscillations are washed out as $g$
is increased and the density of excitations becomes large enough that they
interact, leading to thermalization [34].

The spin-spin correlations also exhibit signatures of confinement, see Fig.
1(d), where we plot a vertical cut of $C_{ij}$ with spin $i$ chosen in the
fourth row. The spread of correlations is strongly suppressed for small $g/J$.
At $g/J=0.5$, we find non-monotonous behaviour with the same frequency as the
oscillations of entanglement. The light-cone component becomes more apparent
as $g$ is increased, and at $g/J=1.5$ it is clearly the dominant contribution.
It should be mentioned that the light-cone never vanishes completely, but is
only suppressed. Its apparent disappearance is an artefact of the choice of a
plotting scale [9, 14].

Figure 2:  Spectrum of $\sigma^{z}$, given by the Fourier transform
$\mathcal{F}(\sigma^{z})$. (a) Spectrum of magnetization for two cases of
$g/J$. Dashed vertical lines correspond to $\omega/J=4,8$, and $12$. (b)
Heatmap of the spectrum for a range of $g/J$. Black dashed lines correspond to
the transition energies obtained by perturbation theory.

Nature of excitations \- The oscillating behaviour of local observables
implies that the dynamics is generated by repeated creation and annihilation
of local excitations. Their nature can be inferred from the spectral density
of local observables [9, 34, 14, 24]. In Fig. 2 we show the results of the
Fourier transform of the magnetization from Fig. 1(b) between $tJ=10$ and
$tJ=80$, denoted by $\mathcal{F}(\sigma^{z})$. We multiply the timeseries by
the Hamming window to exclude possible transient events at the edges of the
interval [35, 36]. An example of the spectral density for two values of $g/J$
is shown in Fig. 2(a), while a heatmap for a broader range is in Fig. 2(b).

The three peaks match the excitation energies of the dominant dynamical
processes. The most prominent peak, starting at $\omega=8J$ for $g=0$,
corresponds to the transition between the initial state and a state with one
flipped spin, $\psi_{1}$. The peak at $12J$ corresponds to a transition
between $\psi_{0}$ and $\psi_{2}$, a state where two neighbouring spins are
flipped, with a total domain wall length of $6$. Finally, the peak at $4J$ is
produced by the transition between $\psi_{1}$ and $\psi_{2}$.

The $g$-dependence of the excitation energies can be captured by perturbative
corrections to the product states. The expressions up to the second order in
$g$ are:

$\begin{split}E_{0}&=-\frac{g^{2}}{8J}N^{2}+\mathcal{O}(g^{4}),\\\
E_{1}&=8J-\frac{g^{2}}{8J}\left(N^{2}+6\right)+\mathcal{O}(g^{4}),\\\
E_{2}&=12J-\frac{g^{2}}{8J}N^{2}+\mathcal{O}(g^{4}).\end{split}$ (2)

Note that $N$ is the linear size of our system, which consists of $N^{2}$
sites in total. For complete calculations see the SM.

The black dashed lines in Fig. 2(b) indicate the transition energies obtained
as differences of these energies, denoted $\Delta_{ij}=E_{j}-E_{i}$. They
match the positions of the spectral peaks to surprisingly large values of
transverse field, up to $g\sim J$. This is somewhat surprising for expressions
which are technically only accurate up to the second order in $g$, and
indicates that the contribution of the higher excited modes is negligible
throughout the confining regime.

Note that aside from the extensive correction common to all states, we do not
find a second-order contribution to $E_{2}$. An intuitive interpretation is
that in second order the corrections that come from coupling $\psi_{2}$ to
$\psi_{1}$, and the ones coming from $\psi_{3}$ (superposition of connected
three-spin clusters, total domain length of 8), originate from the same
physical process, ie. adding or removing a spin from an existing cluster. Both
add two domain walls, and being symmetric in this sense, they cancel out. The
processes that couple $\psi_{2}$ to $\psi_{0}$ only appear at fourth order,
and we expect that accounting for those would reproduce the curvature of the
right and left peaks in Fig. 2(b).

The spectrum of excitations hints at the underlying reasons for quick
thermalization beyond $g/J\sim 2$. The $\Delta_{12}$ peak indicates the cost
of flipping a spin next to an already existing spin flip – growing a spin
cluster by one. Quenching to the point where this process is energetically as
expensive as creating a single spin flip populates not only the first few but
also higher excited states, resulting in the formation of larger clusters and
thus thermalization on a timescale of $1/J$.

Because the dynamics of confined systems is dominated by low energy
excitations, it can be accurately described by a simple effective model [16,
14, 24]. This is typically done by projecting the Hamiltonian to a subspace of
low-energy excitations (equivalents of our $\psi_{0}$, $\psi_{1}$ and
$\psi_{2}$). However, this approach leads to certain inconsistencies. Namely,
the matrix element between the ground and first excited state is extensive
(indeed, in our case we find
$\left\langle\psi_{0}\right|H\left|\psi_{1}\right\rangle=gN$), while others
are independent of $N$. Therefore, an effective model including $\psi_{0}$
predicts a spectral gap (equivalent of $\Delta_{01})$ that diverges with
system size. To accurately reproduce the excitation energies, the initial
vacuum state (equivalent of $\psi_{0}$) should be excluded from the low energy
model, and its mass set exactly to zero.

Perturbative calculations imply that the origin of the problem is in the
arbitrary truncation of the subspace of the effective model (see SM for
details). We find that the extensive corrections come from a state coupling to
an extensive number of excitations with one more and one fewer spin flip – it
is always possible to flip a decoupled spin in approximately $N^{2}$ places
for an energy penalty of $8J$. The contributions coming from the state with
one more and one less spin flip cancel out in the perturbative calculations,
and we only find an overall $g^{2}N^{2}/8J$ correction to all states, so that
energy differences remain independent of system size. However, if the Hilbert
space is first truncated into an effective model, the coupling between the
last retained and the first truncated state is not taken into account, the
extensive contributions do not cancel out, and the gaps grow with system size.
While the two approaches produce a similar energy diagram, an advantage of the
perturbative approach is that it can also provide corrections to the
eigenstates.

Interfaces \- Next, we probe the dynamics of interfaces by investigating
quenches of inhomogeneous initial states. We consider two examples: the
stripe, an initial state where two domains of opposite polarized spins are
separated by a flat interface, and the square, where the flipped spins are
arranged into a square domain.

Confinement importantly affects interface dynamics. Whereas interfaces are
typically expected to diffusively melt, in the confining phase they tend to
remain stable up to very long times. This is because classically resonant
processes dominate the dynamics for small $g/J$, ie. only spins with two
neighbours of each spin species can flip, as the initial and final state have
the same energy. These processes only locally reshape the domain wall, and
thus its total length emerges as a dynamic constraint. The interface dynamics
is well understood in the limit of infinite $J$ where non-resonant processes
are explicitly forbidden [25, 26]. Slow thermalization comes as a consequence
of the subsequent fragmentation of Hilbert space into sectors labeled by the
domain wall length, and further into so-called dynamically decoupled Krylov
subspaces; two domains cannot merge through local reshaping of the interface
if they are separated by multiple lattice sites [37, 38].

Results for the stripe are presented in Fig. 3, with a sketch of the initial
state shown in Fig. 3(a). Recall the periodic boundary conditions, which
ensure that the interface is completely flat with no kinks.

Figure 3:  Spread of correlations near an interface. (a) Sketch of the initial
state, with blue representing $\uparrow$ and red $\downarrow$ spins. Black
lines indicate the cuts shown in (d). (b) Sketch of the resonant process along
the interface. (c) Heatmap of the spectral density of
$\left\langle\sigma^{z}\right\rangle(t)$ for a spin at the interface for a
range of $g/J$. White dashed lines correspond to the transition energies of a
freely propagating edge mode. (d) Horizontal ($x$) and vertical ($y$) cuts of
the connected correlations function $C_{ij}$ with respect to the spin at
$(4,4)$ (white square in (a)). Blue dashed lines are $\pm 4gt$, showing that
the interface mode carries the correlations.

In this case, the lowest excitation above the initial state is a spin flip at
the interface. It creates two additional domain walls and costs $2\times 2J$.
However, once this excitation is present, flipping the neighbouring interface
spin does not come with any additional energy penalty, see Fig. 3(b). Such
excitations are thus allowed to freely propagate along the boundary, with the
expected dispersion $E(k)=4J\left(1+g\cos(k)\right)$ [39].

The spectral density of $\left\langle\sigma^{z}\right\rangle$ for a spin at
the interface is shown in Fig. 3(b). The most prominent feature is the
excitation which starts at $\omega=4J$ for $g=0$ and fans out into a continuum
with increasing $g$, as expected of a freely propagating mode. White dashed
lines correspond to the expected envelopes of the transition energies
$\omega=4J\pm 4g$. The continuum is split into discrete peaks as a consequence
of the finite size of the lattice.

Vertical and horizontal cuts of the connected correlation functions are shown
in Fig. 3(c). Because of the interface, the spread of correlations is not the
same in the horizontal and vertical directions. In the horizontal direction
(orthogonal to the interface) correlations show signatures of confinement
equivalent to results in Fig. 1. Interestingly, the correlation spread is
symmetric, meaning that the confined excitations spread across the interface
into the domain with opposite polarisation in the same way as within the
domain. On the other hand, the vertical (y) cuts show very different
behaviour. The freely propagating interface mode spreads the correlations
along the interface with a velocity proportional to $g$ (blue dashed lines in
Fig. 3(c)), even for $g\ll J$. As $g$ is further increased the dynamics lose
their $x$-$y$ asymmetric nature, and at large $g$ correlations spread equally
in either direction.

Finally, we consider the melting of an $8\times 8$ square of $\uparrow$ spins
embedded in a $16\times 16$ lattice and surrounded by $\downarrow$ spins. Fig.
4(a) shows snapshots of magnetization for small and intermediate $g/J$, while
Fig. 4(b) shows the time dependence of magnetization for a set of spins close
to a corner of the square. For extended plots of magnetization as well as
correlation spread, see the SM.

Figure 4:  Melting of a square. (a) Snapshots of magnetization for small and
intermediate $g/J$. (b,c) Magnetization $\left\langle\sigma^{z}\right\rangle$
for a set of spins indicated by the (b) white and (c) orange squares in the
initial state in (a). The panels with thick dashed frames show the same corner
spin.

As shown in Refs. [25, 26], in the $J\rightarrow\infty$ limit and for an
infinitely large corner the melting process can be understood as a series of
resonant spin-flip processes. Because the corner spin has two neighbours of
each spin species it flips on a timescale set by $1/g$. After the corner spin
is flipped, its two neighbours can resonantly flip as well. This triggers a
propagating edge mode in the same way as in Fig. 3. The process is illustrated
in the magnetization snapshots for $g/J=0.1$ in Fig. 4(a), and results in the
square melting into a cross-like shape.

A more quantitative view is presented in Figs. 4(b,c). At small $g/J$ the edge
mode dominates. For spins along the interface shown in panel (c), the delay of
melting is thus linear in the distance from the corner. Probing deeper into
the square along the diagonal (Fig. 4(b)), we find that spins near the center
stay inert up to relatively long times. Note that here we rescale time with
$g$, and find that the small $g/J$ curves in Fig. 4(b,c) follow the same
trend. This indicates that the smaller energy scale $g$ dominantly drives the
dynamics.

Increasing $g$ causes deviations from this simple picture and a crossover
towards a regime where all spins at the interface melt simultaneously. The
amplitude of magnetization oscillations is gradually reduced, and at $g/J=0.9$
we find a monotonous decrease of $\left\langle\sigma^{z}\right\rangle$ for all
spins of the square simultaneously. The melting produces waves of
magnetization radiating from sides of the square into the environment (see
bottom panel of Fig. 4(c), and SM extended plots). This is a sign that the
dynamics is no longer confined, and that the length of the domain wall is not
conserved.

Realization on a quantum simulator \- The Ising model is naturally realized in
various quantum simulation platforms. However, being two-dimensional, the
presented setup is most directly applicable to neutral atom arrays and within
reach of existing experiments. More concretely, Ref. [40] demonstrates control
over 2D arrays of hundreds of Rydberg atoms realizing a nearest-neighbour
antiferromagnetic Ising model with the interaction on the order of $2\,$MHz
and coherence time of $20\mu$s. When driving the transition between the two
spin states at around $1\,$Mhz, the timescales shown in Fig. 4 correspond to
experimental times of a few microseconds. If one chose to relax the nearest
neighbour constraint the atoms could be moved much closer together, leading to
stronger interactions and thus much faster dynamics. See for example Ref. [41]
for an experimental and Ref. [42] for a similar theoretical proposal.

Outlook \- We used tree tensor networks to study constrained dynamics in the
ordered phase of the two-dimensional quantum Ising model. We showed the
signatures of confinement of quasiparticles and extracted the spectrum of
excitations from oscillations of the magnetization. The spectrum can be
adequately reproduced by perturbation theory up to $g\sim J$, while the
failure at larger $g$ is related to thermalization. Then, we investigated the
melting of interfaces of domains with opposite spins, how it depends on
confinement and how it changes with increasing transverse field. We find a
crossover between the constrained dynamics of classically resonant processes
in the limit of small $g$ and diffusive melting at large $g$.

An important implication of the presented work is the fact that these
classical simulations are possible up to such long times. The underlying
reason is slow entanglement growth and thus slow increase of bond dimension
required to retain the accuracy of tensor network states. The confining regime
is thus perfect for benchmarking quantum machines with tensor network
calculations [24].

We see various other open questions that could be addressed either on a
quantum simulator, or with minimal alterations of our setup. One interesting
avenue is the exploration of Floquet dynamics through quenches with a time-
dependent transverse field. Because relatively long timescales are reachable,
it is possible to investigate the regime of small and intermediate
frequencies, and thus explore the regime beyond the $1/\omega$ corrections of
the Magnus expansion. Floquet dynamics have recently been studied in a related
system, where a discrete time crystal was identified in the prethermal state
[43].

Another option is to investigate false vacuum decay. In two dimensions,
dynamical confinement and the presence of false vacuum are decoupled –
confinement is provided by the lattice, while the false vacuum is invoked by a
longitudinal field. A starting point would be to investigate how the melting
in Fig. 4 depends on longitudinal field, which makes the square either a false
or a true vacuum domain.

Expanding on this idea, it would be interesting to probe the scattering of
false vacuum bubbles. Such calculations have been performed in one dimension
[44], but many more scattering channels are expected to arise in two
dimensions. However, to clearly discern scattering from lattice effects one
would have to significantly increase the system size, possibly beyond the
reach of classical simulations.

###### Acknowledgements.

We acknowledge enlightening discussions with Marcello Dalmonte, Roberto
Verdel, Simone Notarnicola, Jernej Rudi Finžgar, and Sourav Nandy. The authors
acknowledge funding from CN00000013 - Italian Research Center on HPC, Big Data
and Quantum Computing (ICSC), the H2020 projects EuRyQa, QuantERA2017 project
QuantHEP, QuantERA2021 project T-NiSQ, the Quantum Technology Flagship project
PASQuanS2, the Italian Ministry of University and Research (MUR) via PRIN2022
project TANQU, the Departments of Excellence grant 2023-2027 Quantum
Frontiers, the German Federal Ministry of Education and Research (BMBF) the
funding program quantum technologies, project QRydDemo, and the World Class
Research Infrastructure - Quantum Computing and Simulation Center (QCSC) of
Padova University. The simulations were performed using the Quantum Green Tea
software version 0.3.23 and Quantum Tea Leaves version 0.4.46 [45]. The
authors acknowledge computational resources of the HPC center Vega at the
Institute of Information Science (IZUM) in Maribor, Slovenia and HPC Leonardo
by Cineca, Italy.

## References

  * Polkovnikov  _et al._ [2011] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Colloquium: Nonequilibrium dynamics of closed interacting quantum systems, Reviews of Modern Physics 83, 863–883 (2011).
  * De Roeck and Verreet [2019] W. De Roeck and V. Verreet, Very slow heating for weakly driven quantum many-body systems (2019).
  * Yang [1989] C. N. Yang, $\eta$ pairing and off-diagonal long-range order in a hubbard model, Physical Review Letters 63, 2144–2147 (1989).
  * Moudgalya  _et al._ [2020] S. Moudgalya, N. Regnault, and B. A. Bernevig, $\eta$-pairing in hubbard models: From spectrum generating algebras to quantum many-body scars, Physical Review B 102, 10.1103/physrevb.102.085140 (2020).
  * Coleman [1977] S. Coleman, Fate of the false vacuum: Semiclassical theory, Physical Review D 15, 2929–2936 (1977).
  * Greensite [2020] J. Greensite,  _An Introduction to the Confinement Problem_ (Springer International Publishing, 2020).
  * James  _et al._ [2019] A. J. A. James, R. M. Konik, and N. J. Robinson, Nonthermal states arising from confinement in one and two dimensions, Phys. Rev. Lett. 122, 130603 (2019).
  * McCoy and Wu [1978] B. M. McCoy and T. T. Wu, Two-dimensional Ising field theory in a magnetic field: Breakup of the cut in the two-point function, Physical Review D 18, 1259–1267 (1978).
  * Kormos  _et al._ [2016] M. Kormos, M. Collura, G. Takács, and P. Calabrese, Real-time confinement following a quantum quench to a non-integrable model, Nature Physics 13, 246–249 (2016).
  * Mazza  _et al._ [2019] P. P. Mazza, G. Perfetto, A. Lerose, M. Collura, and A. Gambassi, Suppression of transport in nondisordered quantum spin chains due to confined excitations, Physical Review B 99, 10.1103/physrevb.99.180302 (2019).
  * Vovrosh and Knolle [2021] J. Vovrosh and J. Knolle, Confinement and entanglement dynamics on a digital quantum computer, Scientific Reports 11, 10.1038/s41598-021-90849-5 (2021).
  * 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, Domain-wall confinement and dynamics in a quantum simulator, Nature Physics 17, 742–747 (2021).
  * Magoni  _et al._ [2021] M. Magoni, P. P. Mazza, and I. Lesanovsky, Emergent bloch oscillations in a kinetically constrained Rydberg spin lattice, Physical Review Letters 126, 10.1103/physrevlett.126.103002 (2021).
  * Liu  _et al._ [2019] F. Liu, R. Lundgren, P. Titum, G. Pagano, J. Zhang, C. Monroe, and A. V. Gorshkov, Confined quasiparticle dynamics in long-range interacting quantum spin chains, Physical Review Letters 122, 10.1103/physrevlett.122.150601 (2019).
  * Pichler  _et al._ [2016] T. Pichler, M. Dalmonte, E. Rico, P. Zoller, and S. Montangero, Real-time dynamics in U(1) lattice gauge theories with tensor networks, Physical Review X 6, 10.1103/physrevx.6.011023 (2016).
  * Verdel  _et al._ [2020] R. Verdel, F. Liu, S. Whitsitt, A. V. Gorshkov, and M. Heyl, Real-time dynamics of string breaking in quantum spin chains, Physical Review B 102, 10.1103/physrevb.102.014308 (2020).
  * Lerose  _et al._ [2020] A. Lerose, F. M. Surace, P. P. Mazza, G. Perfetto, M. Collura, and A. Gambassi, Quasilocalized dynamics from confinement of quantum excitations, Physical Review B 102, 10.1103/physrevb.102.041118 (2020).
  * Verdel  _et al._ [2023] R. Verdel, G.-Y. Zhu, and M. Heyl, Dynamical localization transition of string breaking in quantum spin chains, Physical Review Letters 131, 10.1103/physrevlett.131.230402 (2023).
  * Lagnese  _et al._ [2021] G. Lagnese, F. M. Surace, M. Kormos, and P. Calabrese, False vacuum decay in quantum spin chains, Physical Review B 104, 10.1103/physrevb.104.l201106 (2021).
  * Pomponio  _et al._ [2022] O. Pomponio, M. A. Werner, G. Zaránd, and G. Takacs, Bloch oscillations and the lack of the decay of the false vacuum in a one-dimensional quantum spin chain, SciPost Physics 12, 10.21468/scipostphys.12.2.061 (2022).
  * Ramos  _et al._ [2020] F. B. Ramos, M. Lencsés, J. C. Xavier, and R. G. Pereira, Confinement and bound states of bound states in a transverse-field two-leg Ising ladder, Physical Review B 102, 10.1103/physrevb.102.014426 (2020).
  * Lagnese  _et al._ [2022] G. Lagnese, F. M. Surace, M. Kormos, and P. Calabrese, Quenches and confinement in a Heisenberg–Ising spin ladder, Journal of Physics A: Mathematical and Theoretical 55, 124003 (2022).
  * Konik and Adamov [2009] R. M. Konik and Y. Adamov, Renormalization group for treating 2D coupled arrays of continuum 1D systems, Physical Review Letters 102, 10.1103/physrevlett.102.097203 (2009).
  * Tindall and Sels [2024] J. Tindall and D. Sels, Confinement in the transverse field Ising model on the heavy hex lattice (2024).
  * Balducci  _et al._ [2022] F. Balducci, A. Gambassi, A. Lerose, A. Scardicchio, and C. Vanoni, Localization and melting of interfaces in the two-dimensional quantum Ising model, Physical Review Letters 129, 10.1103/physrevlett.129.120601 (2022).
  * Balducci  _et al._ [2023] F. Balducci, A. Gambassi, A. Lerose, A. Scardicchio, and C. Vanoni, Interface dynamics in the two-dimensional quantum Ising model, Physical Review B 107, 10.1103/physrevb.107.024306 (2023).
  * Shi  _et al._ [2006] Y.-Y. Shi, L.-M. Duan, and G. Vidal, Classical simulation of quantum many-body systems with a tree tensor network, Physical Review A 74, 10.1103/physreva.74.022320 (2006).
  * Silvi  _et al._ [2019] P. Silvi, F. Tschirsich, M. Gerster, J. Jünemann, D. Jaschke, M. Rizzi, and S. Montangero, The Tensor Networks Anthology: Simulation techniques for many-body quantum lattice systems, SciPost Physics Lecture Notes 10.21468/scipostphyslectnotes.8 (2019).
  * Cataldi  _et al._ [2021] G. Cataldi, A. Abedi, G. Magnifico, S. Notarnicola, N. D. Pozza, V. Giovannetti, and S. Montangero, Hilbert curve vs Hilbert space: exploiting fractal 2D covering to increase tensor network efficiency, Quantum 5, 556 (2021).
  * Bauernfeind and Aichhorn [2020] D. Bauernfeind and M. Aichhorn, Time dependent variational principle for tree tensor networks, SciPost Physics 8, 10.21468/scipostphys.8.2.024 (2020).
  * Haegeman  _et al._ [2016] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Unifying time evolution and optimization with matrix product states, Physical Review B 94, 10.1103/physrevb.94.165116 (2016).
  * Eisert and Osborne [2006] J. Eisert and T. J. Osborne, General entanglement scaling laws from time evolution, Physical Review Letters 97, 10.1103/physrevlett.97.150404 (2006).
  * Hashizume  _et al._ [2022] T. Hashizume, I. P. McCulloch, and J. C. Halimeh, Dynamical phase transitions in the two-dimensional transverse-field Ising model, Physical Review Research 4, 10.1103/physrevresearch.4.013250 (2022).
  * Lin and Motrunich [2017] C.-J. Lin and O. I. Motrunich, Quasiparticle explanation of the weak-thermalization regime under quench in a nonintegrable quantum spin chain, Physical Review A 95, 10.1103/physreva.95.023621 (2017).
  * win [2024] Window function — Wikipedia, the free encyclopedia (2024), [Online; accessed 28-May-2024].
  * Harris [1978] F. Harris, On the use of windows for harmonic analysis with the discrete fourier transform, Proceedings of the IEEE 66, 51 (1978).
  * Yoshinaga  _et al._ [2022] A. Yoshinaga, H. Hakoshima, T. Imoto, Y. Matsuzaki, and R. Hamazaki, Emergence of hilbert space fragmentation in Ising models with a weak transverse field, Physical Review Letters 129, 10.1103/physrevlett.129.090602 (2022).
  * Hart and Nandkishore [2022] O. Hart and R. Nandkishore, Hilbert space shattering and dynamical freezing in the quantum Ising model, Physical Review B 106, 10.1103/physrevb.106.214426 (2022).
  * Mbeng  _et al._ [2020] G. B. Mbeng, A. Russomanno, and G. E. Santoro, The quantum Ising chain for beginners (2020).
  * Scholl  _et al._ [2021] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. Läuchli, and A. Browaeys, Quantum simulation of 2D antiferromagnets with hundreds of Rydberg atoms, Nature 595, 233–238 (2021).
  * Bluvstein  _et al._ [2021] D. Bluvstein, A. Omran, H. Levine, A. Keesling, G. Semeghini, S. Ebadi, T. T. Wang, A. A. Michailidis, N. Maskara, W. W. Ho, S. Choi, M. Serbyn, M. Greiner, V. Vuletić, and M. D. Lukin, Controlling quantum many-body dynamics in driven Rydberg atom arrays, Science 371, 1355–1359 (2021).
  * Samajdar and Huse [2024] R. Samajdar and D. A. Huse, Quantum and classical coarsening and their interplay with the Kibble-Zurek mechanism (2024).
  * Shinjo  _et al._ [2024] K. Shinjo, K. Seki, T. Shirakawa, R.-Y. Sun, and S. Yunoki, Unveiling clean two-dimensional discrete time quasicrystals on a digital quantum computer (2024).
  * Milsted  _et al._ [2022] A. Milsted, J. Liu, J. Preskill, and G. Vidal, Collisions of false-vacuum bubble walls in a quantum spin chain, PRX Quantum 3, 10.1103/prxquantum.3.020316 (2022).
  * Ballarin  _et al._ [2024] M. Ballarin, G. Cataldi, A. Costantini, D. Jaschke, G. Magnifico, S. Montangero, S. Notarnicola, A. Pagano, L. Pavešić, M. Rigobello, N. Reinić, S. Scarlatella, and P. Silvi, Quantum TEA: qtealeaves (2024).

## Supplementary Material

### Numerical details

The numerical calculations presented in this work are performed with tree
tensor networks (TTNs) [27, 28]. In this appendix, we provide qualitative
reasoning for the utility of TTNs compared to other tensor networks and
technical aspects of the algorithm for time evolution that we apply, and
present plots of magnetization and entanglement obtained for a set of maximal
bond dimensions.

#### Tree tensor networks

In TTNs, the state is represented by a tensor network in the shape of a
symmetric binary tree, where the bottom tensors (’leaves’ of the tree) carry
the physical indices. The tree tensor network ansatz is particularly
successful at describing states where distant sites are entangled, such as in
chains with long-range interactions. The reason why can be understood through
a qualitative picture of the way quantum states are encoded in tensor
networks. Qualitatively, entanglement between two correlated physical states
is transferred through bonds of the network that connect them. The more
entanglement there is, the bigger the bond dimension required to accurately
describe the state. Generally, this also means entangled states are better
represented by tensor networks with a bigger number of links. An extreme
example of this is the projected entangled-pair state (PEPS) ansatz, where the
tensor network mimics the structure of the 2D lattice and each tensor is
directly connected to its nearest neighbours. However, this introduces loops
in the network, which leads to technical difficulties. Namely, it is not
possible to define the isometry center of a network with loops. Therefore, the
algorithms for contracting such networks cannot rely on the underlying
symmetries. This leads to poor scaling of calculations of local observables
with bond dimension, and limits PEPS to bond dimensions on the order of
$\chi\sim 10$. The TTN ansatz is in between the MPS and PEPS in the sense that
it is the best connected tensor network which by construction has no loops.
Furthermore, because two physical sites separated by $m$ lattice sites in real
space are only separated by $\sim\ln m$ in the TTN, the entanglement between
them is captured more efficiently.

#### Time evolution

The most common way to evolve tensor network states in time is with the time-
dependent variatonal principle (TDVP) [30, 31]. The most popular is the two-
tensor version of the algorithm (TDVP2), where at each step two neighbouring
tensors are contracted, evolved in time, and then split and truncated via SVD.
This algorithm can grow the bond dimensions between tensors at each time step
and is thus a natural choice for simulations starting from weakly entangled
states, such as the product states we use. A much faster version of the same
concept is the single-tensor TDVP (TDVP1), where each tensor is optimized
without contractions with the neighbours. In this case, the bond dimension
cannot change, which presents a problem if the initial state is a product
state with $\chi=1$.

However, we find that the calculations are much faster if one pads the initial
product state with numerical noise on the order of $10^{-16}$ up to bond
dimension $\chi$, and we perform the evolution with TDVP1, without truncation.
We find a speedup in computational time of about $30\times$ for a $8\times 8$
system at $\chi=100$. We also tried hybrid approaches (like performing TDVP2
every $n-$th time step and TDVP1 otherwise, with $n$ around 10). The advantage
here is that the initial bond dimension is not $\chi$ but rather gradually
increases every $n$-th step. This speeds up the initial time steps, but the
advantage is quickly lost once the bond dimension reaches a sizeable fraction
of the one used for TDVP1. An advantage might be expected for the easiest
cases (small $g/J$ limit) where the bond dimension never reaches the maximal
allowed value. However, these cases can also be accurately simulated by a
TDVP1 simulation with smaller $\chi$, which is anyway much faster.
Furthermore, on a technical level having a single convergence parameter
($\chi$) is more desirable compared to having two ($\chi$ and the minimal
truncation $\epsilon$) with non-trivial interplay.

Another advantage of TDVP1 compared to TDVP2 is that the error of the latter
is non-monotonous in the size of the time step. In TDVP2, an SVD and a
truncation are performed at each time step. This introduces a truncation
error, and thus choosing a time step that is too small might actually increase
the total error. In TDVP1, there is no truncation, and the truncation error is
exactly zero. Consequently, TDVP1 also exactly conserves the energy and norm
of the evolved state.

#### Convergence

In Fig. 5, we show the convergence of (a) magnetization and (b) entanglement
with increasing maximal bond dimension $\chi$. The deviation of the
magnetization at given $\chi$ with respect with the results obtained with
$\chi=160$ are shown in Fig. 5(c).

The convergence is strongly dependent on the underlying physics. Deep in the
confined phase at $g\ll J$, results are well converged already for $\chi=60$
or less. At larger $g$, larger $\chi$ is required to obtain well-converged
results and the simulations become unreliable much sooner in simulated time.

We use $\chi=60$ for simulations of longer times, and find the magnetization
results reliable up to $g/J\approx 1.5$. In simulations where we also want to
ensure accurate convergence in entanglement and correlations, we use
$\chi=100$.

Figure 5:  Convergence with bond dimension $\chi$. (a) Magnetization
$\sigma^{z}$ of a single spin calculated with increasing bond dimension $\chi$
for a set of $g/J$. (b) Entanglement entropy calculated with increasing bond
dimension $\chi$ for a set of $g/J$. (c) Log scale plots of the deviation of
magnetization from the results obtained with $\chi=160$,
$\delta\left\langle\sigma^{z}\right\rangle=\left|\left\langle\sigma^{z}\right\rangle(\chi)-\left\langle\sigma^{z}\right\rangle(\chi=160)\right|$.

### Perturbative corrections

This appendix contains calculations of the perturbative corrections up to the
second order in $g$ for the first two excited states above the completely
polarized ground state. We treat the interaction part as the unperturbed part

$H_{0}=-J\sum_{\langle i,j\rangle}\sigma^{z}_{i}\sigma^{z}_{j},$ (3)

and the transverse field term as the perturbation

$V=-g\sum_{i}\sigma^{x}_{i},$ (4)

with $i$ running across the $N\times N$ spins of the 2D system.

#### General remarks

The first-order perturbative corrections to energy of an unperturbed state
$\psi$ are given by:

$\delta E^{(1)}=-\left\langle\psi\right|V\left|\psi\right\rangle.$ (5)

In our case, the unperturbed states are product states of spins polarized in
$z$, while $V$ only contains spin flipping terms $\sigma^{x}$. To obtain a
non-zero overlap with the initial state, one has to flip at least two spins
(or one spin twice), thus the first-order contributions will always be zero.
The same argument applies to all odd-order corrections.

In the second order, the correction to a state $\psi$ is

$\begin{split}\delta
E^{(2)}&=-\sum_{k\neq\psi}\left\langle\psi\right|V\left|k\right\rangle\frac{1}{E_{k}-E_{0}}\left\langle
k\right|V\left|\psi\right\rangle\\\
&=-\sum_{k\neq\psi}\frac{1}{E_{k}-E_{0}}\left|\left\langle\psi\right|V\left|k\right\rangle\right|^{2},\end{split}$
(6)

where $k$ runs across all unperturbed eigenstates of the system that are not
$\psi$, $E_{k}=\left\langle k\right|H_{0}\left|k\right\rangle$ being their
unperturbed energy and
$E_{0}=\left\langle\psi\right|H_{0}\left|\psi\right\rangle$ the unperturbed
energy of $\psi$. In other words, $\delta E^{(2)}$ is the sum of absolute
values of all possible off-diagonal matrix elements weighted by the difference
of the corresponding state’s energy with respect to $\psi$.

In the case of the Ising model all off-diagonal matrix elements are the same,
$g$. To evaluate $\delta E^{(2)}$, we find all states which are a single spin
flip away from $\psi$, group them according to their energy $E_{k}$, and count
the number of such contributions.

#### Ground state

We are looking for the second order correction:

$\begin{split}\delta
E_{0}^{(2)}&=-\sum_{k}\left\langle\psi_{0}\right|V\left|k\right\rangle\frac{1}{E_{k}-E_{0}}\left\langle
k\right|V\left|\psi_{0}\right\rangle\\\
&=-\sum_{k}\frac{1}{E_{k}-E_{0}}\left|\left\langle\psi_{0}\right|V\left|k\right\rangle\right|^{2}.\end{split}$
(7)

The matrix element $\left\langle\psi_{0}\right|V\left|k\right\rangle$ is non-
zero only for states with a single flipped spin,
$\left|k\right\rangle=\sigma^{x}_{i}\left|\psi_{0}\right\rangle$, with $i$
denoting the coordinates of the flip. Evidently, there are $N^{2}$ such
states, each with the matrix element
$\left\langle\psi_{0}\right|V\left|k\right\rangle=g$ and excitation energy
$E_{k}-E_{0}=4\times 2J$ coming from the four domain walls.

The sum in Eq. (7) thus gives

$\delta E_{0}^{(2)}=-\sum_{i}^{N}\frac{g^{2}}{8J}=-\frac{g^{2}}{8J}N^{2}.$ (8)

#### First excited state

To find the $g$-dependent correction to the excited states, we repeat the
procedure for the state with one flipped spin. By examining the matrix
elements
$\left\langle\psi_{0}\right|H\sigma^{x}_{i}\left|\psi_{0}\right\rangle$, one
notices that the only excited state coupled to $\psi_{0}$ is the zero-momentum
superposition of spin flips

$\left|\psi_{1}\right\rangle=\frac{1}{\sqrt{N^{2}}}\sum_{i}\sigma^{x}_{i}\left|\psi_{0}\right\rangle.$
(9)

It is thus enough to evaluate the perturbative expansion only for this state.
Again the first order correction $\delta E_{1}^{(1)}=0$, while in second order
we have

$\delta
E_{1}^{(2)}=\sum_{k}\frac{1}{E_{k}-E_{1}}\left|\left\langle\psi_{1}\right|V\left|k\right\rangle\right|^{2},$
(10)

with $E_{1}=4\times 2J$.

In this case, the non-zero contributions come from states with two spin flips:

$\left|k\right\rangle=\sigma^{x}_{i}\sigma^{x}_{j}\left|\psi_{0}\right\rangle.$
(11)

Depending on the relative positions of $i$ and $j$, there are three distinct
classes of $\left|k\right\rangle$:

  1. 1.

$i=j$;
$\sigma^{x}_{i}\sigma^{x}_{i}\left|\psi_{0}\right\rangle=\left|\psi_{0}\right\rangle$.
There is exactly one such state, with intermediate energy
$E_{k}-E_{1}=-4\times 2J$.

  2. 2.

$i,j$ neighbours, denoted
$\left|\psi_{2}^{\prime}\right\rangle_{ij}=\sigma^{x}_{i}\sigma^{x}_{j}\left|\psi_{0}\right\rangle$.
There are $\frac{1}{2}4N^{2}$ such states (4 for each lattice site, but
$\frac{1}{2}$ for double counting). Flipping neighbouring spins only creates 6
domain walls, thus the intermediate energy of this type of process is $2\times
2J$.

  3. 3.

$i,j$ neither neighbours nor equal, denoted
$\left|\psi_{2}^{\prime\prime}\right\rangle_{ij}$. There are
$\frac{1}{2}N^{2}(N^{2}-5)$ such states, $\frac{1}{2}$ again for double
counting. The intermediate energy is $4\times 2J$.

Now, we compute the matrix elements for the three classes.

For $\psi_{0}$, we have

$\begin{split}\left\langle\psi_{0}\right|V\left|\psi_{1}\right\rangle&=\left\langle\psi_{0}\right|-g\sum_{i}\sigma^{x}_{i}\frac{1}{\sqrt{N^{2}}}\sum_{j}\sigma^{x}_{j}\left|\psi_{0}\right\rangle\\\
&=-\frac{g}{\sqrt{N^{2}}}\sum_{ij}\left\langle\psi_{0}\right|\sigma^{x}_{i}\sigma^{x}_{j}\left|\psi_{0}\right\rangle\\\
&=-\frac{g}{\sqrt{N^{2}}}\sum_{ij}\delta_{ij}\\\
&=-\frac{g}{\sqrt{N^{2}}}N^{2},\end{split}$ (12)

and the total contribution to $\delta E_{1}^{(2)}$ is

$\frac{1}{E_{0}-E_{1}}\left|\left\langle\psi_{0}\right|V\left|\psi_{1}\right\rangle\right|^{2}=-\frac{1}{8J}\times
g^{2}N^{2}.$ (13)

For the basis states with neighbouring flipped spins $k$ and $l$,
$\left|\psi_{2}^{\prime}\right\rangle_{kl}$, we find

$\begin{split}\left\langle\psi_{2}^{\prime}\right|_{kl}V\left|\psi_{1}\right\rangle&=-\frac{g}{\sqrt{N^{2}}}\left\langle\psi_{2}^{\prime}\right|_{kl}\sum_{ij}\sigma^{x}_{i}\sigma^{x}_{j}\left|\psi_{0}\right\rangle\\\
&=-\frac{g}{\sqrt{N^{2}}}\times 2,\end{split}$ (14)

as the sum over $i$ and $j$ gives two ways to match the two flipped spins in
$\psi_{2}^{\prime}$. The total contribution of all such states to the energy
shift is thus

$2N^{2}\frac{1}{E_{2}^{\prime}-E_{1}}\left|\left\langle\psi_{2}^{\prime}\right|_{kl}V\left|\psi_{1}\right\rangle\right|^{2}=2N^{2}\times\frac{4g^{2}}{N^{2}}\times\frac{1}{4J}=\frac{2g^{2}}{J}.$
(15)

Finally, for $\psi_{2}^{\prime\prime}$ the matrix element is the same as
above:

$\left\langle\psi_{2}^{\prime\prime}\right|_{kl}V\left|\psi_{1}\right\rangle=-\frac{g}{\sqrt{N^{2}}}\times
2,$ (16)

however, the intermediate energy is $8J$, and the number of states is bigger
by a factor of $N^{2}-5$, producing another extensive contribution to the
energy shift:

$\frac{1}{2}N^{2}(N^{2}-5)\times\frac{4g^{2}}{N^{2}}\times\frac{1}{8J}=2\frac{g^{2}}{8J}\left(N^{2}-5\right).$
(17)

The complete second-order energy correction is the sum of these terms:

$\begin{split}\delta
E_{1}^{(2)}&=-\left(-\frac{g^{2}}{8J}N^{2}+\frac{2g^{2}}{J}+\frac{g^{2}}{8J}2\left(N^{2}-5\right)\right)\\\
&=-\frac{g^{2}}{8J}\left(N^{2}+6\right).\end{split}$ (18)

Importantly, the extensive ($\propto N^{2}$) contribution is the same as for
the ground state, see Eq. (8), resulting in a $N$-independent gap.

#### Second excited state

The second excited state is of type $\psi_{2}^{\prime}$, with two neighbouring
spins flipped. Again, it is enough to consider the zero-momentum superposition

$\left|\psi_{2}^{\prime}\right\rangle=\frac{1}{2N^{2}}\sum_{ij}\sigma^{x}_{i}\sigma^{x}_{j}\left|\psi_{0}\right\rangle,$
(19)

with $i$, $j$ neighbours.

The relevant classes of states coupled to
$\left|\psi_{2}^{\prime}\right\rangle$ have three spin flips over the ground
state,
$\sigma^{x}_{i}\sigma^{x}_{j}\sigma^{x}_{k}\left|\psi_{0}\right\rangle$. They
are:

  1. 1.

two of $i$, $j$, $k$ equal, producing a state with a single flipped spin,
$\left|\psi_{1}\right\rangle$. There is $N^{2}$ such states, and the
intermediate energy is $-2\times 2J$.

  2. 2.

$i$, $j$, and $k$ forming a connected three spin cluster, the state denoted
$\left|\psi_{3}^{\prime}\right\rangle=\sigma^{x}_{i}\sigma^{x}_{j}\sigma^{x}_{k}\left|\psi_{0}\right\rangle$.
There are $6N^{2}$ such clusters on a $N\times N$ lattice. The intermediate
energy comes from adding two domain walls, and is thus $2\times 2J$.

  3. 3.

$k$ disconnected from $i$ and $j$, denoted
$\left|\psi_{3}^{\prime\prime}\right\rangle$. There are $2N^{2}(N^{2}-8)$ such
states – for each $i$, $j$ cluster the disconnected spin can be flipped on the
remaining $N^{2}-8$ sites.

For the first case, we find the matrix element

$\left\langle\psi_{1}\right|V\left|\psi_{2}^{\prime}\right\rangle=\left\langle\psi_{0}\right|\sigma^{x}_{i}\sum_{j}g\sigma^{x}_{j}\frac{1}{\sqrt{2N^{2}}}\sum_{kl}\sigma^{x}_{k}\sigma^{x}_{l}\left|\psi_{0}\right\rangle=g\frac{4}{\sqrt{2N^{2}}}.$
(20)

In words, there are four ways to match any one flipped spin by starting from a
superposition of all possible two-spin clusters and flipping one of them. The
$\psi_{1}$ contribution to the correction is thus

$N^{2}\times\frac{-1}{4J}\times\frac{16g^{2}}{2N^{2}}.$ (21)

In the second case, the matrix element is

$\begin{split}\left\langle\psi_{3}^{\prime}\right|V\left|\psi_{2}^{\prime}\right\rangle=g\frac{2}{\sqrt{2N^{2}}}.\end{split}$
(22)

Each connected three-spin cluster can be constructed by flipping a
neighbouring spin of two distinct two-spin clusters. The total contribution is

$6N^{2}\times\frac{1}{4J}\times\frac{4g^{2}}{2N^{2}}.$ (23)

Finally, the case of a decoupled flipped spin is similar to the third case of
the $\psi_{1}$ corrections, but here with a different geometrical factor. The
matrix element is

$\left\langle\psi_{3}^{\prime\prime}\right|V\left|\psi_{2}^{\prime}\right\rangle=g\frac{1}{2N^{2}},$
(24)

and thus the total contribution

$N^{2}\left(N^{2}-8\right)\times\frac{1}{4J}\times\frac{g^{2}}{2N^{2}}.$ (25)

When summing all contributions all $N$-independent terms cancel out, and we
obtain the second order correction to the energy of $\psi_{2}^{\prime}$:

$\delta E_{2}^{(2)}=-g\frac{g^{2}}{8J}N^{2}.$ (26)

See the main text around Fig. 2 for a qualitative interpretation of the
obtained result.

### Additional results for the square

In this section, we provide additional results for the dynamics of the square.
Fig. 6 shows a grid of magnetization for the lower right quadrant of the
square (the panels with blue frames) and two surrounding rows and columns (red
frames). Time is rescaled by $g$, as in Fig. 4.

As discussed in the main text, we observe a crossover between two limiting
behaviours; at small $g/J$ the dynamics consists of a sequence of resonant
spin flips which originate from the corner, while with increasing $g/J$ the
spins at the interface melt simultaneously. Spins within the square remain
inert until the wave of melting reaches them. Interestingly, because the
resonant process is a complete spin flip which propagates ballistically, the
magnetization at intermediate times ($2\sim 6tg$) generally reaches lower
values for small $g/J$, and the square seems to start melting more quickly.
However, while the small-$g$ dynamics leads to multiple oscillations of
magnetization, the large $g/J$ curves monotonously decay to their thermal
value. Longer time simulations are required to accurately determine the nature
of the thermalized state.

The diffusive dynamics at large $g/J$ spreads outside of the initial square,
while it is confined to it at small $g/J$. This is most apparent in the set of
spins that neighbour the square, where we find a large initial jump of
magnetization (note that the red-framed panels have a $y$-scale from $-1$ to
$-0.5$). Its magnitude increases approximately linearly with $g/J$. For
smaller $g/J$ the environment is effectively unperturbed at the scale shown
here.

Figure 6:  Time evolution of magnetization for the lower right quarter of the
square (blue frames) and two rows of surrounding spins (red frames). The
y-range of the blue panels is from $-1$ to $1$, while for red panels it is
$-1$ to $-0.5$.

In Fig. 7, we show the correlations $C_{ij}$ where $i$ is the spin in the
corner of the square, and $j$ are taken along a vertical (top) or a diagonal
(bottom) cut along the system. For $g/J=0.1$ we find the expected linear
correlation spread along the side of the square, while it is somewhat slower
along the diagonal. Interestingly, because the edge mode propagates much
faster, the correlations with the spin in the opposite corner of the square
begin growing approximately at the same time as the ones with the neighbouring
spin along the diagonal.

As $g/J$ is increased, we find the expected light-cone spread, and the spread
of correlations outside of the initial spread. An interesting feature that
appears at larger $g/J$ are areas where $C_{ij}$ decreases (at $tg\sim 1$).
This is probably due to interference of signals propagating along different
paths from the corner spin.

Figure 7:  Correlations $C_{ij}$ with the $i$-th spin at the corner of the
square, at coordinates $(12,12)$ (counted from the upper left corner as
$(1,1)$).