# Uncovering Universal Entanglement Properties of Correlated Quantum Matter
in Quantum Simulation

Manoj K. Joshi  The first three coauthors contributed equally. Institute for
Quantum Optics and Quantum Information, Austrian Academy of Sciences,
Technikerstraße 21a, 6020 Innsbruck, Austria  Institut für Experimentalphysik,
Universität Innsbruck, Technikerstraße 25, 6020 Innsbruck, Austria  Christian
Kokail  The first three coauthors contributed equally. Institute for Quantum
Optics and Quantum Information, Austrian Academy of Sciences, Technikerstraße
21a, 6020 Innsbruck, Austria  Institute for Theoretical Physics,
Technikerstraße 21a, 6020 Innsbruck, Austria  Rick van Bijnen  The first three
coauthors contributed equally. Institute for Quantum Optics and Quantum
Information, Austrian Academy of Sciences, Technikerstraße 21a, 6020
Innsbruck, Austria  Institute for Theoretical Physics, Technikerstraße 21a,
6020 Innsbruck, Austria  Florian Kranzl  Institute for Quantum Optics and
Quantum Information, Austrian Academy of Sciences, Technikerstraße 21a, 6020
Innsbruck, Austria  Institut für Experimentalphysik, Universität Innsbruck,
Technikerstraße 25, 6020 Innsbruck, Austria  Torsten V. Zache  Institute for
Quantum Optics and Quantum Information, Austrian Academy of Sciences,
Technikerstraße 21a, 6020 Innsbruck, Austria  Institute for Theoretical
Physics, Technikerstraße 21a, 6020 Innsbruck, Austria  Rainer Blatt  Institute
for Quantum Optics and Quantum Information, Austrian Academy of Sciences,
Technikerstraße 21a, 6020 Innsbruck, Austria  Institut für Experimentalphysik,
Universität Innsbruck, Technikerstraße 25, 6020 Innsbruck, Austria  Christian
F. Roos  Institute for Quantum Optics and Quantum Information, Austrian
Academy of Sciences, Technikerstraße 21a, 6020 Innsbruck, Austria  Institut
für Experimentalphysik, Universität Innsbruck, Technikerstraße 25, 6020
Innsbruck, Austria  Peter Zoller  Institute for Quantum Optics and Quantum
Information, Austrian Academy of Sciences, Technikerstraße 21a, 6020
Innsbruck, Austria  Institute for Theoretical Physics, Technikerstraße 21a,
6020 Innsbruck, Austria

# Exploring the Entanglement Structure of Many-Body Wavefunctions
on a Quantum Simulator

Manoj K. Joshi  The first three coauthors contributed equally. Institute for
Quantum Optics and Quantum Information, Austrian Academy of Sciences,
Technikerstraße 21a, 6020 Innsbruck, Austria  Institut für Experimentalphysik,
Universität Innsbruck, Technikerstraße 25, 6020 Innsbruck, Austria  Christian
Kokail  The first three coauthors contributed equally. Institute for Quantum
Optics and Quantum Information, Austrian Academy of Sciences, Technikerstraße
21a, 6020 Innsbruck, Austria  Institute for Theoretical Physics,
Technikerstraße 21a, 6020 Innsbruck, Austria  Rick van Bijnen  The first three
coauthors contributed equally. Institute for Quantum Optics and Quantum
Information, Austrian Academy of Sciences, Technikerstraße 21a, 6020
Innsbruck, Austria  Institute for Theoretical Physics, Technikerstraße 21a,
6020 Innsbruck, Austria  Florian Kranzl  Institute for Quantum Optics and
Quantum Information, Austrian Academy of Sciences, Technikerstraße 21a, 6020
Innsbruck, Austria  Institut für Experimentalphysik, Universität Innsbruck,
Technikerstraße 25, 6020 Innsbruck, Austria  Torsten V. Zache  Institute for
Quantum Optics and Quantum Information, Austrian Academy of Sciences,
Technikerstraße 21a, 6020 Innsbruck, Austria  Institute for Theoretical
Physics, Technikerstraße 21a, 6020 Innsbruck, Austria  Rainer Blatt  Institute
for Quantum Optics and Quantum Information, Austrian Academy of Sciences,
Technikerstraße 21a, 6020 Innsbruck, Austria  Institut für Experimentalphysik,
Universität Innsbruck, Technikerstraße 25, 6020 Innsbruck, Austria  Christian
F. Roos  Institute for Quantum Optics and Quantum Information, Austrian
Academy of Sciences, Technikerstraße 21a, 6020 Innsbruck, Austria  Institut
für Experimentalphysik, Universität Innsbruck, Technikerstraße 25, 6020
Innsbruck, Austria  Peter Zoller  Institute for Quantum Optics and Quantum
Information, Austrian Academy of Sciences, Technikerstraße 21a, 6020
Innsbruck, Austria  Institute for Theoretical Physics, Technikerstraße 21a,
6020 Innsbruck, Austria

# Exploring Large-Scale Entanglement in Quantum Simulation

Manoj K. Joshi  The first three coauthors contributed equally. Institute for
Quantum Optics and Quantum Information, Austrian Academy of Sciences,
Technikerstraße 21a, 6020 Innsbruck, Austria  Institut für Experimentalphysik,
Universität Innsbruck, Technikerstraße 25, 6020 Innsbruck, Austria  Christian
Kokail  The first three coauthors contributed equally. Institute for Quantum
Optics and Quantum Information, Austrian Academy of Sciences, Technikerstraße
21a, 6020 Innsbruck, Austria  Institute for Theoretical Physics,
Technikerstraße 21a, 6020 Innsbruck, Austria  Rick van Bijnen  The first three
coauthors contributed equally. Institute for Quantum Optics and Quantum
Information, Austrian Academy of Sciences, Technikerstraße 21a, 6020
Innsbruck, Austria  Institute for Theoretical Physics, Technikerstraße 21a,
6020 Innsbruck, Austria  Florian Kranzl  Institute for Quantum Optics and
Quantum Information, Austrian Academy of Sciences, Technikerstraße 21a, 6020
Innsbruck, Austria  Institut für Experimentalphysik, Universität Innsbruck,
Technikerstraße 25, 6020 Innsbruck, Austria  Torsten V. Zache  Institute for
Quantum Optics and Quantum Information, Austrian Academy of Sciences,
Technikerstraße 21a, 6020 Innsbruck, Austria  Institute for Theoretical
Physics, Technikerstraße 21a, 6020 Innsbruck, Austria  Rainer Blatt  Institute
for Quantum Optics and Quantum Information, Austrian Academy of Sciences,
Technikerstraße 21a, 6020 Innsbruck, Austria  Institut für Experimentalphysik,
Universität Innsbruck, Technikerstraße 25, 6020 Innsbruck, Austria  Christian
F. Roos  Institute for Quantum Optics and Quantum Information, Austrian
Academy of Sciences, Technikerstraße 21a, 6020 Innsbruck, Austria  Institut
für Experimentalphysik, Universität Innsbruck, Technikerstraße 25, 6020
Innsbruck, Austria  Peter Zoller  Institute for Quantum Optics and Quantum
Information, Austrian Academy of Sciences, Technikerstraße 21a, 6020
Innsbruck, Austria  Institute for Theoretical Physics, Technikerstraße 21a,
6020 Innsbruck, Austria

###### Abstract

Entanglement is a distinguishing feature of quantum many-body systems, and
uncovering the entanglement structure for large particle numbers in quantum
simulation experiments is a fundamental challenge in quantum information
science. Here we perform experimental investigations of entanglement based on
the entanglement Hamiltonian, as an effective description of the reduced
density operator for large subsystems. We prepare ground and excited states of
a 1D XXZ Heisenberg chain on a 51-ion programmable quantum simulator and
perform sample-efficient ‘learning’ of the entanglement Hamiltonian for
subsystems of up to 20 lattice sites. Our experiments provide compelling
evidence for a local structure of the entanglement Hamiltonian. This
observation marks the first instance of confirming the fundamental predictions
of quantum field theory by Bisognano and Wichmann, adapted to lattice models
that represent correlated quantum matter. The reduced state takes the form of
a Gibbs ensemble, with a spatially-varying temperature profile as a signature
of entanglement. Our results also show the transition from area to volume-law
scaling of Von Neumann entanglement entropies from ground to excited states.
As we venture towards achieving quantum advantage, we anticipate that our
findings and methods have wide-ranging applicability to revealing and
understanding entanglement in many-body problems with local interactions
including higher spatial dimensions.

##  I Introduction

Figure 1:  Learning the entanglement structure of variationally prepared
ground and heated quantum many-body states. (a) We study the XXZ model with
Hamiltonian $\hat{H}$ given in Eq. (3). This defines our EH ansatz
$\tilde{H}_{A}(\bm{\beta})$ according to Eq. (2) with $\bm{\beta}$ the
entanglement temperature profile on a subsystem of length $L_{A}$ in the
51-ion chain. (b) Experimental procedure for state preparation and data
analysis. A variational quantum circuit (see Appendix G) first prepares
correlated quantum many-body states. In the second step, we collect
frequencies of bit strings sampled in different Pauli bases (see Appendix D).
The data are subjected to an entanglement Hamiltonian tomography (EHT)
procedure which finds the optimal EH $\tilde{H}_{A}(\bm{\beta})$ best
reproducing the experimentally obtained frequencies (see Appendix D). (c)
Entanglement properties of variationally prepared many-body states on 51 ions,
obtained via EHT for different anisotropies $\Delta=1$ and $\Delta=1.7$,
respectively. Blue squares show the results for the VQE ground states, while
red diamonds show data for heated states. The energy spectrum on the left
indicates VQE ground and excited state energies as fractions of the entire
spectral range. The upper panels show the von Neumann entropies as a function
of subsystem size $L_{A}$ which we obtain from the learned EH via
$S_{A}^{\text{VN}}=\braket{\tilde{H}_{A}(\bm{\beta})}+\log[Z_{A}(\bm{\beta})]$.
Dashed lines show the corresponding theoretical curves obtained from fitting
the ansatz $\rho_{A}(\bm{\beta})$ to an MPS simulation of the experiment.
Heated-up states indicate a clear volume-law scaling of the entanglement
entropy $S_{A}^{\text{VN}}\sim L_{A}$ as opposed to area-law scaling for the
VQE ground states $S_{A}^{\text{VN}}\sim\text{const}$. Lower panels depict the
optimal EH parameters for $L_{A}=12$ determined from EHT (see AppendixD).
Transparent lines show the results of all connected 12-site subsystems from
the central 27 ions of the 51-ion chain. Solid marked lines represent the mean
over all these subsystems.

Entanglement is the crucial ingredient that sets apart the quantum world from
its classical counterpart. It is a fundamental concept that has garnered
intense research interest due to its implications for various aspects of
quantum physics, from foundational aspects to quantum computation to condensed
matter systems and quantum chemistry Gross  _et al._ (2023). In quantum many-
body problems, entanglement leads to an exponential scaling of complexity with
system size. While classical simulations struggle to capture this complexity,
quantum simulation experiments have the ability to naturally represent large-
scale entanglement - being quantum systems themselves. Recent years have seen
tremendous progress in large-scale quantum simulation experiments touching
upon the boundaries of what is classically simulatable Joshi  _et al._ (2022);
Zhang  _et al._ (2017); Bohnet  _et al._ (2016); Chen  _et al._ (2023);
Semeghini  _et al._ (2021); Sompet  _et al._ (2022); Léonard  _et al._ (2022);
Zhang  _et al._ (2022). The quantum many-body systems in such experiments are
typically only locally interacting. That is, operators appearing in the system
Hamiltonian act only on local clusters of adjacent particles, with important
consequences for the eigenstates of the system and leading to universal
features of the entanglement structure contained in them.

