# A quantum information perspective on meson melting

Mari Carmen Bañuls<EMAIL_ADDRESS>Max Planck Institute of Quantum
Optics, 85748 Garching bei München, Germany  Munich Center for Quantum Science
and Technology (MCQST), 80799 München, Germany  Michal P. Heller
<EMAIL_ADDRESS>Department of Physics and Astronomy, Ghent
University, 9000 Ghent, Belgium  Karl Jansen<EMAIL_ADDRESS>Deutsches
Elektronen-Synchrotron (DESY), Platanenallee 6, 15738 Zeuthen, Germany
Johannes Knaute<EMAIL_ADDRESS>Racah Institute of Physics,
The Hebrew University of Jerusalem, Jerusalem 91904, Givat Ram, Israel  Viktor
Svensson<EMAIL_ADDRESS>Division of Solid State Physics and
NanoLund, Lund University, S-221 00 Lund, Sweden

###### Abstract

We propose to use quantum information notions to characterize thermally
induced melting of nonperturbative bound states at high temperatures. We apply
tensor networks to investigate this idea in static and dynamical settings
within the Ising quantum field theory, where bound states are confined fermion
pairs – mesons. In equilibrium, we identify the transition from an exponential
decay to a power law scaling with temperature in an efficiently computable
second Rényi entropy of a thermal density matrix as a signature of meson
melting. Out of equilibrium, we identify as the relevant signature the
transition from an oscillatory to a linear growing behavior of reflected
entropy after a thermal quench. These analyses apply more broadly, which
brings new ways of describing in-medium meson phenomena in quantum many-body
and high-energy physics.

###### Contents

  1. I Introduction and motivation
  2. II Preliminaries
    1. II.1 Phenomenology of the meson melting process in QCD
    2. II.2 The 2$D$ Ising QFT
  3. III Methods and concepts
    1. III.1 Tensor network simulations
    2. III.2 Entropies
    3. III.3 Quantum quenches and quasiparticle model
  4. IV Results from thermal equilibrium
  5. V Results from dynamics after thermal quenches
    1. V.1 Setup
    2. V.2 Effective temperature effects
    3. V.3 Quench results in the high temperature scaling limit
  6. VI Discussion and outlook
  7. A Review of QCD and holographic approaches
    1. A.1 Theoretical methods in QCD
    2. A.2 Insights from holographic models
  8. B Entanglement structure and tensor network ansätze
  9. C Details of the transfer operator
  10. D Prony signal analysis

##  I Introduction and motivation

Emergent phenomena of quantum field theories (QFTs) under extreme conditions
can pose significant challenges for their theoretical description both in
condensed matter and particle physics [1, 2, 3, 4]. Many problems of this type
originate from quantum chromodynamics (QCD) and so is the main actor in the
present article: the melting of mesons [5, 6]. With quantum technologies
offering a significant long-term promise for studying fundamental physics
problems, we believe it is the right time to frame the problem of melting of
bound states using their native language of quantum information.

The phenomenon of meson melting, i.e. the breakdown of bound states into their
elementary constituents at high enough temperatures, is of fundamental
interest in QCD, where it is relevant for the physics of the early universe
and for the understanding of nuclear collisions. Whereas a full explanation of
the process from first principles is not available, the current
phenomenological description understands the meson melting as a dynamical
process in which out-of-equilibrium mesons interact with a thermal background
[5].

In this article, we propose a new approach to the meson melting phenomenon
through the study of entanglement measures. As a first step in this direction,
we consider the two-dimensional Ising QFT, which contains both integrable and
nonintegrable parameter regimes exhibiting meson excitations. In particular,
we study static and time-dependent thermal states in different temperature
regimes in the thermodynamic and continuum limit of the quantum Ising model
near its critical point by analyzing several entropic quantities. For this
purpose we employ numerical tensor network (TN) simulations, a tool
originating from quantum information that directly gives access to such
entanglement measures [7, 8, 9]. The dynamical setup is generated by a quantum
quench – a common method, which creates an out-of-equilibrium state through an
(instantaneous) parameter change in the underlying Hamiltonian. This protocol
can be seen as a theoretical method by itself, since it can induce entirely
new dynamical features to the system. We use it to mimic a dynamical situation
as in the modern QCD viewpoint on meson melting. Notice that, recently,
entanglement measures have been used in a similar vein to characterize
($1+1$)-dimensional meson scattering events [10].

Our TN approach is also closely related to the field of quantum simulations
[11, 12, 13, 14], in which quantum hardware is utilized to study many-body
problems via digital or analog implementations. Recently, there has been a
considerable amount of attention to the potential use of these technologies
for studies of meson physics, both in experimental measurements [15, 16, 17,
18] as well as in theoretical proposals [19, 20, 21, 22, 23], which provides
an additional motivation for our study.

This article is organized as follows. Section II is devoted to relevant
physical background information. In particular, in section II.1 we review
physical concepts for a phenomenological description of meson melting in QCD
that motivate our work. Section II.2 introduces the Ising QFT emerging in the
continuum limit of the quantum Ising model as the studied theory in this
article. In section III we describe the relevant concepts of entanglement
properties in quantum many-body systems, which underlie TN ansätze and
numerical algorithms for our simulations. In section IV we then present
simulation results for different scaling behaviors of the Rényi entropy
density in thermal equilibrium states. In section V we analyze entanglement
entropies after thermal quantum quenches at different effective temperatures.
In both cases, we consider mesonic parameter regimes of the Ising QFT and find
signatures that indicate the melting of meson states at high temperatures. The
appendices are devoted to further background information and computational
details.

##  II Preliminaries

###  II.1 Phenomenology of the meson melting process in QCD

Spectral functions are the main quantities studied within QFT to describe the
melting of mesons. In this section we provide an introduction to the known
phenomenology of this effect based on the recent and comprehensive review [5].
For more detailed discussions and further background information we refer to
appendix A.1.

In thermal equilibrium the spectral function $\rho$ is given by the negative
imaginary part of the retarded correlator,

$\rho(\vec{p},\omega)=-\operatorname{Im}[D_{R}(\vec{p},\omega)].$ (1)

Here, the retarded correlator $D_{R}$ is defined w.r.t. a meson operator
$M(\vec{x},t)$ as the thermal expectation value
$D_{R}(\vec{x},t,\vec{x}_{0},t_{0})\equiv\theta(t-t_{0})\operatorname{Tr}\\{\rho_{\beta}[M(\vec{x},t),M^{\dagger}(\vec{x}_{0},t_{0})]\\}$,
where $\rho_{\beta}$ is the thermal density operator at inverse temperature
$\beta$ and $\theta(t)$ is the Heaviside function. In a QCD context, $\rho$ is
typically considered in Fourier space via a Wigner transformation, in which
relative frequency and momentum coordinates $(\vec{p},\omega)$ are used as in
(1). The spectral function in this form probes the particle content of the
physical system. Meson bound states show up as peaks at frequency values given
by their masses.

Within QCD, there are several established methods to calculate spectral
functions of quark bound states. 111We refer to appendix A.1 for a short
overview and discussion of them. They form the basis for our current
understanding, which was initiated by the influential work [25] on suppression
and melting of mesons. Starting from an analogy with Debye screening in an
electromagnetic plasma, the authors argue that the QGP created in a nuclear
collision weakens the binding energies of meson bound states in the thermal
environment, causing a suppression of meson detection rates. Building up on
that, the work [26] put forward the picture of a sequential process. In this
static scenario, an in-medium Hamiltonian describes mesons as its eigenstates.
The potential term in the Hamiltonian gets weaker for increasing temperatures,
such that bound states transform into scattering states at a unique melting
temperature, where a sharp transition takes place. Weakly bound mesons states
melt first, while strongly bound ones survive longer in the thermal medium.

The work [27] opened a new way of thinking about this process by showing that
the quark potential (24) is complex. This implies that a meson bound state is
a dynamical system, which interacts constantly with its thermal environment by
scattering processes. In this modern picture, the melting process is best
understood from the behavior of spectral functions in dependence of the
temperature. With increasing temperature, the meson peaks move towards smaller
frequency values; their magnitude is decreasing while their thermal width is
increasing. Peaks at larger frequencies, i.e. lower binding energies, melt
first and their peak structure dissolves completely into a continuum. The
combination of thermal broadening and decreasing of the binding between the
quarks means that mesons transition into other (unknown) states. The dynamical
origin of this process is not known from this description. As a consequence,
the definition of a melting temperature in this time-dependent, sequential
scenario is not uniquely determinable any more and a less important concept.
The overall intuitive picture drawn in [5] is that the thermal environment
acts as a sieve that filters out weakly bound meson states, which are exposed
to intensifying scattering events as the temperature is raised. Apart from
direct approaches in QCD, also holographic models corroborated this
phenomenological picture of a sequential melting process. Some insights from
these methods are reviewed in appendix A.2.

From these elaborations, it becomes obvious that only a fully real-time
treatment can describe all dynamical processes contributing to the melting of
meson bound states. Recently, open quantum systems approaches to this problem
provided new insights in this direction. They treat the interactions between a
system (quark-antiquark pair) and its environment (thermal medium) in a real-
time Schrödinger formalism. Interestingly, it could link the melting process
to the quantum mechanical effect of decoherence (see e.g. [28, 29] and further
references in [5]).

###  II.2 The 2$D$ Ising QFT

The Ising spin model is defined by the Hamiltonian

$H=-J\,\left(\sum_{j=1}^{N-1}\sigma^{z}_{j}\sigma^{z}_{j+1}+h\sum_{j=1}^{N}\sigma^{x}_{j}+g\sum_{j=1}^{N}\sigma^{z}_{j}\right),$
(2)

where $\sigma_{j}^{\alpha}$ is the corresponding Pauli matrix at position $j$,
and we have considered an open chain of $N$ sites. Throughout our numerical
calculations, we use $J\equiv 1$. This sets an overall energy scale associated
with the lattice. Furthermore, $h$ and $g$ are, respectively, the transverse
and longitudinal fields.

In the pure transverse field Ising model ($g=0$), it is well known that a
quantum phase transition takes place at the critical point $J=h=1$ from a
ferromagnetic (ordered) phase ($h<1$) towards a paramagnetic (disordered)
phase ($h>1$). In the thermodynamic (infinite system size $N\to\infty$) and
continuum limit (vanishing lattice spacing $a\equiv 2/J\to 0$), the infrared
(IR) regime of the Hamiltonian (2) is effectively described by the Majorana
fermion QFT [30]

$H_{\text{IR}}=\int_{-\infty}^{\infty}\mathrm{d}x\,\left\\{\frac{i}{4\pi}\left(\psi\partial_{x}\psi-\bar{\psi}\partial_{x}\bar{\psi}\right)-\frac{iM_{h}}{2\pi}\bar{\psi}\psi+{\cal
C}M_{g}^{15/8}\,\sigma\right\\}.$ (3)

The free fermion mass $M_{h}$ is related to the transverse lattice parameter
$h$ through the relation $M_{h}\equiv 2J|1-h|$, while the longitudinal mass
scale is given by $M_{g}\equiv\,{\cal D}J\,|g|^{8/15}$. Here, ${\cal C}\approx
0.062$ and ${\cal D}\approx 5.416$ are numerical constants [30, 31].

As discussed in detail in [30], the QFT defined by the Hamiltonian (3) appears
in the scaling limit $M_{h}/J\to 0$ for fixed ratio $M_{h}/M_{g}$. It
comprises the following important classes in the vicinity of the critical
point. First, if $M_{h}=M_{g}=0$ (i.e. at the critical point $\\{h=1,g=0\\}$),
the Hamiltonian (3) represents the Ising conformal field theory (CFT), which
is a free Majorana fermion CFT with central charge $c=\frac{1}{2}$. There are
two scalar primary Hermitian operators,

$\epsilon\equiv i\bar{\psi}\psi\sim\sigma_{j}^{x}$ (4)

with scaling dimension $\Delta_{\epsilon}=1$ and $\sigma\sim\sigma_{j}^{z}$
with $\Delta_{\sigma}=\frac{1}{8}$. For $M_{h}\neq 0$, $M_{g}=0$, the
continuum limit of the transverse Ising model in both the ferromagnetic and
paramagnetic phase is integrable in terms of a massive free fermion QFT,
characterized by the fermion mass $M_{h}$. Elementary excitations in the
ferromagnetic phase can be interpreted as domain walls. As soon as a
longitudinal field is turned on ($g\neq 0$), these are non-perturbatively
confined [32]. The two remaining important regimes therefore give rise to
mesonic excitations with the following properties.

  * •

$M_{h}=0$, $M_{g}\neq 0$: This regime describes the integrable and interacting
E8 QFT of Zamolodchikov [33], which is mathematically captured by the
exceptional simple Lie algebra of rank 8. It contains 8 stable mesons as
nonperturbative fermionic bound states in the spectrum. Their masses $M_{n}$
are known and can be expressed as analytical ratios in terms of the lightest
mass $M_{1}\equiv M_{g}$.

  * •

$M_{h}\neq 0$, $M_{g}\neq 0$: This parameter range with both $\epsilon$ and
$\sigma$ perturbations turned on represents an interacting nonintegrable QFT
with stable and unstable meson bound states [33, 34, 35, 36].

We are primarily interested in studying the latter nonintegrable regime at
non-zero temperature, in which TN simulations have the most predictive power
since they provide an ab initio way of capturing the underlying meson physics.
Despite the simplicity of the original lattice model (2), we therefore can
study a highly nontrivial class of QFTs in (3), which we denote as the two-
dimensional (2$D$) Ising QFT in the following discussion.

##  III Methods and concepts

###  III.1 Tensor network simulations

We perform numerical simulations using matrix product states (MPS) and matrix
product operators (MPO) [7, 8, 37, 9]. These are TN ansätze that provide
efficient representations of the state and operators of the quantum many-body
system (see the appendix B for details). More concretely, a MPS parametrizes a
vector of a tensor product vector space by a number (typically one per site)
of tensors of dimension $d\times\chi\times\chi$, where $d$ is the physical
dimension of the degree of freedom associated to a lattice site, and $\chi$,
called bond dimension, determines the size of the variational family, and
bounds the entanglement entropy of the MPS. In the broadest sense, a MPO is
simply a MPS for the space of operators, where the local dimension is $d^{2}$
instead of $d$.

