# Divide-and-conquer embedding for QUBO quantum annealing

Minjae Jo1, Michael Hanks1,† and M. S. Kim1 1QOLS, Blackett Laboratory,
Imperial College London SW7 2AZ, UK<EMAIL_ADDRESS>
###### Abstract

Quantum annealing promises to be an effective heuristic for complex NP-hard
problems. However, clear demonstrations of quantum advantage are wanting,
primarily constrained by the difficulty of embedding the problem into the
quantum hardware. Community detection methods such as the Girvin–Newman
algorithm can provide a divide-and-conquer approach to large problems. Here,
we propose a problem-focused division for embedding, deliberately worsening
typical measures of embedding quality to improve the partial solutions we
obtain. We apply this approach first to the highly irregular graph of an
integer factorisation problem and, passing this initial test, move on to
consider more regular geometrically frustrated systems. Our results show that
a problem-focused approach to embedding can improve performance by orders of
magnitude.

##  1 Introduction

Quantum annealing (QA) is a method for energy minimisation that exploits the
highly non-classical quantum tunnelling effect. Tunnelling can be a more
efficient cooling mechanism than classical thermal fluctuation when the
potential landscape contains tall and narrow barriers, making QA a promising
near-term heuristic for hard optimisation problems [1, 2, 3, 4, 5]. The
difficulty of such problems scales exponentially in the number of input bits,
so that genuine quantum advantage can be expected with problem sizes well
below the 5000 qubit scale of the current D:Wave Advantage system4.1 device.
This device has consequently been a source of much recent interest, with
potential applications ranging across machine learning [6], Kosterlitz-
Thouless model [7], SU(2) lattice gauge theory [8], financial crisis [9], and
portfolio optimisation [10, 11].

For larger-scale problems the D:Wave quantum annealer rapidly loses its
advantage due to the problem of embedding. For example, when such problems are
submitted to the D:Wave quantum hardware, each binary variable may be
represented by a chain of physical qubits, where all physical qubits in such a
chain should be in the same state by the end of the annealing process.
Maintaining a chain with many qubits is not easy and the solution may be less
frequently optimal due to chain breaks and the freeze out phenomenon. The
quantum device therefore displays a problem-hardware trade-off: Accurate
representation of the problem may lead to a less efficient embedding, and
consequently lower device performance.

Divide-and-conquer algorithms, splitting problems into smaller pieces for more
efficient solution, are well-known in classical computing [12]. Examples
include the sorting algorithms [13] or the Strassen algorithm [14] for matrix
multiplication. Divide-and-conquer approaches have also been considered in
certain quantum contexts [15, 16], and could be helpful to achieve advantage
with the D:Wave quantum annealer.

Most existing literature seeking advantage with the D:Wave device relies on
minimising the number of physical qubits in the embedding, without regard for
the problem structure and in particular for the edge weights [17, 18]. Here,
we demonstrate improved performance with the opposite approach, dividing the
problem graph even at the expense of embedding with a larger physical qubit
number. Partitioning is performed using the Girvin–Newman algorithm for
community-detection [19, 20], and we refer to this approach throughout the
text as Problem-Focused Embedding (PFE) to contrast it with hardware-focused
qubit number minimisation.

As an initial proof of principle, we demonstrate improved performance for
problems on irregular graphs where a PFE is most naturally applied. While we
chose integer factorisation [21] for our demonstration, many other examples
such as the anti-ferromagnetic Ising model in scale-free graphs [22] may be
cited. We then look at a class of interesting NP-hard problems on regular
structures: Geometrically frustrated lattice systems. In these systems, ground
state degeneracy triggers novel phases of matter including stripe phases [23,
24] and spin liquids [25, 26, 27, 28] characterised by their long-range
correlation and the absence of ordinary magnetic order. In particular, the
Kagome lattice provides a prototype of difficult frustrated-lattice problems
and can be physically realized in SrCr8-xGa4+xO19 [29]. Thus, finding the
ground state of a Kagome lattice has attracted much attention in condensed
matter physics; however, obtaining high-quality solutions using heuristics
such as simulated annealing is not easy due to many the local minima and
limited gradient information.

Our results show that tight clusters are not required to observe the benefit
of a divide-and-conquer strategy. As a measure for performance we take
$F={-\log_{10}{\left({\epsilon}\right)}}/{\tau}$, where $\tau$ is the mean
time required for the quantum device to return each independent set of results
and $\epsilon$ is the probability that the device fails to return the global
optimum solution (see [30] for a discussion of alternatives). Despite the fact
that a PFE typically requires more qubits (due to the overhead of packing sub-
problems in the device), we see better performance metrics.

The remainder of this paper is organised as follows. In Section 2, we describe
problems in quantum annealing in terms of their graph structure, introducing
the concepts of sub-graphs and partial solutions. Our PFE approach is
presented, alongside a description of the Girvin–Newman algorithm. Next, in
Section 3 experimental results demonstrating advantage in the D:Wave system
for both regular and irregular structures are presented. Finally, in Section 4
we summarise our work and discuss ways to identify the potential advantage of
QA.

##  2 Graph Structure and Quantum Annealing

In this section, we discuss the useful graph description of QA problems and
introduce the PFE. We begin with the notion of QA and a limitation arising
from embedding. To overcome the limitation, the PFE consists of two following
steps: divide the problem graph and then merge the partial solutions. We shall
use the Girvin–Newman algorithm to divide the graph with community structure.
We discuss sub-graph reduction and partial solutions using the example of the
1D Ising chain. Based on these concepts, the PFE is introduced using the
example of integer factorisation.

###  2.1 Viewing Quantum Annealing as a Graph

QA aims to minimise a quadratic function of binary variables in the form of a
ground state search for the Ising-type Hamiltonian [31]

$\displaystyle H_{\text{target}}$
$\displaystyle=\sum_{i}h_{i}\hat{\sigma}_{z,i}+\sum_{i>j}J_{ij}\hat{\sigma}_{z,i}\hat{\sigma}_{z,j}\,,$
(1)

where $\hat{\sigma}_{z,i}$ describes the Pauli-Z operator for the $i$th qubit,
while $h_{i}$ and $J_{ij}$ are the qubit biases and coupling strengths. This
is an NP-Hard problem in general, and its solution is obtained by initialising
the system instead in the (known, separable) ground state of the Hamiltonian

$\displaystyle H_{\text{initial}}$ $\displaystyle=\sum_{i}\hat{\sigma}_{x,i}$
(2)

and then varying the relative magnitudes of $H_{\text{initial}}$ and
$H_{\text{target}}$ according to

$\displaystyle H\left(t\right)$
$\displaystyle=f\left(t\right)H_{\text{initial}}+g\left(t\right)H_{\text{target}},$
$\displaystyle g\left(0\right)$ $\displaystyle=0,$
$\displaystyle\lim_{t\rightarrow\infty}f\left(t\right)$ $\displaystyle=0\,,$
(3)

where $f\left(t\right)$ and $g\left(t\right)$ are real, non-negative,
monotonic functions of time $t$. $f\left(t\right)$ is said to be the strength
of quantum fluctuations, as the $\hat{\sigma}_{x}$ terms in
$H_{\text{initial}}$ facilitate transitions between the eigenstates of
$H_{\text{target}}$. In the adiabatic limit, the target ground state is
obtained deterministically. Due to noise and time constraints, adiabatic
evolution is rarely observed in practice; quantum tunnelling becomes the key
mechanism for computational advantage.

QA is often compared to digital simulated annealing (SA), where the
temperature parameter plays a similar role to the strength of quantum
fluctuations. In SA, the idea is to use thermal fluctuations to allow the
system to escape from local minima of the potential landscape so that the
system reaches the global minimum under an appropriate annealing schedule
(rate of decrease of temperature). For an energy barrier $\Delta E$, the
escape probability is $\sim\exp(-\Delta E/T)$ at temperature $T$. In QA, the
system will likely tunnel through the barrier if it is narrow: The tunnelling
probability is $\sim\exp(-w\sqrt{\Delta E}/\Gamma)$ for tunnelling fluctuation
field $\Gamma$ and barrier width $w$. For high and narrow potential barriers,
QA utilising the tunnelling effect can therefore be more efficient.

Connectivity constraints do not affect $H_{\text{initial}}$, as it contains no
coupling terms. The coupling terms in $H_{\text{target}}$ must however be
mapped to interactions in the physical device. A coupling term in
$H_{\text{target}}$ may be simulated either in the form of direct interactions
between neighbouring physical qubits, or indirectly in the form of a chain of
interacting ancillae. Considering these qubits as nodes of a graph, and
coupling terms or interactions as its edges, the problem of mapping an
optimisation problem to the quantum annealing device corresponds to the graph
problem of minor embedding: A first graph is said to be a minor of a second if
it can be obtained by deleting vertices and contracting edges (identifying
neighbouring vertices) in this second graph. A quadratic unconstrained binary
optimisation (QUBO) problem can be embedded in a quantum annealing device if
the graph of coupling terms in $H_{\text{target}}$ can be expressed as a minor
of the graph of physical qubit interactions in the device. The typical measure
of the quality of such an embedding is the overhead in the number of physical
qubits, which corresponds to the number of edge-contractions required to
obtain the logical coupling graph from the physical interaction graph.

###  2.2 Sub-graphs and Partial Solutions

For many problems, we can identify ahead of time that the solution is
contained in only a small region of the total configuration space. Consider
the $1$D periodic chain as a simple example, as displayed in Fig. 1.

Figure 1:  Example of sub-graph reduction for the 1D Ising chain. Here
$J_{ij}$ terms along the $1$D chain are all equal initially, as are $h_{i}$
terms. (a) The initial Ising chain. Spins have been arranged at a global
optimum to emphasise how this is preserved, though the specific configuration
of the chain will not be known initially. Spins comprising the sub-graph to be
reduced are boxed and separated in the centre. (b) Enumeration of
configurations and energies for the three spins in the sub-graph to be
reduced, and identification (in red) of the four lowest-energy solutions
exhausting all possible boundary-spin configurations. (c) The modified Ising
chain. The modified red component of the graph is boxed and separated in the
centre. Also shown in red are the modified self-energy and interaction terms
for the two spins in this component. These new terms are chosen to preserve
the energy level structure associated with the low-energy configurations of
the original sub-graph.

We begin by looking at only a small connected segment, enumerating all
possible solutions and their energies, which we will refer to as the segment’s
‘internal configurations’. Now, when this segment is incorporated into the
global solution for the total chain, the appropriate internal configuration
will depend only on the states of those bits on the boundary between the
segment and the rest of the chain — all other bits within the segment are used
only to select the internal configuration of minimal energy satisfying these
boundary conditions. This means that we can ignore all but four solutions from
the internal segment, taking the minimal-energy solution for each of the four
boundary-bit configurations. The problem may then be remapped so that the
initial segment is represented by only the two boundary bits, with energies
and an interaction chosen to correspond to the energy gaps for the four
solutions we retain. Applying this method recursively, the $1$D chain (or any
tree structure) can be solved in linear time regardless of the problem
parameters.