Investigations of bipartite entanglement start with considering a partition of
the system of interest into a subsystem $A$ and its complement $\bar{A}$ (see
Fig. 1 (a)). For a system prepared in a many-body state $\ket{\Psi}$,
entanglement between the two can be quantified via the von Neumann
entanglement entropy (EE) $S_{A}^{\text{VN}}=-\text{Tr}(\rho_{A}\log\rho_{A})$
where $\rho_{A}=\text{Tr}_{\bar{A}}\ket{\Psi}\bra{\Psi}$ describes the reduced
density matrix of subsystem A. For many-body ground states of locally
interacting systems, one expects a sub-extensive _area-law_ scaling, where the
EE only grows with the size of the boundary $\partial A$ of the subsystem
Cirac  _et al._ (2021). This area law scaling lies at the heart of efficient
tensor-network approximations in classical simulations of many-body ground
states Cirac  _et al._ (2021); Ganahl  _et al._ (2023). In contrast, generic
excited states with energies well above the ground state will exhibit _volume-
law_ scaling, a growth of the EE with the subsystem size reminiscent of the
extensive behavior of thermodynamic entropy Miao and Barthel (2021). However,
extracting such information about entanglement from large-scale experiments
remains challenging, predominantly because of the difficulty of performing
tomography of $\rho_{A}$ on large subsystems Acharya  _et al._ (2019) for the
purpose of evaluating $S_{A}^{{\rm VN}}$. Moreover, identifying quantifiers
that capture the entanglement pattern, beyond a simple scalar value
$S_{A}^{\rm VN}$, and in relation to the geometry and topology of the
subsystem has proven to be a difficult task Amico  _et al._ (2008).

Here we engage these challenges in an experimental setting by considering the
Entanglement (or Modular) Hamiltonian (EH) $\tilde{H}_{A}$, which describes
the reduced density matrix of a subsystem $A$ through $\rho_{A}\sim
e^{-\tilde{H}_{A}}$. For ground states of many-body Hamiltonians with local
interactions, this entanglement Hamiltonian is conjectured to have a simple
operator structure. According to fundamental predictions in quantum field
theory (QFT) by Bisognano and Wichmann Bisognano and Wichmann (1975, 1976),
the entanglement Hamiltonian can be expressed as a spatial deformation of the
system Hamiltonian, i.e., it is composed of the same local operators appearing
in the translationally invariant system Hamiltonian, but acquiring a spatially
dependent prefactor (see also below, Eq. (2)). Rather than performing full
subsystem tomography of $\rho_{A}$, we instead ‘learn’ a local entanglement
Hamiltonian from experimental data. This allows us to study entanglement
properties of large subsystems, and experimentally investigate predictions
from QFT. Moreover, the entanglement Hamiltonian provides a unique insight
into entanglement patterns by offering an interpretation of $\hat{\rho}_{A}$
as a Gibbs state with a locally varying inverse temperature, or ’entanglement
temperature’, quantifying how subregions of the subsystem are entangled with
the outside world Dalmonte  _et al._ (2018); Giudici  _et al._ (2018); Zhu
_et al._ (2020).

In this work, we prepared ground and excited states of the 1D Heisenberg XXZ
model with $N=51$ spins in a trapped-ion platform, and extracted entanglement
Hamiltonians for subsystem sizes up to $L_{A}=20$ lattice sites. We find the
first experimental evidence for an entanglement Hamiltonian in the form of a
deformation of the system Hamiltonian, in line with a fundamental prediction
by Bisognano and Wichmann (BW) Bisognano and Wichmann (1975, 1976) and its
extension to conformal field theories (CFTs) Casini  _et al._ (2011); Cardy
and Tonni (2016). In addition, our results provide us with the Von Neumann
entanglement entropy displaying the area-to-volume law transition Miao and
Barthel (2021). We verify the accuracy of the fitted local entanglement
Hamiltonians from independent experimental data, i.e. we assign fidelities to
the results without having to resort to theoretical simulations.

##  II Locality of the Entanglement Hamiltonian

Our study primarily focuses on the entanglement Hamiltonian (EH). In the
following discussion, we highlight the remarkable observation that the EH
exhibits a local structure in the ground states of quantum field theories
(QFTs), as demonstrated in Bisognano and Wichmann (1975, 1976); Casini  _et
al._ (2011); Cardy and Tonni (2016); Dalmonte  _et al._ (2022). Moreover, this
fundamental finding has been extended to lattice models Dalmonte  _et al._
(2022), providing insights into the operator content of the EH in locally
interacting quantum many-body systems. This allows for efficient protocols to
comprehend the EH, which serves as the foundation for measuring entanglement
in large subsystems within our quantum simulation experiments Kokail  _et al._
(2021).

Relativistic quantum field theory (RQFT) in $d+1$ dimensional Minkowski space
makes predictions about the entanglement structure of the vacuum (ground)
state $\ket{\Omega}$. The Bisognano and Wichmann theorem Bisognano and
Wichmann (1975, 1976) states that the reduced density operator of the ground
state when partitioned into a semi-infinite half-space
$A=\\{\bm{x}\in\mathbb{R}^{d}\;|\;x_{1}>0\\}$ and its complement $\bar{A}$,
takes the form of a Gibbs state,

$\displaystyle\rho_{A}\equiv{\rm
Tr}_{\bar{A}}\ket{\Omega}\bra{\Omega}=\frac{1}{Z_{A}}e^{-\int_{A}{\sf
d}^{d}x\,\beta(\bm{x}){\mathscr{H}}(\bm{x})}\equiv\frac{1}{Z_{A}}e^{-\tilde{H}_{A}}$
(1)

Here, ${\mathscr{H}}(\bm{x})$ is the Hamiltonian density of the RQFT, and
$\beta(\bm{x})\sim x_{1}$ is the inverse temperature profile, that follows a
linear ramp with $x_{1}$ being the coordinate perpendicular to the plane
$\partial A$ cutting the infinite half-spaces. This result generalizes to
conformal field theories (CFTs) Casini  _et al._ (2011); Cardy and Tonni
(2016), predicting in particular that for a region $A$ as a ball of radius $R$
and radial coordinate $r=|\bm{x}|$, the inverse temperature profile has the
form of a parabola, $\beta(\bm{x})\sim(R^{2}-r^{2})/2R$ (see Appendix for
details).

Notably, Eq. (1) implies that the entanglement Hamiltonian $\tilde{H}_{A}$ is
local and shares the same operator structure as the system Hamiltonian but is
spatially deformed according to the profile $\beta(\bm{x})$ defining an
‘entanglement temperature’ $T(\bm{x})=1/\beta(\bm{x})$. This temperature
decreases with increasing distance from the cut, indicating that the dominant
contributions to the entanglement, due to low-lying eigenfunctions of
$\tilde{H}_{A}$, are supported close to the cut.

By adapting and applying BW arguments to the EH for ground states of strongly
interacting lattice models with local interactions, a conjecture has emerged
that generalizes Eq. (1) Dalmonte  _et al._ (2022). This conjecture has been
supported by both numerical and analytical investigations across various many-
body models. To elaborate, let’s consider a spatial deformation of a lattice
Hamiltonian:

$\displaystyle\hat{H}=\sum_{j}\hat{h}_{j}\quad\overset{\text{deform}}{\xrightarrow{\hskip
23.0pt}}\quad\tilde{H}_{A}=\sum_{j\in A}\beta_{j}\hat{h}_{j}+\dots,$ (2)

where $\hat{h}_{j}$ represents few-body terms that act on a neighborhood of
lattice sites $j$, and $\beta_{j}$ denotes the deformation profile. According
to the conjecture, the ground state of $\hat{H}$ gives rise to a reduced
density operator $\rho_{A}$ for a connected subsystem $A$, which assumes the
form described in Eq. (1), with $\tilde{H}_{A}$ derived from Eq. (2). It is
expected that this statement holds true for states that can be effectively
described by a continuum field theory, to which the BW theorem can be applied,
with minor non-universal corrections indicated by the dots.

In the current context, these arguments suggest a general operator structure
for the entanglement Hamiltonian that can be experimentally explored. By
parametrizing $\tilde{H}_{A}$ as in (2) and testing for potential deviations,
we employ a learning protocol called EHT Kokail  _et al._ (2021); Anshu  _et
al._ (2021), owing to the locality of $\tilde{H}_{A}$. In other words, the
number of samples required to learn the coefficients $\beta_{j}$ within a
given error scales polynomially with the number of terms $\hat{h}_{j}$.

Simultaneously, this procedure enables the direct measurement of the Von
Neumann entanglement entropy (EE), given by $S_{A}^{\rm VN}={\rm
Tr}(\rho_{A}\tilde{H}_{A})+\log Z_{A}$. Furthermore, this methodology can also
be applied to more general excited or thermal states. In the latter case, a
Gibbs state with a flat inverse temperature profile $\beta=\beta_{j}$ is
expected (with boundary corrections Garrison and Grover (2018)). Extracting
the EH for states at various energies, spanning from the lowest to the middle
of the spectrum of $\hat{H}$, facilitates the experimental observation of the
transition from an area law to a volume law for the EE.

##  III Experimental setup and model

A programmable trapped-ion quantum simulator serves as our experimental
platform for studying the entanglement structure in correlated quantum many-
body states. In our setup, a linear chain of $N=51\ ^{40}$Ca+ ions is held in
a linear Paul trap using highly anisotropic confining potentials. The spin
states are encoded into long-lived electronic states
$\ket{\downarrow}=\ket{S_{1/2},m=+1/2}$ and
$\ket{\uparrow}=\ket{D_{5/2},m=+5/2}$, defining the computational basis.
Global entangling operations are realized via quench dynamics of an XY model
with controllable long-range interactions, which is engineered via exploiting
the spin and motional degrees of freedom of the trapped ions (see Appendix B).
Our system allows spatially-resolved addressing and detection, enabling
arbitrary single-qubit rotations and high-fidelity readout (see Appendix C).

To study universal features of the entanglement structure on the trapped-ion
platform, we focus on realizing low-energy states of the Heisenberg XXZ model
with open boundary conditions

$\displaystyle\hat{H}_{\text{XXZ}}=J\sum_{j=1}^{N-1}\left(\hat{S}_{j}^{x}\hat{S}_{j+1}^{x}+\hat{S}_{j}^{y}\hat{S}_{j+1}^{y}+\Delta\hat{S}_{j}^{z}\hat{S}_{j+1}^{z}\right),$
(3)

where $\hat{S}_{j}^{\alpha}$ denote spin-1/2 operators acting on lattice sites
$j$. For an anisotropy parameter $-1<\Delta\leq 1$, the critical regime, this
model is gapless and its low-energy physics is described by a CFT with central
charge $c\\!=\\!1$ Calabrese and Cardy (2004). Preparing low-energy states in
this regime allows us not only to experimentally study lattice analogs of the
BW theorem on half partitions, but also finite bulk intervals. We further
present results analyzing the change of the temperature profile outside the
critical regime $\Delta>1$ as well as for highly excited states of the XXZ
chain.

##  IV Results