The MPS/MPO ansatz allows us to work directly in the infinite spatial volume
(thermodynamic limit), by assuming translational invariance and studying a
unit cell. For a nearest-neighbor model like (2), we use a two-site unit cell
and apply the infinite time-evolving block decimation (iTEBD) algorithm [38]
to simulate both imaginary and real-time evolution. The former allows us to
efficiently approximate thermal states as an MPO [39, 40], while the latter is
employed in the subsequent quench studies. The iTEBD algorithm is based on a
Suzuki-Trotter decomposition (see [41] and references therein) of the time
evolution operator into discrete time steps and a sequential application of
two-body operators acting on the MPS or MPO ansatz. The resulting bond
dimension increases, in general, under this operation, and in practice it is
necessary to truncate it. This is one source of error in the numerical
algorithms. While it also limits the reachable time scales under real-time
evolution, this will prove not to be a significant obstacle for the problem at
hand.

###  III.2 Entropies

In our work several entropic quantities play a pivotal role in unraveling the
meson melting. In particular, we make use of the von Neumann entropy and of
the more general Rényi entropies. The former is defined as

$S(\rho)=-\operatorname{Tr}\left[\rho\log\rho\right],$ (5)

where $\rho$ is the density matrix describing the state of the system. The
Rényi entropies 222Observe that in the limit $\alpha\searrow 1$,
$S_{\alpha}(\rho_{A})$ reduces to $S(\rho_{A})$. In the continuum or in
presence of gauge fields, the definition of reduced density matrices can lead
to mathematical subtleties, which, however, will not be relevant for our
subsequent spin chain calculations. are defined as

$S_{\alpha}(\rho)=\frac{1}{1-\alpha}\log(\operatorname{Tr}\rho^{\alpha}).$ (6)

In the case of pure states, when we consider a bipartition in subsystems A and
B, the corresponding entropies for the reduced density matrix
$\rho_{A}\equiv\operatorname{Tr}_{B}\rho$ of a subsystem measure the
entanglement with respect to this bipartition, and the von Neumann entropy (5)
$S(\rho_{A})=S(\rho_{B})$ is then referred to as the entanglement entropy.

In the case of pure states given by MPS, these entropies can be efficiently
calculated, because it is easy to recover the Schmidt decomposition across
each cut between a pair of sites, which in turn provides the eigenvalues of
the corresponding reduced density matrix for half the chain. In the case of
mixed states given by MPO, like the ones representing thermal states, these
quantities cannot be efficiently recovered, in general. We can however compute
the 2-Rényi entropy (or another low integer order) of the whole state.

Another option utilized by us is to consider a purification of the mixed
state, i.e. embedding a given mixed state as a reduced density matrix of a
pure state from an enlarged Hilbert space. A canonical purification of a mixed
state, demonstrated here on an example of a thermal density matrix, is given
by the thermofield double state

$|\Psi_{\beta}\rangle\propto e^{-\beta H/2}|\Phi\rangle,$ (7)

where $|\Phi\rangle$ is a maximally entangled state in a doubled system. Such
a canonical purification can be constructed also for a general mixed state.
For the canonical purification there is a natural identification of spatial
subsystems in the original system and corresponding subsystems in the copy.
This leads to the notion of the reflected entropy, which is the entropy
associated with tracing out the same spatial subsystems in both Hilbert space
factors in the canonical purification. One can efficiently compute the
reflected entropy using the MPS/MPO techniques and it will play an important
role in the rest of the paper. The reflected entropy was originally introduced
in [43] in a holographic setting as a proposed dual to twice the area of the
entanglement wedge cross section. Because of this, it has received a
considerable amount of attention in the context of the AdS/CFT correspondence
and QFTs, see e.g. [44, 45, 46, 47, 48, 49, 50, 51]. In particular, the works
[52, 53, 54] discussed it also for CFT quenches, whereas [55] analyzed a
variety of related mixed state entanglement measures as a holographic probe of
confinement.

The purification of a mixed state is not unique, since any unitary
transformation on the ancillary copy leave the state of the system unchanged.
This affects entanglement entropies in which bipartition includes parts in the
auxiliary system. One can define a version of the reflected entropy, called
the entanglement of purification [56], as the minimal one over all possible
purifications. The reflected entropy provides an upper bound to such an
optimal quantity, which in general is hard to compute. Note that the area of
the entanglement wedge cross section is conjectured to correspond to the
entanglement of purification [57]. If this conjecture is correct, it implies
that in holographic states optimization allows to lower the reflected entropy
by a factor of two by applying a unitary on a purifying system.

In the present work, we employ the TN ansatz to extract the discussed entropic
quantities for the Ising QFT. Before closing this section, we want to point
out that in a parallel vein, the authors of [58] calculated recently
entanglement measures in (1+1)-dimensional theories using a Hamiltonian
truncation approach [59], which provides an alternative computational ansatz.

###  III.3 Quantum quenches and quasiparticle model

We now turn to the concept of quantum quenches to study a dynamical situation
in the context of meson melting. A quantum quench means a rapid change of a
Hamiltonian parameter governing the time evolution of a physical system.
Initiated primarily by theworks [60, 61, 62], quantum quenches became a
theoretical framework to probe and characterize dynamical properties of
quantum many-body systems [63]. The modification of the underlying Hamiltonian
during time evolution kicks the system out of equilibrium. This drastic change
can induce new processes and effects through dynamical interactions.

Our studies are motivated by the seminal work [64]. Therein, the authors
analyzed the entanglement growth for quantum quenches in the Ising model. In
more detail, iTEBD methods were used to prepare the ground state
$\ket{\psi_{0}}$ of some Hamiltonian $H_{0}$, which then was evolved with a
Hamiltonian $H_{1}$ in the ferromagnetic phase in presence of meson states,
i.e. at finite longitudinal field values. It was observed that the existence
of mesons causes oscillations in the entanglement entropy at frequencies given
by their masses and mass differences. If the initial state is in the
paramagnetic phase (i.e. a quench across the phase boundary), the associated
entanglement growth persists mostly forever, whereas it is always bounded if
$\ket{\psi_{0}}$ is located in the ferromagnetic phase and mesons are produced
at rest. While the original work [64] considered a semi-classical parameter
regime away from the quantum critical point, the recent paper [65] observed
such entanglement oscillations also in the E8 QFT regime of the Ising model.
Similar observations were made for gauge theories in [66, 67, 68] as well as
for quasiparticles in paramagnetic quenches [69, 70]. Further aspects of
entanglement dynamics and thermalization behavior in spin chain models are
discussed, e.g., in [71, 72, 73].

Many properties of the entanglement production and spreading after quantum
quenches can be understood with a quasiparticle model. Originally developed in
[60] for the Ising model, it was later found to hold also for generic
integrable systems (in absence of mesons) [74, 75, 76] and thermofield double
states [77]. A generalization of the model was provided in [78]. In essence,
the model describes the initial state as a source of independent entangled
quasiparticle pairs, which are created at any given point and move after the
quench with opposite momentum and velocity through the quantum many-body
system and therefore spread quantum correlations. In lattice systems, the
maximum speed of propagation is limited by the Lieb–Robinson bound [79]. As a
consequence, the entropy of a finite-sized subregion increases linearly as
long as only one partner of the particle pair leaves it. When also the second
partner leaves the interval, the entanglement entropy saturates. This picture
is consistent with general entanglement scaling laws in the time evolution
found in [80]. For the special case of a semi-infinite spatial bipartition
this implies an indefinite entanglement growth. However, also the boundedness
of entanglement production for the meson case can be explained within the
quasiparticle picture [64, 81]. Namely, the existence of a confining
potential, due to the presence of meson states, induces that quasiparticles
bounce back as they get separated from their partners, causing the oscillatory
and bounded behavior for mesons produced at rest.

In this paper, we consider quenches in which the system is initially prepared
in a thermal equilibrium state, and then is evolved with a different
Hamiltonian. Such finite temperature quenches are much less explored than
those starting from a ground date, see [82, 83, 84, 85, 86, 87] for selected
results. In this scenario, we analyze the time evolution of different entropic
quantities to detect the melting of mesons at high temperatures. Since we are
interested in nonintegrable parameter regimes, we will make use of the
described quasiparticle model interpretations in our subsequent meson analyses
only in a qualitative fashion.

##  IV Results from thermal equilibrium

Figure 1: Left: Temperature dependence of the Rényi entropy density $s_{2}$. The curves are calculated using iMPO simulation at constant ratio $M_{h}/M_{g}\approx 0.09$ in the nonintegrable ferromagnetic phase using the parametrization (10). Right: Prony result of the $n=1$ curve of $s_{2}(\beta J)$ in the complexified $\beta J$ plane. The plots demonstrate that the low temperature behavior of $s_{2}$ is dominated by an exponential decay with frequencies matching the meson masses $M_{1}/J$ and $M_{2}/J$ (extracted from [35]), shown as black circles in the right panel. Table 1: High-temperature scaling exponents $p$ in the functional power-law ansatz (16) for several quantities $f$ and the parametrization (10). The numerical values are obtained from a fit to the iMPO simulations and are in good agreement with the analytical CFT expectation in the scaling limit. $f$ | $n=0$ | $n=1$ | $n=2$ | $n=3$ | exact
---|---|---|---|---|---
$s_{2}$ | 0.926 | 1.073 | 1.106 | 1.078 | 1
$E$ | 1.171 | 1.645 | 1.933 | 2.022 | 2
$\partial E/\partial\beta$ | 1.176 | 2.027 | 2.709 | 2.955 | 3

Figure 2: High-temperature scaling of the second Rényi entropy density $s_{2}$
(top) and thermal energy density $E$ (bottom) in the nonintegrable
ferromagnetic phase at $M_{h}/M_{g}\approx 0.09$ for the parametrization (10).
The grey line at $\beta M_{1}=1$ indicates the high-temperature threshold.
Dashed lines indicate the lattice scale for each curve. When taking the
scaling limit (increasing $n$) the results indicate the emergence of a power-
law behavior matching the CFT expectation, cf. Tab. 1.

In section II.1 we have discussed that in a QCD context meson melting is
indicated by a thermal broadening and movement of peaks in the in-medium
spectral function. For the Ising QFT, our previous work [88] studied the
retarded thermal correlation function w.r.t. the $\epsilon$ operator (4),
which is a closest analogue in this system to the simplest meson operator in
QCD. We did observe sharp meson peaks, whose residues decrease with the
temperature. Further, we observed that growing temperature causes the
appearance of singularities associated to meson mass differences, which we
interpreted as in-medium effects (for detailed discussions of the both effects
see also [89]). However, over the range of the temperatures explored we did
not see definite signatures of the meson broadening. This, of course, does not
mean that the mesons do not melt, but rather that the chosen observable in
conjunction with the range of parameters and the methods adopted are not
suitable to observe such an effect. In particular, our previous studies
required long time evolutions to resolve the peaks associated with the lowest
lying mesons, which became prohibitively costly for TN methods due to the
required size of the bond dimension at large temperatures.

For this reason in the present work we shift our viewpoint to the study of
entropies in the ferromagnetic regime of the Ising QFT. These quantities are
sensitive to classical and quantum correlations in the many-body system, and
can provide a complementary probe of meson features at high temperatures. In
this section, we consider a static setting, namely thermal equilibrium states
of the Ising model (2) at inverse temperature $\beta$, and study their second
Rényi entropy per site, defined as

$s_{2}=-\frac{1}{N}\log\frac{\operatorname{Tr}\rho_{\beta}^{2}}{(\operatorname{Tr}\rho_{\beta})^{2}}$
(8)

for the (unnormalized) thermal state $\rho_{\beta}=\exp[-\beta H]$.

We find a MPS approximation to the purification (7) in the thermodynamic limit
using the iTEBD algorithm to approximate $\exp[-(\beta/2)H]$. The transfer
operator of this MPS determines the norm of the state (and can be chosen to be
normalized by appropriate gauge choice). In order to determine the density
$s_{2}$, we thus just need the dominant eigenvalue $\eta$ of the transfer
operator of $\rho_{\beta}^{2}$, $\mathbb{E}_{\rho_{\beta}^{2}}$,333See
appendix C for details on the transfer operator.

$s_{2}=-\log\eta.$ (9)

The value of $\eta$ can be determined efficiently by considering the action of
$\mathbb{E}_{\rho_{\beta}^{2}}$ on an eigenvector and using iterative
eigensolvers.

For our meson studies we choose exemplarily the Ising QFT parameter point
specified by $M_{h}/M_{g}\approx 0.09$ in the nonintegrable ferromagnetic
phase. Additionally, we vary the individual masses $M_{h}$ and $M_{g}$
according to the parametrization

$M_{h,g}^{(n)}=\frac{M_{h,g}^{(0)}}{2^{n}},$ (10)

where $M_{h}^{(0)}/J=0.125$, $M_{g}^{(0)}/J\approx 1.356$ and the ratio
$M_{h}/M_{g}$ is kept constant. For increasing values of $n$, the masses are
decreasing, i.e. we approach the critical point in the phase diagram. This
means nothing else than taking the scaling limit previously described in
section II.2. The resulting masses $M_{1}/J$ of the first meson, which we read
out from [35], are hence also decreasing. For a fixed range of lattice spacing
values, $\beta J\in[1,16]$, we therefore can probe different physical
temperatures, measured in units $\beta M_{1}$.

Before discussing the results of these simulations, let us consider the
expected scaling behavior in the transverse field Ising model, i.e. without
longitudinal field perturbations causing confinement. In this case the model
can be mapped to a free fermion system (see e.g. [91]) from which the
partition function follows as

$Z_{\beta}=\operatorname{Tr}\left[\operatorname{\mathrm{e}}^{-\beta
H}\right]=\prod_{k}2\cosh\left(\frac{\beta\epsilon_{k}}{2}\right),$ (11)

where the single-particle energies are given by
$\epsilon_{k}=2J\sqrt{1+h^{2}-2h\cos(k)}$ for discrete wave number $k$. The
definition (8) implies immediately that the second Rényi entropy density can
be calculated as $s_{2}=2\beta[f(2\beta)-f(\beta)]$ from the free energy
density $f(\beta)=-\frac{1}{N\beta}\log Z_{\beta}$. In the thermodynamic
limit, this leads to

$s_{2}(\beta)=\log
2-\frac{1}{\pi}\int_{0}^{\pi}\mathrm{d}k\log\left[\frac{\cosh(\beta\epsilon_{k})}{\cosh^{2}(\frac{\beta\epsilon_{k}}{2})}\right].$
(12)

At low temperatures, i.e. in the limit $\beta\to\infty$, and away from
criticality the dominant contribution to $s_{2}$ originates from $k=0$, which
gives

$s_{2}\sim\frac{\operatorname{\mathrm{e}}^{-2\beta
J|h-1|}}{\sqrt{\frac{\pi\beta Jh}{|h-1|}}}.$ (13)

Note that the term $|h-1|$ is proportional to the free fermion mass
$M_{h}=2J|1-h|$. On the other hand, at the critical point ($h=1$), the
integral can be evaluated to give a power-law decay with the inverse
temperature,