Remapping the problem to eliminate variables comes at no cost when the problem
graph is a tree structure. This may be highly significant for scale-free
graphs and others with many nodes of degree one or two. However, for many
problems we expect to be presented with a much more densely connected graph,
with many small but non-zero edges corresponding to measurement samples. For
these more general graphs such elimination introduces new interaction terms of
order up to the degree of the removed node [16]. Conversely, reducing the
order of interactions (or ‘quadratization’) requires the introduction of new
variables [32]. While divide-and-conquer strategies have been applied in the
more general context of quantum approximate optimization algorithms (QAOA)
[33], they have thus far been capable only of limited application in QUBO-
restricted experimental quantum annealing. We next introduce a modified
divide-and-conquer method suitable for quantum annealing.

###  2.3 A Divide-and-Conquer Method for QUBO-Restricted Quantum Annealing

It is intuitive that solutions for distinct clusters should retain their form
in spite of interconnecting edges, when those edges have very small weights.
This intuition is supported by the fact that the eigenvalues $\lambda_{k}$ of
an Hermitian matrix $A$ are well-conditioned under small perturbations $\delta
A$,

$\displaystyle\max_{k}\left\lvert\lambda_{k}\left(A+\delta
A\right)-\lambda_{k}\left(A\right)\right\rvert\leq\left\lVert\delta
A\right\rVert_{2},$ (4)

and is implicit in spectral methods for graph sparsification [34]. Consider
two highly connected clusters connected by a number of weak edges, as
represented diagrammatically in Figure 2.

Figure 2:  Diagrammatic representation of sampling local solutions in a
tightly clustered problem graph. Here the problem graph is grouped into two
tightly connected clusters, with a number of weaker edges between the two
clusters. The value $a$, twice the sum of boundary edge magnitudes, determines
an upper bound on the energies of local configurations. Configurations
exceeding this bound may be discarded, reducing the local search space.
Supposing $N$ valid local configurations exist in the first cluster, and $M$
in the second, the overall problem complexity is then reduced to the sum of
the individual sub-problems, and an additional merging term of order $N\cdot
M$.

While the ground state of each cluster may no longer correspond to the target
solution for the entire graph, we observe that the energy penalty introduced
by the connecting weak edges is bounded above by a value equal to twice the
sum of the edge weight magnitudes (twice because variables take the values
$\pm 1$). Call this bound $a$. When the global solution is projected onto a
local state in each cluster, each resulting local state has an energy at most
$a$ larger than its local ground state; solutions with local energies
exceeding this bound could be improved by simply choosing the local ground
state, regardless of the boundary configuration. Supposing that $N$ of the
total $2^{S}$ internal energies fall within the gap $a$ of the ground state
($S$ nodes in the cluster), the cluster is then effectively reduced to the
smaller $N$-dimensional system. In fact, once these low energy solutions have
been determined, a stricter gap follows using a boundary sum with edges from
only those nodes where each solution differs from the local ground state.

One remaining question is how to sample efficiently from our restricted set of
local solutions, which may differ widely from one another in their particular
configurations. Fortunately, unlike other heuristics such as genetic or Tabu
search [35, 36], simulated and quantum annealing are capable of sampling
states with a wider range of energies in the global configuration space by
manipulating the cooling schedule and final temperature.

###  2.4 A Problem-Focused Embedding

Figure 3:  Example of the PFE for integer factoring problem $143=13\times 11$
with 12 qubits. (a) The problem graph is partitioned as two sub-graphs (the
ten red nodes and two blue nodes) by the Girvin–Newman algorithm. The
thickness of the line represents the edge weight. (b) The sampled energy
distribution for the larger red sub-graph. Only configurations lower than the
blue boundary line are considered in searching for global solutions. The final
temperature (quantum fluctuation magnitude) has been chosen to sample this
region of interest.

To implement our divide-and-conquer strategy, we use a simple partitioning
method based on the edge betweenness centrality, defined as the number of the
shortest paths that go through an edge in a graph [37]. The problem graph is
divided into sub-graphs using the weighted Girvin–Newman algorithm [19, 20].
This algorithm detects communities from the problem graph by gradually
removing bottleneck edges between them. When several edges have the same
betweenness centrality, one is removed at random.

Figure 3 shows a simple application of our divide-and-conquer strategy to the
factorization of $143=13\times 11$. In the QUBO representation, this problem
reduces to a graph with 12 nodes. The problem graph is divided into smaller
sub-graphs, each sub-graph is solved by QA or SA, and we merge the solutions.
Hereafter, we focus only on bipartition, as a proof of concept.

##  3 Experimental results

Next we apply the PFE of Section 2 to several small-scale examples of NP-hard
problems. This approach is naturally adopted for irregular problem graphs such
as that of the integer factoring problem described in Fig. 3, and so we begin
there. We then apply the PFE approach to regular, geometrically-frustrated
systems, to show that tight clusters are not always required.

Our results $F/F_{SA}$ are expressed relative to a standard application of
digital simulated annealing111using a single process on an Apple M1 chip. This
unit of measure is arbitrary and we note that, particularly for irregular
graph structures, SA does not necessarily represent the most efficient
classical method. We choose this form both for succinctness in the tables to
follow and to emphasise the increase quantum advantage, as distinct from any
bare classical advantage arising from our division of the problem.

###  3.1 An Irregular Graph Structure

Problem graphs with irregular structures introduce two complications:

  1. 1.

They have broad structural features, often exploitable via specific classical
strategies.

  2. 2.

They are highly non-planar and usually contain long-range interactions, making
device embedding difficult.

These constraints imply the need for hybrid classical and quantum solutions.

One application that generates irregular problem graph structures is integer
factorisation. Integer factorisation is an NP-hard problem for a classical
computer, and consequently forms the basis of RSA public key cryptography.
While Shor’s algorithm achieves polynomial complexity for a hypothetical
universal quantum computer, the largest number factored using Shor’s algorithm
at present is 21 [38], so there is still a long way to go before practically
meaningful factorisation is achieved in this manner. Recently, QA with the
D:Wave device was proposed as a factorisation heuristic using a transformation
that maps the problem into a QUBO model [21].

Table 1 shows the relative performance $F/F_{\rm SA}$ of standard QA and the
PFE for the factorisation of $1591=43\times 37$.

Table 1: Results of the relative performance metric of factoring $1591=43\times 37$ with 29 qubits using the standard and PFE approaches. QA is performed by D:Wave Advantage_system4.1. The number of samples is $10^{3}$. The number of physical nodes $N_{\rm P}$, the sample estimate for the success probability using Bayes’ inference $p_{\rm Bayes}$, the total run time in microseconds (including overheads) $\tau_{\rm Real}$, and performance metric $F/F_{\rm SA}$ (relative to that of a standard application of SA) are listed.  Method | $N_{\rm P}$ | $p_{\rm Bayes}$ | $\tau_{\rm Real}$ | $F/F_{\rm SA}$
---|---|---|---|---
Standard QA | $73$ | $10^{-3}$ | $170,302$ | $0.26$
PFE (QA$+$QA) |  $66$ & $2$ | $0.17$ | $405,995$ | $19.98$
PFE (SA$+$SA) | $-$ | $0.17$ | $652,691$ | $12.43$

The relative performance metric of the PFE is enhanced approximately $76$
times beyond standard QA (as proposed, for instance, in [21]) and the required
number of physical qubits is decreased. Applying the same divide-and-conquer
strategy in an entirely classical manner to SA, a relative improvement of only
about $12$ times is observed over the standard approach.

###  3.2 A Regular Planar Structure

The regular structures of geometrically frustrated lattices in many-body
physics have two attractive qualities for QA:

  1. 1.

Their low degree and local interactions should ease embedding in the quantum
device.

  2. 2.

The lack of clustering or other broad structure means that heuristics such as
simulated annealing may often be the most fruitful approach.

Despite these advantages, the performance of standard QA rapidly decreases as
we near about one hundred qubits, due to growing embedding constraints. This
rough $100$ qubit limit presents a challenge for the practical application of
QA to many body physics.

A divide-and-conquer approach may help to raise the tolerable qubit number. It
is not immediately clear, however, whether we should expect to observe
advantage in problem division for a highly regular structure. Next we
therefore investigate the performance of a PFE approach to ground state search
on the Kagome lattice, shown in Figure 4. The Kagome lattice is a prototypical
example of geometric frustration in many body physics.

Figure 4:  The (periodic) Kagome lattice, or trihexagonal tiling, forms a
degree-$4$ regular graph, and may be obtained by deleting nodes in a
triangular lattice.

Table 2 shows the relative performance $F/F_{\rm SA}$ of standard QA and the
PFE approach for QA and SA.

Table 2: Results of the relative performance metric of the Kagome lattice with $90$ qubits using the standard QA and the PFE approach to both QA and SA. QA is performed by D:Wave Advantage_system4.1. The number of samples is $10^{3}$ and $h=J$. The number of physical nodes $N_{\rm P}$, the sample estimate for the success probability using Bayes’ inference $p_{\rm Bayes}$, the total run time in microseconds (including overheads) $\tau_{\rm Real}$, and performance metric $F/F_{\rm SA}$ (relative to that of a standard application of SA) are listed. Method | $N_{\rm P}$ | $p_{\rm Bayes}$ | $\tau_{\rm Real}$ | $F/F_{\rm SA}$
---|---|---|---|---
Standard QA | $127$ | $0.02$ | $152,904$ | $2.30$
PFE (QA$+$QA) |  $54$ & $48$ | $0.92$ | $527,265$ | $69.47$
PFE (SA$+$SA) | $-$ | $0.92$ | $1,537,878$ | $23.83$

In comparison to a standard application of QA or SA, the PFE approach
increases the total run time due to the final matching of component solutions,
but increases to an even greater degree the probability of discovering the
global optimum with each sample. We observe that the PFE approach improves the
performance by approximately $30$ times for QA ($69/2.3$), and approximately
$24$ times when applied classically for SA. While the bulk of the improvement
seems therefore classical, there is nonetheless a significant further boost
that is unique to the quantum device.

##  4 Discussion

We have proposed and tested experimentally a problem-focused divide-and-
conquer approach to quantum annealing. Dividing the problem graph improves the
quality of partial solutions by improving the embedding quality for
interactions within each partition, facilitating quantum advantage. We
successfully applied this approach to irregular graph structures, where
advantage is most expected, but also to highly regular geometrically
frustrated lattices.

We expect that the improved performance we observe for frustrated lattices
with many qubits will be applicable in many-body physics: Frustration in the
Kagome lattice can lead to the emergence of a flat band, with strongly-
correlated atomic behaviour determined by interactions and topology. Our
initial results suggest that QA with the D:Wave system could be fruitfully
applied to ground state discovery in strongly-correlated systems, identifying
critical phenomena in frustrated models and predicting experimental results in
strongly-correlated materials.

In contrast to similar recent work [16], our approach requires only quadratic
interaction terms, and is therefore compatible with D:Wave quantum devices.
Further, the residual problem complexity depends in our case on the ratio
between boundary edge weights and the spectral density of low-energy states,
rather than the boundary node number.

For our proof-of-concept results in Section 3, we implemented community
detection using the weighted Girvin–Newman algorithm. We acknowledge however
that, community detection being a well studied problem, alternative methods
for community detection may improve performance. We use a performance metric
making use of knowledge of the global minimum (through $\epsilon$, the
probability of failing to return this minimum value). Presently, D:Wave
devices demonstrate quantum advantage with problems roughly on the order of
$150$ qubits. We expect therefore that problems solved on contemporary quantum
annealing devices will not be those that cannot be solved by a classical
machine, but rather those for which a classical machine may solve the problem
in a less cost-efficient manner. As a result, the use of global information
for test problems of comparable size is justified.