Figure 2: Entanglement temperature profiles for different subsystem sizes at
the critical point $\Delta=1$ (a) Energy distribution of VQE states calculated
from the mean and variance of $\hat{H}$ in the corresponding MPS wave
functions. (b) Local inverse temperatures $\beta_{j}$ for different subsystem
sizes in the bulk and at the boundary of the 51-ion chain obtained from EHT
(see Appendix). The lower panels (blue curves) show the results for the VQE
ground state up to $L_{A}=20$ sites. Temperature profiles up to $L_{A}=12$
sites are obtained from an EH ansatz with local fit parameters $\beta_{j}$
(blue squares) and operator components $\hat{h}_{j}$ as defined in Fig. 1 (a).
For $L_{A}>12$, we describe the temperature profile with a second-order
polynomial $\beta_{j}=q_{0}+q_{1}j+q_{2}j^{2}$, introducing 3 global fit
parameters $\\{q_{m}\\}_{m=0}^{2}$. The orange pentagons plotted in the lower
right panel of (b) shows the Uhlmann fidelity of the learned
$\rho_{A}(\bm{\beta})$ with respect to the corresponding $\rho_{A}$ from the
exact ground state for subsystems of size $L_{A}=7$ that we sweep through the
ion chain. The observed fidelity drop at the edges of the chain causes the
learned coefficients $\beta_{j}$ to deviate from the ones of the exact ground
state (red dash-dotted lines) in a boundary region. The coefficients
$\beta_{j}$ observed in the heated state (red diamonds) are consistent with
the uniform temperature profile of a thermal Gibbs state with a notable
decrease of the local inverse temperature compared to the VQE ground state.
(c) Verification of the EHT procedure via cross-fidelity estimation (see
Appendix F). Reduced density matrices of size $L_{A^{\prime}}=5$ are computed
from the learned $\rho_{A}(\bm{\beta})$ and cross-verified against independent
data taken from the experiment. Circles represent the maximum fidelity
$\mathcal{F}_{\text{max}}^{\text{sim}}$ with respect to theoretical
simulations. Error bars have been obtained via Jack-knifing and are smaller
than symbols if not shown.  Figure 3: Entanglement Hamiltonian of disjoint
5-site subsystems. (a) Local temperature profiles of disconnected subsystems
with separations $d_{12}=1$ (orange) and $d_{12}=8$ (blue). The inset
indicates the increased inverse temperatures $\beta_{j}$ of subsystems with
separations $d_{12}=1$ as opposed to subsystems with a separation $d_{12}=8$,
resulting in the decrease of mutual information as demonstrated in panel (b).
The inset further illustrates an asymmetry of the temperature profiles which
we interpret as a temperature gradient (as numerically observed in Sels and
Wouters (2017) in the context of free-fermion models). (b) Colormaps of the EH
parameters $\beta_{ij}$ for different subsystem separations. (c) Upper panel:
mutual information
$I_{AB}=S_{A}^{\text{VN}}+S_{B}^{\text{VN}}-S_{AB}^{\text{VN}}$ as a function
of subsystem separation $d_{12}$ together with theoretical predictions from
the MPS wave function. Lower panel: maximum fidelity with respect to an
independent data set as a function of subsystem separation $d_{12}$, averaged
over many subsystems. The orange lines depict the fidelity of an ansatz
$\tilde{H}_{A}(\bm{\beta})$ which includes cross-links $\beta_{ij}$ with $i\in
A$ and $j\in B$, while the blue line uses the standard ansatz for two
independent subsystems. For small separations, cross-links $\beta_{ij}$ lead
to a significant boost in fidelity. Error bars are smaller than symbols

Our main experimental procedure consists of variational state preparation
followed by Entanglement Hamiltonian Tomography as outlined in Fig. 1 (b). We
performed experimental preparations of approximate ground and excited states
for an XXZ model (cf. Eq. (3)) by optimizing quantum circuits generated
through quench dynamics in a variational quantum eigensolver (VQE) feedback
loop Kokail  _et al._ (2019). This resulted in variational states that
correspond to superpositions of eigenstates from finite energy windows within
the XXZ model (indicated by the blue and red brackets in the energy bar of
Fig. 1 (c)).

We successfully prepared states with an energy distance from the true ground
state that corresponds to approximately 2% of the entire spectral range (see
Appendix G). Furthermore, states with high average energy, located in the
middle of the spectrum, can be prepared efficiently by quenching the initial
state using the native interaction Hamiltonian of the ion chain (see Appendix
B), and subsequently applying the same VQE circuit utilized for preparing low-
energy state. We will refer to the states prepared using this method as ‘VQE
heated states’.

Subsequently, we analyze the entanglement properties of the prepared many-body
states via EHT, using experimental samples from 243 Pauli bases as base data.
In each of these Pauli bases, we collect 200 samples from quantum projective
measurements. EHT is performed from an ansatz of the EH of the form
$\tilde{H}_{A}(\bm{\beta})=\sum_{j}\beta_{j}\hat{h}_{j}$ with operator
components $\hat{h}_{j}$ defined in Fig. 1 (a). The EHT procedure is
independently verified via cross-fidelity check with respect to independent
experimental data sets and theoretical simulations (see below).

#### Area-law and volume-law scaling of the entanglement entropy.

Our main results are summarized in Fig. 1 (c), in which we analyze the
entanglement structure of VQE ground and heated states for different values of
the anisotropy $\Delta$. For subsystems $A$ in the bulk of the chain, we
observe the anticipated distinct behavior for the low-energy and excited
states. While the former exhibits an approximately constant EE, consistent
with an area law of entanglement, the heated-up states exhibit a
characteristic growth of the EE, which is consistent with a linear, volume
law, scaling $S_{A}\propto L_{A}$ with subsystem size $L_{A}=2,\dots,12$. This
behavior is intimately related to the characteristic shape of the EH
temperature profiles $\beta_{j}$, displayed in the lower panels of Fig. 1 (c).
The parabolic shape of the profiles in the VQE ground state results in spins
near the boundary of $A$ providing the dominant contribution to the
entanglement with the environment $\bar{A}$, thus capturing the essential
feature of area-law entanglement. In contrast, the profiles for the excited
states $\beta_{j}$ exhibit a relatively flat plateau within the bulk of the
subsystem, with only small differences observed between the boundary and bulk
spins. In either case, we find a smooth profile for the learned parameters,
consistent with expectations from CFT.

Furthermore, our results provide indications of a distinctive behavior between
temperature profiles of the VQE ground states in the critical and non-critical
regimes. As can be seen in the lower right panel of Fig. 1 (c) ($\Delta=1.7$),
the flanks of the profile $\beta_{j}$ exhibit an approximately constant slope,
to be contrasted to the parabolic profile in the critical regime ($\Delta=1$).
Our findings are consistent with the analytical results of free-particle
systems Eisler  _et al._ (2020), where EH parameters of non-critical chains
follow a triangular deformation. This suggests that the different shapes
(parabolic vs. triangular) of the EH parameters reflect the distinct
functional behavior in the decay of correlation functions (power-law vs.
exponential).

#### Scaling behavior of entanglement temperature profiles with $L_{A}$.

We now turn to analyzing the entanglement structure in more detail, by
studying the behavior of the EH as a function of subsystem size for bulk and
boundary regions at the critical anisotropy $\Delta=1$. In Fig. 2 (a), we
display normalized variance and mean energies for variationally prepared
ground and the excited states for which we summarise the results below. Fig. 2
(b) summarizes the resulting entanglement temperature profiles $\beta_{j}$,
encoding the information of subsystem density matrices of the 51-ion chain up
to $L_{A}=20$ sites. For the VQE ground state, the profiles exhibit a
parabolic shape for all subsystem sizes $L_{A}$ whose height grows
approximately linearly with $L_{A}$, i.e. individual spins that are close to
the subsystem’s edges contribute dominantly to the entanglement. Although the
original BW predictions are made for ground states, we find that even for the
approximate ground states, which are superpositions of states spanning a
finite energy range of the spectrum, the parabolic profile remains robust. In
Appendix H.2, we show numerically for the lowest 200 states that each
individual eigenstate exhibits a parabolically deformed EH. The heated-up
states however, are much higher in energy and here the inverse entanglement
temperature profile flattens and forms a plateau, whose height stays
approximately constant as a function of subsystem size, resulting in a linear
scaling of entanglement entropy (see Fig. 1 (c)) reminiscent of a locally
thermalized state.

The BW theorem predicts a linear slope of the temperature profile
$\beta_{j}\sim j$ for a bi-partition of an infinite system into two halves,
and we expect this profile to bend over in the presence of a boundary Kokail
_et al._ (2021). These expectations are confirmed by our EH learning procedure
for the VQE ground state as illustrated for subsystems at the edge, mimicking
a half-infinite subsystem, of the ion chain in Fig. 2 (b), in agreement with
exact ground state simulations (red dashed-dotted line). The profiles close to
the boundary of the ion chain, however, reveal that VQE state preparation is
less accurate there, which we attribute to finite-size artifacts due to the
low depth of our variational circuit. In order to quantify this effect, we
compute the average fidelity of 7-qubit reduced density matrices compared to
the corresponding subsystems in the exact ground state (see lower right panel
of Fig. 2 (b)), which shows a clear deviation from the exact ground state in
the boundary region. We note that the primary contribution to the density
matrix stems from spins located in close proximity to the entanglement cut,
which clarifies the slight drop in geometric mean fidelity (see Appendix F) at
the boundary.

#### Verification.

We perform direct experimental verification of the reconstructed density
matrices, by employing protocols similar to those described in Refs. Elben
_et al._ (2020, 2023). The procedure we employ works as follows. Having
obtained $\rho_{A}(\bm{\beta})$ for a given subsystem $A$, we further split
the subsystem $A$ into 2 subintervals $A_{1}$ and $A_{2}$, and compute reduced
density matrices for the subinterval $A_{1}$ via
$\rho_{A_{1}}(\bm{\beta})=\text{Tr}_{A_{2}}\left[\rho_{A}(\bm{\beta})\right].$
We then cross-verify the model $\rho_{A_{1}}(\bm{\beta})$ against an
independently taken data set via a Hilbert-Schmidt fidelity estimation (see
Elben  _et al._ (2020) and Appendix F for details) and average the resulting
fidelity over all connected subintervals $A_{1}$ of $A$. Fig. 2 (c),
summarizes the results of this verification procedure up to $L_{A}=12$ sites,
for subinterval sizes $N_{A_{1}}=5$. For both VQE ground state and heated
state, the geometric mean fidelity $\mathcal{F}_{\text{mean}}$ is
significantly above 90 %, while the max fidelity $\mathcal{F}_{\text{max}}$ in
the heated state approaches $\sim$85% (as defined in Appendix F).

Moreover, we perform the same verification procedure not only with respect to
independent data sets from the experiment but also to data from theoretical
simulations. To this end, we simulate the variational circuits using a time-
dependent variational principle with Matrix Product States (MPS) on 51 spins
and compute reduced density matrices of subintervals
$\rho_{A_{1}}^{\text{MPS}}$ from the MPS wave function, and with these, we
perform direct fidelity estimation with the experimental data. The results,
denoted as $\mathcal{F}_{\text{sim}}$ in Fig. 2 (c), demonstrate that the
maximum fidelities between $\rho_{A_{1}}(\bm{\beta})$ and
$\rho_{A_{1}}^{\text{MPS}}$ are consistent with the experimental fidelities
within the error bars.

#### Entanglement structure of disjoint subsystems.

So far, we have focused on a single connected subsystem and demonstrated the
universal applicability of the BW arguments. We now investigate _disconnected_
subsystems $A\cup B$ with regions $A$ and $B$ that are separated by a distance
$d_{12}$, where no universal predictions are available. We find that the
reduced density operator is well captured by an EH of the form