$s_{2}\sim\frac{\pi}{16}\left(\frac{1}{\beta J}+\frac{1}{16(\beta
J)^{3}}+\ldots\right).$ (14)

The $1/\beta$ behavior is the (free) CFT result and the $1/\beta^{3}$
correction is a lattice effect.

The exponential scaling at low temperatures in the massive free fermion regime
is clearly distinct from the power-law decay at high temperatures at
criticality. While we do not have a prediction for the scaling in the
nonintegrable regime, it is tempting to assume that the temperature dependence
could be equally determined by the existence of a mass gap in the spectrum. To
test this hypothesis, we need to analyze both the low-temperature regime,
given for $\beta M_{1}\gg 1$, and the high-temperature regime, given by $\beta
M_{1}\lesssim 1$. In the latter case, we additionally need to ensure that the
result is not dominated by lattice excitations, i.e. $\beta J\gtrsim 1$.

Fig. 1 (left panel) shows the behavior of $s_{2}$ as a function of the inverse
temperature $\beta M_{1}$. In all cases, we observe an exponential decay at
low temperatures. In the right panel, we analyze this regime quantitatively
(exemplified for $n=1$) by applying a signal analysis technique based on the
Prony method onto $s_{2}$ as a function of $\beta J$ and decomposing it into
harmonic contributions with complex parameters. 444We refer to appendix D for
more details on the method. Two clear poles on the imaginary axis are visible.
They quantify the exponential decay at low temperatures. The values agree with
the masses of the first two meson states (shown as black circles). Based on
these findings, we postulate the following low-temperature scaling of the
Rényi entropy density

$s_{2}\sim\sum_{i}c_{i}\operatorname{\mathrm{e}}^{-\beta M_{i}}$ (15)

for some coefficients $c_{i}$ and meson masses $M_{i}$, whereby the lowest
mass $M_{1}$ is obviously dominating. This form of the asymptotic scaling is
the natural generalization of the massive free fermion result
$s_{2}\sim\operatorname{\mathrm{e}}^{-\beta M_{h}}$ in (13).

Fig. 2 (top panel) shows again $s_{2}$ as a function of $\beta M_{1}$, but
emphasizing the high-temperature regime in a double-logarithmic plot. To the
left of the grey line the high temperature regime $\beta M_{1}\lesssim 1$
starts, while the dashed lines denote the estimated lattice scale for each
curve. As $n$ is increasing, i.e. the individual masses $M_{h,g}$ are
decreasing, one can observe the emergence of a linear scaling in between these
two curves, indicating a power-law behavior. Assuming the general functional
dependence

$f(\beta)\sim\beta^{-p},$ (16)

we extract in Tab. 1 the values of the power-law exponents $p$ from a linear
fit to the numerical data in between these two scales. For $f=s_{2}$, the
value $p=1$ is obtained with a very good accuracy of only a few percent. This
implies that the high-temperature scaling behavior matches the CFT expectation
from (14). Based on the clear observation from fig. 1 of meson states in the
entropy scaling at low temperatures and their absence at high temperatures, we
interpret this signature as strong evidence for meson melting. In this regime,
which is not yet influenced by the lattice scale, melted meson states do not
leave an imprint on the entropy scaling any more, which is why the scaling is
identical to a critical system. To further corroborate this finding, we
additionally analyse also the scaling behavior of the thermal energy density
(bottom panel in Fig. 2) and its derivative. Both these thermodynamic
quantities show a similar power-law scaling at high temperatures. The
corresponding exponents are tabulated in Tab. 1. At $n=3$ they agree well with
the analytical expectation, which is set by the scaling dimension $\Delta=2$
of the associated CFT operator. Interestingly, the Rényi entropy $s_{2}$ seems
to allow an even better estimation of the scaling exponent for all values of
$n$, indicating the potential of entropic quantities to identify the meson
melting process.

In summary, we have shown that the scaling behavior of the second Rényi
entropy density provides a clear signature of the asymptotic high and low-
temperature regime in a nonintegrable Ising QFT with meson bound states. At
low temperatures, we found that the behavior is dominated by the exponential
damping of the form (15). At high temperatures, the power-law dependence (16)
matches the CFT behavior. We interpret the latter scaling as signaling the
fact that the meson states have been melted. In fact, a similar analysis for
the E8 regime yields the same qualitative picture, providing evidence that the
results are independent from integrability properties of the underlying
system. Naively, it might look trivial that a massive QFT matches at high
temperatures the CFT behavior due to the fact that the temperature then
provides effectively the only scale in the system. (A similar situation for
the massive free fermion regime in the Ising model is discussed for example in
[91].) However, in our present studies we do have nonperturbative meson bound
states in the spectrum, which are “visible” at low temperatures, and we
observe a smooth transition to the high-temperature regime, which we can
disentangle from the UV lattice scale. We therefore see the scaling properties
of the second Rényi entropy as a clear witness of the meson melting process in
a static situation. While we have focused on the second order Rényi entropy,
it is possible that this conclusion holds also for higher-order ones.

##  V Results from dynamics after thermal quenches

###  V.1 Setup

We consider quenches from the thermal equilibrium state of a pre-quench
Hamiltonian $H_{0}$. At time $t=0$ an instantaneous global quench is applied
and the system evolves with Hamiltonian $H_{1}$. In particular, we fix for
$H_{1}$ the parameters that correspond to $M_{h}^{(0)}$ and $M_{g}^{(0)}$ from
the parametrization (10) (namely $h=0.9375$ and $g\approx 0.0746$), as this
corresponds to the smallest time period $2\pi J/M_{1}\approx 4.2$ (for
$M_{1}/J\approx 1.5$) of the first meson, and we know, from our earlier
equilibrium studies [88], that in this regime, MPO simulations and Prony
signal analysis [93] allowed a reliable identification of the QFT meson
masses.

The previously considered global Rényi entropy density $s_{2}$ in eq. (9) is
constant under unitary time evolution, and computing its corresponding value
for a finite subsystem is computationally too demanding for large bond
dimensions. But the MPS approximation of the purification gives us access to
the reflected entropies, which correspond to partitioning the purification in
half. If $\\{\lambda_{\alpha}\\}$ are the Schmidt values across a cut of the
purified state, the von Neumann and 2-Rényi reflected entropies are given by

$\displaystyle\tilde{s}_{1}$
$\displaystyle=-\sum_{\alpha}\lambda_{\alpha}^{2}\log(\lambda_{\alpha}^{2}),$
(17) $\displaystyle\tilde{s}_{2}$
$\displaystyle=-\log\left(\sum_{\alpha}\lambda_{\alpha}^{4}\right).$ (18)

In this section, we study the growth and oscillatory behavior of these
quantities after a thermal quench.

Table 2:  Comparison of the pre-quench temperatures $\beta J$ of the initial thermal state (w.r.t. $H_{0}$) and the resulting post-quench effective temperatures $\beta^{*}J$ and $\beta^{*}M_{1}$ (w.r.t. $H_{1}$) for the different quench types shown in Fig. 3. initial temperature $\beta J$ | effective temperature $\beta^{*}J$, ($\beta^{*}M_{1}$)
---|---
| 1 | 2 | 3 | 4 | 5 | 6
16 | 6.10 | 3.55 | 1.55 | 0.91 | 3.55 | 0.93
| (9.3) | (5.4) | (2.4) | (1.4) | (5.4) | (1.4)
2 | 2.00 | 1.92 | 1.34 | - | - | -
| (3.0) | (2.9) | (2.0) |  |  |
0.97 | - | - | 0.91 | - | - | -
|  |  | (1.4) |  |  |
0.5 | 0.50 | - | - | - | - | -
| (0.8) |  |  |  |  |
Figure 3: Overview of the different quench protocols in the transverse ($h$)
vs. longitudinal ($g$) field plane. Thermal initial states are prepared at
inverse temperature $\beta$ for the parameter locations marked by green dots.
By quenching to the final nonintegrable ferromagnetic point, denoted by red
arrows towards the black cross, an effective temperature $\beta^{*}$ as listed
in Tab. 2 is induced by the thermal quench. The E8 regime is marked by the
grey dotted line.

Figure 4: Time dependence of the first (dashed curves) and second (solid
curves) reflected entropies for the quench protocols (1)-(4) (as defined in
Fig. 3) with an initial temperature $\beta J=16$. With increasing effective
temperature (in order of the quench labels) the entanglement growth is getting
enhanced and suppresses the entanglement oscillations. The plateau arising in
quench type (4) at late times is an unphysical effect caused by the truncated
bond dimension (see main text for further discussions).

After the quench, the state is not in thermal equilibrium for $H_{1}$.
However, we can define an effective temperature, which corresponds to the
thermal equilibrium of $H_{1}$ at the same (conserved) energy of the system.
The effective inverse temperature $\beta^{*}$ is then given by the condition
555Note that we always consider densities here since we work with the unit
cell in the iMPO formalism. For our global quench setup, the relation (19) is
then identically also valid for the whole system.

$\frac{\operatorname{Tr}\left[H_{1}\operatorname{\mathrm{e}}^{-\beta
H_{0}}\right]}{\operatorname{Tr}\left[\operatorname{\mathrm{e}}^{-\beta
H_{0}}\right]}=\frac{\operatorname{Tr}\left[H_{1}\operatorname{\mathrm{e}}^{-\beta^{*}H_{1}}\right]}{\operatorname{Tr}\left[\operatorname{\mathrm{e}}^{-\beta^{*}H_{1}}\right]}.$
(19)

Since we keep the final Hamiltonian $H_{1}$ fixed, the effective temperature
can be increased either by varying the initial temperature $\beta$ w.r.t. a
fixed pre-quench Hamiltonian $H_{0}$ or by modifying the parameters of $H_{0}$
for fixed $\beta$. We will employ both methods in our studies. In the latter
scenario, we choose the initial points in the phase diagram as shown in Fig.
3. The different quench protocols, labelled by the numbers in Fig. 3, result
in the effective temperatures $\beta^{*}J$ (measured w.r.t. $H_{1}$)
summarized in Tab. 2 (determined up to two digits) for several initial
temperatures $\beta J$ (measured w.r.t. $H_{0}$). The quench Hamiltonian
$H_{1}$ corresponds to the fixed parameters $h=0.9375$ and $g\approx 0.0746$.
In the quench protocol 1, we have $h=0.93$, and for 5 $h=0.8732$. In type 6,
the longitudinal field is chosen asymptotically large as $g=100$. The
remaining initial parameters are identifiable from Fig. 3. Note in particular
that protocol 4 starts in the classical regime, whereas 2 is in the E8 phase
and 3 at the critical point. The different quench types encompass pure
transverse quenches ( 1, 2, 4, 5) and longitudinal quenches 6. Otherwise, both
fields are quenched. We consider the case of an initial low temperature (i.e.
close to the ground state) at $\beta J=16$ as well as higher initial
temperatures, cf. Tab. 2. Note that protocol 5 is tuned such that is has the
same effective temperature as 2 for $\beta J=16$. Similarly, type 4 and 6 at
$\beta J=16$ and 3 at $\beta J=0.97$ result in (nearly) the same effective
temperature. Overall, we have adjusted our setup to achieve both low physical
temperatures $\beta^{*}M_{1}\gg 1$ as well as high temperatures
$\beta^{*}M_{1}\lesssim 1$, while avoiding the lattice regime $\beta^{*}J\ll
1$. 666Note that the effective temperatures $\beta^{*}J=2.00$ and
$\beta^{*}J=0.50$ for type (1) mean that the difference to the initial $\beta
J$ is smaller than our considered precision of two digits.

###  V.2 Effective temperature effects

We start with the analysis of quenches from an initial state close to the
ground state of $H_{0}$, prepared at low temperature $\beta J=16$. Fig. 4
shows the results for $\tilde{s}_{1}$ (dashed curves) and $\tilde{s}_{2}$
(solid curves) as a function of time (in units of $J$) for the quench
protocols 1\- 4. As listed in Tab. 2, the effective temperatures are in
increasing order from $\beta^{*}M_{1}\approx 9.3$ to $\beta^{*}M_{1}\approx
1.4$. As can be appreciated from Fig. 4, the reflected entropy grows faster
for larger effective temperature. For our discussion, we limit ourselves to
regimes which are numerically achievable with our MPO simulations. When the
entanglement grows fast, this is possible only up to a certain time, after
which the observed entropy saturates due to a finite bond dimension effect
(see e.g. Fig. 4 after $tJ\approx 20$). Thus our analysis is limited to early
and intermediate time scales.

In the quench types 1 and 2, which are shown in detail in the left column,
$\tilde{s}_{1}$ and $\tilde{s}_{2}$ exhibit an oscillatory behavior, similar
to the ground state case first discussed in [64]. As we will analyze in detail
below, these entanglement oscillations are caused by the QFT meson states at
characteristic frequencies corresponding to their masses. Type (2) (orange
curves) shows on top of these entanglement oscillations an overall increase.
This can be explained in the quasiparticle model either through the production
of mesons with finite velocities, or the fact that the arising envelope
frequencies (that bound the entanglement growth when mesons are produced at
rest) would require (much) longer time scales. The latter case is plausible
when keeping in mind that we have chosen a parameter point in the QFT regime
of the Ising model close to its critical point. The meson masses, e.g. in
absolute units of $J$, are smaller as compared to a semi-classical regime and
hence result in much longer time recurrences. From Fig. 4 it becomes apparent
that the entanglement oscillations are heavily suppressed as the effective
temperature is raised from type 1 to 4. This indicates that the meson states
do not dominate the entanglement growth at their characteristic mass
frequencies anymore. In particular, the first reflected entropy of quench
protocol 4 (red dashed curve in the right panel) grows fast at early times,
with an apparently linear behavior. Below, we will study each of these effects
in detail.

1

2

3

Figure 5: Results of the Prony signal analysis (top) and Fourier spectra
(bottom) of $\tilde{s}_{2}$ for the quench protocols (1)-(3) (from left to
right) with an initial temperature $\beta J=16$ and effective temperatures
$\beta^{*}M_{1}$ as shown. Solid vertical lines represent meson masses from
ref. [35] and dashed lines mark the continuum threshold of $2M_{1}$. Green
dash-dotted lines mark all possible mass differences of the first 4 mesons.
Both methods allow an identification of meson states as poles in the complex
frequency plane or peaks in the Fourier spectrum. When the effective
temperature is raised, meson poles become fuzzier, while the spectrum
increases and flattens out.