We identify two interesting questions arising from our results, which could
motivate future work: Firstly, though we restrict ourselves to a bipartition
for proof-of-concept, three or perhaps multiple sub-graphs could become
optimal as the problem size increases. How might the density of low-energy
states be estimated to assess the trade-off between a more efficient embedding
as a result of further sub-division and the cost of merging partial solutions?
Secondly, rather than solving each partition separately, it may prove fruitful
to develop a whole-problem embedding, giving priority to intra-community
links. At present, strong coupling in ancillary chain structures interferes
with the solution in each partition. Is an alternative construction for
ancillary chains possible, that would allow communities to interact coherently
without such strong interference?

## Acknowledgments

This research was supported by the education and training program of the
Quantum Information Research Support Center, funded through the National
research foundation of Korea (NRF) by the Ministry of science and ICT (MSIT)
of the Korean government, No. 2021M3H3A103657313 (MJ). We acknowledge the
Samsung GRC project and the UK Hub in Quantum Computing and Simulation, part
of the UK National Quantum Technologies Programme with funding from UKRI EPSRC
grant EP/T001062/1. This work was supported by the Engineering and Physical
Sciences Research Council of the UK (EPSRC) Grant number EP/W032643/1.

## References

## References

  * [1] Hernandez M and Aramon M 2017 Quantum Information Processing 16 133
  * [2] Babbush R, Love P J and Aspuru-Guzik A 2014 Scientific Reports 4 6603
  * [3] Farhi E, Goldstone J, Gutmann S, Lapan J, Lundgren A and Preda D 2001 Science 292 472–475 (Preprint https://www.science.org/doi/pdf/10.1126/science.1057726) URL https://www.science.org/doi/abs/10.1126/science.1057726
  * [4] Willsch D, Willsch M, Gonzalez Calaza C D, Jin F, De Raedt H, Svensson M and Michielsen K 2022 Quantum Information Processing 21 141
  * [5] Yaacoby R, Schaar N, Kellerhals L, Raz O, Hermelin D and Pugatch R 2022 Phys. Rev. E 105(3) 035305 URL https://link.aps.org/doi/10.1103/PhysRevE.105.035305
  * [6] Benedetti M, Realpe-Gómez J, Biswas R and Perdomo-Ortiz A 2017 Phys. Rev. X 7(4) 041052 URL https://link.aps.org/doi/10.1103/PhysRevX.7.041052
  * [7] King A D, Carrasquilla J, Raymond J, Ozfidan I, Andriyash E, Berkley A, Reis M, Lanting T, Harris R, Altomare F, Boothby K, Bunyk P I, Enderud C, Fréchette A, Hoskinson E, Ladizinsky N, Oh T, Poulin-Lamarre G, Rich C, Sato Y, Smirnov A Y, Swenson L J, Volkmann M H, Whittaker J, Yao J, Ladizinsky E, Johnson M W, Hilton J and Amin M H 2018 Nature 560 456–460 ISSN 0028-0836, 1476-4687 URL http://www.nature.com/articles/s41586-018-0410-x
  * [8] A Rahman S, Lewis R, Mendicelli E and Powell S 2021 Physical Review D 104 034501 ISSN 2470-0010, 2470-0029 URL https://link.aps.org/doi/10.1103/PhysRevD.104.034501
  * [9] Ding Y, Gonzalez-Conde J, Lamata L, Martín-Guerrero J D, Lizaso E, Mugel S, Chen X, Orús R, Solano E and Sanz M 2021 arXiv:1904.05808 [cond-mat, physics:quant-ph] URL http://arxiv.org/abs/1904.05808
  * [10] Palmer S, Sahin S, Hernandez R, Mugel S and Orus R 2021 arXiv:2106.06735 [quant-ph, q-fin] URL http://arxiv.org/abs/2106.06735
  * [11] Phillipson F and Bhatia H S 2020 arXiv:2012.01121 [quant-ph, q-fin] URL http://arxiv.org/abs/2012.01121
  * [12] Rosenberg A L and Heath L S 2005 Graph separators, with applications (New York: Kluwer Academic Publishers) ISBN 978-0-306-46977-0 oCLC: 757106468
  * [13] Sedgewick R 1978 Commun. ACM 21 847–857 ISSN 0001-0782 URL https://doi.org/10.1145/359619.359631
  * [14] Strassen V 1969 Numerische Mathematik 13 354–356
  * [15] Song W, Lim Y, Jeong K, Ji Y S, Lee J, Kim J, Kim M S and Bang J 2022 Quantum Science and Technology 7 025009 ISSN 2058-9565 URL https://iopscience.iop.org/article/10.1088/2058-9565/ac51b0
  * [16] Guerreschi G G 2021 Solving Quadratic Unconstrained Binary Optimization with divide-and-conquer and quantum algorithms arXiv:2101.07813 [quant-ph] URL http://arxiv.org/abs/2101.07813
  * [17] Klymko C, Sullivan B D and Humble T S 2014 Quantum Information Processing 13 709–729 URL https://doi.org/10.1007/s11128-013-0683-9
  * [18] Choi V 2008 arXiv:0804.4884 [quant-ph] URL http://arxiv.org/abs/0804.4884
  * [19] Girvan M and Newman M E J 2002 Proceedings of the National Academy of Sciences 99 7821–7826 (Preprint https://www.pnas.org/doi/pdf/10.1073/pnas.122653799) URL https://www.pnas.org/doi/abs/10.1073/pnas.122653799
  * [20] Newman M E J 2004 Phys. Rev. E 70(5) 056131 URL https://link.aps.org/doi/10.1103/PhysRevE.70.056131
  * [21] Jiang S, Britt K A, McCaskey A J, Humble T S and Kais S 2018 Scientific Reports 8 17667 ISSN 2045-2322 URL http://www.nature.com/articles/s41598-018-36058-z
  * [22] Bartolozzi M, Surungan T, Leinweber D B and Williams A G 2006 Phys. Rev. B 73(22) 224419 URL https://link.aps.org/doi/10.1103/PhysRevB.73.224419
  * [23] Bramwell S T and Gingras M J P 2001 Science 294 1495–1501 (Preprint https://www.science.org/doi/pdf/10.1126/science.1064761) URL https://www.science.org/doi/abs/10.1126/science.1064761
  * [24] Zhao K, Deng H, Chen H, Ross K A, Petříček V, Günther G, Russina M, Hutanu V and Gegenwart P 2020 Science 367 1218–1223 (Preprint %****␣ms.bbl␣Line␣125␣****https://www.science.org/doi/pdf/10.1126/science.aaw1666) URL https://www.science.org/doi/abs/10.1126/science.aaw1666
  * [25] Broholm C, Aeppli G, Espinosa G P and Cooper A S 1990 Phys. Rev. Lett. 65(25) 3173–3176 URL https://link.aps.org/doi/10.1103/PhysRevLett.65.3173
  * [26] Balents L 2010 Nature 464 199–208 URL https://doi.org/10.1038/nature08917
  * [27] Cépas O and Canals B 2012 Phys. Rev. B 86(2) 024434 URL https://link.aps.org/doi/10.1103/PhysRevB.86.024434
  * [28] Lin T, Ke X, Thesberg M, Schiffer P, Melko R G and Gingras M J P 2014 Phys. Rev. B 90(21) 214433 URL https://link.aps.org/doi/10.1103/PhysRevB.90.214433
  * [29] Ramirez A P, Espinosa G P and Cooper A S 1990 Phys. Rev. Lett. 64(17) 2070–2073 URL https://link.aps.org/doi/10.1103/PhysRevLett.64.2070
  * [30] King J, Yarkoni S, Nevisi M M, Hilton J P and McGeoch C C 2015 arXiv:1508.05087 [quant-ph] URL http://arxiv.org/abs/1508.05087
  * [31] Kadowaki T and Nishimori H 1998 Physical Review E 58 5355–5363 ISSN 1063-651X, 1095-3787 URL https://link.aps.org/doi/10.1103/PhysRevE.58.5355
  * [32] Dattani N 2019 Quadratization in discrete optimization and quantum mechanics arXiv:1901.04405 [physics, physics:quant-ph] URL http://arxiv.org/abs/1901.04405
  * [33] Li J, Alam M and Ghosh S 2022 IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 1–1 ISSN 0278-0070, 1937-4151 URL https://ieeexplore.ieee.org/document/9911690/
  * [34] Batson J, Spielman D A, Srivastava N and Teng S H 2013 Communications of the ACM 56 87–94 ISSN 0001-0782, 1557-7317 URL https://dl.acm.org/doi/10.1145/2492007.2492029
  * [35] Boros E, Hammer P L and Tavares G 2007 Journal of Heuristics 13 99–132 ISSN 1381-1231, 1572-9397 URL http://link.springer.com/10.1007/s10732-007-9009-3
  * [36] Kochenberger G, Hao J K, Glover F, Lewis M, Lü Z, Wang H and Wang Y 2014 Journal of Combinatorial Optimization 28 58–81 ISSN 1382-6905, 1573-2886 URL http://link.springer.com/10.1007/s10878-014-9734-0
  * [37] Freeman L C 1977 Sociometry 40 35 ISSN 00380431 URL https://www.jstor.org/stable/3033543?origin=crossref
  * [38] Martín-López E, Laing A, Lawson T, Alvarez R, Zhou X Q and O’Brien J L 2012 Nature Photonics 6 773–776 URL https://doi.org/10.1038/nphoton.2012.259


# Experimental verification of the area law of mutual information in a quantum
field simulator

Mohammadamin Tajik<EMAIL_ADDRESS>Vienna Center for Quantum
Science and Technology, Atominstitut, TU Wien, Stadionallee 2, 1020 Vienna,
Austria  Ivan Kukuljan  Max-Planck Institute of Quantum Optics, Hans-
Kopfermann-Str. 1, 85748 Garching, Germany  Munich Center for Quantum Science
and Technology (MCQST), Schellingstr. 4, D-80799 München  Spyros Sotiriadis
Dahlem Center for Complex Quantum Systems, Freie Universität Berlin, Berlin,
Germany  Bernhard Rauer  Vienna Center for Quantum Science and Technology,
Atominstitut, TU Wien, Stadionallee 2, 1020 Vienna, Austria  Laboratoire
Kastler Brossel, ENS-Université PSL, CNRS, Sorbonne Université, Collége de
France, 24 rue Lhomond, 75005 Paris, France  Thomas Schweigler  Vienna Center
for Quantum Science and Technology, Atominstitut, TU Wien, Stadionallee 2,
1020 Vienna, Austria  Federica Cataldini  Vienna Center for Quantum Science
and Technology, Atominstitut, TU Wien, Stadionallee 2, 1020 Vienna, Austria
João Sabino  Vienna Center for Quantum Science and Technology, Atominstitut,
TU Wien, Stadionallee 2, 1020 Vienna, Austria  Instituto de Telecomunicações,
Physics of Information and Quantum Technologies Group, Av. Rovisco Pais 1,
1049-001, Lisbon, Portugal  Instituto Superior Técnico, Universidade de
Lisboa, Av. Rovisco Pais 1, 1049-001, Lisbon, Portugal  Frederik Møller
Philipp Schüttelkopf  Si-Cong Ji  Vienna Center for Quantum Science and
Technology, Atominstitut, TU Wien, Stadionallee 2, 1020 Vienna, Austria  Dries
Sels  Department of physics, New York University, New York City, NY 10003, USA
Center for Computational Quantum Physics, Flatiron Institute, 162 5th Ave, New
York, NY 10010, USA  Eugene Demler  Institute for Theoretical Physics, ETH
Zurich, Wolfgang-Pauli-Str. 27, 8093 Zurich, Switzerland  Jörg Schmiedmayer
<EMAIL_ADDRESS>Vienna Center for Quantum Science and Technology,
Atominstitut, TU Wien, Stadionallee 2, 1020 Vienna, Austria

###### Abstract

Theoretical understanding of the scaling of entropies and the mutual
information has led to significant advances in the research of correlated
states of matter, quantum field theory, and gravity. Measuring von Neumann
entropy in quantum many-body systems is challenging as it requires complete
knowledge of the density matrix. In this work, we measure the von Neumann
entropy of spatially extended subsystems in an ultra-cold atom simulator of
one-dimensional quantum field theories. We experimentally verify one of the
fundamental properties of equilibrium states of gapped quantum many-body
systems, the area law of quantum mutual information. We also study the
dependence of mutual information on temperature and the separation between the
subsystems. Our work is a crucial step toward employing ultra-cold atom
simulators to probe entanglement in quantum field theories.

The study of quantum information measures is central to a wide range of areas
in physics, from condensed matter and atomic physics to high energy physics
and gravity [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. Some of the most commonly
studied quantities in quantum information are the (entanglement) entropy and
the quantum mutual information. If a system described by the density matrix
$\varrho$ is composed of subsystems $A$ and $B$, then the von Neumann (vN)
entropy of subsystem $A$ is defined as

$S_{A}=-\text{Tr}\left(\varrho_{A}\ln\left(\varrho_{A}\right)\right)\,,$ (1)

where $\varrho_{A}=\text{Tr}_{B}\left(\varrho\right)$ is the reduced density
matrix of subsystem $A$. If the state $\varrho$ is pure
($\text{Tr}\left(\varrho^{2}\right)=1$), then the vN entropy is a measure of
entanglement between $A$ and $B$, thus called entanglement entropy, where
$S_{A}=S_{B}$. For mixed states ($\text{Tr}\left(\varrho^{2}\right)<1$), the
vN entropy captures both classical and quantum correlations, and it is no
longer a good measure of entanglement. For mixed states, several other
measures and witnesses of entanglement have been studied, with the positivity
of partial transpose criterion  and quantum discord being prominent examples
[1, 13]. In cases where, rather than entanglement, the shared amount of
information between two subsystems $A$ and $B$ is of interest, the quantum
mutual information (MI),

$I(A:B)=S_{A}+S_{B}-S_{A\cup B}\,,$ (2)

is a central object of study. It measures the total amount of correlation
between the two subsystems, including all higher-order correlations for both
pure and mixed states. For pure states, the value of the MI is equal to twice
the entanglement entropy of one of the subsystems.

Figure 1: Schematic of the experimental protocol. a, The experimental protocol
starts by cooling down a pair of tunnelling coupled superfluids in a double-
well potential with a finite single particle tunnelling rate $J$, typically
$\approx 2\pi\times 1\,$\mathrm{Hz}$$. The radial trapping frequency is
$\omega_{\perp}=2\pi\times 1.4\,$\mathrm{kHz}$$ and typical linear densities
are $n_{\mathrm{1D}}\approx 70\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{-1}$$.
b, Right after the cooling process, the tunnelling rate is changed to zero in
$2\,$\mathrm{ms}$$. At $t=0$, the two condensates are already uncoupled and
their independent evolution following the Tomonaga-Luttinger liquid
Hamiltonian starts. The measured phase correlations at different times are
used to fit the full covariance matrix $\Gamma$. c, To calculate the vN
Entropy for a subsystem, $S(\Gamma_{A})$, we use the correlations in that
specific region (yellow shaded area).  Note that $Q$, $R$, and $P$ have
different units and the range of the color axis are different. All tomography
results are presented and discussed in SM.

Information theory measures reveal one of the fundamental properties of
quantum many-body systems, the area laws [3, 14]. It was first noticed in
gravitational physics, that surprisingly, the entropy of a black hole is
proportional to the surface area of its event horizon and not its volume [11].
Interestingly, a similar property was found in quantum many-body systems: The
vN entropy of ground states of systems with a gapped Hamiltonian scales with
the surface area of the subsystem and not its volume [14]. Even more general,
thermal states of systems with a gapped Hamiltonian exhibit an area law of
mutual information [15]. This means that the information shared between parts
of a quantum many-body system is only considerable over a short distance, set
by the correlation length. In particular, such a bound on the required
information to model a many-body system provides the foundations for the
overwhelming success of tensor network-based methods [2]. In contrast, it is
known that critical systems described by conformal field theory exhibit a
logarithmic scaling instead of the area law [9]. In thermal states, while the
mutual information has an area or a log law, the vN entropy will exhibit a
volume law where it is proportional to the volume of the subsystem.

Extracting quantum information measures in quantum many-body systems has been
the aim of several experiments [16, 17, 18, 19, 20]. Calculating the vN
entropy requires access to the density matrix of the full system, $\varrho$,
which usually  requires a full state tomography in different experimental
platforms. Instead of $\varrho$, several techniques are developed to measure
the purity, $\text{Tr}\left(\varrho^{2}\right)$, which enables the calculation
of second-order Rényi entropy,
$S_{2}=-\ln(\text{Tr}\left(\varrho^{2}\right))$. These methods can be based on
the interference of two identical copies of a quantum system [17] or
randomised measurements on a single copy [20]. Note that, in these examples,
the purity is directly measured but the full state is not reconstructed, hence
the calculation of vN entropy is not feasible.

Figure 2: Area law of MI and volume law of vN entropy. a, A system of size
$L=49\,$\mathrm{\SIUnitSymbolMicro m}$$ is divided into a subsystem $A$ of
length $l$ and its complement subsystem of length $(L-l)$. b, The experimental
results for $I(A:A^{\mathrm{C}})$, $S_{A}$, and $S_{A^{\mathrm{C}}}$,
calculated based on $N=7$ modes, are plotted as bullets with error bars
representing the $95\%$ confidence intervals obtained via bootstrapping [21].
The shaded areas show the $95\%$ confidence interval for the theory
predictions, considering the uncertainty in the estimated temperature and
tunnelling rate $J$. c, Close up of the measured MI.

Nevertheless, several optical lattice setups are able to measure the vN
entropy, which is obtained either from a single site reduced density matrix or
as a classical thermodynamic entropy of the whole system [16, 19]. In special
cases, when the system is diagonal enough, even many-body vN entropy can be
accessed [22]. However, the measurement of vN entropy between extended spatial
subsystems has so far remained elusive, as well as the verification of the
predicted area law scaling of the MI.

In this work, we address these challenges and study the scaling of the vN
entropy and the MI with subsystem size in a continuous quantum many-body
system. Our setup is a pair of tunnelling-coupled quasi-one-dimensional (1D)
ultracold Bose gases (see Fig. 1a), cooled down and trapped below an atom chip
[23]. Along the longitudinal axis, $z$, the clouds are confined in box-like
potentials with hard walls, created by superposing magnetic and optical dipole
potentials [24]. In one of the transverse directions, the atoms are trapped in
a double-well potential, created by dressing with radiofrequency fields. The
single-particle tunnelling rate $J$ between the two condensates is adjusted by
changing the amplitude of the radiofrequency fields [25].

Figure 3: Shared information content between two spatially separated
subsystems. Here, we calculate the MI of two disjoint subsystems with the same
length, $l/L=0.15$, as a function of the distance between them, $d$. As
demonstrated in the inset, while the subsystem $A$ is kept fixed on the left
edge of the system, $B$ is shifted away to the other edge. Bullets and
shadings represent experimental and theoretical data for $I(A:B)$ respectively
(see Fig. 2 for more details on the error bars and shaded area). The solid
black line is an exponential fit with
$l_{\mathrm{fit}}=5.1\,$\mathrm{\SIUnitSymbolMicro m}$$. The finite bias $b$
results from estimating a positive quantity, MI, using finite statistics.

To prepare a state in thermal equilibrium, we directly cool down an atomic
cloud of 87Rb in a strongly coupled double-well, using standard techniques of
laser cooling and evaporative cooling [26]. The quantum fields describing each
condensate can be written in phase-density representation as
$\psi_{\mathrm{n}}(z)=\exp(i\theta_{\mathrm{n}})\sqrt{\rho_{\mathrm{n}}}$,
with $\mathrm{n}=1,2$. The spatially resolved relative phase between the two
condensates, $\varphi(z)=\theta_{1}(z)-\theta_{2}(z)$, is extracted from
interference images, taken $15.6\,$\mathrm{ms}$$ after releasing the atoms
from the trap and letting them fall freely. In the limit of low energy
excitations, the correlations of the relative phase are well described by the
sine-Gordon Hamiltonian [27]. Expanding the interaction term in the case of
high tunnelling rates leads to the massive Klein-Gordon Hamiltonian,

$H_{\mathrm{KG}}=\int_{0}^{L}\mathrm{d}z\biggl{[}g_{\mathrm{1D}}\delta\rho^{2}(z)+\frac{\hbar^{2}n_{\mathrm{1D}}}{4m}\left(\partial_{z}\varphi(z)\right)^{2}\\\
+\hbar Jn_{\mathrm{1D}}\varphi^{2}(z)\biggr{]}\,,$ (3)

where $L$ is the length of a system with uniform averaged 1D density
$n_{\mathrm{1D}}$, $g_{\mathrm{1D}}$ is the 1D interaction strength, $m$ is
the mass of an atom, and the relative density,
$\delta\rho(z)=\left[\rho_{1}(z)-\rho_{2}(z)\right]/2$, is the conjugate field
of the relative phase, fulfilling
$[\varphi(z),\delta\rho(z^{\prime})]=i\mathrm{\delta}(z-z^{\prime})$ . A
direct measurement of $\delta\rho(z)$ is unfeasible in our current
experimental setup. Hence, to reconstruct the full covariance matrix of the
initial state,

$\Gamma=\matrixquantity[Q&R\\\ R^{\mathrm{T}}&P]\,,$ (4)

we use a tomography procedure [28]. Here,
$Q_{ij}=\left<\varphi(z_{i})\varphi(z_{j})\right>$,
$P_{ij}=\left<\delta\rho(z_{i})\delta\rho(z_{j})\right>$,
$R_{ij}=\left<\frac{1}{2}\left\\{\varphi(z_{i}),\delta\rho(z_{j})\right\\}\right>$
with $i,j\in\\{1,\ldots,N\\}$, where $z_{i}$ denotes different points on a
discrete grid with $N$ points. The upper limit for $N$ is given by the
resolution of the imaging system, which limits our access to higher momentum
modes (larger than $N$) and enforces an ultra-violet (UV) cut-off.

To apply this tomographic method, we ramp up the barrier between the two
strongly coupled condensates in $2\,$\mathrm{ms}$$, and let them evolve
independently in the uncoupled double-well ($J=0$).  We then directly measure
the phase-phase correlations for different evolution times after the quench.
The post-quench dynamics follow the Tomonaga-Luttinger liquid Hamiltonian.
Over time, the initial eigenmodes of the relative density rotate into the
phase quadrature  and vice versa, enabling us to access the information about
these eigenmodes by fitting the initial second-order correlation functions of
phase-density and density-density to the measured evolution of the phase-phase
correlations in the momentum space. A thorough explanation of the
reconstruction process and its results is given in SM.

In the experiment, we prepare initial states that are thermal equilibrium
states of Klein-Gordon Hamiltonian. The quadratic form of this Hamiltonian, as
shown in (3), implies that the prepared initial states are Gaussian. After a
fast quench of $J$ to zero, the initial state evolves under another quadratic
Hamiltonian, ensuring that the state remains Gaussian under evolution [29]. To
confirm  Gaussianity, we measure the normalized averaged connected fourth-
order correlation function, $M^{(4)}$ , and show that the higher order
correlations are negligible. Note that the tomography process does not include
any restrictions regarding the Gaussianity of the reconstructed state.

Having Gaussian states significantly simplifies the calculation of the quantum
information measures - an otherwise highly non-trivial task in quantum field
theory [30]. Gaussian states are fully described by their covariance matrix
$\Gamma$. Its symplectic spectrum is obtained by diagonalizing
$i\mathcal{J}\Gamma$ where $\mathcal{J}=\matrixquantity[0&I\\\ -I&0]$ is the
symplectic unit. The symplectic spectrum consists of pairs $\pm\gamma_{n}$,
$n\in\\{1,\ldots,N\\}$. It encodes the complete information that is contained
in the covariance matrix  [14]. Consequently, it can be used to reconstruct
the full density matrix of the state and the measures of quantum information.
In particular, the vN entropy is given by

$S(\Gamma)=\sum_{n=1}^{N}\left[\left(\gamma_{n}+\tfrac{1}{2}\right)\ln\left(\gamma_{n}+\tfrac{1}{2}\right)\right.\\\
\left.-\left(\gamma_{n}-\tfrac{1}{2}\right)\ln\left(\gamma_{n}-\tfrac{1}{2}\right)\right]$
(5)

For non-Gaussian states, neglecting higher-order correlations and estimating
the vN entropy based on the covariance matrix gives a lower bound to the
actual entropy [31].

Figure 4: Temperature dependence of MI and vN entropy. The bullets represent
experimental data for $S_{A}$ (left axis) and $I(A:A^{\mathrm{C}})$ (right
axis) with $N=7$ for measurements with three different temperatures, but
otherwise similar parameters. The size of the subsystem $A$ is $l/L=0.4$ (see
Fig. 2a). The error bars show the $95\%$ confidence intervals achieved via
bootstrapping. The shaded areas represent the calculations based on the
theoretical model for different tunnelling rates, between $J=2\pi\times
0.1\,$\mathrm{Hz}$$ (dashed lines) and $2\pi\times 1\,$\mathrm{Hz}$$ (solid
lines). The inset shows the theoretical predictions for $S_{A}$ (left axis)
and $I(A:A^{\mathrm{C}})$ (right axis) for three different values of $N$, for
$J=2\pi\times 0.4\,$\mathrm{Hz}$$ without considering the effect of finite
optical resolution (see SM).

Having the reconstructed initial covariance matrix at hand, we use expression
(5) to calculate the vN entropy of any subsystem $A$, $S_{A}(\Gamma_{A})$ (see
Fig. 1c). Using equation (2), the MI between two subsystems $A$ and $B$ is
calculated. To observe the scaling of the vN entropy and the MI, we calculate
the vN entropy of subsystems with different lengths and consecutively the MI
with their complement subsystems as illustrated in Fig. 2a. As expected for
thermal states, the vN entropy is in the volume law regime, depending linearly
on the size of the subsystem (Fig. 2b). Measuring the vN entropy allows us to
study the scaling of the mutual information. We find an area law for MI, with
a plateau forming in the bulk of the system (Fig. 2c). Our results represent
an experimental verification of one of the elementary features of quantum
many-body systems [15].

We continue by studying the dependence of MI on the distance between two
subsystems. In this case, we calculate the MI of two subsystems $A$ and $B$,
separated by a gap of length $d$. The results are presented in Fig. 3. As
expected, the MI decreases as the two subsystems get further apart. We can
extract a decay length by fitting an exponential function to the experimental
data. The fitted decay length,
$l_{\mathrm{fit}}=5.1\,(3.7,\,8)\,$\mathrm{\SIUnitSymbolMicro m}$$ agrees with
the correlation length
$l_{\mathrm{C}}=6.8\,(6,\,7.7)\,$\mathrm{\SIUnitSymbolMicro m}$$ calculated
based on the experimental parameters (see SM). The intervals represent the
$95\%$ confidence intervals obtained via bootstrapping.

For typical temperatures of our experiment, the vN entropy depends linearly on
the temperature. The linear dependence changes to a logarithmic dependence by
introducing a finite UV cut-off, as shown in the inset of Fig. 4 (dashed
lines). As temperature increases, the symplectic eigenvalues grow, and the
calculated entropy using a finite number of modes saturates. The MI, however,
regardless of the UV cut-off, reaches a finite asymptotic value, given by the
classical correlations [32]. Reducing the number of modes reduces the
asymptotic value due to the limited available information in the modes taken
into account, as presented in the inset of Fig. 4 (solid lines). In Fig. 4,
measured vN entropy and MI for three different temperatures are illustrated
along with the theoretical predictions. The measurements agree with the theory
calculated for the extracted parameters and $N=7$ lowest modes.

It is important to stress that the only assumption we make to calculate the vN
entropy and the MI is that the post-quench dynamics follow a Tomonaga-
Luttinger liquid Hamiltonian, which has been confirmed in the previous work
[29]. Our measurements do not rely on any assumption related to the
Gaussianity of the initial state. We rather confirm that the initial state is
Gaussian by measuring the higher-order correlations. Even for non-Gaussian
initial states, our results would represent a lower bound to the entropy of
the full state.

The results presented here are a step towards the more ambitious goal of
measuring many-body entanglement in a continuous 1D quantum system.
Reconstruction of the full covariance matrix enables us to calculate any
entanglement measure applicable to both pure and mixed states, such as
logarithmic negativity. However, there are two main limiting factors
preventing us from detecting the entanglement. The entanglement can only be
detected if a sufficient number of momentum modes are measured whose mode
occupation numbers are close to the ground state. In our current experiments,
the nonzero temperature ($10\,$-$\,100\,$\mathrm{nK}$$) of the Bose gas keeps
the occupation numbers of the lower momentum modes too high above the value of
the ground state. At the same time, the finite optical resolution introduces a
soft cut-off which not only prevents us from measuring higher momentum modes
but also modifies the lower modes that can be measured. Improving any of these
two aspects would make the measurement of the entanglement possible.

Another promising direction for future work is to go beyond quadratic models
and detect entanglement in an interacting model. It has already been
demonstrated that the atom chip experiments can successfully simulate the
sine-Gordon model, and higher-order correlation functions can be measured
[27]. Developing a tomography procedure for this setting would give us access
to entanglement properties in interacting many-body quantum systems [33, 34].

Figure 5: Additional results for area law of MI and volume law of vN entropy.
a, Experimental results for $I(A:A^{\mathrm{C}})$, $S_{A}$, and
$S_{A^{\mathrm{C}}}$, calculated based on $N=7$ modes (left) and $N=8$ modes
(right) (see Fig. 2 for detailed explanation). The extracted parameters
including their $95\%$ confidence intervals in parenthesis are given above. b,
Close up of the measured MI.

###### Acknowledgements.

We want to thank M. Gluza, J. Eisert, I. Cirac, I. Mazets and S. Erne for
their helpful discussions. This work is supported by the DFG/FWF Research Unit
FOR 2724 ‘Thermal machines in the quantum world’, the FQXi program on
‘Information as fuel’. The Flatiron Institute is a division of the Simons
Foundation. ED acknowledges support from the ARO grant number
W911NF-20-1-0163. The work of I.K. was supported by the Max-Planck-Harvard
Research Center for Quantum Optics (MPHQ). F.C., F.M., and J. Sabino
acknowledge support from the Austrian Science Fund (FWF) in the framework of
the Doctoral School on Complex Quantum Systems (CoQuS). J. Sabino acknowledges
support by the Fundação para a Ciência e Tecnologia, Portugal
(PD/BD/128641/2017). T.S. acknowledges support from the Max Kade Foundation
through a postdoctoral fellowship. B.R. acknowledges support by the European
Union’s Horizon 2020 research and innovation programme under the Marie
Skłodowska-Curie grant agreement No. 888707. D.S. is partially supported by
AFOSR: Grant FA9550-21-1-0236. S.S. acknowledges support by the European
Union’s Horizon 2020 research and innovation programme under the Marie
Skłodowska-Curie grant agreement No. 101030988 and by the Slovenian Research
Agency (ARRS) under the grant QTE (Grant No. N1-0109).

## Supplemental Material

### Experimental realization and measurements

We realise a pair of strongly tunnelling coupled 1D superfluids by cooling
down ${}^{87}\mathrm{Rb}$ atoms in a double-well potential in an atom chip
setup. The initial state is prepared by cooling the atoms directly into a
double-well potential. The initial state in this case is a thermal equilibrium
state with typical temperatures of $30$-$120\,$\mathrm{nK}$$ and linear atomic
densities of $\approx 70\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{-1}$$ in a
box-like potential with length $\approx 50\,$\mathrm{\SIUnitSymbolMicro m}$$.

To achieve lower effective temperatures, we first cool down the atoms in a
dressed single-well potential, where the first excited state is in the vacuum
state. Slowly splitting the cloud into two parts, maps the ground state and
the first excited state to two states with a smaller energy gap and symmetric
and anti-symmetric wave functions. The resulting prethermalized state has a
lower effective temperature in the anti-symmetric modes (relative degrees of
freedom) [35]. This method has been deployed in the measurement with effective
temperature $T_{\mathrm{eff}}=14\,$\mathrm{nK}$$ presented in Fig. 4 and 5.

To probe the system, we turn off all the traps, and we let the atoms fall
freely for $15.6\,$\mathrm{ms}$$. We measure the projected 2D atomic density
distributions via absorption imaging, from which we extract the relative phase
between the two condensates for different points along the 1D direction $z$.
Due to the destructive nature of the imaging process, we repeat the
measurement hundreds of times to accumulate statistics. Thus, all the
expectation values calculated are obtained through ensemble averaging.

### Quantum field simulation using coupled quasi-one-dimensional superfluids

As it has been discussed in several earlier works (see, e.g., Refs. [27, 29,
36]), low energy excitations of coupled parallel one-dimensional gases of
weakly interacting atoms can be utilised as quantum field simulator of the
sine-Gordon Hamiltonian,

$H_{\mathrm{sG}}=\int_{0}^{L}\mathrm{d}z\biggl{[}g_{\mathrm{1D}}\delta\rho^{2}+\frac{\hbar^{2}n_{\mathrm{1D}}}{4m}\left(\partial_{z}\varphi\right)^{2}\\\
-2\hbar Jn_{\mathrm{1D}}\cos(\varphi)\biggr{]}\,.$ (6)

This model describes the relative phase, $\varphi$, and the relative density
fluctuations of two superfluids (see Fig. 1a). These two fields are canonical
conjugate of each other, i.e.
$[\varphi(z),\delta\rho(z^{\prime})]=-\mathrm{i}\delta(z-z^{\prime})$. In (6),
$m$ is the atomic mass, $n_{\mathrm{1D}}$ the uniform atomic density of the
condensates, $g_{\mathrm{1D}}$ the inter-atomic interaction, $J$ the single-
particle tunnelling rate.

For sufficiently cold gases in the strong coupling regime, i.e., when the
phase coherence length

$\displaystyle\lambda_{T}=\frac{2\hbar^{2}n_{\mathrm{1D}}}{mk_{B}T}$ (7)

is larger than the healing length of the relative phase (correlation length)

$\displaystyle l_{C}=\sqrt{\frac{\hbar}{4mJ}}\ ,$ (8)

the cosine term in (6) can be expanded to second order and be approximated by
the quadratic Klein-Gordon model. Introducing the sound velocity $c$,
Luttinger parameter $K$ and the Klein-Gordon quasi-particle mass, $M$, given
in terms of the microscopic parameters by

$\displaystyle c$
$\displaystyle=\sqrt{\frac{g_{\mathrm{1D}}n_{\mathrm{1D}}}{m}}\ ,$ (9)
$\displaystyle K$
$\displaystyle=\frac{\hbar\pi}{2}\sqrt{\frac{n_{\mathrm{1D}}}{mg_{\mathrm{1D}}}}\
,$ (10) $\displaystyle M$ $\displaystyle=2m\sqrt{\frac{\hbar
J}{g_{\mathrm{1D}}n_{\mathrm{1D}}}}\ ,$ (11)

the Klein-Gordon (KG) Hamiltonian can be written as

$H_{\mathrm{KG}}=\frac{\hbar
c}{2}\int_{0}^{L}\mathrm{d}z\,\biggl{[}\frac{\pi}{K}\delta\rho^{2}(z)+\frac{K}{\pi}\left(\partial_{z}\varphi(z)\right)^{2}\biggr{]}\\\
+\frac{M^{2}c^{4}}{2\hbar
c}\int_{0}^{L}\mathrm{d}z\;\frac{K}{\pi}\varphi^{2}(z)\ .$ (12)

Note that the first two terms in (12) are the Tomonaga-Luttinger liquid (TLL)
Hamiltonian. In table 1, the relevant parameters for each measurement are
listed, where $\omega_{M}=Mc^{2}/\hbar$ is the KG mass in units of angular
frequency.

### Reconstruction of the initial full covariance matrix

To extract the full covariance matrix that characterises the state of the
system, we use the quantum tomography method developed in [28]. Given that
only one of the two canonical variables (the phase) is accessible through
experimental measurements, its canonically conjugate variable (the density
fluctuations) can be accessed indirectly by letting the system evolve under a
harmonic Hamiltonian with known mode frequencies and measuring the phase at
different times. For each of the harmonic modes, the dynamics correspond to a
rotation in phase space, so that over time, the initial density variance turns
into phase variance and vice versa. Thus, we can fully reconstruct the initial
covariance matrix from phase measurements at a sufficiently large number of
different times. This reconstruction is done by first going to Fourier space,
where the modes evolve independently, and then fitting the data for the time
evolution of the phase covariance of each mode to the known functions
expressing this data in terms of the initial phase and density correlations.

However, one mode needs special treatment: The zero mode, i.e. the mode
corresponding to zero momentum, therefore zero energy. This mode does not
rotate in phase space; instead, it moves at a constant velocity. This means
that the zero mode phase variance does not oscillate in time but grows as a
quadratic function of time instead, an effect known as phase diffusion [37].
Moreover, because of the compactified nature of the phase-field, which means
that phases differing by $2\pi$ should be considered identical, even though
the phase grows in time with no bound, measurements can only observe its
growth within the interval $[-\pi,+\pi]$. Therefore, the estimation of the
initial phase and density correlations of the zero mode is different from that
of the others. For the zero mode, we fit a quadratic function instead of an
oscillatory function of time, and we restrict the fit to times before reaching
the upper bound due to compactification.

More specifically, the dynamics are chosen to follow the Tomonaga-Luttinger
liquid model with Hamiltonian

$H_{\mathrm{TLL}}=\int_{0}^{L}dz\left[g_{\mathrm{1D}}(\delta\rho(z))^{2}+\frac{\hbar^{2}n_{1D}}{4m}\left(\partial_{z}\varphi(z)\right)^{2}\right]$
(13)

For a hard-wall box trap, the vanishing of the particle current at the edges
of the system means that the effective boundary conditions are of Neumann type
$\partial_{z}\varphi(x=0)=\partial_{z}\varphi(x=L)=0$. In this case, using
cosine eigenfunctions

$\displaystyle f_{n}^{\varphi}(z)$
$\displaystyle=\begin{cases}2\left(n\hbar\pi\sqrt{\frac{n_{1D}}{g_{\mathrm{1D}}m}}\right)^{-1/2}\cos(n\frac{\pi}{L}z)\,,&n>0\\\
1\,,&n=0\end{cases}$ $\displaystyle f_{n}^{\delta\rho}(z)$
$\displaystyle=\begin{cases}-\frac{1}{L}\left(n\hbar\pi\sqrt{\frac{n_{1D}}{g_{\mathrm{1D}}m}}\right)^{1/2}\cos(n\frac{\pi}{L}z)\,,&n>0\\\
-\frac{1}{L}\,,&n=0\end{cases}\,,$ (14)

the Hamiltonian can be diagonalised in terms of cosine Fourier modes:

$\displaystyle H$ $\displaystyle=\frac{\hbar
u}{2}\delta\rho_{0}^{2}+\sum_{n=1}^{\infty}\frac{\hbar\omega_{n}}{2}\left[\delta\rho_{n}^{2}+\varphi_{n}^{2}\right]\,,$
(15)

with

$\displaystyle u$ $\displaystyle=2g_{\mathrm{1D}}/\hbar L$ (16)
$\displaystyle\omega_{n}$ $\displaystyle=ck_{n}\,.$ (17)

Furthermore:

$\displaystyle k_{n}$ $\displaystyle=n\pi/L$ (18) $\displaystyle c$
$\displaystyle=\sqrt{g_{\mathrm{1D}}n_{\mathrm{1D}}/m}$ (19) $\displaystyle
g_{\mathrm{1D}}$
$\displaystyle=\hbar\omega_{\perp}a_{\mathrm{s}}\frac{2+3a_{\mathrm{s}}n_{\mathrm{1D}}}{\left(1+2a_{\mathrm{s}}n_{\mathrm{1D}}\right)}\,.$
(20)

Here, $c$ is the speed of sound and $g_{\mathrm{1D}}$ the density-broadened 1D
interaction strength [38], $a_{\mathrm{s}}=5.2\,$\mathrm{nm}$$ is the three-
dimensional scattering length [39] and $m=1.44\times 10^{-25}\,$\mathrm{kg}$$
the mass of a ${}^{87}\mathrm{Rb}$ atom.

As mentioned before and reflected in the TLL Hamiltonian in Fourier space (Eq.
15), we have to treat the zero mode separately. We begin with the harmonic
part of the Hamiltonian ($n>0$) where the time evolution of the modes is given
by

$\displaystyle\delta\rho_{n}(t)$
$\displaystyle=\delta\rho_{n}(0)\cos(\omega_{n}t)+\varphi_{0}(0)\sin(\omega_{n}t)$
(21) $\displaystyle\varphi_{n}(t)$
$\displaystyle=\varphi_{0}(0)\cos(\omega_{n}t)-\delta\rho_{n}(0)\sin(\omega_{n}t)\,,$
(22)

In the experiment, we use matter interferometry to measure spatially resolved
relative phase between two superfluids, from which, a referenced second-order
correlation is calculated for each time step:

$\Phi^{2}_{ab}(t)=\left<\bigl{(}\varphi(z_{a},t)-\varphi(z_{0},t)\bigr{)}\bigl{(}\varphi(z_{b},t)-\varphi(z_{0},t)\bigr{)}\right>\,.$
(23)

Note that subtracting the phase of an arbitrary reference position $z_{0}$
will only remove the zero mode and does not affect any of the higher modes.

Expanding $\Phi^{2}$ with eigenfunctions (Reconstruction of the initial full
covariance matrix), gives

$\Phi^{2}_{ab}(t)=\sum^{N}_{j,k=1}f^{a,b}_{j,k}\left<\varphi_{j}(t)\varphi_{k}(t)\right>$
(24)

where

$f^{a,b}_{j,k}=\bigl{(}f^{\phi}_{j}(z_{a})-f^{\phi}_{j}(z_{0})\bigr{)}\bigl{(}f^{\phi}_{k}(z_{b})-f^{\phi}_{k}(z_{0})\bigr{)}\,.$
(25)

Using the equation of motion (22), and defining
$\tilde{Q}_{jk}=\left<\varphi_{j}(0)\varphi_{k}(0)\right>$,
$\tilde{R}_{jk}=\left<\frac{1}{2}\left\\{\varphi_{j}(0),\delta\rho_{k}(0)\right\\}\right>$,
and $\tilde{P}_{jk}=\left<\delta\rho_{j}(0)\delta\rho_{k}(0)\right>$, we
obtain

$\begin{split}\Phi^{2}_{ab}(t)&=\sum^{N}_{j,k=1}f^{a,b}_{j,k}\cos(\omega_{j}t)\cos(\omega_{k}t)\tilde{Q}_{jk}\\\
&+\sum^{N}_{j,k=1}(-f^{a,b}_{j,k}-f^{a,b}_{j,k})\cos(\omega_{j}t)\sin(\omega_{k}t)\tilde{R}_{jk}\\\
&+\sum^{N}_{j,k=1}f^{a,b}_{j,k}\sin(\omega_{j}t)\sin(\omega_{k}t)\tilde{P}_{jk}\,.\end{split}$
(26)

Equation (26) stands in the heart of the tomography process: The goal is to
find the elements of $\tilde{Q}$, $\tilde{R}$, and $\tilde{P}$ using an
optimization process. Note that the left hand side is calculated using
measured phase profiles in the experiment and in the right hand side,
$f^{a,b}_{j,k}$ and $\omega_{j}$ can be calculated from the experimental
parameters, as shown in (Reconstruction of the initial full covariance matrix)
and (17) respectively. Results are presented in SM. For a more detailed
explanation of the tomography process, please refer to [28].

Unlike all higher modes, the zero mode, which corresponds to the constant in
space eigenfunction, is not a harmonic oscillator mode of the Tomonaga-
Luttinger liquid Hamiltonian (equation (13)). This is because only one of the
canonical variables, $\delta\rho_{0}$, is present in the Hamiltonian for the
zero mode. As a result, the time evolution of the zero mode is given by

$\displaystyle\delta\rho_{0}(t)$
$\displaystyle=\delta\rho_{0}(0)=\mathrm{const},$ (27)
$\displaystyle\varphi_{0}(t)$
$\displaystyle=-u\delta\rho_{0}(0)t+\varphi_{0}(0)\,,$ (28)

which means that the phase variance grows with time as

$\left\langle\varphi_{0}(t)^{2}\right\rangle=\left\langle\varphi_{0}^{2}\right\rangle_{t=0}-\left\langle\left\\{\varphi_{0},\delta\rho_{0}\right\\}\right\rangle_{t=0}ut+\left\langle\delta\rho_{0}^{2}\right\rangle_{t=0}u^{2}t^{2}$
(29)

However, because of the compactified nature of the phase-field, its zero mode
component $\varphi_{0}$ is not a well-defined, measurable operator. Only
imaginary exponentials of the form $e^{{\rm i}n\varphi_{0}}$ for integer $n$
are well-defined. Nevertheless, under the assumption the initial state is
Gaussian in terms of the zero mode (as also for all other modes too) and given
that it remains Gaussian under the dynamics following from $H_{\mathrm{TLL}}$,
we can derive the zero mode variance from the mean value of $e^{{\rm
i}n\varphi_{0}}$ using the cumulant expansion formula for the special case of
Gaussian random variables

$\left\langle\exp\left(\mathrm{\mathrm{i}}\varphi_{0}\right)\right\rangle=\exp\left(\mathrm{\mathrm{i}}\left\langle\varphi_{0}\right\rangle-\left\langle\varphi_{0}^{2}\right\rangle/2\right)$
(30)

From the above, we find

$\left\langle\varphi_{0}^{2}\right\rangle=-2\log\lvert\left\langle\exp\left(\mathrm{\mathrm{i}}\varphi_{0}\right)\right\rangle\rvert$
(31)

Therefore, in order to extract the zero mode part of the covariance matrix
$\left\langle\varphi_{0}^{2}\right\rangle_{t=0},\left\langle\delta\rho_{0}^{2}\right\rangle_{t=0}$
and
$\left\langle\left\\{\varphi_{0},\delta\rho_{0}\right\\}\right\rangle_{t=0}$
in the initial state we calculate the zero mode variance of the phase at each
time from (30) and fit this with the theoretical formula (31).

Having the covariance matrix in the Fourier space for the first $N$ modes, we
use a discrete Fourier transformation based on the eigenfunctions (equation
(Reconstruction of the initial full covariance matrix)) to calculate the
covariance matrix $\Gamma$ in the real space. We chose the cut-off based on
the reconstructed occupation numbers. We only take into account modes with
physical (positive) occupation numbers. This covariance matrix is used to
calculate vN entropy and MI, as discussed in the main text and Fig. 1C.

### Covariance matrix of Klein-Gordon model in thermal equilibrium

The theory predictions in our work are calculated based on the covariance
matrix of the thermal equilibrium states of Klein-Gordon (KG) model (equation
3), which is given by [40]:

$\displaystyle Q_{ij}$ $\displaystyle=\frac{\pi}{2KL}\frac{\hbar
c}{Mc^{2}}\coth\left(\frac{Mc^{2}}{2k_{B}T}\right)$
$\displaystyle+\frac{\pi}{KL}\sum_{n=1}^{N}\frac{\hbar
c}{\epsilon_{n}}\coth\left(\frac{\epsilon_{n}}{2k_{B}T}\right)\cos(k_{n}z_{i})\cos(k_{n}z_{j})$
(32) $\displaystyle R_{ij}$ $\displaystyle=0$ (33) $\displaystyle P_{ij}$
$\displaystyle=\frac{K}{2\pi L}\frac{Mc^{2}}{\hbar
c}\coth\left(\frac{Mc^{2}}{2k_{B}T}\right)$ $\displaystyle+$
$\displaystyle\frac{K}{\pi L}\sum_{n=1}^{N}\frac{\epsilon_{n}}{\hbar
c}\coth\left(\frac{\epsilon_{n}}{2k_{B}T}\right)\cos(k_{n}z_{i})\cos(k_{n}z_{j})\,,$
(34)

with the dispersion relation
$\epsilon_{n}=\sqrt{\hbar^{2}k_{n}^{2}c^{2}+M^{2}c^{4}}$. Here, $M$ is the KG
mass, $L$ the system size, $K$ the Luttinger parameter, $T$ the temperature,
$N$ the UV cut-off and the rest of the parameters are defined before. The next
section will explain how $M$ and $T$ are estimated based on the measured data.
To include the effect of the finite imaging resolution, we convolve the
theoretical calculations with a Gaussian point-spread function with a standard
deviation $\sigma_{\mathrm{PSF}}\approx 3\,$\mathrm{\SIUnitSymbolMicro m}$$
[38] (see next section).

### Estimation of the temperature and Klein-Gordon mass

In order to compute theoretical predictions for the mutual information in the
initial state based on the assumption that these are the thermal states of the
KG model, we need to estimate two effective parameters, the mass and the
temperature. We do this by fitting the results of the tomographic
reconstruction for the mode variances to those corresponding to KG thermal
states. Given that the modes are decoupled from each other both initially and
throughout the tomography dynamics, the estimation of the KG mode frequency
can be done independently for each mode. Having estimated the mode
frequencies, we can then verify if they follow the theoretical dispersion
relation of the KG model and extract the corresponding mass parameter by a
fit.

The relation between the post-quench quadratures and the initial (pre-quench)
KG state with mode occupation number $N_{0n}$ is given by

$\displaystyle\left\langle\varphi_{n}^{2}\right\rangle$
$\displaystyle=\frac{\omega_{n}}{\omega_{0n}}\left(N_{0n}+\frac{1}{2}\right)$
(35) $\displaystyle\left\langle\delta\rho_{n}^{2}\right\rangle$
$\displaystyle=\frac{\omega_{0n}}{\omega_{n}}\left(N_{0n}+\frac{1}{2}\right)$
(36)

where $\omega_{0n}$ and $\omega_{n}$ are the pre- and post-quench mode
frequencies respectively. Both $\left\langle\varphi_{n}^{2}\right\rangle$ and
$\left\langle\delta\rho_{n}^{2}\right\rangle$ are achieved via tomography and
$\omega_{n}$ is given by the equation (17). From the equations (35) and (36)
we can calculate both $\omega_{0n}$ and $N_{0n}$:

$\displaystyle N_{0n}$
$\displaystyle=\sqrt{\left\langle\varphi_{n}^{2}\right\rangle\left\langle\delta\rho_{n}^{2}\right\rangle}-\frac{1}{2}$
(37) $\displaystyle\omega_{0n}$
$\displaystyle=\omega_{n}\sqrt{\frac{\left\langle\delta\rho_{n}^{2}\right\rangle}{\left\langle\varphi_{n}^{2}\right\rangle}}\,.$
(38)

Assuming the initial state is thermal, we use the following fit function to
extract the temperature, $T=(\beta k_{\mathrm{B}})^{-1}$:

$\displaystyle N_{0n}^{\mathrm{fit}}$
$\displaystyle=\mathrm{exp}\left({-k_{n}^{2}\sigma_{\mathrm{PSF}}^{2}/2}\right)(\frac{1}{\mathrm{exp}(\hbar\omega_{0n}\beta)-1}+\frac{1}{2})-\frac{1}{2}\,.$
(39)

Note that we also introduce the effect of the imaging system by multiplying
the modes with a Gaussian point spread function with width
$\sigma_{\mathrm{PSF}}$ which corresponds to a convolution with a Gaussian
point-spread function in the real space [38]. To extract the KG mass, we fit
the KG dispersion relation,

$\omega_{0n}^{\mathrm{fit}}=\sqrt{c^{2}k_{n}^{2}+M^{2}c^{4}/\hbar^{2}}\,,$
(40)

to the calculated $\omega_{0n}$.

### Tomography results

Figure 6: Tomography results for three measurements from lowest to highest temperature (M1 to M3). For each measurement, the elements of the covariance matrix in Fourier and position space are presented. On the right, the diagonal elements of the covariance matrix in the Fourier space are plotted. The error bars represent 95% confidence intervals obtained via bootstrapping. Solid lines are the theory calculations corresponding to the extracted temperature and KG mass for each measurement. Measurement | $T\,($\mathrm{n}\mathrm{K}$)$ | $\lambda_{T}\,($\mathrm{\SIUnitSymbolMicro m}$)$ | $\omega_{M}/2\pi\,($\mathrm{Hz}$)$ | $l_{\mathrm{C}}\,($\mathrm{\SIUnitSymbolMicro m}$)$ | $q=\lambda_{T}/\ell_{J}$ | $r=\hbar\omega_{M}/k_{\mathrm{B}}T$
---|---|---|---|---|---|---
M1 | 14 | 37 | 34 | 7.3 | 5 | 0.11
M2 | 48 | 18 | 43 | 6.8 | 2.6 | 0.043
M3 | 94 | 9.5 | 40 | 7.4 | 1.3 | 0.021
Table 1: List of relevant parameters for the three measurements. $T$ is the
temperature, $\lambda_{T}$ the thermal coherence length,
$\omega_{M}=Mc^{2}/\hbar$ the KG mass in angular frequency units, and $l_{C}$
the healing length of the phase.  Figure 7: Extraction of temperature and KG
mass from the mode occupation number and dispersion relation a, For each
measurement $N_{0n}$ is plotted along with a fit curve (solid line),
$N_{0n}^{\mathrm{fit}}$. b, Dispersion relation for three different
measurements are plotted. For every quantity, the error bars represent 95%
confidence intervals obtained via bootstrapping.

The results of the reconstruction process for the three measurements are
presented in Fig. 6 and 7 and the corresponding parameters are written in
Table 1. We show the covariance matrix elements obtained after the tomography
in the Fourier space, as well as the elements in the position space after a
discrete inverse Fourier transform. Note that the only constraint on the
elements of $\tilde{\Gamma}$ is that the resulting state is physical i.e. the
Heisenberg uncertainty relation holds [28]. We do not restrict the state by
any further assumption. For example, we do not assume that the state is
Gaussian (see next section) or $\tilde{R}$ is small. The only assumption in
the tomography process is that the post-quench dynamics follow the Tomonaga-
Luttinger liquid Hamiltonian for short times as mentioned before.

### Fourth order correlation

In order to show that the reconstructed state is Gaussian, we measure the
normalised averaged connected fourth-order correlation function, $M^{(4)}$.
Please refer to previous extensive works by our group [29, 41, 38]. The
results for the four measurements are shown in Fig. 8.

Figure 8: Fourth order correlations The error bars represent 95% confidence
intervals obtained via bootstrapping. The shaded area is the typically
calculated bias [41].

### Effect of finite imaging resolution

Figure 9: Imaging effect on mutual information. Mutual information for $l=L/2$
for different cut-off, $N$. The error bars represent 95% confidence interval.
The shaded area is the theoretical calculation and the solid black line is the
residual MI due to the effect of finite imaging resolution. Figure 10:
Dependence of $l_{\mathrm{fit}}$ on the total number of modes, $N$. a, For
different $N$, the same quantities are plotted as in Fig. 3 in the main text
but here, we shifted the theory calculation by $b$, to represent better the
agreement between theory, experimental data and the fit. Note that changing
the total number of modes consequently changes the size of the sub-region. b,
The resulting decay length, $l_{\mathrm{fit}}$ is plotted for different $N$.
the solid part of the error bars represents 68% confidence interval (standard
error) and the extended transparent part 95% confidence interval. The solid
orange line and the shaded area show the $l_{\mathrm{C}}=6.8(6.0,7.7)$ and its
95% confidence interval.

We investigate the effect of the finite imaging resolution and the ultra-
violet cut-off on the measurement of mutual information (MI). The Gaussian
point-spread function of the imaging system acts as a soft cut-off. This cut-
off not only limits the number of observable modes but also affects the lower
modes that can be measured [38]. In length scales close to and smaller than
the imaging resolution, the blurring effect creates correlations, producing
unphysical MI. We theoretically quantify this ”fake” MI by considering a state
with high temperature, and short correlation length (much shorter than the
imaging resolution), whose MI is practically zero ($\sim 10^{-6}$). For this
state, we calculate the MI of the half-system, $l=L/2$, for different numbers
of modes, $N$. The result is presented in Fig. 9 with solid black lines along
with the experimental data. The error bars are $95\%$ confidence intervals
obtained via bootstrapping and the shaded area is the $95\%$ confidence
interval for the theoretical calculation. Note that the calculated residual MI
only depends on the imaging resolution. The figure suggests that our measured
MI is significantly higher than what is calculated to be the residual MI due
to the imaging effect.

We further investigate the effect of the hard cut-off on the decay of MI as
discussed in the main text (Fig. 3). In Fig. 10a, for each $N$, the
information similar to Fig. 2 in the main text is shown. In Fig. 10b, the
dependence of $l_{\mathrm{fit}}$ on the total number of modes, $N$, is
plotted. We observe that for all values of $N$, experimental data and the
theoretical calculations are in good agreement. As the number of modes
decreases, sub-region size and spacing increase, which lowers the quality of
fit, resulting in large error bars. Nevertheless, the decay length stays
consistent for different cut-offs. Note that for $N<5$ the sub-regions become
too large and reliable fitting is no longer possible.

## References

  * Horodecki  _et al._ [2009] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Rev. Mod. Phys. 81, 865 (2009).
  * Orús [2014] R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Annals of Physics 349, 117 (2014).
  * Schuch  _et al._ [2008] N. Schuch, M. M. Wolf, F. Verstraete, and J. I. Cirac, Entropy scaling and simulability by matrix product states, Phys. Rev. Lett. 100, 030504 (2008).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 91, 021001 (2019).
  * Senthil [2015] T. Senthil, Symmetry-protected topological phases of quantum matter, Annual Review of Condensed Matter Physics 6, 299 (2015).
  * Callan and Wilczek [1994] C. Callan and F. Wilczek, On geometric entropy, Physics Letters B 333, 55 (1994).
  * Holzhey  _et al._ [1994] C. Holzhey, F. Larsen, and F. Wilczek, Geometric and renormalized entropy in conformal field theory, Nuclear Physics B 424, 443 (1994).
  * Srednicki [1993] M. Srednicki, Entropy and area, Physical Review Letters 71, 666 (1993).
  * Calabrese and Cardy [2004] P. Calabrese and J. Cardy, Entanglement entropy and quantum field theory, Journal of Statistical Mechanics: Theory and Experiment 2004, P06002 (2004).
  * Calabrese and Cardy [2009] P. Calabrese and J. Cardy, Entanglement entropy and conformal field theory, Journal of Physics A: Mathematical and Theoretical 42, 504005 (2009).
  * Polchinski [2017] J. Polchinski, The Black Hole Information Problem, in  _New Frontiers in Fields and Strings_ (WORLD SCIENTIFIC, Boulder, Colorado, 2017) pp. 353–397.
  * Ryu and Takayanagi [2006] S. Ryu and T. Takayanagi, Holographic Derivation of Entanglement Entropy from the anti–de Sitter Space/Conformal Field Theory Correspondence, Physical Review Letters 96, 181602 (2006).
  * Modi  _et al._ [2012] K. Modi, A. Brodutch, H. Cable, T. Paterek, and V. Vedral, The classical-quantum boundary for correlations: Discord and related measures, Rev. Mod. Phys. 84, 1655 (2012).
  * Eisert  _et al._ [2010] J. Eisert, M. Cramer, and M. B. Plenio, Colloquium: Area laws for the entanglement entropy, Reviews of Modern Physics 82, 277 (2010).
  * Wolf  _et al._ [2008] M. M. Wolf, F. Verstraete, M. B. Hastings, and J. I. Cirac, Area laws in quantum systems: Mutual information and correlations, Phys. Rev. Lett. 100, 070502 (2008).
  * Luo  _et al._ [2007] L. Luo, B. Clancy, J. Joseph, J. Kinast, and J. E. Thomas, Measurement of the entropy and critical temperature of a strongly interacting Fermi gas, Phys. Rev. Lett. 98, 080402 (2007).
  * Islam  _et al._ [2015] R. Islam, R. Ma, P. M. Preiss, M. Eric Tai, A. Lukin, M. Rispoli, and M. Greiner, Measuring entanglement entropy in a quantum many-body system, Nature 528, 77 (2015).
  * Omran  _et al._ [2015] A. Omran, M. Boll, T. A. Hilker, K. Kleinlein, G. Salomon, I. Bloch, and C. Gross, Microscopic observation of Pauli blocking in degenerate fermionic lattice gases, Phys. Rev. Lett. 115, 263001 (2015).
  * Cocchi  _et al._ [2017] E. Cocchi, L. A. Miller, J. H. Drewes, C. F. Chan, D. Pertot, F. Brennecke, and M. Köhl, Measuring entropy and short-range correlations in the two-dimensional Hubbard model, Phys. Rev. X 7, 031025 (2017).
  * 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).
  * Efron and Tibshirani [1986] B. Efron and R. Tibshirani, Bootstrap methods for standard errors, confidence intervals, and other measures of statistical accuracy, Statistical Science , 54 (1986).
  * Evrard  _et al._ [2021] B. Evrard, A. Qu, J. Dalibard, and F. Gerbier, Observation of fragmentation of a spinor Bose-Einstein condensate, Science 373, 1340 (2021).
  * Folman  _et al._ [2000] R. Folman, P. Krüger, D. Cassettari, B. Hessmo, T. Maier, and J. Schmiedmayer, Controlling Cold Atoms using Nanofabricated Surfaces: Atom Chips, Phys. Rev. Lett. 84, 4749 (2000).
  * Tajik  _et al._ [2019] M. Tajik, B. Rauer, T. Schweigler, F. Cataldini, J. Sabino, F. S. Møller, S.-C. Ji, I. E. Mazets, and J. Schmiedmayer, Designing arbitrary one-dimensional potentials on an atom chip, Optics Exp. 27, 33474 (2019).
  * Schumm  _et al._ [2005] T. Schumm, S. Hofferberth, L. M. Andersson, S. Wildermuth, S. Groth, I. Bar-Joseph, J. Schmiedmayer, and P. Krüger, Matter-wave interferometry in a double well on an atom chip, Nature Physics 1, 57 (2005).
  * Rauer  _et al._ [2016] B. Rauer, P. Grišins, I. E. Mazets, T. Schweigler, W. Rohringer, R. Geiger, T. Langen, and J. Schmiedmayer, Cooling of a One-Dimensional Bose Gas, Phys. Rev. Lett. 116, 030402 (2016).
  * Schweigler  _et al._ [2017] T. Schweigler, V. Kasper, S. Erne, I. Mazets, B. Rauer, F. Cataldini, T. Langen, T. Gasenzer, J. Berges, and J. Schmiedmayer, Experimental characterization of a quantum many-body system via higher-order correlations, Nature 545, 323 (2017).
  * Gluza  _et al._ [2020] M. Gluza, T. Schweigler, B. Rauer, C. Krumnow, J. Schmiedmayer, and J. Eisert, Quantum read-out for cold atomic quantum simulators, Comm. Phys. 3, 12 (2020).
  * Schweigler  _et al._ [2021] T. Schweigler, M. Gluza, M. Tajik, S. Sotiriadis, F. Cataldini, S.-C. Ji, F. S. Møller, J. Sabino, B. Rauer, J. Eisert, and J. Schmiedmayer, Decay and recurrence of non-Gaussian correlations in a quantum many-body system, Nature Phys. 17, 559 (2021).
  * Serafini [2017] A. Serafini,  _Quantum Continuous Variables: A Primer of Theoretical Methods_ (CRC Press, Taylor & Francis Group, Boca Raton, 2017).
  * Eisert and Wolf [2007] J. Eisert and M. M. Wolf, Gaussian quantum channels, in  _Quantum Information with Continuous Variables of Atoms and Light_ (PUBLISHED BY IMPERIAL COLLEGE PRESS AND DISTRIBUTED BY WORLD SCIENTIFIC PUBLISHING CO., 2007) pp. 23–42.
  * Katsinis and Pastras [2020] D. Katsinis and G. Pastras, An inverse mass expansion for the mutual information in free scalar qft at finite temperature, Journal of High Energy Physics 2020, 91 (2020).
  * Kukuljan  _et al._ [2018] I. Kukuljan, S. Sotiriadis, and G. Takacs, Correlation functions of the quantum sine-Gordon model in and out of equilibrium, Phys. Rev. Lett. 121, 110402 (2018).
  * Emonts and Kukuljan [2022] P. Emonts and I. Kukuljan, Reduced density matrix and entanglement of interacting quantum field theories with hamiltonian truncation (2022).
  * Gring  _et al._ [2012] M. Gring, M. Kuhnert, T. Langen, T. Kitagawa, B. Rauer, M. Schreitl, I. Mazets, D. A. Smith, E. Demler, and J. Schmiedmayer, Relaxation and prethermalization in an isolated quantum system, Science 337, 1318 (2012).
  * Rauer  _et al._ [2018] B. Rauer, S. Erne, T. Schweigler, F. Cataldini, M. Tajik, and J. Schmiedmayer, Recurrences in an isolated quantum many-body system, Science 360, 307 (2018).
  * Berrada  _et al._ [2013] T. Berrada, S. van Frank, R. Bücker, T. Schumm, J.-F. Schaff, and J. Schmiedmayer, Integrated Mach–Zehnder interferometer for Bose–Einstein condensates, Nature Communications 4 (2013).
  * Schweigler [2019] T. Schweigler,  _Correlations and dynamics of tunnel-coupled one-dimensional Bose gases_ , Ph.D. thesis, Techniche Universität Wien (2019).
  * van Kempen  _et al._ [2002] E. G. van Kempen, S. J. Kokkelmans, D. J. Heinzen, and B. J. Verhaar, Interisotope Determination of Ultracold Rubidium Interactions from Three High-Precision Experiments, Physical Review Letters 88, 093201 (2002).
  * Le Bellac [1996] M. Le Bellac,  _Thermal field theory_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 1996).
  * Gluza  _et al._ [2022] M. Gluza, T. Schweigler, M. Tajik, J. Sabino, F. Cataldini, F. S. Møller, S.-C. Ji, B. Rauer, J. Schmiedmayer, J. Eisert, and S. Sotiriadis, Mechanisms for the emergence of Gaussian correlations, SciPost Phys. 12, 113 (2022).