$\displaystyle\tilde{H}_{A\cup B}=\sum_{ij\in A\cup
B}\beta_{ij}\;\hat{\bm{S}}_{i}\cdot\hat{\bm{S}}_{j}\;,$ (4)

where
$\hat{\bm{S}}=\left(\hat{S}^{x},\hat{S}^{y},\sqrt{\Delta}\hat{S}^{z}\right)^{T}$
and $\beta_{ij}=\delta_{i,j-1}\beta_{i}$ whenever $i$ and $j$ are within the
same sub-subsystem.

We analyze disconnected subsystems for the $\Delta=1$ dataset. As shown in
Fig. 3 (a), the intra-subsystem profiles $\beta_{j}$ approach the expected
parabolic shape for large separations, indicating that $A$ and $B$ become
statistically independent. For short distances, these profiles are modified
and acquire an asymmetry, in qualitative agreement with analytic predictions
for specific CFTs Casini and Huerta (2009); Eisler  _et al._ (2022) and
reminiscent of an entropic force Sels and Wouters (2017). We further quantify
this effect by calculating the mutual information $I_{AB}$ shown in Fig. 3
(c), which increases for small $d_{12}$ and approaches a small constant value
for large $d_{12}$, in agreement with our theoretical simulations. The
correctness of the ansatz (4) is verified by again computing the fidelity
$\mathcal{F}_{\text{max}}$, which exceeds $90\%$, for all values of $d_{12}$.
Omitting the additional terms in the ansatz that connect the subsystems A and
B leads to markedly lower fidelities (see Fig. 3 (c). The values of all
necessary fit parameters $\beta_{ij}$ are shown in Fig. 3 (b), where the
emergence (vanishing) of inter-subsystem coupling with decreasing (increasing)
distance becomes apparent.

To summarize, despite the small size of the subsystems studied here, our
findings provide the first experimental evidence in favor of bi-local
generalization of Eq. (2), compatible with a CFT calculation for a massless
Dirac field Casini and Huerta (2009); Eisler  _et al._ (2022). We again
emphasize that a general prediction for the EH for disconnected subsystems is
presently not available. Applying our approach to different models with a
known effective CFT description can therefore help to improve our
understanding of entanglement properties for general CFTs.

##  V Conclusions and Outlook

The entanglement Hamiltonian provides a powerful tool to study entanglement in
correlated quantum matter governed by local Hamiltonians. In case the EH is
local, it is not only efficiently learnable from experimental data, but it
also provides a readily interpretable ’entanglement temperature’ profile
providing insights into the entanglement structure of the underlying quantum
many-body state. Our work presents the first experimental observation of a
local EH in strongly interacting lattice models, as an extension of
predictions from BW originally made in the context of ground states of QFTs.
Interestingly, we observe that the local structure of the EH is robust and
persists over a large range of low-energy states. We have studied low and high
energy states of the Heisenberg model with 51 spins in a trapped ion quantum
simulator. The local operator structure of the learned EHs is verified by
direct fidelity estimation from independent data, as well as numerical
simulations providing excellent agreement. Our methods also enable clear
observation of the transition from an area law of entanglement to a volume law
in excited states, with entanglement temperature profiles transitioning from
strongly deformed to near-uniform distributions.

We anticipate that the entanglement characteristics of ground states explained
by BW arguments apply to a broad class of many-body systems with local
Hamiltonians including higher spatial dimensions and fermionic systems.
Furthermore, the toolset used for measuring the operator structure of the
entanglement Hamiltonian can be applied to all present-day programmable
quantum simulation platforms. These advancements provide a framework for
exploring entanglement-related phenomena in experiments systematically. For
instance, the approach can be used to recognize topologically ordered phases
of matter through entanglement spectroscopy Zache  _et al._ (2022); Li and
Haldane (2008); Chandran  _et al._ (2014), or it can be used to test new
concepts providing further insights into entanglement structure Murciano  _et
al._ (2022); Kim  _et al._ (2022). Furthermore, the Ryu-Takayanagi conjecture
Ryu and Takayanagi (2006) quantitatively relates entanglement properties of
CFTs to the geometry of a dual gravitational theory, enabling the indirect
study of gravity through the holographic principle, where recent quantum
simulation experiments Periwal  _et al._ (2021) provide the necessary
programmability of interactions. A shift of focus to the EH as the central
object of study in investigations of entanglement in many-body systems thus
opens the door to a broad class of new physics to be explored on programmable
quantum simulators.

## Acknowledgements

CK and PZ thank Dries Sels for discussions. This project has received funding
from the European Union’s Horizon 2020 research and innovation programme under
grant agreement No 101113690 (PASQuanS2.1). CK, RvB, TVZ, and PZ were
supported by the US Air Force Office of Scientific Research (AFOSR) via IOE
Grant No. FA9550-19-1-7044 LASCEM, the Austrian Research Promotion Agency
(FFG) contract 884471 (ELQO, RvB), and by the Simons Collaboration on Ultra-
Quantum Matter, which is a grant from the Simons Foundation (651440, PZ). MJ,
FK, CFR, and RB acknowledge the financial support for the experiment from the
Austrian Science Fund through the SFB BeyondC (F7110), and the Institut für
Quanteninformation GmbH. This research was supported in part by the National
Science Foundation under Grant No. NSF PHY-1748958. The computational results
presented have been achieved (in part) using the HPC infrastructure LEO of the
University of Innsbruck. Simulations were performed using iTensor Fishman  _et
al._ (2022).

## Author contribution

MJ and FK developed and conducted the experiment under the guidance of RB and
CR. CK, RvB, TVZ, and PZ proposed the research and developed the quantum
protocols. CK, RvB, and MJ performed the data analysis. CK, RvB, TVZ, and PZ
wrote the manuscript, and MJ contributed texts on experimental setups. All
authors contributed to the discussion of the results.

##  Appendix A Characterizing bi-partite entanglement

The main text considers bipartite entanglement properties of the ground and
excited states. For a quantum system in a pure state $\ket{\Psi}$, the
entanglement properties with respect to a bipartition $A\\!:\\!\bar{A}$ are
specified by the Schmidt decomposition,
$\ket{\Psi}=\sum_{\alpha=1}^{\chi_{A}}\lambda_{\alpha}\ket{\Phi^{\alpha}_{A}}\otimes\ket{\Phi^{\alpha}_{\bar{A}}}$.
Here, $\lambda_{\alpha}$ are Schmidt coefficients, and the Schmidt rank
$\chi_{A}$ serves as a proxy of entanglement. Defining a reduced density
matrix ${\rho}_{A}={\rm Tr}\left[\ket{\Psi}\bra{\Psi}\right]$, and an EH
$\tilde{H}_{A}$ as in Eq. (1), we identify the Schmidt vectors
$\ket{\Phi^{\alpha}_{A}}$ with the eigenvectors of $\rho_{A}$ and
$\tilde{H}_{A}$. The entanglement spectrum (ES) is defined via
$\lambda_{A}^{2}=e^{-\xi_{\alpha}}$. Thus, knowledge of $\rho_{A}$, or
equivalently $\tilde{H}_{A}$ – as provided by sample-efficient learning of the
entanglement Hamiltonian in Appendix D – fully specifies bipartite
entanglement.

Area and volume law scaling are defined via the bipartite Von Neumann
entanglement entropy defined as $S^{\rm VN}_{A}\leavevmode\nobreak\
=\leavevmode\nobreak\ -{\rm Tr}\left(\rho_{A}\log\rho_{A}\right)$, or in terms
of the EH as $S_{A}^{\rm VN}\leavevmode\nobreak\ =\leavevmode\nobreak\ {\rm
Tr}(\rho_{A}\tilde{H}_{A})+\log Z_{A}$. Area law behavior, as is
characteristic for many-body ground states, is identified as the scaling
$S^{\rm VN}_{A}\propto L_{A}^{d-1}$ (or $\log L_{A}$ for $d=1$ at critical
points), where $L_{A}$ is the linear extent of the subsystem $A$ and $d$
denotes the number of spatial dimensions. In contrast, volume law scaling is
given by $S_{A}\propto V_{A}=L_{A}^{d}$, as expected, e.g. for a thermodynamic
entropy. For the $d=1$ Heisenberg model, area and volume law scaling are
observed for (approximate) ground and excited states, respectively, in Fig. 1
(c).

##  Appendix B Approximated power-law interactions

The experimental platform is discussed in previous references Kranzl  _et al._
(2022); Joshi  _et al._ (2022). In our trapped-ion quantum simulator, an
approximated power-law type spin-spin interaction is engineered by laser
fields driving the electronic and transverse motional degrees of freedom of
the trapped ions, allowing us to perform entangling steps in the variational
optimization loop. The interaction is described by

$\displaystyle\hat{H}_{\text{XY}}=\sum_{i,j>i}J_{ij}(\hat{\sigma}^{+}_{i}\hat{\sigma}^{-}_{j}+\hat{\sigma}^{-}_{i}\hat{\sigma}^{+}_{j}),$
(5)

with $J_{ij}\simeq J_{0}/|i-j|^{\alpha}$. The power-law exponent $\alpha$ is
controlled by adjusting the laser frequency of a three-tone laser field that
manipulates electronic and transverse motional degrees of freedom of the ion
chain Jurcevic  _et al._ (2014). For the current study, the exponent is tuned
to $\alpha\approx 0.82$.

In the numerical simulations, we used a spin-spin coupling matrix constructed
from the experimental measurements to compare the experiment and theory. In
the previous studies Joshi  _et al._ (2022); Kranzl  _et al._ (2022), the
influence of the third tone, which is used for compensating the light-induced
shifts on the spectator electronic levels, on the spin-spin coupling for long
ion strings was neglected. Here, we revise the calculation of the $J_{ij}$
matrix and compare the experimentally measured spin-spin coupling matrix with
the calculated one. For this calculation, we use the following expression
after taking all motional modes ($2N$) and their respective coupling to each
frequency component of the laser beam into account. Here, the coupling between
the $i^{\text{th}}$ and $j^{\text{th}}$ ion is described by

$J_{ij}\\!=\\!\frac{\hbar
k^{2}}{4m}\sum_{n=1}^{2N}M_{i}^{(n)}M_{j}^{(n)}\\!\left[\frac{\Omega_{i}^{\text{(B/R)}}\Omega_{j}^{\text{(B/R)}}}{\omega_{\text{B/R}}^{2}-\omega_{n}^{2}}+\frac{1}{2}\frac{\Omega_{i}^{{(\text{C})}}\Omega_{j}^{(\text{C})}}{\omega_{\text{C}}^{2}-\omega_{n}^{2}}\right],$
(6)

where $M_{i}^{(n)}$ is the normalised mode amplitude of the $i^{\text{th}}$
ion and the $n^{\text{th}}$ motional mode with mode frequency $\omega_{n}$
Porras and Cirac (2004). $\omega_{\text{R/B/C}}$ and $\Omega_{\text{R/B/C}}$
are the detunings and Rabi frequencies of the laser beam from the two-level
atomic transition. Here, subscripts R, B, and C denote three frequencies of
the laser beam that contains red-detuned, blue-detuned, and compensation
beams. $k$ is the wavenumber of the 729 nm laser beam and $m$ is the mass of
the calcium ion.

In our experiments, the laser is detuned by
$\pm(\omega_{\text{COM}}+2\pi\times 25\leavevmode\nobreak\ \text{kHz})$ from
the carrier transition, while the center-of-mass mode frequency is
$\omega_{\text{COM}}=2\pi\times 2.93$ MHz. Experimentally measured nearest-
neighbor terms of the engineered spin-spin interaction are shown in Fig. 4
(a), where solid lines are theory results. The engineered spin-spin
interaction is then used to experimentally examine the quench dynamics of a
single spin initialized to spin up in the middle of the ion chain while having
all others in the spin-down state. The flip-flop interaction coherently drives
the excitation to other locations of the ion chain while keeping the total
magnetization conserved. The experimental results are shown in Fig. 4 (b),
where solid lines are numerical results. The full numerically simulated
$J_{ij}$ matrix for the experimentally measured parameters is shown in Fig. 4
(c).

Figure 4:  Effective interactions and spin dynamics for a 51-ion chain (a)
Experimentally measured nearest-neighbor interaction terms $J_{i,i+1}$,
compared to theoretical calculations (solid line). (b) The quench dynamics of
a single spin initialized to a spin-up state in the middle of the ion chain,
while other spins initialized to a spin-down state, under the engineered flip-
flop type interaction plotted in discs and solid lines are numerical results.
(c) Theoretically calculated interaction matrix for the experimental
parameters.

##  Appendix C System characterization

In our experimental platform, the input state is prepared with a global
$X(\pi/2)$ gate over all ions, followed by an operation with a far-detuned
laser beam that is tightly addressed over even ion sites out of all 51 ions
while performing light shift gates. The unaddressed ions are prepared to spin
down, while addressed ions are prepared to spin up after applying another
global $X(-\pi/2)$ gate. The state preparation is subjected to inevitable
dephasing processes, thus a spin-echo sequence is employed while splitting the
sequential addressing operation into two parts to improve the state
preparation. The state preparation fidelity
$\mathcal{F}=|\braket{\downarrow\uparrow\downarrow\uparrow...}{\psi_{0}}|^{2}$
is measured to be 0.75(7) for the whole ion chain, which corresponds to a
single particle state preparation fidelity of $0.994(2)$. In the fidelity
estimation, the error bars account for fluctuations over different days of
measurement outcomes. An important remark: here, we perform direct
fluorescence measurements (i.e. the measurements on the $z$ basis) and detect
the individual ion on the EMCCD and compare the measured bit-string to the
ideal state to calculate the fidelity. From our independent measurements, the
detection error is measured to be smaller than $10^{-3}$ per particle, thus
the drop in the present fidelity can be assigned to the state preparation.

In addition to the state preparation error discussed above, there are also
measurement errors in our experimental system: here, we discuss the results of
two-qubit tomography performed for all nearest-neighbor pairs of the Néel
state (the input state for the variational optimization). The average fidelity
of the reconstructed two-qubit state is estimated to be
$\mathcal{F_{\text{two-qubit}}}=0.980(5)$, i.e. $\mathcal{F_{\text{single-
qubit}}}=0.989(5)$. This estimated fidelity also accounts for the state
preparation discussed previously. To examine the major source of error in the
measurements, we perform direct single-qubit tomography of the same input
state while utilizing only global $x$, $y$, and $z$ base measurements, in
contrast to the two-qubit case where local rotations are performed to account
for all 9 bases measurements. Here, the average single-qubit state fidelity is
estimated to be $0.994(3)$, which is higher than the estimates from the two-
qubit tomography reconstruction. In summary, these analyses imply that in our
system the leading errors, while carrying out measurements into different
bases, arise from the local rotations which are performed with the help of a
tightly focused beam.

##  Appendix D Entanglement Hamiltonian tomography and error mitigation

To reconstruct a reduced density matrix $\hat{\rho}_{A}$ from experimental
data, and to gain insight into its entanglement structure, we use the
procedure of entanglement Hamiltonian tomography (EHT) as introduced in Ref.
Kokail  _et al._ (2021). Here, the reduced density matrix is assumed to be of
the form

$\displaystyle\rho_{A}(\bm{\beta})=\exp\left(-\sum_{j\in
A}\beta_{j}\hat{h}_{j}\right)/Z_{A}(\bm{\beta}).$ (7)

The parameters $\beta_{j}$ are free variables to be fitted to the data, and
$Z(\bm{\beta})=\text{Tr}[\exp(-\sum_{j\in A}\beta_{j}\hat{h}_{j})]$ is a
constant that ensures trace normalization. The operators

$\displaystyle\hat{h}_{j}=\frac{J}{2}(\hat{S}^{+}_{j}\hat{S}^{-}_{j+1}+\mathrm{H.c.})+\Delta\hat{S}^{z}_{j}\hat{S}^{z}_{j+1}$
(8)

are associated with a link between two adjacent sites $j,j+1\in A$, and form a
decomposition of the Heisenberg Hamiltonian, i.e., such that
$\hat{H}=\sum_{j}\hat{h}_{j}$.

State preparation and measurement errors in the experiment can be accounted
for by applying a quantum operation consisting of a depolarizing map with rate
$p_{1}$, and spontaneous emission from $\ket{\uparrow}$ to $\ket{\downarrow}$
with rate $p_{2}$, of the form

$\displaystyle\mathcal{D}_{p}[\rho_{A}(\bm{\beta})]=\prod_{i=1}^{N_{A}}\mathcal{D}_{p}^{(i)}\rho_{A}(\bm{\beta}),$
(9)

where $\mathcal{D}_{p}^{(i)}$ is a quantum operation, applied to particle $i$
of a density matrix $\hat{\rho}$:

$\displaystyle\mathcal{D}_{p}^{(i)}\hat{\rho}=\sum_{k}E_{k}\hat{\rho}E_{k}^{\dagger},$
(10)

with the following $E_{k}$ acting on particle $i$:

$\displaystyle E_{0}$
$\displaystyle=\sqrt{1-3p_{1}/4}\ket{\downarrow}\bra{\downarrow}+\sqrt{1-3p_{1}/4-p_{2}}\ket{\uparrow}\bra{\uparrow}$
(11) $\displaystyle E_{1}$
$\displaystyle=\sqrt{\frac{p_{1}}{4}}\hat{\sigma}^{x},\ \
E_{2}=\sqrt{\frac{p_{1}}{4}}\hat{\sigma}^{y},\ \
E_{3}=\sqrt{\frac{p_{1}}{4}}\hat{\sigma}^{z},$ (12) $\displaystyle E_{4}$
$\displaystyle=\sqrt{p_{2}}\hat{\sigma}^{-}=\sqrt{p_{2}}\ket{\downarrow}\bra{\uparrow}.$
(13)

The rates $p_{1}$ and $p_{2}$ are calibrated through an analysis of the total
magnetization of the system after applying the state preparation circuit.
Since the circuit conserves the magnetization, any deviations from the
intended initial state magnetization can be attributed to the decoherence
channel, and the values of $p_{1}$ and $p_{2}$ can be determined from the
variance and mean of the magnetization.

Experimental data is collected as bit strings measured in the computational
(Z) basis after applying unitary basis rotations

$\displaystyle U^{(\bm{\alpha})}=\bigotimes_{j\in A}u_{j}^{(\alpha_{j})}$ (14)

where $\hat{u}_{j}^{(\alpha_{j})}$ are single particle basis rotations applied
at site $j$, rotating from the $\alpha_{j}=x,y,z$ basis to the $Z$ basis. We
use $3^{5}$ different measurement settings $\bm{\alpha}$, corresponding to a
tomographically complete basis set for all contiguous 5-site subsystems of the
51 ion chain (and in particular of subsystem A). While these measurements are
not tomographically complete for larger subsystems, they still yield enough
information to reliably fit the restricted ansatz Eq. (7). In each basis, we
take $10^{2}$ measurements, each measurement comprising a single bit string in
the computational basis.

The free parameters $\beta_{j}$ of the ansatz (7) are subsequently fitted to
the acquired data in a least-squares sense, i.e., minimizing a cost function

$\displaystyle\chi^{2}\\!=\\!\sum_{{\bm{\alpha}}}\sum_{\bm{s}}\left[P^{(\bm{\alpha})}_{\bm{s}}\\!-\\!\text{Tr}\left(\mathcal{D}_{p}[\rho_{A}(\bm{\beta})]U^{({\bm{\alpha}})}\ket{\bm{s}}\bra{\bm{s}}U^{({\bm{\alpha}})\dagger}\right)\right]^{2},$
(15)

where $P^{(\bm{\alpha})}_{\bm{s}}$ is the experimentally observed probability
for measuring bit string $\bm{s}$ after applying the basis transformation
${\bm{\alpha}}$. This cost function represents the difference between the
observed frequencies of occurrence of the bit strings, and the corresponding
expectation values from the density matrix ansatz
$\mathcal{D}_{p}[\rho_{A}(\bm{\beta})]$ for a particular choice of parameters
$\bm{\beta}$.

Finally, we note that the EHT procedure effectively provides a built-in method
for error mitigation, i.e., it is possible to filter out the effects of
decoherence. Namely, since the contributions of decoherence are explicitly
fitted in the form of the quantum operation $\mathcal{D}_{p}$, we are thus
able to isolate the coherent part $\rho_{A}({\bm{\beta}})$ from the
experimental data.

##  Appendix E Data post-processing

The XXZ-model studied in the main text exhibits a global
$\mathbb{Z}_{2}$-symmetry

$\displaystyle[\hat{H},\hat{\mathcal{P}}]=0\ \text{with}\
\hat{\mathcal{P}}=\bigotimes_{i}\hat{S}_{i}^{x}$ (16)

resulting in the situation that eigenstates with opposite total magnetization
$\ket{\Phi_{n}}$ and $\hat{\mathcal{P}}\ket{\Phi_{n}}$ are degenerate. The
variational circuits used in the experiment (see Appendix G) conserve total
magnetization, hence starting with an initial state $\ket{\Psi_{0}}$ of given
magnetization, the circuit is only able to prepare approximations to one of
the states $\ket{\Psi_{G}}$ or $\hat{\mathcal{P}}\ket{\Psi_{G}}$.

The original prediction of BW, however, is only valid for systems with a
_unique_ ground state. The ground state (GS) degeneracy of the XXZ model, on
the other hand, depends on the total number of sites $N$. For even $N$ the GS
is unique with magnetization $M=0$, while for odd $N$ the GS is two-fold
degenerate with $M=\pm 1$. In the limit $N\rightarrow\infty$, this distinction
vanishes in the sense that the corresponding reduced density matrices
$\rho_{A}$ converge to a single result. In order to test the BW prediction for
the experimental system with an _odd_ number of sites, we approximate a pure
ground state with mean magnetization given by the superposition
$|\Psi_{G}\rangle+\hat{P}|\Psi_{G}\rangle$, where $|\Psi_{G}\rangle$ is one of
the $M=\pm 1$ ground states. Explicitly, we calculate observables instead for
the mixture
$\ket{\Psi_{G}}\bra{\Psi_{G}}+\hat{\mathcal{P}}\ket{\Psi_{G}}\bra{\Psi_{G}}\hat{\mathcal{P}}$.
Independent of the experimental analysis, we have numerically confirmed that
this procedure converges to the correct expectation values in the bulk in the
limit $N\rightarrow\infty$.

##  Appendix F Verification

To verify the learning procedure, we determine a (mixed-state) fidelity
between the experimental quantum state under study, described by the density
matrix $\rho_{\rm exp}\equiv\rho_{1}$, and the reconstructed density matrix
from EHT, $\rho_{A}(\bm{\beta})\equiv\rho_{2}$. In particular, we analyze the
reconstructed density matrix in terms of two different Hilbert-Schmidt
fidelities defined in Liang  _et al._ (2019), given by the _maximum_ fidelity

$\displaystyle\mathcal{F}_{\textrm{max}}(\rho_{1},\rho_{2})=\frac{\text{Tr}(\rho_{1}\rho_{2})}{\max\\{\text{Tr}(\rho_{1}^{2}),\text{Tr}(\rho_{2}^{2})\\}},$
(17)

and the _geometric mean_ fidelity

$\displaystyle\mathcal{F}_{\textrm{mean}}(\rho_{1},\rho_{2})=\frac{\text{Tr}(\rho_{1}\rho_{2})}{\sqrt{\text{Tr}(\rho_{1}^{2})\text{Tr}(\rho_{2}^{2})}},$
(18)

which both measure the overlap between $\rho_{1}$ and $\rho_{2}$, normalized
by their purities.

As shown in Ref. Elben  _et al._ (2020), terms of the form
$\text{Tr}(\rho_{i}\rho_{j})$ for $i,j=1,2$, as occurring in Eqs. (17)-(18),
can be evaluated from outcomes of measurements performed in sufficiently many
measurement bases. Whereas Ref. Elben  _et al._ (2020) suggests the use of
randomized measurement bases Elben  _et al._ (2023), we use here the set of
tomographically complete Pauli measurements for all contiguous subsystems of
size 5, i.e. the same measurement bases used for the EHT protocol described
above.

Specifically, we denote by $P^{(1)}_{\bm{\alpha}}(\bm{s})$ the frequency of
having observed a particular bitstring $\bm{s}$ in the experiment (where
$\rho_{\rm exp}$ is realized and Pauli basis rotation $\hat{U}^{(\alpha)}$ has
been applied) and
$P^{(2)}_{\bm{\alpha}}(\bm{s})=\text{Tr}\left(\rho_{A}(\bm{\beta})U^{(\alpha)}\ket{\bm{s}}\bra{\bm{s}}U^{(\alpha)\dagger}\right)$,
i.e. the expectation value of observing bitstring $\bm{s}$ in the
reconstructed state $\rho_{A}(\bm{\beta})$. Then, we obtain the overlap
$\text{Tr}(\rho_{i}\rho_{j})$ for $i=1,j=2$ and purities
$\text{Tr}(\rho_{i}\rho_{j})$ for $i=j=1,2$ via Elben  _et al._ (2020)

$\displaystyle\text{Tr}(\rho_{i}\rho_{j})=\frac{2^{N_{A}}}{N_{\bm{\alpha}}}\sum_{\bm{\alpha}}\sum_{\bm{s},\bm{s}^{\prime}}(-2)^{-\mathcal{D}[\bm{s},\bm{s}^{\prime}]}{P_{{\bm{\alpha}}}^{(i)}(\bm{s})P_{{\bm{\alpha}}}^{(j)}(\bm{s}^{\prime})},$
(19)

where the Hamming distance $\mathcal{D}[\bm{s},\bm{s}^{\prime}]$ between two
strings $\bm{s}$ and $\bm{s}^{\prime}$ is defined as the number of entries
where $s_{k}\neq{s}^{\prime}_{k}$, i.e.
$\mathcal{D}[\bm{s},{\bm{s}}^{\prime}]\equiv\\#\left\\{k\in\\{1,\dots,N_{A}\\}\,|\,s_{k}\neq{s}^{\prime}_{k}\right\\}$
(see also Refs. Flammia and Liu (2011); da Silva  _et al._ (2011)).

Eq. (19) provides a direct experimental verification of the reconstructed
density matrix via the Hilbert-Schmidt fidelity, requiring no further theory
input such as simulations, and can be evaluated from the same type of
measurements employed for EHT. Importantly, however, the measurements used for
fidelity estimation should be independent from those used in EHT, to avoid
false correlations and biasing. We, therefore, split the total dataset for
each quantum state into two, where one-half of the data is used to reconstruct
$\rho_{A}(\bm{\beta})$, and the other half is subsequently used to evaluate
the fidelity (17). To evaluate the purity $\text{Tr}(\rho_{1}\rho_{1})$, the
verification data set is split once more into two sets to perform an unbiased
estimation. That is, in Eq. (19), we evaluate
$P_{{\bm{\alpha}}}^{(j)}(\bm{s})$ from one dataset, and
$P_{{\bm{\alpha}}}^{(j)}(\bm{s}^{\prime})$ from the other.

Finally, since the measurement basis set is only tomographically complete for
contiguous subsystems of size 5, and since we have taken only a limited number
of measurements, we have found that the fidelity estimation becomes inaccurate
for subsystems larger than 5. For subsystems of size $L_{A}>5$ we therefore
compute the averaged 5-site fidelity for all contiguous sub-subsystems
contained in the subsystem, as a measure of the fidelity of the total
subsystem.

##  Appendix G Variational state preparation

The experiment prepares variational quantum states Cerezo  _et al._ (2021) by
alternatingly applying two types of unitaries. The first type of operation
consists of an entangling operation of the form
$\hat{U}_{XY}(\theta)=\exp(-\mathrm{i}\theta\hat{H}_{XY})$, which applies for
a variable duration $\theta$ the native interaction Hamiltonian described in
Eq. (5). The second type of operation consists of single particle rotations,
applied to each second site, of the form
$\hat{U}_{Z}(\theta)=\exp\left(-\mathrm{i}\sum_{k=1}^{\lfloor
N/2\rfloor}\hat{\sigma}^{z}_{2k}\theta/2\right)$. Starting from an initial
Nèel state $\ket{\psi_{0}}=\ket{\downarrow\uparrow\downarrow\ldots}$, the
variational quantum states are thus of the form

$\displaystyle\ket{\Psi(\bm{\theta})}=\hat{U}_{Z}(\theta_{M})\hat{U}_{XY}(\theta_{M-1})\cdots\hat{U}_{Z}(\theta_{2})\hat{U}_{XY}(\theta_{1})\ket{\psi_{0}}.$
(20)

The parameters $\bm{\theta}$ of the variational state are optimized in a
feedback loop with a classical computer running an optimization algorithm that
attempts to minimize the expectation value of the energy of the state under
the Heisenberg Hamiltonian $\hat{H}$. From 30 measurements in each of the
$X,Y,Z$ bases, the quantity
$\bra{\Psi(\bm{\theta})}\hat{H}\ket{\Psi(\bm{\theta})}$ is estimated, serving
as a cost function for the classical optimizer. We use a variant of the SPSA
algorithm Kandala  _et al._ (2017), enhanced with a Gaussian process surrogate
model Rasmussen and Williams (2005). The search is warm-started with an
initial guess consisting of optimal parameters from a numerical exact
optimization for 13 particles. The variational optimization on the experiment
serves to refine these parameters. We reach energies, normalized to the total
spectral range, of $0.048(2)$ for $\Delta=1$ and $0.056(8)$ for $\Delta=1.7$.

For the ‘heated states’ discussed in the main text, we apply one additional
quench with the native entangling operation $\hat{H}_{XY}$. This quench is
applied to the initial state, before running the variational circuit with
optimal parameters. We found that quenches to the initial state are more
effective at heating the state than similar length quenches applied after
executing the circuit. The duration of the heating quenches is $1.87$ ms and
$1.5$ ms for $\Delta=1$ and for $\Delta=1.7$, respectively. For these
quenches, the measured energies on the quantum simulator are $0.37(1)$ and
$0.197(4)$.

We note that there also exist deterministic state preparation protocols for
eigenstates of the XXZ chain, making use of the integrability of the model, in
terms of so-called algebraic Bethe circuits Sopena  _et al._ (2022). Given the
requirement of a universal gate set to implement these circuits, we have
employed simpler short-depth variational circuits here.

##  Appendix H Bisognano-Wichmann theorem and extensions

###  H.1 The Entanglement Hamiltonian in Quantum Field Theories

Entanglement properties of quantum many-body systems can rarely be described
analytically. A notable exception is given by the Bisognano-Wichmann (BW)
theorem Bisognano and Wichmann (1975, 1976); Witten (2018), which applies to
the ground state $|\Omega\rangle$ of any relativistic quantum field theory
(RQFT) in $(d+1)$-dimensional Minkowski space. Given the underlying
Hamiltonian $\hat{H}=\int{\sf d}^{d}x\,\mathscr{H}(x)$ where
$\mathscr{H}(x)=T^{00}(x)$ is determined by the energy-momentum tensor
$T^{\mu\nu}(x)$, the reduced density operator as given in Eq. (1) of the main
text can be calculated exactly for the special case of a bi-partition of space
into two halves, $\mathbb{R}^{d}=A\cup B$ with
$A=\\{x=(x_{1},\dots,x_{d})\in\mathbb{R}^{d}|x_{1}>0\\}$. The corresponding EH
takes the form of a ‘deformation’ of $\hat{H}$,

$\displaystyle\tilde{H}_{A}=\int_{x\in A}{\sf
d}^{d}x\,\beta(x)\mathscr{H}(x)-F\;,$ (21)

with a linearly increasing ‘local inverse temperature’

$\displaystyle\beta(x)=\beta^{\text{BW}}(x)=2\pi x_{1},$ (22)

and a normalization constant $F$.

The classic result of BW can be extended when the QFT exhibits conformal
invariance Casini  _et al._ (2011); Cardy and Tonni (2016), which allows to
map the half-space into other regions $A\subset\mathbb{R}^{d}$. In particular,
for a solid sphere of radius $R$, i.e. $A=\\{x\in\mathbb{R}^{d}|x^{2}\leq
R^{2}\\}$, the EH is again local with

$\displaystyle\beta(x)=\beta^{\text{CFT}}(x)=2\pi\frac{(R^{2}-x^{2})}{2R}\;.$
(23)

We emphasize that the predictions in Eqs. (22) and (23) are valid in arbitrary
spatial dimensions $d$ and only require Lorentz or conformal invariance of the
QFT, respectively.

Translating these analytical results from the continuum to a spatial lattice
suggests approximate linear and parabolic deformations such as

$\displaystyle\beta^{\text{BW}}(x)\quad\rightarrow$
$\displaystyle\quad\beta^{\text{BW}}_{n}\sim n\;,$ $\displaystyle
n=0,1,2,\dots\;,$ (24) $\displaystyle\beta^{\text{BW}}(x)\quad\rightarrow$
$\displaystyle\quad\beta^{\text{CFT}}_{n}\sim n(N-n)\;,$ $\displaystyle
n=-N,\dots N$ (25)

written here for ${d}=1$ for simplicity [see also Eq. (2)].

For disconnected intervals as studied in the main text, no generally
applicable extension of the BW theorem is presently known. For reference, we
state here the result for a massless Dirac field on a line with $A=A_{+}\cup
A_{-}\subset\mathbb{R}$ and two intervals $A_{\pm}=(\pm a,\pm b)$ Casini and
Huerta (2009); Eisler  _et al._ (2022)

$\displaystyle\tilde{H}_{A}=\int_{x\in
A}dx\left[\beta_{\text{loc.}}(x)\mathscr{H}(x)+\beta_{\text{bi-
loc.}}(x)\mathscr{H}_{\text{bi-loc.}}(x,x_{c}(x))\right]\;.$ (26)

Here, $\mathscr{H}_{\text{bi-loc.}}(x,x_{c})$ is a bi-local operator that
connects every $x\in A_{\pm}$ with a conjugate point $x_{c}(x)=-ab/x\in
A_{\mp}$, and the spatial deformations take the form

$\displaystyle\beta_{\text{loc.}}(x)$
$\displaystyle=\frac{(b^{2}-x^{2})(x^{2}-a^{2})}{2(b-a)(ab+x^{2})}\;,$ (27)
$\displaystyle\beta_{\text{bi-loc.}}(x)$
$\displaystyle=\frac{ab}{x(ab+x^{2})}\beta_{\text{loc.}}(x)\;.$ (28)

Note that the local function $\beta_{\text{loc.}}(x)$ interpolates between the
expected parabolic shapes for a single subsystem when $A_{+}$ and $A_{-}$
touch ($b=R$ with $a\rightarrow 0$) and two independent subsystems when
$A_{+}$ and $A_{-}$ are far away ($b=a+2R$ with $a\rightarrow\infty$), while
the bi-local function $\beta_{\text{bi-loc.}}(x)$ vanishes in both limits.

###  H.2 Applicability to general low-energy states

Figure 5:  Numerical analysis of the Entanglement Hamiltonian of the lowest
200 eigenstates of a 51-site XXZ chain (a) Local temperature profiles of the
EH for subsystems of $L_{A}=9$ sites in the center of the spin chain for the
individual eigenstates. The inverse temperature profiles are colored according
to their eigenenergies, which are measured in units of the total spectral
range. (b) Verification of the local structure of the EH via computing Uhlmann
fidelities between the density matrices $\rho_{A}(\bm{\beta})$, reconstructed
from EHT, and the exact density matrices $\rho_{A}$ obtained via performing a
Schmidt decomposition on the individual eigenstates $\ket{\Phi_{k}}$ with
$k=1\dots 200$. GS denotes the temperature profile as well as the fidelity for
the ground state wave function.

Predictions about the spatial structure of the Entanglement Hamiltonian (EH)
from RQFT and Conformal Field Theory (CFT) apply to _ground states_ of local
theories (see main text). In the present paper, we study the EH for subsystems
of variationally prepared quantum many-body states on a trapped-ion quantum
simulator. We find that these states can be represented as superpositions of a
finite number of low-lying eigenstates of the target model. While numerical
verification of the local structure of the EH has been established for ground
states of quantum lattice models Dalmonte  _et al._ (2022), it is less clear
if such findings hold true for low-lying excited states or their
superpositions. In the following, we study the EH’s spatial profile in low-
lying eigenstates of a 51-site XXZ chain and verify its local structure using
fidelity estimations with respect to the exact density matrices.

Specifically, we compute excited states of the 51-site XXZ model using the
Density Matrix Renormalization Group (DMRG) with matrix product states (MPS).
In particular, for calculating the $k^{\text{th}}$ excited state, we modify
the XXZ Hamiltonian $\hat{H}$ by adding projectors on the $k-1$ pre-computed
eigenstates $\ket{\phi_{k}}$, with a constant weight factor $w$:

$\displaystyle\hat{H}\rightarrow\hat{H}+w\sum_{q<k}\ket{\phi_{q}}\bra{\phi_{q}}.$
(29)

Within this scheme, DMRG minimizes the energy of $\hat{H}$ while
simultaneously minimizing the overlap with all previously computed excited
states. Using this strategy, we compute MPS representations for the lowest 200
excited states of the XXZ model.

In Fig. 5 (a) we performed EHT for subsystems of $L_{A}=9$ sites for the
lowest 200 excited states. As can be seen, each of the individual eigenstates
exhibits a parabolically deformed EH, where the temperature profile for
higher-lying states acquires a flat plateau in the bulk of the subsystem.

In order to verify the validity of the local temperature profiles, in Fig. 5
(b) we compute the Uhlmann fidelity
$\mathcal{F}(\rho_{1},\rho_{2})=\left(\text{Tr}\sqrt{\sqrt{\rho_{1}}\rho_{2}\sqrt{\rho_{1}}}\right)^{2}$
between the density matrices $\rho_{A}(\bm{\beta})$, obtained from EHT, and
the exact density matrices
$\rho_{A}=\text{Tr}_{\bar{A}}\left(\ket{\Phi_{k}}\bra{\Phi_{k}}\right)$. While
the overall fidelities consistently exceed $\sim\,94\%$ for the lowest 200
eigenstates, Fig. 5 (b) reveals that the spread in fidelity increases for
higher-lying excited states. However, even at the highest energies
investigated, eigenstates can be found where the reduced density matrix can be
described by a local EH with up to $\sim$ 99% fidelity.

##  Appendix I Numerical justification of the data post-processing procedure

Figure 6:  Uhlmann fidelity $\mathcal{F}$, maximum fidelity
$\mathcal{F}_{\text{max}}$ and geometric-mean fidelity
$\mathcal{F}_{\text{mean}}$ (see Appendix F) of the ‘correct’ density matrix
$\rho_{A}^{\text{sym}}$ with respect to the modified density matrix
$\rho_{A}^{\text{mod}}$ as a function of total system size $L$, for a 5-site
subsystem from the center of the spin chain. The plot illustrates that
contributions to the density matrix $\rho_{A}^{\text{sym}}$ originating from
off-diagonal coherences in the state $\ket{\Psi^{\text{sym}}}$ vanish in the
limit $L\rightarrow\infty$.

As discussed in the main text, in our experiment we prepare approximate ground
states $\ket{\Psi_{G}}$ of the XXZ model for $N=51$ sites in different
parameter regimes. These states have good quantum numbers with respect to
total magnetization $\hat{M}=\sum_{j}\hat{S}_{j}^{z}$ but are not eigenstates
of the global $\mathbb{Z}_{2}$ symmetry operator
$\hat{\mathcal{P}}=\bigotimes_{j}\hat{S}_{j}^{x}$. This is because our
variational circuits prepare states with a fixed finite magnetization $M\neq
0$, while the global $\mathbb{Z}_{2}$ operation changes the magnetization $M$
to $-M$. In Appendix H, we discuss a data post-processing technique that
effectively simulates the presence of a state of the form
$\ket{\Psi_{G}^{\text{sym}}}=\frac{1}{\sqrt{2}}\left(\ket{\Psi_{G}}+\hat{\mathcal{P}}\ket{\Psi_{G}}\right)$
in experimental measurements, restoring the global $\mathbb{Z}_{2}$ symmetry
and facilitating the investigation of the BW predictions through an analysis
of a unique ground state (see Methods). This is achieved via an approximation,
modifying the experimentally obtained sample data to make them appear as if
they originate from a state of the form

$\displaystyle\rho_{A}^{\text{mod}}=\frac{1}{2}\left(\rho_{A}+\hat{\mathcal{P}}_{A}\rho_{A}\hat{\mathcal{P}}_{A}\right),$
(30)

with $\hat{\mathcal{P}}_{A}=\bigotimes_{j\in A}\hat{S}_{j}^{x}$ and $\rho_{A}$
the reduced density matrix of a subsystem of $\ket{\Psi_{G}}$. In the
following, we justify this procedure via numerical simulations and show that
it converges to the ‘correct’ density matrix

$\displaystyle\rho_{A}^{\text{sym}}=\text{Tr}_{\bar{A}}\left(\ket{\Psi_{G}^{\text{sym}}}\bra{\Psi_{G}^{\text{sym}}}\right),$
(31)

in the limit $L\rightarrow\infty$, where $L$ denotes the total number of
particles.

To verify the procedure, we compute $\rho_{A}^{\text{sym}}$ of a bulk region
of the superposition state $\ket{\Psi_{G}^{\text{sym}}}$ and compare it to the
modified density matrix $\rho_{A}^{\text{mod}}$ of Eq. (30). Expanding the
reduced density matrices of a subsystem $A$ in a basis of Pauli strings
$\hat{\sigma}_{j}^{\alpha_{j}}\hat{\sigma}_{j+1}^{\alpha_{j+1}}\dots\hat{\sigma}_{j+L_{A}-1}^{\alpha_{j+L_{A}-1}}$,
we can find an analytical expression for the difference between
$\rho_{A}^{\text{sym}}$ and $\rho_{A}^{\text{mod}}$, given by

$\displaystyle\begin{split}\rho_{A}^{\text{sym}}-\rho_{A}^{\text{mod}}&=\sum_{\\{\alpha_{j}\\}}\bra{\Psi_{G}}\\{\hat{\sigma}_{j}^{\alpha_{j}}\hat{\sigma}_{j+1}^{\alpha_{j+1}}\dots\hat{\sigma}_{j+L_{A}-1}^{\alpha_{j+L_{A}-1}},\hat{\mathcal{P}}\\}\ket{\Psi_{G}}\\\
&\times\hat{\sigma}_{j}^{\alpha_{j}}\hat{\sigma}_{j+1}^{\alpha_{j+1}}\dots\hat{\sigma}_{j+L_{A}-1}^{\alpha_{j+L_{A}-1}}\end{split}$
(32)

where $\\{\cdot,\cdot\\}$ denotes the anticommutator and
$\hat{\sigma}_{j}^{\alpha_{j}}\in[\hat{\mathbb{1}}_{j},\hat{\sigma}_{j}^{x},\hat{\sigma}_{j}^{y},\hat{\sigma}_{j}^{z}]$.
The expectation values appearing in Eq. (32) contain Pauli strings of the form
$\hat{\sigma}_{1}^{x}\cdots\hat{\sigma}_{j-1}^{x}\hat{\sigma}_{j}^{\alpha_{j}}\hat{\sigma}_{j+1}^{\alpha_{j+1}}\dots\hat{\sigma}_{j+L_{A}-1}^{\alpha_{j+L_{A}-1}}\hat{\sigma}_{j+L_{A}}^{x}\cdots\hat{\sigma}_{L}^{x}$
which span the whole spin chain. Intuitively, expectation values of large
Pauli strings vanish for $L\gg$, which we numerically verify in Fig. 6.

In particular, Fig. 6 illustrates different density matrix fidelities (see
Appendix F) between $\rho_{A}^{\text{sym}}$ and $\rho_{A}^{\text{mod}}$ for
density matrices computed for a subsystem of size $L_{A}=5$ from the center of
the spin chain as a function of the total system size $L$. As can be seen,
$\rho_{A}^{\text{sym}}$ and $\rho_{A}^{\text{mod}}$ become approximately
identical at $L\sim 10^{3}$, indicating the vanishing of expectation values of
large Pauli stings as $L\rightarrow\infty$. Contributions to the density
matrix $\rho_{A}^{\text{sym}}$ that stem from off-diagonal coherences in
$\ket{\Psi^{\text{sym}}}$ thus vanish for $L\rightarrow\infty$. This allows us
to replace the coherent superposition state $\ket{\Psi^{\text{sym}}}$ by an
incoherent mixture
$\ket{\Psi_{G}}\bra{\Psi_{G}}+\hat{\mathcal{P}}\ket{\Psi_{G}}\bra{\Psi_{G}}\hat{\mathcal{P}}$
for large system sizes, for which reduced density matrices take the form of
Eq. (30).

##  Appendix J Spatial structure of entanglement eigenstates

Figure 7: Energy density of entanglement eigenstates in comparison to the
ground state. (a) Energy density in the ground state of a 100-site XXZ model
at the isotropic point $J=\Delta=1$. (b) Absolute difference between energy
densities of the ground state wave function and entanglement eigenstates of a
10-site subsystem from the bulk of the chain. Note that the excess energy in
the dominant Schmidt vectors $\ket{\Phi_{A}^{\alpha}}$ with $\alpha=0,1,2$
with respect to the ground state $\ket{\Psi_{G}}$ is concentrated close to the
entanglement cut.

In this section we provide further evidence for the interpretation of the
deformation parameters $\beta_{j}$ in the EH as a “local inverse
(entanglement) temperature”. To be explicit, we consider again the EH
$\tilde{H}_{A}=\sum_{j\in A}\beta_{j}h_{j}$ with a finite subsystem $A$ in the
bulk for the ground state $\ket{\Psi_{G}}$ of
$\hat{H}_{\text{XXZ}}=\sum_{j}h_{j}$ as in the main text. The fact that
$\rho_{A}\propto e^{-\tilde{H}_{A}}$ is a Gibbs state with $\beta_{j}$ small
close to the entanglement cut and larger towards the middle of the subsystem
intuitively suggests that the “dominant contributions to the entanglement live
close to the boundary”.

To make this intuition more precise, we write

$\displaystyle\rho_{A}=\frac{1}{Z_{A}}e^{-\tilde{H}_{A}}=\sum_{\alpha}e^{-\xi_{\alpha}}\ket{\Phi_{A}^{\alpha}}\bra{\Phi_{A}^{\alpha}}$
(33)

with $\xi_{\alpha}$ and $\ket{\Phi_{A}^{\alpha}}$ the eigenvalues and -states
of $\tilde{H}_{A}$, i.e. the entanglement spectrum and the Schmidt vectors,
respectively. The dominant contribution to the entanglement is thus attributed
to the Schmidt vectors with the smallest eigenvalues. We expect that these
vectors carry excitations that are dominantly supported close to the
entanglement cut.

We reveal this anticipated spatial structure by calculating the average energy
density $\langle\hat{h}_{j}\rangle$ and $j\in A$ (see Fig. 7). In Fig. 7a, we
show the energy density for $\ket{\Psi_{G}}$, which exhibits a homogeneous
profile (up to a 2-site unit cell effect) as expected for a translationally
invariant system. In contrast, the energy density in the low-lying Schmidt
vectors is strongly inhomogeneous as demonstrated in Fig. 7 (b), where we plot
the absolute difference in energy density in $\ket{\Phi_{A}^{\alpha}}$ w.r.t.
$\ket{\Psi_{G}}$. We interpret the fact that this difference is largest close
the entanglement cut as localized excitations in the dominant Schmidt vectors
which are supported in this boundary region.

This spatial “localization” of entanglement finds its most prominent
manifestation in the context of topological order. As originally pointed out
by Li and Haldane Li and Haldane (2008), the low-lying entanglement spectrum
for a subsystem of topologically ordered state carries a fingerprint of the
structure associated to the edge state CFT, which is again due to dominant
Schmidt vectors contributing excitations that are mainly supported close to
the entanglement cut. Our findings suggest that such interpretations also
apply in a more general context.

## References

  * Gross  _et al._ (2023) D. Gross, A. Sevrin, and P. Zoller,  _The Physics of Quantum Information_ (WORLD SCIENTIFIC, 2023) https://www.worldscientific.com/doi/pdf/10.1142/13369 .
  * Joshi  _et al._ (2022) M. K. Joshi, F. Kranzl, A. Schuckert, I. Lovas, C. Maier, R. Blatt, M. Knap, and C. F. Roos, Science 376, 720 (2022).
  * Zhang  _et al._ (2017) J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z. X. Gong, and C. Monroe, Nature 551, 601 (2017).
  * Bohnet  _et al._ (2016) J. G. Bohnet, B. C. Sawyer, J. W. Britton, M. L. Wall, A. M. Rey, M. Foss-Feig, and J. J. Bollinger, Science 352, 1297 (2016).
  * Chen  _et al._ (2023) C. Chen, G. Bornet, M. Bintz, G. Emperauger, L. Leclerc, V. S. Liu, P. Scholl, D. Barredo, J. Hauschild, S. Chatterjee, _et al._ , Nature 616, 691 (2023).
  * Semeghini  _et al._ (2021) G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, _et al._ , Science 374, 1242 (2021).
  * Sompet  _et al._ (2022) P. Sompet, S. Hirthe, D. Bourgund, T. Chalopin, J. Bibo, J. Koepsell, P. Bojović, R. Verresen, F. Pollmann, G. Salomon, _et al._ , Nature 606, 484 (2022).
  * Léonard  _et al._ (2022) J. Léonard, S. Kim, J. Kwan, P. Segura, F. Grusdt, C. Repellin, N. Goldman, and M. Greiner, arXiv preprint arXiv:2210.10919 (2022).
  * Zhang  _et al._ (2022) W.-Y. Zhang, M.-G. He, H. Sun, Y.-G. Zheng, Y. Liu, A. Luo, H.-Y. Wang, Z.-H. Zhu, P.-Y. Qiu, Y.-C. Shen, _et al._ , arXiv preprint arXiv:2210.02936 (2022).
  * Cirac  _et al._ (2021) J. I. Cirac, D. Perez-Garcia, N. Schuch, and F. Verstraete, Reviews of Modern Physics 93, 045003 (2021).
  * Ganahl  _et al._ (2023) M. Ganahl, J. Beall, M. Hauru, A. G. Lewis, T. Wojno, J. H. Yoo, Y. Zou, and G. Vidal, PRX Quantum 4, 010317 (2023).
  * Miao and Barthel (2021) Q. Miao and T. Barthel, Phys. Rev. Lett. 127, 040603 (2021).
  * Acharya  _et al._ (2019) A. Acharya, T. Kypraios, and M. Guţǎ, Journal of Physics A: Mathematical and Theoretical 52, 234001 (2019).
  * Amico  _et al._ (2008) L. Amico, R. Fazio, A. Osterloh, and V. Vedral, Rev. Mod. Phys. 80, 517 (2008).
  * Bisognano and Wichmann (1975) J. J. Bisognano and E. H. Wichmann, J. Math. Phys. 16, 985 (1975).
  * Bisognano and Wichmann (1976) J. J. Bisognano and E. H. Wichmann, J. Math. Phys. 17, 303 (1976).
  * Dalmonte  _et al._ (2018) M. Dalmonte, B. Vermersch, and P. Zoller, Nature Physics 14, 827 (2018).
  * Giudici  _et al._ (2018) G. Giudici, T. Mendes-Santos, P. Calabrese, and M. Dalmonte, Phys. Rev. B 98, 134403 (2018).
  * Zhu  _et al._ (2020) W. Zhu, Z. Huang, Y.-C. He, and X. Wen, Phys. Rev. Lett. 124, 100605 (2020).
  * Casini  _et al._ (2011) H. Casini, M. Huerta, and R. C. Myers, Journal of High Energy Physics 2011, 1 (2011).
  * Cardy and Tonni (2016) J. Cardy and E. Tonni, Journal of Statistical Mechanics: Theory and Experiment 2016, 123103 (2016).
  * Dalmonte  _et al._ (2022) M. Dalmonte, V. Eisler, M. Falconi, and B. Vermersch, Annalen der Physik 534, 2200064 (2022).
  * Kokail  _et al._ (2021) C. Kokail, R. van Bijnen, A. Elben, B. Vermersch, and P. Zoller, Nature Physics 17, 936 (2021).
  * Anshu  _et al._ (2021) A. Anshu, S. Arunachalam, T. Kuwahara, and M. Soleimanifar, Nature Physics 17, 931 (2021).
  * Garrison and Grover (2018) J. R. Garrison and T. Grover, Phys. Rev. X 8, 021026 (2018).
  * Calabrese and Cardy (2004) P. Calabrese and J. Cardy, Journal of statistical mechanics: theory and experiment 2004, P06002 (2004).
  * Sels and Wouters (2017) D. Sels and M. Wouters, in  _Journal of Physics: Conference Series_ , Vol. 880 (IOP Publishing, 2017) p. 012047.
  * Kokail  _et al._ (2019) C. Kokail, C. Maier, R. van Bijnen, T. Brydges, M. K. Joshi, P. Jurcevic, C. A. Muschik, P. Silvi, R. Blatt, C. F. Roos, _et al._ , Nature 569, 355 (2019).
  * Eisler  _et al._ (2020) V. Eisler, G. D. Giulio, E. Tonni, and I. Peschel, Journal of Statistical Mechanics: Theory and Experiment 2020, 103102 (2020).
  * Elben  _et al._ (2020) A. Elben, B. Vermersch, R. van Bijnen, C. Kokail, T. Brydges, C. Maier, M. K. Joshi, R. Blatt, C. F. Roos, and P. Zoller, Phys. Rev. Lett. 124, 010504 (2020).
  * Elben  _et al._ (2023) A. Elben, S. T. Flammia, H.-Y. Huang, R. Kueng, J. Preskill, B. Vermersch, and P. Zoller, Nature Reviews Physics 5, 9 (2023).
  * Casini and Huerta (2009) H. Casini and M. Huerta, Classical and quantum gravity 26, 185005 (2009).
  * Eisler  _et al._ (2022) V. Eisler, E. Tonni, and I. Peschel, Journal of Statistical Mechanics: Theory and Experiment 2022, 083101 (2022).
  * Zache  _et al._ (2022) T. V. Zache, C. Kokail, B. Sundar, and P. Zoller, Quantum 6, 702 (2022).
  * Li and Haldane (2008) H. Li and F. D. M. Haldane, Phys. Rev. Lett. 101, 010504 (2008).
  * Chandran  _et al._ (2014) A. Chandran, V. Khemani, and S. L. Sondhi, Phys. Rev. Lett. 113, 060501 (2014).
  * Murciano  _et al._ (2022) S. Murciano, V. Vitale, M. Dalmonte, and P. Calabrese, Phys. Rev. Lett. 128, 140502 (2022).
  * Kim  _et al._ (2022) I. H. Kim, B. Shi, K. Kato, and V. V. Albert, Phys. Rev. B 106, 075147 (2022).
  * Ryu and Takayanagi (2006) S. Ryu and T. Takayanagi, Phys. Rev. Lett. 96, 181602 (2006).
  * Periwal  _et al._ (2021) A. Periwal, E. S. Cooper, P. Kunkel, J. F. Wienand, E. J. Davis, and M. Schleier-Smith, Nature 600, 630 (2021).
  * Fishman  _et al._ (2022) M. Fishman, S. R. White, and E. M. Stoudenmire, SciPost Phys. Codebases , 4 (2022).
  * Kranzl  _et al._ (2022) F. Kranzl, M. K. Joshi, C. Maier, T. Brydges, J. Franke, R. Blatt, and C. F. Roos, Phys. Rev. A 105, 052426 (2022).
  * Jurcevic  _et al._ (2014) P. Jurcevic, B. P. Lanyon, P. Hauke, C. Hempel, P. Zoller, R. Blatt, and C. F. Roos, Nature 511, 202 (2014).
  * Porras and Cirac (2004) D. Porras and J. I. Cirac, Phys. Rev. Lett. 92, 207901 (2004).
  * Liang  _et al._ (2019) Y.-C. Liang, Y.-H. Yeh, P. E. M. F. Mendonça, R. Y. Teh, M. D. Reid, and P. D. Drummond, Rep. Prog. Phys. 82, 076001 (2019).
  * Flammia and Liu (2011) S. T. Flammia and Y.-K. Liu, Phys. Rev. Lett. 106, 230501 (2011).
  * da Silva  _et al._ (2011) M. P. da Silva, O. Landon-Cardinal, and D. Poulin, Phys. Rev. Lett. 107, 210404 (2011).
  * Cerezo  _et al._ (2021) M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, _et al._ , Nature Reviews Physics 3, 625 (2021).
  * Kandala  _et al._ (2017) A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * Rasmussen and Williams (2005) C. E. Rasmussen and C. K. I. Williams,  _Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)_ (The MIT Press, 2005).
  * Sopena  _et al._ (2022) A. Sopena, M. H. Gordon, D. García-Martín, G. Sierra, and E. López, Quantum 6, 796 (2022).
  * Witten (2018) E. Witten, Reviews of Modern Physics 90, 045003 (2018).