These findings are in general valid for both $\tilde{s}_{1}$ and
$\tilde{s}_{2}$ as visible in Fig. 4. Observe, however, that the growth of
$\tilde{s}_{2}$ is suppressed as compared to $\tilde{s}_{1}$. For this reason,
it turns out that the second reflected entropy allows for a more precise
identification of meson masses from the entanglement oscillations when
evaluating both the Prony signal analysis method as well as the Fourier
spectra. That is why we consider $\tilde{s}_{2}$ for a detailed analysis in
the following discussion. The results are presented in Fig. 5 for protocols
1\- 3 (from left to right). At the lowest effective temperature
$\beta^{*}M_{1}\approx 9.3$ for quench type 1 (left column), the Prony
analysis allows a clear identification of the first 5 meson poles, which are
in good agreement with their QFT mass values when compared to the predictions
in [35] (solid vertical lines). Additional features are visible also at
$M_{6}$ to $M_{8}$, which, however, become fuzzier and hence more uncertain.
In addition, the continuum threshold at $2M_{1}$ is identifiable as a vertical
line of poles, indicating a branch cut (shown as the dashed vertical line). In
the corresponding Fourier spectra, the meson poles translate into peaks at
their frequency values. The peak of the first meson is largely dominating over
the other ones. Overall, the Fourier spectrum is decreasing over several
orders of magnitude towards larger frequencies. The green dash-dotted lines in
the Fourier spectra mark all 6 possible mass differences between the first 4
meson states. At their respective values, kinks appear in the spectrum,
corresponding to additional poles in the Prony plot. They lie at frequencies
smaller than the first meson mass as well as between $M_{1}$ and $M_{2}$.

As the effective temperature is raised to $\beta^{*}M_{1}\approx 5.4$ in
protocol 2 (middle column) and $\beta^{*}M_{1}\approx 2.4$ in protocol 3
(right column), the meson pole identifications become fuzzier in the Prony
analysis. Notably, the absolute values of the corresponding Fourier spectra
are raised by 3 and 6 orders of magnitude, respectively. For protocol 3, the
Fourier spectrum is flattening out (relative to the overall magnitude). Since
in all cases the post-quench nonintegrable QFT regime is identical, we
attribute this effect to the increasing effective temperature. This hypothesis
is corroborated by data from a different quench with the same effective
temperature (see Fig. 6). When comparing the meson peak heights to the overall
magnitude of the Fourier spectrum, the impact of the meson mass frequencies
seems to get suppressed at high temperatures. We see this thermal suppression
effect on the Fourier spectrum as a first hint of how the melting of meson
states at even higher temperatures can be detected. A sequential process in
similarity to QCD (cf. the discussion in section II.1) could show up in this
description through a stepwise decrease and disappearance of individual meson
peaks in the Fourier spectrum of the reflected entropies. The available data
in Fig. 5 are, however, inconclusive about whether this process is consistent
with this picture.

Figure 6: Comparison of the quench protocols (2) (blue curves) and (5) (yellow
curves) for an initial temperature $\beta J=16$ and the same effective
temperature $\beta^{*}M_{1}\approx 5.4$. Top: Time dependence of the first
(dashed) and second (solid) reflected entropy in both models. Bottom:
Corresponding Fourier spectra of the second reflected entropy. The different
regime of the initial state (E8 vs. nonintegrable ferromagnetic) causes only
mild differences in the behavior of entanglement growth and oscillations.

To analyze the influence of the initial state in more detail, we compare
quenches from different initial states at the same effective temperature. Fig.
6 shows the time evolution of the reflected entropies and resulting Fourier
spectra of quench protocols 2 and 5 starting from a pre-quench thermal state
at $\beta J=16$. Both states have the same effective temperature
$\beta^{*}J\approx 3.55$ or $\beta^{*}M_{1}\approx 5.4$, but are initially
prepared in the E8 regime versus the nonintegrable ferromagnetic phase. The
figure shows that $\tilde{s}_{1}$ (dashed curves) and $\tilde{s}_{2}$ (solid
curves) differ only marginally in the two quench protocols and exhibit a
slight phase shift. The first meson peak dominates the Fourier spectrum
equally in both models. For the available frequency resolution and range, both
Fourier spectra are in close correspondence with each other. These results
indicate that the described features, at least for this particular case, are
primarily driven by thermal effects, set by the effective temperature, while
the influence of the initial state seems subdominant. Below, we will find a
similar situation also for other models at even higher temperatures, providing
evidence for the robustness of the physical interpretation about thermal
effects from our findings.

1

2

3

Figure 7: Time dependence of the reflected entropies for the quench protocols
(1)-(3) (from left to right) with an initial temperature $\beta J=2$ and
effective temperatures $\beta^{*}M_{1}$ as shown. The first (dashed) and
second (solid) reflected entropy are shown for several bond dimensions (blue
vs. orange curves), indicating that the results have not yet converged at late
times.

1

2

3

Figure 8: Results of the Prony signal analysis (top) and Fourier spectra
(bottom) of $\tilde{s}_{2}$ for the quench protocols (1)-(3) (from left to
right) with an initial temperature $\beta J=2$ and effective temperatures
$\beta^{*}M_{1}$ as shown. Results are shown for the $\chi=600$ curves in Fig.
7. Background lines are as in Fig. 5. The plots exemplify the effect of
increased effective temperature through raising the initial temperature as
compared to Fig. 5: Except the E8 initial regime, meson states do not leave a
considerable impact on the entanglement growth through entanglement
oscillations.

We now study the situation of raised effective temperature by starting from an
initial thermal state at a higher temperature $\beta J=2$. In this case, the
entanglement growth is larger and it becomes important to monitor the behavior
of the iTEBD simulation in dependence of the bond dimension. The obtained
reflected entropies are shown in Fig. 7 for the quench protocols 1\- 3 (from
left to right) for two different bond dimensions, $\chi$=400 (blue) and 600
(orange). The effective temperatures for protocol 1 and 2 are very close,
$\beta^{*}M_{1}\approx 3.0$ and $\beta^{*}M_{1}\approx 2.9$, resulting in
similar curves for the reflected entropies. For quench type 3, the effective
temperature $\beta^{*}M_{1}\approx 2.0$ causes an entanglement growth about
twice as fast. At long times we observe the results for both bond dimensions
differ, indicating that they are not fully converged. This effect is milder
for $\tilde{s}_{2}$ (both on an absolute and relative scale), which is most
likely because $\tilde{s}_{2}$ is more local than $\tilde{s}_{1}$. Due to the
larger entanglement growth, the discrepancy is larger for quench type 3.

Being aware of these limitations, we compare in Fig. 8 the Prony results and
Fourier spectra of $\tilde{s}_{2}$ for the different quench protocols (for the
highest bond dimension) for an overall qualitative assessment. In all
examples, features in the complex frequency become fuzzier in the Prony
analysis when compared to the quench results in Fig. 5 at lower temperatures.
Some of the meson poles fully disappear, in particular for type 1 and 3 (top
left and right panel). Observe also that the branch cuts indicating the
continuum threshold disappeared in all types. The corresponding Fourier
spectra are increased by several orders of magnitude and now flattened out in
all cases. Although type 1 and 2 have a comparable effective temperature, the
Prony result for the latter (top middle panel) still allows a much clearer
identification of some meson states in the Prony analysis. This difference
seems to be rooted in the fact that initial states were prepared in distinct
phases with different properties.

Figure 9: Time dependence of the first (dashed) and second (solid) reflected
entropy for quench protocols (3), (4) and (6) with an effective temperature
$\beta^{*}M_{1}\approx 1.4$. The vertical dashed grey line indicates
convergence in the bond dimensions up to $tJ\lesssim 9$. The curves differ
only mildly although initialized in different phases and temperatures,
implying that the described features are induced by increasing the effective
temperature. Figure 10: Time dependence of $\tilde{s}_{1}$ (dashed blue) and
$\tilde{s}_{2}$ (solid blue) for quench protocol (1) with an effective
temperature $\beta^{*}M_{1}\approx 0.8$. The vertical dashed grey line
indicates convergence in the bond dimensions up to $tJ\lesssim 7.5$. The
results for $\tilde{s}_{1,2}$ in the nonintegrable ferromagnetic phase (blue
curves) are compared to a massive free fermion counterpart (orange curves) at
the same effective temperature in units of the lattice spacing. The black and
grey dashed and solid lines mark linear fit functions for $\tilde{s}_{1}$ and
$\tilde{s}_{2}$ at early times, respectively. The results demonstrate the
emergence of an overall linear growth behavior at high effective temperatures
in the nonintegrable ferromagnetic phase, which is consistent with the
interpretation of melted meson states which do not dominate the entanglement
growth anymore. This behavior is analogous to the free fermion regime, in
which no meson bound states are present.

The overall phenomenological picture that arises from these analyses is that
an increase of the effective temperature leads to an enhancement of the
entanglement growth, which suppresses the impact of entanglement oscillations
caused by meson states in the thermal bath. In absence of mesons, the
quasiparticle model predicts a fully linear growth of entanglement entropies.
We therefore aim in the following discussion to identify a linear growth of
the reflected entropies as a signature of meson melting at the highest
attainable temperatures. For this line of reasoning, we consider in Fig. 9 the
time evolution of reflected entropies for the quench types 3, 4 and 6 with an
(nearly) identical effective temperature $\beta^{*}M_{1}\approx 1.4$. Due to
the large entanglement growth, numerical results are converged in bond
dimension only up to times $tJ\lesssim 9$ (shown by the grey dashed line). At
later times the entropies saturate (cf. also the red curves in Fig. 4 for type
4), a signature of the truncation error. Altogether, there are only marginal
differences for $\tilde{s}_{1}$ between all quench protocols and a small
scaling difference for $\tilde{s}_{2}$ in type 3 (blue solid curve) compared
to the others. Since the initial states are in very different regimes
(classical, ferromagnetic and critical), we interpret the overall agreement as
a confirmation that the discussed effects on the entanglement growth are
indeed thermally induced. Fig. 9 therefore exhibits directly also a linear
function, which is obtained by fitting $\tilde{s}_{1}$ in the time interval
$2\leq tJ\leq 9$ as $0.09+0.42(tJ)$. The data for $\tilde{s}_{1}$ are bounded
by such a linear growth, but do not yet follow exactly this scaling behavior.
On the other hand, the curves for $\tilde{s}_{2}$ grow less strongly in a
nonlinear fashion. In this case, the time ranges for which the data are
converged are too short to enable a robust Prony or Fourier signal analysis.

The highest effective temperature $\beta^{*}M_{1}\approx 0.8$ was achieved for
quench protocol 1. The results, displayed by blue curves (dashed for
$\tilde{s}_{1}$, solid for $\tilde{s}_{2}$) in Fig. 10, converged in bond
dimension up to $tJ\approx 7.5$ (indicated by the vertical dashed grey line).
One can observe a strong growth of both reflected entropies. In the time
interval $1.5\leq tJ\leq 7.5$, we can fit the functions $0.04+0.67(tJ)$ and
$0.18+0.53(tJ)$ to $\tilde{s}_{1}$ and $\tilde{s}_{2}$, respectively. In
contrast to the simulations at intermediate temperatures, both reflected
entropies now seem to follow this linear growth behavior after a brief initial
time period with irregular growth induced through the quench itself. The same
phenomenological picture arises independently in quench type 2 (not shown
here), i.e. with an initial state in the E8 regime, at the same temperature.
As alluded before, the existence of mesons necessarily implies the appearance
of entanglement oscillations in the time evolution of entropic quantities,
consistent with a quasiparticle model interpretation, if applicable in the
nonintegrable range. Since in the present case we have $\beta^{*}M_{1}<1$ in
the physical high-temperature regime, we explain the linear behavior through
the fact that the mesons are melted, i.e. do not contribute entanglement
oscillations to the overall growth of entanglement entropies anymore.

This interpretation is corroborated through a quench study in the (massive)
free fermion regime. Complementary to protocol 1, we choose for that case a
transverse quench using the parameters
$(h_{0},g_{0},\beta)=(0.93,0,0.5)\rightarrow(h_{1},g_{1},\beta^{*})=(0.9375,0,0.50)$.
The resulting reflected entropies (orange curves) are compared in Fig. 10 to
the previous nonintegrable counterpart (blue curves). Although the effective
temperatures are not directly comparable anymore due to the different mass
scales, the differences in the entanglement growth between these two quench
types are relatively small. In particular, the free fermion data confirm the
linear growth of the reflected entropies after the initial quench period. We
can respectively fit the functions $0.01+0.68(tJ)$ and $0.08+0.59(tJ)$ to
$\tilde{s}_{1}$ and $\tilde{s}_{2}$, which are shown as grey curves and can be
compared to the black ones for the previous nonintegrable regime. These
results confirm that a linear growth of reflected entropies is a property of
the free fermion regime, which is in accordance with the quasiparticle model
description. It therefore strengthens our interpretation on meson melting for
the nonintegrable ferromagnetic case, since this qualitative behavior can be
explained when meson states are not present in the physical system.

###  V.3 Quench results in the high temperature scaling limit

Figure 11: Dependence of $\tilde{s}_{1}$ (dashed) and $\tilde{s}_{2}$ (solid)
on the physical time $tM_{1}$ for quenches of type (1) in the high temperature
scaling limit, when mesons are expected to be melted (left panel). The post
quench masses are varied according to the parametrization (10), resulting in
higher effective temperatures for increasing $n$. Only the time intervals of
convergence are shown. The reflected entropies are growing more rapidly when
approaching the scaling limit for increasing $n$ and can be fitted by linear
functions of the form (20), whose inverse slopes $r^{-1}$ are plotted in the
right panel in dependence on the inverse effective temperature
$\beta^{*}M_{1}$ (blue squares for $\tilde{s}_{1}$, orange circles for
$\tilde{s}_{2}$). The inverse slopes themselves approach the infinite
temperature limit in a linear fashion (solid blue and dashed orange fit
curves), implying a high temperature dependence of the form (21) for the
reflected entropies.

The results of the previous section have shown that we can extract meson
states from entanglement oscillations of reflected entropies at low and
intermediate temperatures in the nonintegrable ferromagnetic phase. Their
masses match quantitative expectations of the Ising QFT. When the regime of
high effective temperatures is reached, the linear entropy growth suppresses
all oscillations. In this case, which we interpreted through the melting of
mesons, the entanglement growth is very large and its QFT origin is not
immediately visible anymore. Similarly to the simulations in thermal
equilibrium, it then makes sense to consider the quench scenario also in the
scaling limit in order to approach the QFT regime rigorously.

We therefore now focus on the scaling limit, defined as $M_{h}/J\to 0$ for
fixed ratio $M_{h}/M_{g}$, for quench protocol 1 at the highest effective
temperatures. This setup leaves some degrees of freedom, since we can
parametrize the pre-quench ($M_{h}^{pre},M_{g}^{pre}$) and post-quench
($M_{h}^{post},M_{g}^{post}$) masses as well as the initial $(\beta)$ and
effective temperatures $(\beta^{*})$. Here, we keep the ratios
$M_{h}^{pre}/M_{g}^{pre}\approx 0.10$ and $M_{h}^{post}/M_{g}^{post}\approx
0.09$ constant as in the original quench setup, and vary the post-quench
masses again according to the parametrization (10). Since we perform a
transverse quench, we have $M_{g}^{pre}=M_{g}^{post}$ and hence
$M_{h}^{pre}/M_{g}^{pre}=M_{h}^{pre}/M_{g}^{post}$.

In the first set of simulations, we choose, identically to the previous case,
the initial temperature $\beta J=0.5$, which results in different physical
effective temperatures $\beta^{*}M_{1}$. The results for the reflected
entropies are shown in the left panel of Fig. 11. The time dependencies of
$\tilde{s}_{1}$ (dashed curves) and $\tilde{s}_{2}$ (solid curves) are shown
in physical units, i.e. $tM_{1}$, where $M_{1}$ is the first meson mass read
out from [35]. We plot only the early time window, in which convergence is
achieved. As we increase $n$, i.e. as we approach the critical point in the
scaling limit $M_{h}/J\to 0$, the reflected entropies grow faster. In all
cases, a linear growth behavior of the form

$\tilde{s}_{1,2}(t)=r\,t\,M_{1}+s_{1,2}(t=0)$ (20)

can be deduced from the results after a very short initial quench phase. Here,
$r$ denotes the slope of the linear dependence. The numerical fit results for
$r^{-1}$ are shown in the right panel of Fig. 11 as a function of the
(inverse) effective temperature $\beta^{*}M_{1}$. Both the results for
$\tilde{s}_{1}$ (blue squares) and $\tilde{s}_{2}$ (orange circles) can be
well fitted by linear functions through the origin as
$r^{-1}_{\tilde{s}_{1}}=(2.93\pm 0.02)\beta^{*}M_{1}$
$r^{-1}_{\tilde{s}_{2}}=(3.60\pm 0.10)\beta^{*}M_{1}$. This, unsurprisingly,
tells us that the slope diverges when reaching the infinite temperature limit.
It, however, shows that this limit is approached with a proper scaling
behavior, set only by the temperature. This means that the reflected entropies
follow in the high temperature regime, up to a numerical constant and shift,
the time dependence

$\tilde{s}_{1,2}(t)\sim\frac{t}{\beta^{*}}.$ (21)

This behavior is consistent with the interpretation that mesons are melted,
since their mass scale does not leave an imprint on the time evolution
anymore.

Figure 12: Dependence of $\tilde{s}_{1}$ (dashed) and $\tilde{s}_{2}$ (solid)
on the physical time $tM_{1}$ for quenches of type (1) in the high temperature
scaling limit. The post quench masses are varied according to the
parametrization (10) at fixed physical temperature $\beta^{*}M_{1}\approx
0.76$. The entropy curves converge for increasing $n$ and confirm an overall
linear increase associated to mesons being melted.

As elaborated before, the scaling limit, in which the Ising QFT emerges, can
be taken also in a different way. While we have so far considered a situation
in which the effective temperature is raised, we now keep
$\beta^{*}M_{1}\approx 0.76$ fixed. 777This parameter value corresponds to the
blue $n=0$ curves in Fig. 11. Using the same parametrization (10) for the
post-quench parameters, this implies the series $\beta^{*}J=\\{0.5,1,2,4\\}$
for the effective temperatures in transverse quenches of type 1. 888The
initial temperature $\beta J$ follows from demanding (19), and $M_{h}^{pre}$
is uniquely determined from the specified ratio $M_{h}^{pre}/M_{g}^{pre}$.
Fig. 12 shows the resulting time dependencies of $\tilde{s}_{1}$ (dashed
curves) and $\tilde{s}_{2}$ (solid curves). It becomes discernible that the
first reflected entropy grows faster when $n$ is increasing from $n=0$ to
$n=2$. The results for $\tilde{s}_{1}$ at $n=2$ (green dashed curve) and $n=3$
(red dashed curve) are coinciding, indicating that convergence is achieved.
For $\tilde{s}_{2}$, convergence is approached slightly slower in the
available time window (solid orange, green and red curves). These findings
confirm the overall qualitative picture of an linear increase in absence of
meson oscillations, which persists in the high-temperature scaling limit of
the Ising QFT.

##  VI Discussion and outlook

In this article we initiated and performed an in-depth study of the meson
melting phenomenon from a quantum information perspective. We consider the
2$D$ Ising QFT emerging in the continuum limit of the quantum Ising model as a
test bench for a novel entanglement viewpoint on the topic. This is possible
because the quantum Ising model contains nonperturbative meson bound states in
integrable and nonintegrable ferromagnetic parameter regimes, which resembles
the relevant confinement property with QCD. It therefore can provide important
insights into phenomena manifested also by much more complicate theories.
While this study is motivated by the analogous QCD effect, we however
emphasize that in the present analysis we do not even attempt to draw any
conclusions for QCD itself.

From a modern QCD viewpoint, the meson melting mechanism is interpreted as a
sequential process that strictly requires a dynamical treatment. While both
QCD and holographic approaches indicate the melting of individual meson states
through a thermal broadening of their corresponding peaks in the spectral
function, in our previous work [88] could not observe an analogous feature in
the complex frequency structure of retarded thermal equilibrium correlators in
ferromagnetic parameter regimes of the 2$D$ Ising QFT. We therefore motivated
in this paper a new paradigm to capture and describe meson melting
characteristics via entropic quantities. For this purpose we employed TN
simulations, which enable the approximation of thermal states as well as their
subsequent real-time evolution. Moreover, TNs easily allow calculations of
several entropies directly in the thermodynamic limit.

In particular, for a thermal equilibrium setting, we showed in section IV that
the scaling of the second Rényi entropy density is determined at low
temperatures through an exponential damping (as a function of the inverse
temperature $\beta$) set by the meson masses. In contrast, at high
temperatures, the behavior is identical to that expected in a CFT (power-law
dependence), meaning that meson states have been melted, i.e. the temperature
sets the only scale in the physical system.

An independent, truly out-of-equilibrium probe is provided by the real time
evolution of reflected entropies after thermal quantum quenches (cf. section
V). In this setting, we prepared initial thermal states in the ferromagnetic
phase of the Ising model and implemented a quantum quench by simulating their
evolution under a fixed post-quench Hamiltonian in the nonintegrable
ferromagnetic phase. Whereas the initial state is not in equilibrium, we can
identify an effective temperature corresponding to the energy density of the
state. We could vary the effective temperature of this setup either through
different parameter regimes in the pre-quench state or by modifying its
initial temperature. While at low effective temperatures mesons give rise to
entanglement oscillations (in line with the vacuum quenches in [64]), we
observed a linear entanglement growth at the highest achieved temperatures,
meaning that the imprint of mesons is fully suppressed. The latter observation
is consistent with a quasiparticle model interpretation that excludes the
presence of meson states, and with MPO simulations in the (massive) free
fermion regime. We similarly infer from these dynamical properties that the
meson melting process took place. For the particular theory under
consideration, the 2$D$ Ising QFT, we do not see conclusive indications to
decide whether this process is consistent with a sequential picture. Such a
feature in principle could be discernible from the temperature dependent meson
content in Fourier spectra of reflected entropies.

Our analyses comprehensively show that Rényi or reflected entropies can serve
as witnesses of the meson melting process both in static and dynamical
settings at finite temperature. While we focused on the relativistic Ising QFT
emerging in the IR limit of the corresponding spin chain model, it would be of
course highly desirable to analyse a similar setting also in more complicated
and higher-dimensional gauge theories. In view of the rapid progress in using
both TN and quantum simulations for this purpose (see [98, 99, 100, 101] for
an overview of current developments and future prospects for high-energy
physics), one should remain optimistic about prospects of studying meson
melting in these richer systems along the lines we developed in the present
article. In particular, some of our considered entropic quantities, like Rényi
entropies, became already measurable in quantum simulation experiments [102].
This opens the opportunity to study similar physical systems and test our
predictions beyond the classical computational realm. In a parallel vein, the
recent work [10] demonstrated the role of entanglement generation in
(1+1)-dimensional QED scattering events using TNs. Following these lines of
research, we are convinced that a quantum information based perspective can
provide relevant insights into many open questions related to mesonic systems
in particle and condensed matter physics.

###### Acknowledgements.

We would like to thank Burkhard Kämpfer and Rico Zöllner for useful
discussions and correspondence. We thank Alvise Bastianello for comments on
the manuscript and discussions on the topic. This work was partly supported by
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy – EXC-2111 – 390814868, and by the EU-QUANTERA
project TNiSQ (BA 6059/1-1). The work of JK was supported in part by a
fellowship from the Studienstiftung des deutschen Volkes (German Academic
Scholarship Foundation). JK is supported by the Israel Academy of Sciences and
Humanities & Council for Higher Education Excellence Fellowship Program for
International Postdoctoral Researchers.

##  Appendix A Review of QCD and holographic approaches

###  A.1 Theoretical methods in QCD

The meson operator entering formula (1) of the spectral function is
constructed as
$M(\vec{x},t)\equiv\bar{q}_{f}(\vec{x},t)\,\Gamma\,q_{f}(\vec{x},t)$ from the
Dirac spinor $q_{f}$ of a quark with flavor $f$. Here, $\Gamma$ is the vertex
operator - a combination of gamma matrices, which selects the spin and angular
momentum. See, e.g., [5] for a detailed selection of possible choices. In the
simplest cases, the identity $\mathds{1}$ or gamma matrix $\gamma^{\mu}$ can
be chosen. Real-time properties such as strength and form of field
fluctuations are captured by the time-ordered two-point function
$D(\vec{x},t,\vec{x}_{0},t_{0})\equiv\langle\mathcal{T}M(\vec{x},t)M^{\dagger}(\vec{x}_{0},t_{0})\rangle_{\beta}=\operatorname{Tr}\\{\rho_{\beta}\mathcal{T}M(\vec{x},t)M^{\dagger}(\vec{x}_{0},t_{0})\\}$,
where $\mathcal{T}$ is the the time-ordering operator and $\rho_{\beta}$ the
thermal density operator. This quantity can be expressed as a path-integral
over the Schwinger–Keldysh contour that contains both the initial conditions
(encoded in thermal density matrix elements) and the quantum dynamics
(enforced by the QCD action). The time-ordered two-point function can be
decomposed as follows

$\displaystyle D(\vec{x},\vec{x}_{0},t,t_{0})$
$\displaystyle=\frac{1}{2}\langle\\{M(\vec{x},t),M^{\dagger}(\vec{x}_{0},t_{0})\\}\rangle_{\beta}$
$\displaystyle\phantom{=}+\frac{1}{2}\operatorname{sign}(t-t_{0})\langle[M(\vec{x},t),M^{\dagger}(\vec{x}_{0},t_{0})]\rangle_{\beta}$
$\displaystyle\equiv
F(\vec{x},t,\vec{x}_{0},t_{0})-\frac{i}{2}\operatorname{sign}(t-t_{0})\rho(\vec{x},t,\vec{x}_{0},t_{0}).$
(22)

Here, we defined the statistical function $F$, which measures the population
of states, and the spectral function $\rho$. Meson peaks in $\rho$ depend on
the momentum $\vec{p}$, following the (relativistic) dispersion relation of
the QFT. Above the continuum threshold, given by twice the mass of the
lightest particle, unbound particle pair peaks can appear within a continuous
background structure. The thermal width of any such peak is proportional to
the inverse lifetime of the corresponding (quasi)particle. Another important
quantity in this context is the binding energy of a bound state, which is
given by the difference of its mass peak to the continuum threshold.

QCD based methods for the calculations of spectral functions have intrinsic
limitations, as we will concisely review here following the detailed
discussion in [5]. First, effective field theory approaches make use of the
fact that meson flavors of large mass can be integrated out to get a
nonrelativistic approximation of QCD. The resulting expansion in energy scales
is, however, truncated, limiting the physical content of the description.
Furthermore, lattice QCD allows a nonperturbative treatment on a discretized
spacetime lattice based on Monte Carlo sampling. However, the computations in
taking the continuum limit are numerically extremely expensive. Both these
methods give access to the Euclidean correlator $D_{E}(\vec{x},\tau)=\langle
M(\vec{x},-i\tau)M^{\dagger}(\vec{0},0)\rangle_{\beta}$, which is related in
Fourier space to the spectral function via a integral transformation. The
extraction of $\rho(\vec{p},\omega)$ from this relation for only a limited
number of discrete, uncertainty-prone lattice data points limits the
quantitative robustness of the results. An alternative approach is based on a
potential method. To define it, consider the real-time Wilson loop

$W_{\Box}(r,t)=\mathcal{P}\exp\left[ig\oint_{r\times
t}A^{\mu}\mathrm{d}z_{\mu}\right]$ (23)

along a rectangular path of size $r$ times $t$ in the spatial and time
direction. The static interquark potential $V_{s}$ is then defined as

$V_{s}(r)=\lim_{t\to\infty}\frac{i\partial_{t}W_{\Box}(r,t)}{W_{\Box}(r,t)}.$
(24)

This potential together with a mass and kinetic term defines a Hamiltonian,
which governs the time evolution of the two-point meson correlator via the
Schrödinger equation. The imaginary part of its solution in Fourier space
yields the spectral function. The accuracy of this method is limited through
the omittance of finite velocity corrections.

Overall, the presented methods of determining in-medium spectral functions for
quark bound states at finite temperature provide the phenomenological picture
of the meson melting process discussed in section II.1 in the main text.

###  A.2 Insights from holographic models

Holographic approaches based on the gauge/gravity duality provide a framework
to study high-temperature properties of quarkonium and to describe the melting
of mesons. The underlying principle is to study a gravitational theory in a
hyperbolic AdS spacetime with an additional holographic coordinate as the dual
to a QFT in Minkowski spacetime. Since this duality translates the strong
coupling regime of the QFT into a weakly-coupled gravity theory, the
considered problem can become computationally more tractable. Thermodynamic
and entanglement properties of gravitational solutions admitting black hole
horizons can then be identified with those of the gauge theory. (For general
introductions into this broad field and further applications we refer to,
e.g., [103, 104, 105, 106].) We briefly want to point out some major lessons
regarding meson melting from these explorations, based on the recent series of
papers [107, 108, 109, 110, 111, 112, 113, 114, 115], in which this topic was
thoroughly discussed.

The basic ingredient is the description of meson states as test particles in a
fixed asymptotic AdS spacetime. 999As test particles, different meson flavors
do not react back between themselves nor with the metric and the dilaton. The
underlying action

$S=\int\mathrm{d}^{4}x\mathrm{d}z\,\sqrt{-g}\operatorname{\mathrm{e}}^{-\phi}G(\phi)F_{mn}F^{mn}$
(25)

contains a dilaton field $\phi$, a flavor coupling function $G(\phi)$, and a
$U(1)$ gauge field $V_{m}$ in the definition
$F_{mn}=\partial_{m}V_{n}-\partial_{n}V_{m}$ of the field strength tensor. The
dilaton breaks conformal symmetry and models the running coupling in QCD. Upon
a field redefinition $\psi$ in Fourier space, the equation of motion of the
gauge field takes in tortoise coordinates 101010The tortoise coordinate $\xi$
follows from the blackening function $f$ and the holographic coordinate $z$ as
$\operatorname{\mathrm{d}\\!}\xi=\operatorname{\mathrm{d}\\!}z/f$. The mass
$m$ is the zero component $p_{0}$ of the four-momentum. the Schrödinger form

$\left(\partial_{\xi}^{2}-[U_{T}-m^{2}]\right)\psi=0$ (26)

with a temperature-dependent potential $U_{T}$. Normalizable solutions of this
equation yield a discrete mass spectrum $m_{n}^{2}$ ($n\in\mathds{N}$). The
original soft-wall model of [118] uses a quadratic ansatz (w.r.t. $z$) for the
dilaton field and the simplest choice $G\equiv 1$ to enforce a linear spectrum
of the form $m_{n}^{2}=c_{0}+c_{1}n$, which matches experimental data of
radial meson excitation spectra at $T=0$. At finite temperature, the potential
well in $U_{T}(z)$ is decreasing, implying the existence of lesser bound
states. In other words, the model is capable of reproducing the sequential
melting of meson states. Further studies in [119, 120, 121, 122, 123, 124,
125, 126] on refinements of this model found, however, unrealistic melting
temperatures much below the QCD deconfinement scale
$\mathcal{O}(150\,\mathrm{MeV})$. The authors of [107, 108, 109, 110, 111,
112, 113, 114] resolved that problem by finding one-parameter extensions of
the metric and dilaton as well a construction principle for the blackening
function that result in consistent melting temperatures, and, at the same
time, exhibit thermodynamic properties consistent with the equation of state
from lattice QCD. While this approach is based on ad-hoc ansätze for the
background, the latter can be determined also self-consistently, for example,
in a Einstein–dilaton model, defined by the action

$S=\frac{1}{16\pi
G_{\text{N}}^{(5)}}\int\mathrm{d}^{4}x\mathrm{d}z\,\left(R-\frac{1}{2}(\partial\phi)^{2}-V(\phi)\right),$
(27)

where the dilaton potential $V(\phi)$ governs the thermodynamic system
properties. It was found, that then only non-trivial flavor functions
$G(\phi)$ allow a consistent description of both meson trajectories (governed
by $U_{T}$) and thermodynamic functions (as a crossover deconfinement
transition). Based on a linear response framework for holographic settings,
one can calculate the spectral function according to eq. (1) from the boundary
asymptotics of the bulk field. In [114], this lead to the observation that the
temperature, at which a peak forms in the spectrum might be much larger than
the melting temperature following from the Schrödinger potential.

Overall, one can fairly conclude that holographic methods do not yet reach the
level to provide robust quantitative predictions of all QCD relevant features
(meson trajectories, thermodynamic quantities, spectral functions). As noted
also in [5], some of the holographic models, e.g. in the recent work [127],
predict the melting of spectral functions towards larger frequencies at high
temperatures, which is in direct contradiction to QCD results.

In contrast to QCD based methods, holographic approaches provide, however, the
advantage of giving access to entanglement measures. Starting with the work
[128], it was observed that holographic entanglement entropy, which according
to the seminal work [129, 130] is given as a minimal bulk surface area, can
serve as a probe (i.e. order parameter) of the confinement-deconfinement
transition. Further explorations of this quantity for different holographic
models mimicking QCD properties can be found in [131, 132, 133, 134, 135, 136,
137, 138, 139, 140, 141, 142]. However, no such studies so far analyzed this
phenomenon explicitly related to the meson melting process. Our TN simulations
in this article therefore provide the genuinely first study of entanglement
properties in this physical context.

##  Appendix B Entanglement structure and tensor network ansätze

The starting point of our considerations is a generic quantum many-body system
with $N$ elementary constituents of local physical dimension $d$ living on
discrete lattice sites. An arbitrary pure state vector $\ket{\Psi}$ in the
tensor product Hilbert space
$\mathcal{H}=\mathcal{H}_{1}\otimes\mathcal{H}_{1}\otimes\ldots\otimes\mathcal{H}_{N}$
of such a system can be written as

$\ket{\Psi}=\sum_{i_{1},i_{2},\ldots,i_{N}=1}^{d}\psi_{i_{1},i_{2},\ldots,i_{N}}\ket{i_{1}}\otimes\ket{i_{2}}\otimes\cdots\otimes\ket{i_{N}},$
(28)

where $\psi_{i_{1},i_{2},\ldots,i_{N}}$ are complex coefficients w.r.t. some
basis vectors $\ket{i_{1}},\ket{i_{2}},\ldots,\ket{i_{N}}$. The exponentially
large set of $d^{N}$ wave function coefficients
$\psi_{i_{1},i_{2},\ldots,i_{N}}$ easily exceeds classical computational
resources for an exact representation, which is thus limited to very small
systems. TN ansätze provide families of states with an efficient
representation, in which the coefficients of the wave function are
parametrized by a polynomial number of tensors (see e.g. [7, 8, 9, 143, 144]
for introductions to the topic).

In this work we make use of the matrix product state (MPS) and matrix product
operator (MPO) ansatz. For a system with $N$ sites of physical dimension $d$,
a MPO is an operator with the form

$\displaystyle O=\sum_{\begin{subarray}{c}i_{1},i_{2},\ldots,i_{N}=1\\\
j_{1},j_{2},\ldots,j_{N}=1\end{subarray}}^{d}\operatorname{Tr}$
$\displaystyle(A_{i_{1}j_{1}}^{(1)}A_{i_{2}j_{2}}^{(2)}\cdots
A_{i_{N}j_{N}}^{(N)})$ (29)
$\displaystyle\ket{i_{1}}\bra{j_{1}}\otimes\ket{i_{2}}\bra{j_{2}}\otimes\cdots\otimes\ket{i_{N}}\bra{j_{N}}.$

Here, the local rank-4 tensors $A_{\alpha\beta i_{k}j_{k}}^{(k)}$
$(k=1,\ldots,N)$ have two physical indices $i_{k},j_{k}$ (corresponding to the
bra and ket operation) and virtual indices $\alpha,\beta$ with bond dimension
$\chi$. For open boundary conditions, the edge tensors ($k=1,\,N$) reduce to
rank-3 tensors, instead. For fixed physical indices, the coefficient is a
product of matrices, hence the nomenclature. The ansatz (29) can be used to
represent density operators, but the positivity condition cannot be imposed
directly at the level of the individual tensors. It is however possible to
further restrict the form of the tensors to enforce positivity, by using a
purification. In particular, thermal states of local Hamiltonians admit an
efficient approximation of this form, as discussed in section III.

##  Appendix C Details of the transfer operator

The MPO approximation of the unit cell of a thermal state can be represented
graphically as a TN diagram as follows

$\vbox{\hbox{\includegraphics[width=164.77771pt]{figures/rho_beta.png}}}.$
(30)

Here, the two columns (blue and orange) indicate the repeating sites in the
2-site unit cell of the translational invariant chain. Each box represents a
tensor $A_{\alpha\beta i_{k}j_{k}}^{(k)}$ as in the MPO definition (29).
Virtual bond indices are drawn horizontally, physical ones vertically. Solid
boxes symbolize $\sqrt{\rho_{\beta}}$ and striped boxes respectively the
adjoint $\sqrt{\rho_{\beta}}^{\dagger}$. As discussed before, this
construction ensures positivity of the state. The thermal transfer operator is
then defined as

$\vbox{\hbox{\includegraphics[width=156.10345pt]{figures/transferop_col.png}}},$
(31)

where an overall trace over the physical indices is implied.

##  Appendix D Prony signal analysis

The basic idea of the Prony method [93] is to reconstruct a function $f(t)$ as
a sum of complex exponentials,

$f(t)=\sum_{k=1}^{K}c_{k}\operatorname{\mathrm{e}}^{-i\,\omega_{k}\,t},$ (32)

where $c_{k}$ and $\omega_{k}$ are respectively complex coefficients and
frequencies, and $K$ is the (variable) total number of modes. Prony based
methods make use of the linearity of the ansatz (32) to determine $c_{k}$ and
$\omega_{k}$ independently. In contrast to standard Fourier methods, this
provides the advantage of extracting also complex valued frequencies. Observe
in particular that real frequencies correspond to oscillations, while negative
(positive) imaginary contributions give rise to an exponential decay (growth).
In our previous work [88] we developed a signal analysis technique, which
performs a Prony analysis on a finite subrange of the (time) domain of the
variable $t$ and subsequently shifts the analysis window towards later times.
Each identified mode is then visualized as in Fig. 1 with a unique color in
the complex frequency plane (on a scale from blue to red). Discrete
frequencies will appear as stable modes across different analysis windows,
whereas streaks (a sequence of poles in different time windows) indicate
branch cuts.

## References

  * Coleman [2015] P. Coleman,  _Introduction to Many-Body Physics_ (Cambridge University Press, 2015)
  * Greiner, Schramm, and Stein [2007] W. Greiner, S. Schramm, and E. Stein,  _Quantum Chromodynamics_ (Springer-Verlag Berlin Heidelberg, 2007)
  * Busza, Rajagopal, and van der Schee [2018] W. Busza, K. Rajagopal, and W. van der Schee,  _Heavy Ion Collisions: The Big Picture, and the Big Questions,_ Ann. Rev. Nucl. Part. Sci. 68, 339 (2018), arXiv:1802.04801 [hep-ph]
  * Berges  _et al._ [2021] J. Berges, M. P. Heller, A. Mazeliauskas, and R. Venugopalan,  _QCD thermalization: Ab initio approaches and interdisciplinary connections,_ Rev. Mod. Phys. 93, 035003 (2021), arXiv:2005.12299 [hep-th]
  * Rothkopf [2020] A. Rothkopf,  _Heavy Quarkonium in Extreme Conditions,_ Phys. Rept. 858, 1 (2020), arXiv:1912.02253 [hep-ph]
  * Friman  _et al._ [2011] B. Friman, C. Hohne, J. Knoll, S. Leupold, J. Randrup, R. Rapp, and P. Senger,  _The CBM Physics Book: Compressed Baryonic Matter in Laboratory Experiments,_ Lect. Notes Phys. 814, 1 (2011)
  * Verstraete, Murg, and Cirac [2008] F. Verstraete, V. Murg, and J. Cirac,  _Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems,_ Advances in Physics 57, 143 (2008), arXiv:0907.2796 [quant-ph]
  * Schollwöck [2011] U. Schollwöck,  _The density-matrix renormalization group in the age of matrix product states,_ Annals of Physics 326, 96 (2011), arXiv:1008.3477 [cond-mat]
  * Orus [2014] R. Orus,  _A Practical Introduction to Tensor Networks: Matrix Product States and Projected Entangled Pair States,_ Annals Phys. 349, 117 (2014), arXiv:1306.2164 [cond-mat.str-el]
  * Rigobello  _et al._ [2021] M. Rigobello, S. Notarnicola, G. Magnifico, and S. Montangero,  _Entanglement generation in (1+1)D QED scattering processes,_ Phys. Rev. D 104, 114501 (2021), arXiv:2105.03445 [hep-lat]
  * Hauke  _et al._ [2012] P. Hauke, F. M. Cucchietti, L. Tagliacozzo, I. Deutsch, and M. Lewenstein,  _Can one trust quantum simulators?_ Rep. Prog. Phys. 75, 082401 (2012)
  * Cirac and Zoller [2012] J. I. Cirac and P. Zoller,  _Goals and opportunities in quantum simulation,_ Nat. Phys. 8, 264 (2012)
  * Alexeev  _et al._ [2021] Y. Alexeev _et al._ ,  _Quantum Computer Systems for Scientific Discovery,_ P. R. X. Quantum. 2, 017001 (2021), arXiv:1912.07577 [quant-ph]
  * Fraxanet, Salamon, and Lewenstein [2022] J. Fraxanet, T. Salamon, and M. Lewenstein,  _The Coming Decades of Quantum Simulation,_ (2022), arXiv:2204.08905 [quant-ph]
  * Tan  _et al._ [2021] W. L. Tan _et al._ ,  _Domain-wall confinement and dynamics in a quantum simulator,_ Nature Phys. 17, 742 (2021), arXiv:1912.11117 [quant-ph]
  * Vovrosh and Knolle [2021] J. Vovrosh and J. Knolle,  _Confinement and entanglement dynamics on a digital quantum computer,_ Scientific Reports 11, 11577 (2021), arXiv:2001.03044 [cond-mat.str-el]
  * Vovrosh  _et al._ [2021] J. Vovrosh, K. E. Khosla, S. Greenaway, C. Self, M. Kim, and J. Knolle,  _Simple mitigation of global depolarizing errors in quantum simulations,_ Phys. Rev. E 104, 035309 (2021), arXiv:2101.01690 [quant-ph]
  * Mildenberger  _et al._ [2022] J. Mildenberger, W. Mruczkiewicz, J. C. Halimeh, Z. Jiang, and P. Hauke,  _Probing confinement in a $\mathbb{Z}_{2}$ lattice gauge theory on a quantum computer,_ (2022), arXiv:2203.08905 [quant-ph]
  * 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,_ Phys. Rev. B 102, 014308 (2020), arXiv:1911.11382 [cond-mat.stat-mech]
  * Surace and Lerose [2021] F. M. Surace and A. Lerose,  _Scattering of mesons in quantum simulators,_ New J. Phys. 23, 062001 (2021), arXiv:2011.10583 [cond-mat.quant-gas]
  * 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 [cond-mat.quant-gas]
  * Knaute and Hauke [2022] J. Knaute and P. Hauke,  _Relativistic meson spectra on ion-trap quantum simulators,_ Phys. Rev. A 105, 022616 (2022), arXiv:2107.09071 [cond-mat.str-el]
  * Vovrosh  _et al._ [2022] J. Vovrosh, R. Mukherjee, A. Bastianello, and J. Knolle,  _Dynamical hadron formation in long-range interacting quantum spin chains,_ (2022), arXiv:2204.05641 [cond-mat.quant-gas]
  * Note [1] We refer to appendix A.1 for a short overview and discussion of them.
  * Matsui and Satz [1986] T. Matsui and H. Satz,  _$J/\psi$ Suppression by Quark-Gluon Plasma Formation,_ Phys. Lett. B 178, 416 (1986)
  * Karsch, Kharzeev, and Satz [2006] F. Karsch, D. Kharzeev, and H. Satz,  _Sequential charmonium dissociation,_ Phys. Lett. B 637, 75 (2006), arXiv:hep-ph/0512239
  * Laine  _et al._ [2007] M. Laine, O. Philipsen, P. Romatschke, and M. Tassler,  _Real-time static potential in hot QCD,_ JHEP 03, 054 (2007), arXiv:hep-ph/0611300
  * Kajimoto  _et al._ [2018] S. Kajimoto, Y. Akamatsu, M. Asakawa, and A. Rothkopf,  _Dynamical dissociation of quarkonia by wave function decoherence,_ Phys. Rev. D 97, 014003 (2018), arXiv:1705.03365 [nucl-th]
  * De Boni [2017] D. De Boni,  _Fate of in-medium heavy quarks via a Lindblad equation,_ JHEP 08, 064 (2017), arXiv:1705.03567 [hep-ph]
  * Rakovszky  _et al._ [2016] T. Rakovszky, M. Mestyán, M. Collura, M. Kormos, and G. Takács,  _Hamiltonian truncation approach to quenches in the Ising field theory,_ Nucl. Phys. B911, 805 (2016), arXiv:1607.01068 [cond-mat.stat-mech]
  * Hódsági, Kormos, and Takács [2018] K. Hódsági, M. Kormos, and G. Takács,  _Quench dynamics of the Ising field theory in a magnetic field,_ SciPost Phys. 5, 027 (2018), arXiv:1803.01158 [cond-mat.stat-mech]
  * 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,_ Phys. Rev. D 18, 1259 (1978)
  * Zamolodchikov [1989] A. B. Zamolodchikov,  _Integrals of Motion and S-Matrix of the (Scaled) T=T c Ising Model with Magnetic Field,_ Int. J. Mod. Phys. A4, 4235 (1989)
  * Delfino, Grinza, and Mussardo [2005] G. Delfino, P. Grinza, and G. Mussardo,  _Decay of particles above threshold in the Ising field theory with magnetic field,_ arXiv (2005), 10.1016/j.nuclphysb.2005.12.024, hep-th/0507133
  * Fonseca and Zamolodchikov [2006] P. Fonseca and A. Zamolodchikov,  _Ising spectroscopy. I. Mesons at $T<T_{c}$,_ (2006), arXiv:hep-th/0612304 [hep-th]
  * Zamolodchikov [2013] A. Zamolodchikov,  _Ising Spectroscopy II: Particles and poles at $T>T_{c}$,_ (2013), arXiv:1310.4821 [hep-th]
  * Cirac and Verstraete [2009] J. I. Cirac and F. Verstraete,  _Renormalization and tensor product states in spin chains and lattices,_ Journal of Physics A: Mathematical and Theoretical 42, 504004 (2009)
  * Vidal [2007] G. Vidal,  _Classical simulation of infinite-size quantum lattice systems in one spatial dimension,_ Phys. Rev. Lett. 98, 070201 (2007), arXiv:cond-mat/0605597
  * Hastings [2006] M. B. Hastings,  _Solving gapped hamiltonians locally,_ Phys. Rev. B 73, 085115 (2006)
  * Molnar  _et al._ [2015] A. Molnar, N. Schuch, F. Verstraete, and J. I. Cirac,  _Approximating gibbs states of local hamiltonians efficiently with projected entangled pair states,_ Phys. Rev. B 91, 045138 (2015)
  * Hatano and Suzuki [2005] N. Hatano and M. Suzuki,  _Finding Exponential Product Formulas of Higher Orders,_ Lect. Notes Phys. 679, 37 (2005), arXiv:math-ph/0506007
  * Note [2] Observe that in the limit $\alpha\searrow 1$, $S_{\alpha}(\rho_{A})$ reduces to $S(\rho_{A})$. In the continuum or in presence of gauge fields, the definition of reduced density matrices can lead to mathematical subtleties, which, however, will not be relevant for our subsequent spin chain calculations.
  * Dutta and Faulkner [2021] S. Dutta and T. Faulkner,  _A canonical purification for the entanglement wedge cross-section,_ JHEP 03, 178 (2021), arXiv:1905.00577 [hep-th]
  * Jeong, Kim, and Nishida [2019] H.-S. Jeong, K.-Y. Kim, and M. Nishida,  _Reflected Entropy and Entanglement Wedge Cross Section with the First Order Correction,_ JHEP 12, 170 (2019), arXiv:1909.02806 [hep-th]
  * Bao and Cheng [2019] N. Bao and N. Cheng,  _Multipartite Reflected Entropy,_ JHEP 10, 102 (2019), arXiv:1909.03154 [hep-th]
  * Chu, Qi, and Zhou [2020] J. Chu, R. Qi, and Y. Zhou,  _Generalizations of Reflected Entropy and the Holographic Dual,_ JHEP 03, 151 (2020), arXiv:1909.10456 [hep-th]
  * Bueno and Casini [2020a] P. Bueno and H. Casini,  _Reflected entropy, symmetries and free fermions,_ JHEP 05, 103 (2020a), arXiv:2003.09546 [hep-th]
  * Li, Chu, and Zhou [2020] T. Li, J. Chu, and Y. Zhou,  _Reflected Entropy for an Evaporating Black Hole,_ JHEP 11, 155 (2020), arXiv:2006.10846 [hep-th]
  * Bueno and Casini [2020b] P. Bueno and H. Casini,  _Reflected entropy for free scalars,_ JHEP 11, 148 (2020b), arXiv:2008.11373 [hep-th]
  * Camargo  _et al._ [2021] H. A. Camargo, L. Hackl, M. P. Heller, A. Jahn, and B. Windt,  _Long Distance Entanglement of Purification and Reflected Entropy in Conformal Field Theory,_ Phys. Rev. Lett. 127, 141604 (2021), arXiv:2102.00013 [hep-th]
  * Hayden, Parrikar, and Sorce [2021] P. Hayden, O. Parrikar, and J. Sorce,  _The Markov gap for geometric reflected entropy,_ JHEP 10, 047 (2021), arXiv:2107.00009 [hep-th]
  * Kudler-Flam, Kusuki, and Ryu [2020] J. Kudler-Flam, Y. Kusuki, and S. Ryu,  _Correlation measures and the entanglement wedge cross-section after quantum quenches in two-dimensional conformal field theories,_ JHEP 04, 074 (2020), arXiv:2001.05501 [hep-th]
  * Moosa [2020] M. Moosa,  _Time dependence of reflected entropy in rational and holographic conformal field theories,_ JHEP 05, 082 (2020), arXiv:2001.05969 [hep-th]
  * Kudler-Flam, Kusuki, and Ryu [2021] J. Kudler-Flam, Y. Kusuki, and S. Ryu,  _The quasi-particle picture and its breakdown after local quenches: mutual information, negativity, and reflected entropy,_ JHEP 03, 146 (2021), arXiv:2008.11266 [hep-th]
  * Jain and Mahapatra [2020] P. Jain and S. Mahapatra,  _Mixed state entanglement measures as probe for confinement,_ Phys. Rev. D 102, 126022 (2020), arXiv:2010.07702 [hep-th]
  * Terhal  _et al._ [2002] B. M. Terhal, M. Horodecki, D. W. Leung, and D. P. DiVincenzo,  _The entanglement of purification,_ Journal of Mathematical Physics 43, 4286–4298 (2002)
  * Takayanagi and Umemoto [2018] T. Takayanagi and K. Umemoto,  _Entanglement of purification through holographic duality,_ Nature Phys. 14, 573 (2018), arXiv:1708.09393 [hep-th]
  * Emonts and Kukuljan [2022] P. Emonts and I. Kukuljan,  _Reduced density matrix and entanglement in interacting quantum field theory with Hamiltonian truncation,_ (2022), arXiv:2202.11113 [quant-ph]
  * Yurov and Zamolodchikov [1991] V. P. Yurov and A. B. Zamolodchikov,  _Truncated fermionic space approach to the critical 2-D Ising model with magnetic field,_ Int. J. Mod. Phys. A 6, 4557 (1991)
  * Calabrese and Cardy [2005] P. Calabrese and J. L. Cardy,  _Evolution of entanglement entropy in one-dimensional systems,_ J. Stat. Mech. 0504, P04010 (2005), arXiv:cond-mat/0503393
  * Calabrese and Cardy [2006] P. Calabrese and J. L. Cardy,  _Time-dependence of correlation functions following a quantum quench,_ Phys. Rev. Lett. 96, 136801 (2006), arXiv:cond-mat/0601225
  * Calabrese and Cardy [2007] P. Calabrese and J. Cardy,  _Quantum Quenches in Extended Systems,_ J. Stat. Mech. 0706, P06008 (2007), arXiv:0704.1880 [cond-mat.stat-mech]
  * Mitra [2018] A. Mitra,  _Quantum quench dynamics,_ Annual Review of Condensed Matter Physics 9, 245 (2018), https://doi.org/10.1146/annurev-conmatphys-031016-025451
  * Kormos  _et al._ [2017] 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 (2017), arXiv:1604.03571 [cond-mat.stat-mech]
  * Castro-Alvaredo  _et al._ [2020] O. A. Castro-Alvaredo, M. Lencsés, I. M. Szécsényi, and J. Viti,  _Entanglement Oscillations near a Quantum Critical Point,_ Phys. Rev. Lett. 124, 230601 (2020), arXiv:2001.10007 [cond-mat.stat-mech]
  * Magnifico  _et al._ [2020] G. Magnifico, M. Dalmonte, P. Facchi, S. Pascazio, F. V. Pepe, and E. Ercolessi,  _Real Time Dynamics and Confinement in the $\mathbb{Z}_{n}$ Schwinger-Weyl lattice model for 1+1 QED,_ Quantum 4, 281 (2020), arXiv:1909.04821 [quant-ph]
  * Chanda  _et al._ [2020] T. Chanda, J. Zakrzewski, M. Lewenstein, and L. Tagliacozzo,  _Confinement and lack of thermalization after quenches in the bosonic Schwinger model,_ Phys. Rev. Lett. 124, 180602 (2020), arXiv:1909.12657 [cond-mat.stat-mech]
  * Lerose  _et al._ [2019] A. Lerose, F. M. Surace, P. P. Mazza, G. Perfetto, M. Collura, and A. Gambassi,  _Quasilocalized dynamics from confinement of quantum excitations,_ Phys. Rev. B 102, 041118 (2019), arXiv:1911.07877 [cond-mat.stat-mech]
  * Collura, Kormos, and TakÃ¡cs [2018] M. Collura, M. Kormos, and G. TakÃ¡cs,  _Dynamical manifestation of the Gibbs paradox after a quantum quench,_ Physical Review A 98 (2018), 10.1103/physreva.98.053610, arXiv:1801.05817 [cond-mat.stat-mech]
  * Pomponio, PristyÃ¡k, and TakÃ¡cs [2019] O. Pomponio, L. PristyÃ¡k, and G. TakÃ¡cs,  _Quasi-particle spectrum and entanglement generation after a quench in the quantum potts spin chain,_ Journal of Statistical Mechanics: Theory and Experiment 2019, 013104 (2019), arXiv:1810.05539 [cond-mat.stat-mech]
  * Leviatan  _et al._ [2017] E. Leviatan, F. Pollmann, J. H. Bardarson, and E. Altman,  _Quantum thermalization dynamics with Matrix-Product States,_ (2017), arXiv:1702.08894 [cond-mat.stat-mech]
  * Lencses, Pomponio, and Takacs [2020] M. Lencses, O. Pomponio, and G. Takacs,  _Relaxation and entropy generation after quenching quantum spin chains,_ SciPost Phys. 9, 011 (2020), arXiv:2004.09550 [cond-mat.stat-mech]
  * Birnkammer, Bastianello, and Knap [2022] S. Birnkammer, A. Bastianello, and M. Knap,  _Prethermalization in confined spin chains,_ (2022), arXiv:2202.12908 [cond-mat.stat-mech]
  * Alba and Calabrese [2017] V. Alba and P. Calabrese,  _Entanglement and thermodynamics after a quantum quench in integrable systems,_ Proceedings of the National Academy of Sciences 114, 7947 (2017), arXiv:1608.00614 [cond-mat.str-el]
  * Alba and Calabrese [2018] V. Alba and P. Calabrese,  _Entanglement dynamics after quantum quenches in generic integrable systems,_ SciPost Phys. 4, 017 (2018), arXiv:1712.07529 [cond-mat.stat-mech]
  * De Luca and Bastianello [2020] A. De Luca and A. Bastianello,  _Entanglement front generated by an impurity traveling in an isolated many-body quantum system,_ Phys. Rev. B 101, 085139 (2020), arXiv:1909.07341 [cond-mat.stat-mech]
  * Chapman  _et al._ [2019] S. Chapman, J. Eisert, L. Hackl, M. P. Heller, R. Jefferson, H. Marrochio, and R. C. Myers,  _Complexity and entanglement for thermofield double states,_ SciPost Phys. 6, 034 (2019), arXiv:1810.05151 [hep-th]
  * Bastianello and Collura [2020] A. Bastianello and M. Collura,  _Entanglement spreading and quasiparticle picture beyond the pair structure,_ SciPost Physics 8, 045 (2020), arXiv:2001.01671 [cond-mat.stat-mech]
  * Lieb and Robinson [1972] E. H. Lieb and D. W. Robinson,  _The finite group velocity of quantum spin systems,_ Commun. Math. Phys. 28, 251 (1972)
  * Eisert and Osborne [2006] J. Eisert and T. J. Osborne,  _General Entanglement Scaling Laws from Time Evolution,_ Phys. Rev. Lett. 97, 150404 (2006), arXiv:quant-ph/0603114
  * Scopa, Calabrese, and Bastianello [2021] S. Scopa, P. Calabrese, and A. Bastianello,  _Entanglement dynamics in confining spin chains,_ (2021), arXiv:2111.11483 [cond-mat.stat-mech]
  * Sotiriadis, Calabrese, and Cardy [2009] S. Sotiriadis, P. Calabrese, and J. Cardy,  _Quantum quench from a thermal initial state,_ EPL (Europhysics Letters) 87, 20002 (2009), arXiv:0903.0895 [cond-mat.stat-mech]
  * Zhang, Shen, and Liu [2011] J. M. Zhang, C. Shen, and W. M. Liu,  _Quantum quench dynamics of the bose-hubbard model at finite temperatures,_ Phys. Rev. A 83, 063622 (2011), arXiv:1103.1539 [cond-mat.stat-mech]
  * Bácsi and Dóra [2013] A. Bácsi and B. Dóra,  _Quantum quench in the luttinger model with finite temperature initial state,_ Phys. Rev. B 88, 155115 (2013), arXiv:1307.7582 [cond-mat.str-el]
  * Collura and Karevski [2014] M. Collura and D. Karevski,  _Quantum quench from a thermal tensor state: Boundary effects and generalized Gibbs ensemble,_ Phys. Rev. B 89, 214308 (2014), arXiv:1402.1944 [cond-mat.stat-mech]
  * Granet, Fagotti, and Essler [2020] E. Granet, M. Fagotti, and F. H. L. Essler,  _Finite temperature and quench dynamics in the Transverse Field Ising Model from form factor expansions,_ SciPost Phys. 9, 033 (2020), arXiv:2003.09014 [cond-mat.stat-mech]
  * Kohn and Santoro [2021] L. Kohn and G. E. Santoro,  _Quenching the anderson impurity model at finite temperature: Entanglement and bath dynamics using matrix product states,_ (2021), arXiv:2107.02807 [cond-mat.str-el]
  * Banuls  _et al._ [2020] M. C. Banuls, M. P. Heller, K. Jansen, J. Knaute, and V. Svensson,  _From spin chains to real-time thermal field theory using tensor networks,_ Phys. Rev. Res. 2, 033301 (2020), arXiv:1912.08836 [hep-th]
  * Knaute [2022] J. Knaute,  _Tensor Networks: From Holography to Quantum Field Theory_ , Ph.D. thesis, Freie U., Berlin (2022)
  * Note [3] See appendix C for details on the transfer operator.
  * Sachdev [2011] S. Sachdev,  _Quantum Phase Transitions_, 2nd ed. (Cambridge University Press, 2011)
  * Note [4] We refer to appendix D for more details on the method.
  * Peter [2014] T. Peter,  _Generalized Prony Method_ , Ph.D. thesis, University of Göttingen (2014), as of 08/03/2022 available online at na.math.uni-goettingen.de/pdf/Pe13.pdf
  * Note [5] Note that we always consider densities here since we work with the unit cell in the iMPO formalism. For our global quench setup, the relation (19) is then identically also valid for the whole system.
  * Note [6] Note that the effective temperatures $\beta^{*}J=2.00$ and $\beta^{*}J=0.50$ for type (1) mean that the difference to the initial $\beta J$ is smaller than our considered precision of two digits.
  * Note [7] This parameter value corresponds to the blue $n=0$ curves in Fig.+.1667em11.
  * Note [8] The initial temperature $\beta J$ follows from demanding (19), and $M_{h}^{pre}$ is uniquely determined from the specified ratio $M_{h}^{pre}/M_{g}^{pre}$.
  * Bañuls and Cichy [2020] M. C. Bañuls and K. Cichy,  _Review on novel methods for lattice gauge theories,_ Reports on Progress in Physics 83, 024401 (2020)
  * Bañuls  _et al._ [2020] M. C. Bañuls _et al._ ,  _Simulating Lattice Gauge Theories within Quantum Technologies,_ Eur. Phys. J. D 74, 165 (2020), arXiv:1911.00003 [quant-ph]
  * Bass and Zohar [2021] S. D. Bass and E. Zohar,  _Quantum technologies in particle physics,_ Phil. Trans. A. Math. Phys. Eng. Sci. 380, 20210072 (2021)
  * Davoudi  _et al._ [2022] C. W. B. Z. Davoudi _et al._ ,  _Quantum Simulation for High Energy Physics,_ (2022), arXiv:2204.03381 [quant-ph]
  * 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,  _Probing Rényi entanglement entropy via randomized measurements,_ Science 364, 260 (2019)
  * Papantonopoulos [2011] E. Papantonopoulos,  _From Gravity to Thermal Gauge Theories: The AdS/CFT Correspondence,_ Lect. Notes Phys. 828, 1 (2011)
  * Ammon and Erdmenger [2015] M. Ammon and J. Erdmenger,  _Gauge/Gravity Duality_ (Cambridge University Press, 2015)
  * Nastase [2015] H. Nastase,  _Introduction to the AdS/CFT Correspondence_ (Cambridge University Press, 2015)
  * Casalderrey-Solana  _et al._ [2014] J. Casalderrey-Solana, H. Liu, D. Mateos, K. Rajagopal, and U. A. Wiedemann,  _Gauge/String Duality, Hot QCD and Heavy Ion Collisions_ (Cambridge University Press, 2014)
  * Zöllner and Kämpfer [2016] R. Zöllner and B. Kämpfer,  _Holographically emulating sequential versus instantaneous disappearance of vector mesons in a hot environment,_ Phys. Rev. C 94, 045205 (2016), arXiv:1607.01512 [hep-ph]
  * Zöllner, Wunderlich, and Kämpfer [2017] R. Zöllner, F. Wunderlich, and B. Kämpfer,  _Extended soft-wall model for the QCD phase diagram,_ Acta Phys. Polon. Supp. B 10, 789 (2017), arXiv:1611.04124 [hep-th]
  * Zöllner and Kämpfer [2017a] R. Zöllner and B. Kämpfer,  _Extended soft wall model with background related to features of QCD thermodynamics,_ Eur. Phys. J. A 53, 139 (2017a), arXiv:1701.01398 [hep-ph]
  * Zöllner and Kämpfer [2017b] R. Zöllner and B. Kämpfer,  _Holography at QCD- $T_{\rm c}$,_ J. Phys. Conf. Ser. 878, 012023 (2017b), arXiv:1703.02958 [hep-ph]
  * Zöllner and Kämpfer [2018] R. Zöllner and B. Kämpfer,  _Holographic vector mesons in a dilaton background,_ J. Phys. Conf. Ser. 1024, 012003 (2018), arXiv:1708.05833 [hep-th]
  * Zöllner and Kämpfer [2020a] R. Zöllner and B. Kämpfer,  _Phase structures emerging from holography with Einstein gravity – dilaton models at finite temperature,_ Eur. Phys. J. Plus 135, 304 (2020a), arXiv:1807.04260 [hep-th]
  * Zöllner and Kämpfer [2020b] R. Zöllner and B. Kämpfer,  _Holographic vector meson melting in a thermal gravity-dilaton background related to QCD,_ Eur. Phys. J. ST 229, 3585 (2020b), arXiv:2002.07200 [hep-ph]
  * Zöllner and Kämpfer [2021a] R. Zöllner and B. Kämpfer,  _Quarkonia Formation in a Holographic Gravity–Dilaton Background Describing QCD Thermodynamics,_ Particles 4, 159 (2021a), arXiv:2007.14287 [hep-ph]
  * Zöllner and Kämpfer [2021b] R. Zöllner and B. Kämpfer,  _Holographic bottomonium formation in a cooling strong-interaction medium at finite baryon density,_ Phys. Rev. D 104, 106005 (2021b), arXiv:2109.05824 [hep-th]
  * Note [9] As test particles, different meson flavors do not react back between themselves nor with the metric and the dilaton.
  * Note [10] The tortoise coordinate $\xi$ follows from the blackening function $f$ and the holographic coordinate $z$ as $\operatorname{\mathrm{d}\\!}\xi=\operatorname{\mathrm{d}\\!}z/f$. The mass $m$ is the zero component $p_{0}$ of the four-momentum.
  * Karch  _et al._ [2006] A. Karch, E. Katz, D. T. Son, and M. A. Stephanov,  _Linear confinement and AdS/QCD,_ Phys. Rev. D 74, 015005 (2006), arXiv:hep-ph/0602229
  * Ghoroku and Yahiro [2006] K. Ghoroku and M. Yahiro,  _Holographic model for mesons at finite temperature,_ Phys. Rev. D 73, 125010 (2006), arXiv:hep-ph/0512289
  * Colangelo  _et al._ [2007] P. Colangelo, F. De Fazio, F. Jugeau, and S. Nicotri,  _On the light glueball spectrum in a holographic description of QCD,_ Phys. Lett. B 652, 73 (2007), arXiv:hep-ph/0703316
  * Colangelo  _et al._ [2008] P. Colangelo, F. De Fazio, F. Giannuzzi, F. Jugeau, and S. Nicotri,  _Light scalar mesons in the soft-wall model of AdS/QCD,_ Phys. Rev. D 78, 055009 (2008), arXiv:0807.1054 [hep-ph]
  * Fujita  _et al._ [2009] M. Fujita, K. Fukushima, T. Misumi, and M. Murata,  _Finite-temperature spectral function of the vector mesons in an AdS/QCD model,_ Phys. Rev. D 80, 035001 (2009), arXiv:0903.2316 [hep-ph]
  * Colangelo, Giannuzzi, and Nicotri [2009] P. Colangelo, F. Giannuzzi, and S. Nicotri,  _Holographic Approach to Finite Temperature QCD: The Case of Scalar Glueballs and Scalar Mesons,_ Phys. Rev. D 80, 094019 (2009), arXiv:0909.1534 [hep-ph]
  * Sui  _et al._ [2010] Y.-Q. Sui, Y.-L. Wu, Z.-F. Xie, and Y.-B. Yang,  _Prediction for the Mass Spectra of Resonance Mesons in the Soft-Wall AdS/QCD with a Modified 5D Metric,_ Phys. Rev. D 81, 014024 (2010), arXiv:0909.3887 [hep-ph]
  * Cui, Takeuchi, and Wu [2012] L.-X. Cui, S. Takeuchi, and Y.-L. Wu,  _Thermal Mass Spectra of Vector and Axial-Vector Mesons in Predictive Soft-Wall AdS/QCD Model,_ JHEP 04, 144 (2012), arXiv:1112.5923 [hep-ph]
  * Colangelo, Giannuzzi, and Nicotri [2012] P. Colangelo, F. Giannuzzi, and S. Nicotri,  _In-medium hadronic spectral functions through the soft-wall holographic model of QCD,_ JHEP 05, 076 (2012), arXiv:1201.1564 [hep-ph]
  * Martin Contreras, Diles, and Vega [2021] M. A. Martin Contreras, S. Diles, and A. Vega,  _Heavy quarkonia spectroscopy at zero and finite temperature in bottom-up AdS/QCD,_ Phys. Rev. D 103, 086008 (2021), arXiv:2101.06212 [hep-ph]
  * Klebanov, Kutasov, and Murugan [2008] I. R. Klebanov, D. Kutasov, and A. Murugan,  _Entanglement as a probe of confinement,_ Nucl. Phys. B 796, 274 (2008), arXiv:0709.2140 [hep-th]
  * Ryu and Takayanagi [2006a] S. Ryu and T. Takayanagi,  _Holographic Derivation of Entanglement Entropy from the anti-de Sitter Space/Conformal Field Theory Correspondence,_ Phys. Rev. Lett. 96, 181602 (2006a), arXiv:hep-th/0603001 [hep-th]
  * Ryu and Takayanagi [2006b] S. Ryu and T. Takayanagi,  _Aspects of holographic entanglement entropy,_ JHEP 08, 045 (2006b), arXiv:hep-th/0605073 [hep-th]
  * Lewkowycz [2012] A. Lewkowycz,  _Holographic entanglement entropy and confinement,_ JHEP 05, 032 (2012), arXiv:1204.0588 [hep-th]
  * Kim [2013] N. Kim,  _Holographic entanglement entropy of confining gauge theories with flavor,_ Phys. Lett. B 720, 232 (2013)
  * Kol  _et al._ [2014] U. Kol, C. Núñez, D. Schofield, J. Sonnenschein, and M. Warschawski,  _Confinement, phase transitions and non-locality in the entanglement entropy,_ JHEP 06, 005 (2014), arXiv:1403.2721 [hep-th]
  * Ghodrati [2015] M. Ghodrati,  _Schwinger effect and entanglement entropy in confining geometries,_ Phys. Rev. D 92, 065015 (2015), arXiv:1506.08557 [hep-th]
  * Zhang [2017] S.-J. Zhang,  _Holographic entanglement entropy close to crossover/phase transition in strongly coupled systems,_ Nucl. Phys. B 916, 304 (2017), arXiv:1608.03072 [hep-th]
  * Dudal and Mahapatra [2017] D. Dudal and S. Mahapatra,  _Confining gauge theories and holographic entanglement entropy with a magnetic field,_ JHEP 04, 031 (2017), arXiv:1612.06248 [hep-th]
  * Knaute and Kämpfer [2017] J. Knaute and B. Kämpfer,  _Holographic entanglement entropy in the QCD phase diagram with a critical point,_ Phys. Rev. D 96, 106003 (2017), arXiv:1706.02647 [hep-ph]
  * Ali-Akbari and Lezgi [2017] M. Ali-Akbari and M. Lezgi,  _Holographic QCD, entanglement entropy, and critical temperature,_ Phys. Rev. D 96, 086014 (2017), arXiv:1706.04335 [hep-th]
  * Dudal and Mahapatra [2018] D. Dudal and S. Mahapatra,  _Interplay between the holographic QCD phase diagram and entanglement entropy,_ JHEP 07, 120 (2018), arXiv:1805.02938 [hep-th]
  * Li, Xu, and Huang [2021] Z. Li, K. Xu, and M. Huang,  _The entanglement properties of holographic QCD model with a critical end point,_ Chin. Phys. C 45, 013116 (2021), arXiv:2002.08650 [hep-th]
  * Aref’eva, Patrushev, and Slepov [2020] I. Y. Aref’eva, A. Patrushev, and P. Slepov,  _Holographic entanglement entropy in anisotropic background with confinement-deconfinement phase transition,_ JHEP 07, 043 (2020), arXiv:2003.05847 [hep-th]
  * Lezgi, Ali-Akbari, and Asadi [2021] M. Lezgi, M. Ali-Akbari, and M. Asadi,  _Nonconformality, subregion complexity, and meson binding,_ Phys. Rev. D 104, 026001 (2021), arXiv:2011.11625 [hep-th]
  * Okunishi, Nishino, and Ueda [2022] K. Okunishi, T. Nishino, and H. Ueda,  _Developments in the Tensor Network – from Statistical Mechanics to Quantum Entanglement,_ J. Phys. Soc. Jap. 91, 062001 (2022), arXiv:2111.12223 [cond-mat.stat-mech]
  * Bañuls [2022] M. C. Bañuls,  _Tensor Network Algorithms: a Route Map,_ (2022), arXiv:2205.10345 [quant-ph]