# Generation of scalable many-body Bell correlations in spin chains
with short-range two-body interactions

Marcin Płodzień  ICFO - Institut de Ciencies Fotoniques, The Barcelona
Institute of Science and Technology, 08860 Castelldefels, Barcelona, Spain
<EMAIL_ADDRESS>Tomasz Wasak  Institute of Physics, Faculty of
Physics, Astronomy and Informatics, Nicolaus Copernicus University in Toruń,
Grudzia̧dzka 5, 87-100 Toruń, Poland  Emilia Witkowska  Institute of Physics
PAS, Aleja Lotnikow 32/46, 02-668 Warszawa, Poland  Maciej Lewenstein  ICFO -
Institut de Ciencies Fotoniques, The Barcelona Institute of Science and
Technology, 08860 Castelldefels, Barcelona, Spain  ICREA, Pg. Lluís Companys
23, 08010 Barcelona, Spain  Jan Chwedeńczuk  Faculty of Physics, University of
Warsaw, ul. Pasteura 5, PL-02-093 Warsaw, Poland<EMAIL_ADDRESS>
###### Abstract

Dynamical generation of strong and scalable quantum resources, like many-body
entanglement and Bell correlations, in spin-$1/2$ chains is possible with all-
to-all interactions, either for constant interaction strength realizing one-
axis twisting protocol or for power-law decaying potentials. We show, however,
that such quantum resources can also be dynamically generated with a finite
range of interactions. We identify a necessary critical range and indicate a
critical time when scalable quantum correlations appear. Finally, we show that
the certification of generated states is accessible in the modern quantum
simulator platforms.

Introduction— The potential for future quantum technologies is fuelled by
quantum resources which are the quantum coherence and the triad of many-body
non-classical correlations: entanglement [horodecki2009entanglement],
Einstein-Podolsky-Rosen (EPR) steering [uola2020steering] and the Bell
nonlocality [brunner2014bell]. These are crucial for quantum computing,
cryptography, communication and metrology—the key pillars of the second
quantum revolution [Acin_2018, Eisert2020, Fraxanet2022]. As such, the main
goal of quantum technologies in the next years is generation,
characterization, storage, and certification of many-body quantum states
[Frerot_2023].

The archetypal model allowing generation of scalable quantum resources is a
spin-$1/2$ chain undergoing one-axis twisting (OAT) [kitagawa1993squeezed,
wineland1994squeezed]. In the OAT protocol, the dynamics is governed by a non-
linear Hamiltonian with all-to-all infinite range spin couplings, like
$\hat{S}^{2}_{z}$, where $\hat{S}_{z}$ is a collective spin operator along
$z$-axis, while initially spins are polarized along a direction perpendicular
to the $z$-axis. The OAT protocol dynamically generates quantum states such as
spin-squeezed states useful for high-precision metrology limit
[RevModPhys.90.035005, PhysRevLett.105.053601, PhysRevA.50.67,
MllerRigat2023], many-body entangled and the many-body Bell correlated states
[Tura1256, schmied2016bell, Aloy2019, Baccari2019, Tura2019,
PRXQuantum.2.030329, zukowski2002bell, cavalcanti2007bell, he2011entanglement,
cavalcanti2011unified, spiny.milosz, PhysRevLett.126.210506,
10.21468/SciPostPhysCore.5.2.025, PRXQuantum.2.030329, PhysRevA.46.R6797,
PhysRevA.47.5138, PhysRevA.92.043622, PhysRevLett.100.210401, Li2009,
PhysRevA.96.013823, Kajtoch-sc-2018, Schulte2020, PhysRevLett.126.160402,
PhysRevA.105.022625].

The OAT dynamics can be realized with a variety of quantum simulating
platforms, like ultracold systems utilizing atom-atom collisions
[Treutlein2010, Oberthaler2010, Chapman2012, PhysRevLett.125.033401], atom-
light interactions [PhysRevLett.104.073602, PhysRevLett.105.080403], in
Rydberg atoms [doi:10.1126/science.aax9743, bornet2023scalable, Eckner2023],
array of trapped ions [Bohnet1297, franke2023quantumenhanced], and
superconducting qubits [PhysRevLett.119.180511, doi:10.1126/science.aay0600,
doi:10.1126/sciadv.aba4935]. The recent theoretical proposals for the OAT
simulation with ultracold atoms in optical lattices effectively simulate
Hubbard and Heisenberg models [Kajtoch2018, PhysRevResearch.1.033075,
Plodzien2020, Plodzien2022, Plodzie2023generation, PhysRevResearch.3.013178,
Hernandez2022, Dziurawiec2023, yanes2023spin]. All the above mentioned
proposals are based on effectively induced all-to-all infinite-range
couplings, either with constant amplitude as required by OAT, or with power-
law decaying amplitude like dipolar or van der Waals interactions.

Here, we show that with short-range interactions, scalable many-body
entanglement and Bell correlations can be generated with spin chains with
Ising-type Hamiltonian. In contrast to the long-range interactions and all-to-
all couplings, in our system, the number of connected spins by the
interaction, i.e., the number of two-qubit entangling operations, increases
linearly with the system size for a fixed interaction range, which is an
advantage compared to the former systems, which require quadratic scaling of
the number of operations for the protocol.

Figure 1: Visualization of the spin-connectivity in a chain of $N$ spins for
different interaction ranges $r$. The three panels show nearest-neighbors
interaction with the range $r=1$ (left), finite-range interaction with $r=4$
(middle), and the all-to-all type of coupling with $r=N-1$, present in one-
axis twisting-like protocols (right).

We find a critical range allowing many-body entanglement generation,
independent of the number of spins, characterize the depth of Bell
correlations as a function of the interaction range, and determine the
critical time at which these correlations emerge. Consequently, we show that
scalable quantum resources can be generated with effective short-range
interactions for arbitrary large system size $N\gg 1$. Finally, we show that
the dynamically generated quantum states possessing many-body Bell
correlations can be certified with the help of present experimental
techniques.

Spin chains and Bell correlations— Let us consider a one-dimensional chain of
$N$ spins-$1/2$ described by the following Hamiltonian

$\displaystyle\hat{H}=\sum_{k,l=1}^{N}J_{kl}\hat{\sigma}_{z}^{(k)}\hat{\sigma}_{z}^{(l)},$
(1)

where $J_{kl}$ is the coupling strength (in units of $\hbar$) of the $k$-th
spin interacting with its $l$–th partner, $\hat{\sigma}_{z}^{(k)}$ are Pauli
$z$-operators. Here, we take the simplest finite-range interaction potential
governed by the rectangular function $J_{kl}=1$ for $0<|k-l|\leqslant r$, and
$0$ otherwise. The interaction range $r$ changes from $r=1$ (nearest-neighbor
couplings) to $r=N-1$ (all-to-all couplings), see Fig. 1. A recent experiment
[PhysRevLett.128.113602] showed a high degree of control over the distance-
selective interactions $J_{kl}$.

Let us start with the uncorrelated spins polarized along $x$-axis, orthogonal
to the $z$-axis distinguished by the Hamiltonian $\hat{H}$, and evolve the
system with the time-evolution operator $\hat{U}(\tau)=e^{-i\tau\hat{H}}$,
i.e.,

$|\psi(\tau)\rangle=\hat{U}(\tau)|1\rangle^{\otimes N}_{x},$ (2)

and let the interactions correlate spins. Here, $|1\rangle_{x}$ denotes an
eigenstate of the $x$-axis Pauli operator with positive eigenvalue. In the
limiting case of the all-to-all interactions, $r=N-1$, system dynamics is
equivalent to OAT, where arbitrary entanglement-depth many-body entanglement
and Bell-correlated states can be generated [Plodzien2022]. There arise
natural questions about the critical range $r$, at which quantum features
emerge, and about the relation between the range and the pace at which these
correlations are generated.

A well-suited tool to address these questions is the following correlator

$\displaystyle\mathcal{E}_{N}(\tau)=\left|\langle\hat{\sigma}_{+}^{(1)}\ldots\hat{\sigma}_{+}^{(N)}\rangle\right|^{2}\equiv
2^{Q_{N}-N},$ (3)

where the rising operators are taken along the $x$-axis, i.e.,
$\hat{\sigma}_{+}^{(k)}=\frac{1}{2}(\hat{\sigma}_{y}^{(k)}+i\hat{\sigma}_{z}^{(k)})$.
This choice of orientation is dictated by the observation that the OAT
dynamics generates GHZ-type superpositions of $|1\rangle_{x}^{\otimes N}$ and
$|-1\rangle_{x}^{\otimes N}$ states at time $\tau=\pi/4$.

To see how $\mathcal{E}_{N}$ relates to the local hidden variable (LHV)
theory, replace each $\hat{\sigma}_{x/y}^{(k)}$ with a $c$-number quantity
$\sigma_{x/y}^{(k)}$, that can take binary ($\pm 1$) values, depending on the
random hidden variable $\lambda$. The locality resides in the assumption that
$\sigma$’s depend only on a single label $k$, while realism implies the
presence of a probability distribution $p(\lambda)$ for the hidden variable.
Hence, a correlator $\mathcal{E}_{N}$ consistent with the LHV theory can be
expressed as

$\displaystyle\mathcal{E}_{N}=\left|\int\\!\\!d\lambda\,p(\lambda)\sigma_{+}^{(1)}(\lambda)\ldots\sigma_{+}^{(N)}(\lambda)\right|^{2}.$
(4)

Employing now the Cauchy-Schwarz inequality and the fact that for binary
outcomes $\left|\sigma_{+}^{(k)}(\lambda)\right|^{2}=1/2$, we obtain the bound
$\mathcal{E}_{N}\leqslant 2^{-N}$, which is the $N$-body Bell inequality
[PhysRevLett.126.210506]. Note that in Eq. (4) we did not assume the $N$
subsystems are described by quantum mechanics. If that were the case, the
restriction that the quantum spin lies within the Bloch sphere would yield
$\left|\langle\hat{\sigma}_{+}^{(k)}(\lambda)\rangle\right|^{2}\leqslant 1/4$,
where the mean denotes the trace with the single-qubit density matrix. Taking
all subsystems as quantum-mechanical qubits, we obtain the inequality
$\mathcal{E}_{N}\leqslant 4^{-N}$. Its violation implies the presence of
entanglement in the system. Among these two limiting cases of the structure of
the bounds ($\mathcal{E}_{N}\leqslant 4^{-N}$ and $\mathcal{E}_{N}\leqslant
2^{-N}$), there is a large set of inequalities obtained under assumptions of
quantum-mechanical restrictions for some subset of $N$ parties. Violation of
such inequalities implies the EPR steering [uola2020steering]. Note that the
Bell inequality $\mathcal{E}_{N}>2^{-N}$, or equivalently $Q_{N}>0$ [cf. Eq.
(3)], which is the most difficult to violate, ensures entanglement and the EPR
steering in the system.

Furthermore, the larger the inequality violation, the stronger correlations
are quantified by the correlator, which are related to the concept of
nonlocality depth. To see this, take the maximally nonclassical state, the
Greenberg-Horne-Zeilinger (GHZ) state

$\displaystyle|\psi_{\rm
GHZ}\rangle=\frac{1}{\sqrt{2}}\left(|1\rangle^{\otimes
N}_{x}+|-1\rangle^{\otimes N}_{x}\right),$ (5)

which gives the largest value of the Bell correlator (3), i.e.,
$\mathcal{E}_{N}=\frac{1}{4}$, or more conveniently for this purpose,
$Q_{N}=N-2$ [see the r.h.s. of Eq. (3)]. If a single spin is not Bell-
correlated in the chain, i.e., it correlates with the other $N-1$ spins
according to Eq. (4), we obtain the maximal value of $Q_{N}=N-3$. Hence,
$Q_{N}>N-3$ implies that the nonlocality encompasses all the spins, and we say
that the nonlocality depth, denoted with $\nu$, is $\nu=N$. These
considerations can be generalized to other integer values of $Q_{N}$, and when

$\displaystyle\nu-3<Q_{N}\leqslant\nu-2,$ (6)

then up to $\nu$ spins are Bell-correlated. For further analysis, we introduce
the fraction of the Bell correlated particles defined as $\beta=Q_{N}/N$. For
more on nonlocality depth and the related depth of entanglement, see Refs.
[zukowski2002bell, cavalcanti2007bell, he2011entanglement, he2010bell,
cavalcanti2011unified, spiny.milosz, PhysRevLett.126.210506].

Universal critical interaction range—

Figure 2: Panels (a)-(d) The Bell correlator $Q_{N}(\tau)$ for
$N=64,128,256,300$ spins as a function of time $\tau$ for the interaction
range: $r=1$ (thin-solid-black), $r=2$ (dotted-blue) $r=3$ (dashed-dotted-
green), $r=4$ (thick-solid-red), $r=5$ (dashed-orange) and $r=N-1$ (solid-
dark-blue). The values $Q_{N}>0$ mark the region where the many-body Bell
inequality is violated. (e) The first maximum of $Q_{N}$ with respect to time,
$Q^{\mathrm{max}}_{N}$, as a function of the total number of spins $N$ for
various $r\in[1,5]$. (f) The exponent of the Bell correlator in the scaling
with $N$, approximated as $Q_{N}^{\mathrm{max}}\approx\gamma
N+\mathrm{const}$, or, equivalently,
$2^{N}\max_{\tau}\mathcal{E}_{N}(\tau)\propto 2^{\gamma N}$. The growth of
scalable many-body Bell correlations with $N$ is manifested by positive
$\gamma$ when $r\geqslant 4$.

We search for the critical interaction range $r$ that allows the dynamical
generation of scalable quantum resources for the arbitrary number of spin $N$,
witnessed by Eq. (3), during time-evolution Eq. (2). In principle, for a given
system of size $N$, the minimal value of the range $r$ for which the state
exhibits Bell correlations may depend on $N$. Remarkably, as we show below, we
find a universal value of this range independent of the length of the chain.

First, we note that in general, the exact numerical calculation of the
correlator in Eq. (3) is exponentially hard as it requires working in the full
many-body basis of dimension $2^{N}$ with access to $N$-body operator
$\hat{\sigma}_{+}^{(1)}\ldots\hat{\sigma}_{+}^{(N)}$. However, in the
considered scenario, where $[\hat{H},\hat{\sigma}^{(k)}_{z}]=0$, the Eq. (3),
as we have found, can be expressed as

$\displaystyle\mathcal{E}_{N}(\tau)=\big{|}2^{-N}\sum_{\vec{s},\vec{s}\,^{\prime}}e^{-i(H_{\vec{s}}-H_{\vec{s}^{\prime}})\tau}s_{1}\ldots
s_{N}\big{|}^{2},$ (7)

where $H_{\vec{s}}$ is the $c$-number counterpart of Eq. (1) with operators
$\hat{\sigma}_{z}^{(k)}$ replaced by $s_{k}$ [supp]. The sums over $\vec{s}$
and $\vec{s}\,^{\prime}$ run over $2^{N}$ combinations of $s_{1},\ldots,s_{N}$
and $s^{\prime}_{1},\ldots,s^{\prime}_{N}$ which take values $\pm 1$. The
formula in Eq. (7) allows for the exact calculation of the Bell correlator for
the high number of qubits (here up to $N=300$) since the exponentially growing
dimensionality of the Hilbert space is no longer a computational constraint
111The dynamically generated highly entangled states of $N$ spins-$1/2$ in
OAT-like protocols have volume law scaling of the entangled entropy $S$, i.e.
$S\propto N$. To capture the entangled entropy $S$ the state-of-the-art tensor
network methods for studies of one-dimensional spin chains, i.e. Matrix
Product States, the bond dimension scales as $S\propto\ln{D}$. As such, to
fully capture many-body entangled and Bell correlated states dynamics,
including the most entangled many-body GHZ state, the bond dimension should
scale exponentially with the number of spins $D\propto e^{N}$, which is
impractical for a large number of spins, here up to $N=80$. Our analytical
approach is free from MPS restrictions and from the course of dimensionality
in the considered problem. (for derivation details see [supp]).

Let us consider first the marginal cases. For the aforementioned $r=1$ case,
the Bell correlator is given by the analytical formula
$\mathcal{E}_{N}(\tau)=\sin^{N}(\tau)\cos^{3N-4}(\tau)$. This expression has a
maximum ${\rm max}_{\tau}(\mathcal{E}_{N})\approx 2^{-1.6N+1}$ at
$\tau\simeq\pi/6$, which is exponentially smaller than the Bell limit
$2^{-N}$. Equivalently, $Q_{N}=-0.6N+1$ and deviation from the Bell limit
$Q_{N}=0$ becomes larger with $N$. Therefore, the nearest-neighbors
interactions cannot generate many-body Bell correlations in our case. The
other extreme case is the all-to-all interactions with $r=N-1$, which was
considered in detail in [PhysRevLett.129.250402]. This case realizes the one-
axis twisting protocol and the Bell correlations are present starting from
$\tau\simeq 1.5/N$. The correlator in Eq. (3) reaches its maximal value
$\mathcal{E}_{N}=1/4$ at $\tau=\pi/4$, when the $N$-body GHZ state is formed
[PhysRevA.102.013328, PhysRevLett.129.090403, PhysRevA.107.013311,
https://doi.org/10.48550/arxiv.2302.09829], which is exponentially larger than
the bound in the $N$-body Bell inequality.

In Fig. 2 in panels (a)—(d) we present time evolution of the correlator
$Q_{N}(\tau)$, calculated using Eq. (7), for different intermediate ranges
$1\leqslant r\leqslant 5$ and for all-to-all couplings $r=N-1$ for various
spin numbers $N=64$ (panel a), 128 (b), 256 (c) and 300 (d). We observe that
the correlator only breaks the Bell limit $Q_{N}>0$ when $r\geqslant 4$ even
for large $N$. This is the first indication of the presence of a critical
range. However, it must be verified if the value of $r=4$ is universal, i.e.,
independent of the system size. To this end, we focus on

$\displaystyle Q_{N}^{\rm max}\equiv\max_{\tau}Q_{N}(\tau),$ (8)

i.e., maximized $Q_{N}$ with respect to $\tau$ calculated for different $N$’s.
In Fig. 2(e) we show $Q_{N}^{\rm max}$, as a function of the number of spins
$N$ and for various values of $1\leqslant r\leqslant 5$. For a given $r$,
these maximal values lie on a straight line determining an exponent $\gamma$
in the exponential scaling of $\mathcal{E}_{N}\propto 2^{(\gamma-1)N}$. In
Fig. 2(f), we present the value of $\gamma$ as a function of the interaction
range. For $r\geqslant 4$ we find positive exponents indicating that the
degree of the Bell inequality violation becomes higher with growing $N$ in
contrast to the case of $r\leqslant 3$. Therefore, the change in the sign of
$\gamma$ indicates qualitatively different scaling regimes of the Bell
correlations.

The observed behaviour of $Q_{N}^{\rm max}$ and $\gamma(r)>0$ for $r\geqslant
4$ demonstrates that $r=4$ is indeed a critical range, at which the Bell
correlations are detected. Crucially, $Q_{N}^{\rm max}$ increases with growing
$N$, indicating scalable Bell correlations in this system
[PhysRevLett.126.210506]. Remarkably, although $r=4$ in the limit $N\gg 1$ is
a short-range interaction encompassing an intensive number of neighboring
spins, the system still exhibits an increasing degree of violation of the Bell
inequality with increasing the system size. The invariance of the critical
range with respect to the number of spins proves its universality.

Qualitatively, the strengthening of Bell correlations with growing $r$ can be
explained as follows. To generate the GHZ-like coherence, which is the witness
of Bell correlations [see Eq. (3)] one needs to flip all the spins, hence, to
act with an $N$-body operator
$\hat{\mathcal{B}}_{N}=\hat{\sigma}_{z}^{(1)}\ldots\hat{\sigma}_{z}^{(N)}$ on
the input state, since $\hat{\sigma}_{z}|1\rangle_{x}=|-1\rangle_{x}$. Note
that the evolution operator coupling any pair of spins is

$\displaystyle
e^{-i\tau\hat{\sigma}_{z}^{(k)}\hat{\sigma}_{z}^{(l)}}=\cos(\tau)-i\sin(\tau)\hat{\sigma}_{z}^{(k)}\hat{\sigma}_{z}^{(l)}.$
(9)

As $r$ grows, each spin couples to more neighbors, and the number of possible
combinations of interacting terms giving the operator $\mathcal{B}_{N}$ grows
exponentially, increasing, in consequence, $\mathcal{E}_{N}(\tau)$. This
growth is, however, balanced by the exponential decay of the amplitude of the
initial state. The eventual observation of the Bell inequality violation is,
thus, the effect of the competition between the two indicated mechanisms. We
formalize these intuitions with the help of the developed spin-inversion
asymptotic expansion theory; for quantitative explanations using diagrammatic
approach 222For intermediate values of $r$, in our open-access repository
[code_repo], we provide examples of the time evolution of
$\mathcal{E}_{N}(\tau)$ obtained analytically for concrete values of $N$ and
$r$. and the details on the asymptotic theory, see [supp].

We note here that identifying relevant mechanisms helps clarify the role of
the boundary conditions. The periodic ones provide some additional coupling
between the spins. Indeed, numerical results confirm that in such a case,
$\mathcal{E}_{N}(\tau)$ is, in general, larger as compared to the open
boundary conditions considered here. Our decision to focus on the open ones is
motivated by the geometry of most of the experimental setups.

Figure 3: (a) The critical time $\tau_{c}$ at which the Bell correlator passes
the Bell limit, $Q_{N}=0$, as a function of interaction range $r$. The dashed
curves correspond to the limit $r\lesssim N$, while solid curves correspond to
the limit $r\ll N$. The latter collapses onto the same line. (b) The maximal
fraction of correlated spins $\beta=\nu/N$ as a function of range $r$, cf. Eq.
(6).

Critical time— As evident from Fig. 2(a)—(d), the interaction range $r$
determines the time, denoted with $\tau_{c}$, at which Bell correlations
emerge. This critical time, at which the correlator $Q_{N}(\tau)$ crosses the
Bell limit, depends on $r$, and the generation of correlations accelerates for
larger $r$’s. This dependence can be extracted from Eq. (7) as follows. For an
arbitrary internal spin coupled with $2r$ neighbours, the sum over $s_{k}=\pm
1$ vanishes unless the corresponding phase-term $e^{-i\tau
s_{k}(s_{k-r}+\ldots+s_{k+r})}$ oscillates quickly enough. The sum in the
parenthesis is at most equal to $2r$, hence the phase factor will vary
significantly between $s_{k}=-1$ and $s_{k}=1$ if $\tau\cdot 2r\gtrsim 1$.
Thus, the Bell correlator becomes significantly non-zero if $\tau_{c}\gtrsim
a/r$, where $a$ is some constant. The scaling of $\tau_{c}$, which is
inversely proportional to the interaction range, is confirmed by exact
solution of the dynamics generated by Eq. (1). In Fig. 3(a), we identify the
critical time $\tau_{c}$ when $Q_{N}$ surpasses the Bell limit for different
values of large $N$=80, 100, 120, 140, for which $r\ll N$ and small $N$=8, 10,
12, 14, where $r\lesssim N$. The observed behavior shows that for large
systems, $\tau_{c}$ becomes independent of the system size. By changing the
range $r$ from the critical value $r=4$ to $r=N=14$ we recover the scaling
$\tau_{c}\propto 1/N$ characteristic to the all-to-all coupling case.

Fraction of Bell correlations spins— In Fig. 3(b), we present the dependence
of the fraction of Bell correlated spins $\beta=Q_{N}/N$ as a function of the
range $r$, deduced from $Q_{N}^{\rm max}$, see Eq. (6). In the limit of $N\gg
r$, the results follow the universal curve. This can be understood as follows.
Using the asymptotic exponential scaling for large $N$, i.e.,
$Q_{N}\approx\gamma(r)N$ [see the discussion below Eq. (8)], we can estimate
the nonlocality depth by $\nu_{N}\approx Q_{N}$ up to corrections on the order
of $1/N$. Then, the fraction $\beta\approx\gamma(r)$, which shows that for
large $N$, $\beta$ becomes independent of $N$ as is observed in the figure.
This relation uncovers the physical role of the exponent of $\mathcal{E}_{N}$
as a concentration of Bell correlated particles. Therefore, we find the range
not only controls the critical time $\tau_{c}$, but also the number of Bell-
correlated clusters in the system.

Figure 4:  The dynamics of the estimated many-body Bell correlator $Q_{N}$ for
(a) $N=4$, (b) $N=6$, (c) $N=8$, and (d) $N=10$ spins with interaction range
$r=4$ (lower curves), and all-to-all connections $r=N-1$ (upper curves). Solid
blue lines present exact results, while red circles represent the value of the
reconstructed many-body Bell correlator from the classical shadows tomography.
Standard deviation is marked as a shaded area.

Many-body Bell correlation certification— Measuring high-order correlation
functions is challenging. Recent advances in the control of many-body systems
allowed for measurements of correlation functions up to the $6$-th order
[Dall2013], $2$-nd Rényi entropy for $N=4$ [Islam2015], and $N=5$ particles
[Linke2018] via extraction characteristic of a quantum state using a
controlled-swap gate acting on two copies of the state [Eckert2002], while for
$N=10$ particles [Brydges2019] using the randomized measurements technique
[Vermersch2018, Elben2019, Elben2020, Elben2020_PRL, Rath2021, Elben2022]. The
many-body correlations considered here, however, can be certified by
inspecting only one element of the density matrix, which couples the state
with all spins up with the state with all spins down in the $x$-basis, cf. Eq.
(3). As such, the problem of many-body Bell correlations measurement can be
cast as a quantum state tomography task [PhysRevLett.74.4101,
PhysRevLett.83.3103, PhysRevLett.92.220402, Haffner2005, 5714248, Toth2010,
Moroder_2012, Cramer2010, PhysRevLett.111.020401, Lanyon2017], which recently
has been enhanced by deep neural networks techniques [adri1, Pan2022,
koutny2022PRA, Shahnawaz2021, ma2023attentionbased, palmieri2023enhancing,
dawid2022modern].

To this end, we simulate tomographical reconstruction of the density matrix
with classical shadows tomography (CST) [10.1145/3313276.3316378, Huang2020,
altepeter20044, 7956181, 10.1145/2897518.2897544, Koh2022classicalshadows,
aaronson2018shadow, PhysRevLett.125.200501, Huang2020, PhysRevLett.127.030503,
Elben2023], which was employed for trapped ions [mcginley2022shadow,
PhysRevX.13.011049, doi:10.1126/science.1208001], and Rydberg atom arrays
[Notarnicola_2023]. We prepared $10$ reconstructions of the target density
matrix $\hat{\varrho}(\tau)$ at a given time $\tau$, where each reconstruction
consists of $M$ classical shadows. From each reconstructed
$\hat{\varrho}(\tau)$, we extract the Bell correlator, and based on the
generated collection, we estimated the mean and standard deviation of the
reconstruction ${\cal E}^{*}(\tau)\equiv 2^{Q_{N}^{*}(\tau)-N}$ [supp]. In
Fig. 4, we demonstrate that Bell correlations quantified by $Q_{N}(\tau)$
(solid blue lines) can be successfully certified by the tomographically
reconstructed mean value $Q_{N}^{*}(\tau)$ (red circles) with standard
deviation (shaded areas) for $N=4,6,8,10$ spins with the interaction range
$r=4$, and the all-to-all couplings when $r=N-1$. We prepared $M=10^{5}N$
classical shadow for $N=10$ and $r=4$, and $M=10^{4}N$ otherwise.

An alternative method to full quantum state tomography relies on the concept
of multiple quantum coherences (MQC), which provide extensive information
about the structure of a many-body state and allows for the relation of the
quantum correlators with other physical quantities like the out-of-time-order
correlations [garttner2017measuring, PhysRevLett.120.040402, supp]. In SM we
provide details how $\mathcal{E}_{N}$ can be determined by means of high-order
multiple-quantum intensity in MQC which is feasible experimentally.

Discussion and Conclusions— In this work, we showed that scalable many-body
entanglement and Bell correlations can be generated with the Ising-type
Hamiltonian with finite-range two-body interactions, in contrast to the
broadly considered all-to-all couplings present in the standard one-axis
twisting protocols. We found a critical interaction range that is independent
of the system size, which means that for large system size $N\gg 1$
interactions are short-ranged. Yet, the degree of violating the many-body Bell
inequality grows with the number of spins. Moreover, the studied quantum
resources can be experimentally certified with the classical shadows quantum
state tomography and multiple-quantum coherences measurements.

Our study addresses the constant need for protocols to generate and
certificate $N$-body entanglement in quantum circuits. First, our results have
important consequences for future digital quantum computers operating on large
numbers of qubits, aiming to generate scalable quantum resources dynamically.
In a single trotterization step of the time-evolution in the one-axis twisting
protocol, there is a need for the all-to-all connections between qubits; thus,
the number of two-qubit entangling gates scales quadratically with the number
of qubits [PhysRevLett.119.180511, doi:10.1126/science.aay0600,
doi:10.1126/sciadv.aba4935, Li2018], while for the finite-critical-range
interactions considered here, it scales linearly which is a significant
reduction of the problem complexity. Next, our analytical results allow us to
extract many-body quantum correlations, in principle, for an arbitrary number
of qubits, which can be used to certify and benchmark quantum computers
operating on hundreds of qubits [wurtz2023aquila, Kim2023].

Author contributions— M.P. indicated the existence of the critical interaction
range, performed many-body simulations, analytical calculations in Ref.
[code_repo], and quantum state tomography results. J.Ch. prepared analytical
expressions for the Bell correlator, critical range, and critical time, and
performed many-body calculations. T.W. developed the spin-inversion asymptotic
expansion theory and performed many-body calculations. All the authors
contributed to discussing the results and the manuscript preparation and
revision.

Ackonwledgements— We acknowledge the discussion with Arghavan Safavi-Naini,
Rene Gerritsma, Rainer Blatt and Tomasso Roscilde. J.Ch. was funded by the
National Science Centre, Poland, within the QuantERA II Programme that has
received funding from the European Union’s Horizon 2020 research and
innovation programme under Grant Agreement No 101017733, Project No.
2021/03/Y/ST2/00195. E.W. acknowledges support from the DAINA project of the
Polish National Science Center DEC-2020/38/L/ST2/00375. This research is part
of the project No. 2021/43/P/ST2/02911 co-funded by the National Science
Centre and the European Union’s Horizon 2020 research and innovation programme
under the Marie Skłodowska-Curie grant agreement no. 945339. For the purpose
of Open Access, the author has applied a CC-BY public copyright licence to any
Author Accepted Manuscript (AAM) version arising from this submission. ICFO
group acknowledges support from: ERC AdG NOQIA; Ministerio de Ciencia y
Innovation Agencia Estatal de Investigaciones
(PGC2018-097027-B-I00/10.13039/501100011033,
CEX2019-000910-S/10.13039/501100011033, Plan National FIDEUA
PID2019-106901GB-I00, FPI, QUANTERA MAQS PCI2019-111828-2, QUANTERA DYNAMITE
PCI2022-132919, Proyectos de I+D+I “Retos Colaboración” QUSPIN
RTC2019-007196-7); MICIIN with funding from European Union
NextGenerationEU(PRTR-C17.I1) and by Generalitat de Catalunya; Fundació
Cellex; Fundació Mir-Puig; Generalitat de Catalunya (European Social Fund
FEDER and CERCA program, AGAUR Grant No. 2021 SGR 01452, QuantumCAT
U16-011424, co-funded by ERDF Operational Program of Catalonia 2014-2020);
Barcelona Supercomputing Center MareNostrum (FI-2022-1-0042); EU (PASQuanS2.1,
101113690); EU Horizon 2020 FET-OPEN OPTOlogic (Grant No 899794); EU Horizon
Europe Program (Grant Agreement 101080086 — NeQST), National Science Centre,
Poland (Symfonia Grant No. 2016/20/W/ST4/00314); ICFO Internal “QuantumGaudi”
project; European Union’s Horizon 2020 research and innovation program under
the Marie-Skłodowska-Curie grant agreement No 101029393 (STREDCH). M.P.
acknowledges the support of the Polish National Agency for Academic Exchange,
the Bekker programme no: PPN/BEK/2020/1/00317. Views and opinions expressed
are, however, those of the author(s) only and do not necessarily reflect those
of the European Union, European Commission, European Climate, Infrastructure
and Environment Executive Agency (CINEA), nor any other granting authority.
Neither the European Union nor any granting authority can be held responsible
for them.

## References

  * Horodecki  _et al._ [2009] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Rev. Mod. Phys. 81, 865 (2009).
  * Uola  _et al._ [2020] R. Uola, A. C. S. Costa, H. C. Nguyen, and O. Gühne, Quantum steering, Rev. Mod. Phys. 92, 015001 (2020).
  * Brunner  _et al._ [2014] N. Brunner, D. Cavalcanti, S. Pironio, V. Scarani, and S. Wehner, Bell nonlocality, Rev. Mod. Phys. 86, 419 (2014).
  * Acín  _et al._ [2018] A. Acín, I. Bloch, H. Buhrman, T. Calarco, C. Eichler, J. Eisert, D. Esteve, N. Gisin, S. J. Glaser, F. Jelezko, S. Kuhr, M. Lewenstein, M. F. Riedel, P. O. Schmidt, R. Thew, A. Wallraff, I. Walmsley, and F. K. Wilhelm, The quantum technologies roadmap: a european community view, New Journal of Physics 20, 080201 (2018).
  * Eisert  _et al._ [2020] J. Eisert, D. Hangleiter, N. Walk, I. Roth, D. Markham, R. Parekh, U. Chabaud, and E. Kashefi, Quantum certification and benchmarking, Nature Reviews Physics 2, 382 (2020).
  * Fraxanet  _et al._ [2023] J. Fraxanet, T. Salamon, and M. Lewenstein, The coming decades of quantum simulation, in  _Sketches of Physics_ (Springer International Publishing, 2023) pp. 85–125.
  * Frérot  _et al._ [2023] I. Frérot, M. Fadel, and M. Lewenstein, Probing quantum correlations in many-body systems: a review of scalable methods, Reports on Progress in Physics 86, 114001 (2023).
  * Kitagawa and Ueda [1993a] M. Kitagawa and M. Ueda, Squeezed spin states, Phys. Rev. A 47, 5138 (1993a).
  * Wineland  _et al._ [1994a] D. Wineland, J. Bollinger, W. Itano, and D. Heinzen, Squeezed atomic states and projection noise in spectroscopy, Phys. Rev. A 50, 67 (1994a).
  * Pezzè  _et al._ [2018] L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, Quantum metrology with nonclassical states of atomic ensembles, Rev. Mod. Phys. 90, 035005 (2018).
  * Wolfgramm  _et al._ [2010] F. Wolfgramm, A. Cerè, F. A. Beduini, A. Predojević, M. Koschorreck, and M. W. Mitchell, Squeezed-light optical magnetometry, Phys. Rev. Lett. 105, 053601 (2010).
  * Wineland  _et al._ [1994b] D. J. Wineland, J. J. Bollinger, W. M. Itano, and D. J. Heinzen, Squeezed atomic states and projection noise in spectroscopy, Phys. Rev. A 50, 67 (1994b).
  * Müller-Rigat  _et al._ [2023] G. Müller-Rigat, A. K. Srivastava, S. Kurdziałek, G. Rajchel-Mieldzioć, M. Lewenstein, and I. Frérot, Certifying the quantum fisher information from a given set of mean values: a semidefinite programming approach (2023).
  * Tura  _et al._ [2014] J. Tura, R. Augusiak, A. B. Sainz, T. Vértesi, M. Lewenstein, and A. Acín, Detecting nonlocality in many-body quantum states, Science 344, 1256 (2014).
  * Schmied  _et al._ [2016] R. Schmied, J.-D. Bancal, B. Allard, M. Fadel, V. Scarani, P. Treutlein, and N. Sangouard, Bell correlations in a bose-einstein condensate, Science 352, 441 (2016).
  * Aloy  _et al._ [2019] A. Aloy, J. Tura, F. Baccari, A. Acín, M. Lewenstein, and R. Augusiak, Device-independent witnesses of entanglement depth from two-body correlators, Phys. Rev. Lett. 123, 100507 (2019).
  * Baccari  _et al._ [2019] F. Baccari, J. Tura, M. Fadel, A. Aloy, J.-D. Bancal, N. Sangouard, M. Lewenstein, A. Acín, and R. Augusiak, Bell correlation depth in many-body systems, Phys. Rev. A 100, 022121 (2019).
  * Tura  _et al._ [2019] J. Tura, A. Aloy, F. Baccari, A. Acín, M. Lewenstein, and R. Augusiak, Optimization of device-independent witnesses of entanglement depth from two-body correlators, Phys. Rev. A 100, 032307 (2019).
  * Müller-Rigat  _et al._ [2021] G. Müller-Rigat, A. Aloy, M. Lewenstein, and I. Frérot, Inferring nonlinear many-body bell inequalities from average two-body correlations: Systematic approach for arbitrary spin-$j$ ensembles, PRX Quantum 2, 030329 (2021).
  * Żukowski and Brukner [2002] M. Żukowski and Č. Brukner, Bell’s theorem for general n-qubit states, Phys. Rev. Lett. 88, 210401 (2002).
  * Cavalcanti  _et al._ [2007] E. G. Cavalcanti, C. J. Foster, M. D. Reid, and P. D. Drummond, Bell inequalities for continuous-variable correlations, Phys. Rev. Lett. 99, 210405 (2007).
  * He  _et al._ [2011] Q. He, P. Drummond, and M. Reid, Entanglement, epr steering, and bell-nonlocality criteria for multipartite higher-spin systems, Phys. Rev. A 83, 032120 (2011).
  * Cavalcanti  _et al._ [2011] E. Cavalcanti, Q. He, M. Reid, and H. Wiseman, Unified criteria for multipartite quantum nonlocality, Phys. Rev. A 84, 032115 (2011).
  * Niezgoda  _et al._ [2020] A. Niezgoda, M. Panfil, and J. Chwedeńczuk, Quantum correlations in spin chains, Phys. Rev. A 102, 042206 (2020).
  * Niezgoda and Chwedeńczuk [2021] A. Niezgoda and J. Chwedeńczuk, Many-body nonlocality as a resource for quantum-enhanced metrology, Phys. Rev. Lett. 126, 210506 (2021).
  * Chwedenczuk [2022] J. Chwedenczuk, Many-body Bell inequalities for bosonic qubits, SciPost Phys. Core 5, 025 (2022).
  * Wineland  _et al._ [1992] D. J. Wineland, J. J. Bollinger, W. M. Itano, F. L. Moore, and D. J. Heinzen, Spin squeezing and reduced quantum noise in spectroscopy, Phys. Rev. A 46, R6797 (1992).
  * Kitagawa and Ueda [1993b] M. Kitagawa and M. Ueda, Squeezed spin states, Phys. Rev. A 47, 5138 (1993b).
  * Gietka  _et al._ [2015] K. Gietka, P. Szańkowski, T. Wasak, and J. Chwedeńczuk, Quantum-enhanced interferometry and the structure of twisted states, Phys. Rev. A 92, 043622 (2015).
  * Li  _et al._ [2008] Y. Li, Y. Castin, and A. Sinatra, Optimum spin squeezing in bose-einstein condensates with particle losses, Phys. Rev. Lett. 100, 210401 (2008).
  * Li  _et al._ [2009] Y. Li, P. Treutlein, J. Reichel, and A. Sinatra, Spin squeezing in a bimodal condensate: spatial dynamics and particle losses, The European Physical Journal B 68, 365 (2009).
  * Wang  _et al._ [2017] M. Wang, W. Qu, P. Li, H. Bao, V. Vuletić, and Y. Xiao, Two-axis-twisting spin squeezing by multipass quantum erasure, Phys. Rev. A 96, 013823 (2017).
  * Kajtoch  _et al._ [2018a] D. Kajtoch, E. Witkowska, and A. Sinatra, Spin-squeezed atomic crystal, EPL (Europhysics Letters) 123, 20012 (2018a).
  * Schulte  _et al._ [2020] M. Schulte, C. Lisdat, P. O. Schmidt, U. Sterr, and K. Hammerer, Prospects and challenges for squeezing-enhanced optical atomic clocks, Nature Communications 11, 5955 (2020).
  * Gietka  _et al._ [2021] K. Gietka, A. Usui, J. Deng, and T. Busch, Simulating the same physics with two distinct hamiltonians, Phys. Rev. Lett. 126, 160402 (2021).
  * Comparin  _et al._ [2022] T. Comparin, F. Mezzacapo, and T. Roscilde, Robust spin squeezing from the tower of states of u(1)-symmetric spin hamiltonians, Phys. Rev. A 105, 022625 (2022).
  * Riedel  _et al._ [2010] M. F. Riedel, P. Böhi, Y. Li, T. W. Hänsch, A. Sinatra, and P. Treutlein, Atom-chip-based generation of entanglement for quantum metrology, Nature 464, 1170–1173 (2010).
  * Gross  _et al._ [2010] C. Gross, T. Zibold, E. Nicklas, J. Estève, and M. K. Oberthaler, Nonlinear atom interferometer surpasses classical precision limit, Nature 464, 1165–1169 (2010).
  * Hamley  _et al._ [2012] C. D. Hamley, C. S. Gerving, T. M. Hoang, E. M. Bookjans, and M. S. Chapman, Spin-nematic squeezed vacuum in a quantum gas, Nature Physics 8, 305–308 (2012).
  * Qu  _et al._ [2020] A. Qu, B. Evrard, J. Dalibard, and F. Gerbier, Probing spin correlations in a bose-einstein condensate near the single-atom level, Phys. Rev. Lett. 125, 033401 (2020).
  * Leroux  _et al._ [2010] I. D. Leroux, M. H. Schleier-Smith, and V. Vuletić, Implementation of cavity squeezing of a collective atomic spin, Phys. Rev. Lett. 104, 073602 (2010).
  * Maussang  _et al._ [2010] K. Maussang, G. E. Marti, T. Schneider, P. Treutlein, Y. Li, A. Sinatra, R. Long, J. Estève, and J. Reichel, Enhanced and reduced atom number fluctuations in a bec splitter, Phys. Rev. Lett. 105, 080403 (2010).
  * Omran  _et al._ [2019] A. Omran, H. Levine, A. Keesling, G. Semeghini, T. T. Wang, S. Ebadi, H. Bernien, A. S. Zibrov, H. Pichler, S. Choi, J. Cui, M. Rossignolo, P. Rembold, S. Montangero, T. Calarco, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Generation and manipulation of schrödinger cat states in rydberg atom arrays, Science 365, 570 (2019), https://www.science.org/doi/pdf/10.1126/science.aax9743 .
  * Bornet  _et al._ [2023] G. Bornet, G. Emperauger, C. Chen, B. Ye, M. Block, M. Bintz, J. A. Boyd, D. Barredo, T. Comparin, F. Mezzacapo, T. Roscilde, T. Lahaye, N. Y. Yao, and A. Browaeys, Scalable spin squeezing in a dipolar rydberg atom array, Nature 621, 728 (2023).
  * Eckner  _et al._ [2023] W. J. Eckner, N. Darkwah Oppong, A. Cao, A. W. Young, W. R. Milner, J. M. Robinson, J. Ye, and A. M. Kaufman, Realizing spin squeezing with rydberg interactions in an optical clock, Nature 621, 734 (2023).
  * 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, Quantum spin dynamics and entanglement generation with hundreds of trapped ions, Science 352, 1297 (2016).
  * Franke  _et al._ [2023] J. Franke, S. R. Muleady, R. Kaubruegger, F. Kranzl, R. Blatt, A. M. Rey, M. K. Joshi, and C. F. Roos, Quantum-enhanced sensing on optical transitions through finite-range interactions, Nature 621, 740 (2023).
  * Song  _et al._ [2017] C. Song, K. Xu, W. Liu, C.-p. Yang, S.-B. Zheng, H. Deng, Q. Xie, K. Huang, Q. Guo, L. Zhang, P. Zhang, D. Xu, D. Zheng, X. Zhu, H. Wang, Y.-A. Chen, C.-Y. Lu, S. Han, and J.-W. Pan, 10-qubit entanglement and parallel logic operations with a superconducting circuit, Phys. Rev. Lett. 119, 180511 (2017).
  * Song  _et al._ [2019] C. Song, K. Xu, H. Li, Y.-R. Zhang, X. Zhang, W. Liu, Q. Guo, Z. Wang, W. Ren, J. Hao, H. Feng, H. Fan, D. Zheng, D.-W. Wang, H. Wang, and S.-Y. Zhu, Generation of multicomponent atomic schrödinger cat states of up to 20 qubits, Science 365, 574 (2019), https://www.science.org/doi/pdf/10.1126/science.aay0600 .
  * Xu  _et al._ [2020] K. Xu, Z.-H. Sun, W. Liu, Y.-R. Zhang, H. Li, H. Dong, W. Ren, P. Zhang, F. Nori, D. Zheng, H. Fan, and H. Wang, Probing dynamical phase transitions with a superconducting quantum simulator, Science Advances 6, eaba4935 (2020), https://www.science.org/doi/pdf/10.1126/sciadv.aba4935 .
  * Kajtoch  _et al._ [2018b] D. Kajtoch, E. Witkowska, and A. Sinatra, Spin-squeezed atomic crystal, EPL (Europhysics Letters) 123, 20012 (2018b).
  * He  _et al._ [2019] P. He, M. A. Perlin, S. R. Muleady, R. J. Lewis-Swan, R. B. Hutson, J. Ye, and A. M. Rey, Engineering spin squeezing in a 3d optical lattice with interacting spin-orbit-coupled fermions, Phys. Rev. Research 1, 033075 (2019).
  * Płodzień  _et al._ [2020a] M. Płodzień, M. Kościelski, E. Witkowska, and A. Sinatra, Producing and storing spin-squeezed states and greenberger-horne-zeilinger states in a one-dimensional optical lattice, Phys. Rev. A 102, 013328 (2020a).
  * Płodzień  _et al._ [2022a] M. Płodzień, M. Lewenstein, E. Witkowska, and J. Chwedeńczuk, One-axis twisting as a method of generating many-body bell correlations, Phys. Rev. Lett. 129, 250402 (2022a).
  * Płodzień  _et al._ [2023] M. Płodzień, T. Wasak, E. Witkowska, M. Lewenstein, and J. Chwedeńczuk, Generation of scalable many-body bell correlations in spin chains with short-range two-body interactions (2023), arXiv:2306.06173 [quant-ph] .
  * Mamaev  _et al._ [2021] M. Mamaev, I. Kimchi, R. M. Nandkishore, and A. M. Rey, Tunable-spin-model generation with spin-orbit-coupled fermions in optical lattices, Phys. Rev. Research 3, 013178 (2021).
  * Hernández Yanes  _et al._ [2022a] T. Hernández Yanes, M. Płodzień, M. Mackoit Sinkevičienė, G. Žlabys, G. Juzeliūnas, and E. Witkowska, One- and two-axis squeezing via laser coupling in an atomic fermi-hubbard model, Phys. Rev. Lett. 129, 090403 (2022a).
  * Dziurawiec  _et al._ [2023a] M. Dziurawiec, T. H. Yanes, M. Płodzień, M. Gajda, M. Lewenstein, and E. Witkowska, Accelerating many-body entanglement generation by dipolar interactions in the bose-hubbard model, Physical Review A 107, 10.1103/physreva.107.013311 (2023a).
  * Hernández Yanes  _et al._ [2023] T. Hernández Yanes, G. Žlabys, M. Płodzień, D. Burba, M. M. Sinkevičienė, E. Witkowska, and G. Juzeliūnas, Spin squeezing in open heisenberg spin chains, Phys. Rev. B 108, 104301 (2023).
  * Hollerith  _et al._ [2022] S. Hollerith, K. Srakaew, D. Wei, A. Rubio-Abadal, D. Adler, P. Weckesser, A. Kruckenhauser, V. Walther, R. van Bijnen, J. Rui, C. Gross, I. Bloch, and J. Zeiher, Realizing distance-selective interactions in a rydberg-dressed atom array, Phys. Rev. Lett. 128, 113602 (2022).
  * He  _et al._ [2010] Q. He, E. Cavalcanti, M. Reid, and P. Drummond, Bell inequalities for continuous-variable measurements, Phys. Rev. A 81, 062106 (2010).
  * [62] See Supplementary materials for detailed analytical calculations .
  * Note [1] The dynamically generated highly entangled states of $N$ spins-$1/2$ in OAT-like protocols have volume law scaling of the entangled entropy $S$, i.e. $S\propto N$. To capture the entangled entropy $S$ the state-of-the-art tensor network methods for studies of one-dimensional spin chains, i.e. Matrix Product States, the bond dimension scales as $S\propto\mathop{ln}\nolimits{D}$. As such, to fully capture many-body entangled and Bell correlated states dynamics, including the most entangled many-body GHZ state, the bond dimension should scale exponentially with the number of spins $D\propto e^{N}$, which is impractical for a large number of spins, here up to $N=80$. Our analytical approach is free from MPS restrictions and from the course of dimensionality in the considered problem.
  * Płodzień  _et al._ [2022b] M. Płodzień, M. Lewenstein, E. Witkowska, and J. Chwedeńczuk, One-Axis Twisting as a Method of Generating Many-Body Bell Correlations, Phys. Rev. Lett. 129, 250402 (2022b).
  * Płodzień  _et al._ [2020b] M. Płodzień, M. Kościelski, E. Witkowska, and A. Sinatra, Producing and storing spin-squeezed states and greenberger-horne-zeilinger states in a one-dimensional optical lattice, Phys. Rev. A 102, 013328 (2020b).
  * Hernández Yanes  _et al._ [2022b] T. Hernández Yanes, M. Płodzień, M. Mackoit Sinkevičienė, G. Žlabys, G. Juzeliūnas, and E. Witkowska, One- and two-axis squeezing via laser coupling in an atomic fermi-hubbard model, Phys. Rev. Lett. 129, 090403 (2022b).
  * Dziurawiec  _et al._ [2023b] M. Dziurawiec, T. Hernández Yanes, M. Płodzień, M. Gajda, M. Lewenstein, and E. Witkowska, Accelerating many-body entanglement generation by dipolar interactions in the bose-hubbard model, Phys. Rev. A 107, 013311 (2023b).
  * Yanes  _et al._ [2023] T. H. Yanes, G. Žlabys, M. Płodzień, D. Burba, M. M. Sinkevičienė, E. Witkowska, and G. Juzeliūnas, Spin squeezing in open heisenberg spin chains (2023).
  * Note [2] For intermediate values of $r$, in our open-access repository [code_repo], we provide examples of the time evolution of $\mathcal{E}_{N}(\tau)$ obtained analytically for concrete values of $N$ and $r$.
  * Dall  _et al._ [2013] R. G. Dall, A. G. Manning, S. S. Hodgman, W. RuGway, K. V. Kheruntsyan, and A. G. Truscott, Ideal n-body correlations with massive particles, Nature Physics 9, 341 (2013).
  * 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).
  * Linke  _et al._ [2018] N. M. Linke, S. Johri, C. Figgatt, K. A. Landsman, A. Y. Matsuura, and C. Monroe, Measuring the rényi entropy of a two-site fermi-hubbard model on a trapped ion quantum computer, Phys. Rev. A 98, 052334 (2018).
  * Ekert  _et al._ [2002] A. K. Ekert, C. M. Alves, D. K. L. Oi, M. Horodecki, P. Horodecki, and L. C. Kwek, Direct estimations of linear and nonlinear functionals of a quantum state, Phys. Rev. Lett. 88, 217901 (2002).
  * 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).
  * Vermersch  _et al._ [2018] B. Vermersch, A. Elben, M. Dalmonte, J. I. Cirac, and P. Zoller, Unitary $n$-designs via random quenches in atomic hubbard and spin models: Application to the measurement of rényi entropies, Phys. Rev. A 97, 023604 (2018).
  * Elben  _et al._ [2019] A. Elben, B. Vermersch, C. F. Roos, and P. Zoller, Statistical correlations between locally randomized measurements: A toolbox for probing entanglement in many-body quantum states, Phys. Rev. A 99, 052323 (2019).
  * Elben  _et al._ [2020a] A. Elben, J. Yu, G. Zhu, M. Hafezi, F. Pollmann, P. Zoller, and B. Vermersch, Many-body topological invariants from randomized measurements in synthetic quantum matter, Science Advances 6, eaaz3666 (2020a).
  * Elben  _et al._ [2020b] A. Elben, R. Kueng, H.-Y. R. Huang, R. van Bijnen, C. Kokail, M. Dalmonte, P. Calabrese, B. Kraus, J. Preskill, P. Zoller, and B. Vermersch, Mixed-state entanglement from local randomized measurements, Phys. Rev. Lett. 125, 200501 (2020b).
  * Rath  _et al._ [2021] A. Rath, C. Branciard, A. Minguzzi, and B. Vermersch, Quantum fisher information from randomized measurements, Phys. Rev. Lett. 127, 260501 (2021).
  * Elben  _et al._ [2022] A. Elben, S. T. Flammia, H.-Y. Huang, R. Kueng, J. Preskill, B. Vermersch, and P. Zoller, The randomized measurement toolbox (2022).
  * Leonhardt [1995] U. Leonhardt, Quantum-state tomography and discrete wigner function, Phys. Rev. Lett. 74, 4101 (1995).
  * White  _et al._ [1999] A. G. White, D. F. V. James, P. H. Eberhard, and P. G. Kwiat, Nonmaximally entangled states: Production, characterization, and utilization, Phys. Rev. Lett. 83, 3103 (1999).
  * Roos  _et al._ [2004] C. F. Roos, G. P. T. Lancaster, M. Riebe, H. Häffner, W. Hänsel, S. Gulde, C. Becher, J. Eschner, F. Schmidt-Kaler, and R. Blatt, Bell states of atoms with ultralong lifetimes and their tomographic state analysis, Phys. Rev. Lett. 92, 220402 (2004).
  * Häffner  _et al._ [2005] H. Häffner, W. Hänsel, C. F. Roos, J. Benhelm, D. Chek-al kar, M. Chwalla, T. Körber, U. D. Rapol, M. Riebe, P. O. Schmidt, C. Becher, O. Gühne, W. Dür, and R. Blatt, Scalable multiparticle entanglement of trapped ions, Nature 438, 643 (2005).
  * Gross [2011] D. Gross, Recovering low-rank matrices from few coefficients in any basis, IEEE Transactions on Information Theory 57, 1548 (2011).
  * Tóth  _et al._ [2010] G. Tóth, W. Wieczorek, D. Gross, R. Krischek, C. Schwemmer, and H. Weinfurter, Permutationally invariant quantum tomography, Phys. Rev. Lett. 105, 250403 (2010).
  * Moroder  _et al._ [2012] T. Moroder, P. Hyllus, G. Tóth, C. Schwemmer, A. Niggebaum, S. Gaile, O. Gühne, and H. Weinfurter, Permutationally invariant state reconstruction, New Journal of Physics 14, 105001 (2012).
  * Cramer  _et al._ [2010] M. Cramer, M. B. Plenio, S. T. Flammia, R. Somma, D. Gross, S. D. Bartlett, O. Landon-Cardinal, D. Poulin, and Y.-K. Liu, Efficient quantum state tomography, Nature Communications 1, 10.1038/ncomms1147 (2010).
  * Baumgratz  _et al._ [2013] T. Baumgratz, D. Gross, M. Cramer, and M. B. Plenio, Scalable reconstruction of density matrices, Phys. Rev. Lett. 111, 020401 (2013).
  * Lanyon  _et al._ [2017] B. P. Lanyon, C. Maier, M. Holzäpfel, T. Baumgratz, C. Hempel, P. Jurcevic, I. Dhand, A. S. Buyskikh, A. J. Daley, M. Cramer, M. B. Plenio, R. Blatt, and C. F. Roos, Efficient tomography of a quantum many-body system, Nature Physics 13, 1158 (2017).
  * Palmieri  _et al._ [2020] A. Palmieri, E. Kovlakov, F. Bianchi, D. Yudin, S. Straupe, J. D. Biamonte, and S. Kulik, Experimental neural network enhanced quantum tomography, npj Quantum Inf 6, 20 (2020).
  * Pan and Zhang [2022] C. Pan and J. Zhang, Deep learning-based quantum state tomography with imperfect measurement, International Journal of Theoretical Physics 61, 10.1007/s10773-022-05209-4 (2022).
  * Koutný  _et al._ [2022] D. Koutný, L. Motka, Z. c. v. Hradil, J. Řeháček, and L. L. Sánchez-Soto, Neural-network quantum state tomography, Phys. Rev. A 106, 012409 (2022).
  * Ahmed  _et al._ [2021] S. Ahmed, C. Sánchez Muñoz, F. Nori, and A. F. Kockum, Quantum state tomography with conditional generative adversarial networks, Phys. Rev. Lett. 127, 140502 (2021).
  * Ma  _et al._ [2023] H. Ma, Z. Sun, D. Dong, C. Chen, and H. Rabitz, Attention-based transformer networks for quantum state tomography (2023), arXiv:2305.05433 [quant-ph] .
  * Palmieri  _et al._ [2023] A. M. Palmieri, G. Müller-Rigat, A. K. Srivastava, M. Lewenstein, G. Rajchel-Mieldzioć, and M. Płodzień, Enhancing quantum state tomography via resource-efficient attention-based neural networks (2023), 2309.10616 [quant-ph] .
  * Dawid  _et al._ [2022] A. Dawid, J. Arnold, B. Requena, A. Gresch, M. Płodzień, K. Donatella, K. A. Nicoli, P. Stornati, R. Koch, M. Büttner, R. Okuła, G. Muñoz-Gil, R. A. Vargas-Hernández, A. Cervera-Lierta, J. Carrasquilla, V. Dunjko, M. Gabrié, P. Huembeli, E. van Nieuwenburg, F. Vicentini, L. Wang, S. J. Wetzel, G. Carleo, E. Greplová, R. Krems, F. Marquardt, M. Tomza, M. Lewenstein, and A. Dauphin, Modern applications of machine learning in quantum sciences (2022), arXiv:2204.04198 [quant-ph] .
  * Aaronson and Rothblum [2019] S. Aaronson and G. N. Rothblum, Gentle measurement of quantum states and differential privacy, in  _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_, STOC 2019 (Association for Computing Machinery, New York, NY, USA, 2019) p. 322–333.
  * Huang  _et al._ [2020] H.-Y. Huang, R. Kueng, and J. Preskill, Predicting many properties of a quantum system from very few measurements, Nature Physics 16, 1050 (2020).
  * Altepeter  _et al._ [2004] J. B. Altepeter, D. F. James, and P. G. Kwiat, 4 qubit quantum state tomography, Quantum state estimation , 113 (2004).
  * Haah  _et al._ [2017] J. Haah, A. W. Harrow, Z. Ji, X. Wu, and N. Yu, Sample-optimal tomography of quantum states, IEEE Transactions on Information Theory 63, 5628 (2017).
  * O’Donnell and Wright [2016] R. O’Donnell and J. Wright, Efficient quantum tomography, in  _Proceedings of the Forty-Eighth Annual ACM Symposium on Theory of Computing_, STOC ’16 (Association for Computing Machinery, New York, NY, USA, 2016) p. 899–912.
  * Koh and Grewal [2022] D. E. Koh and S. Grewal, Classical Shadows With Noise, Quantum 6, 776 (2022).
  * Aaronson [2018] S. Aaronson, Shadow tomography of quantum states (2018), arXiv:1711.01053 [quant-ph] .
  * Elben  _et al._ [2020c] A. Elben, R. Kueng, H.-Y. R. Huang, R. van Bijnen, C. Kokail, M. Dalmonte, P. Calabrese, B. Kraus, J. Preskill, P. Zoller, and B. Vermersch, Mixed-state entanglement from local randomized measurements, Phys. Rev. Lett. 125, 200501 (2020c).
  * Huang  _et al._ [2021] H.-Y. Huang, R. Kueng, and J. Preskill, Efficient estimation of pauli observables by derandomization, Phys. Rev. Lett. 127, 030503 (2021).
  * Elben  _et al._ [2023] A. Elben, S. T. Flammia, H.-Y. Huang, R. Kueng, J. Preskill, B. Vermersch, and P. Zoller, The randomized measurement toolbox, Nature Reviews Physics 5, 9 (2023).
  * McGinley and Fava [2022] M. McGinley and M. Fava, Shadow tomography from emergent state designs in analog quantum simulators (2022), arXiv:2212.02543 [quant-ph] .
  * Tran  _et al._ [2023] M. C. Tran, D. K. Mark, W. W. Ho, and S. Choi, Measuring arbitrary physical properties in analog quantum simulation, Phys. Rev. X 13, 011049 (2023).
  * Lanyon  _et al._ [2011] B. P. Lanyon, C. Hempel, D. Nigg, M. Müller, R. Gerritsma, F. Zähringer, P. Schindler, J. T. Barreiro, M. Rambach, G. Kirchmair, M. Hennrich, P. Zoller, R. Blatt, and C. F. Roos, Universal digital quantum simulation with trapped ions, Science 334, 57 (2011), https://www.science.org/doi/pdf/10.1126/science.1208001 .
  * Notarnicola  _et al._ [2023] S. Notarnicola, A. Elben, T. Lahaye, A. Browaeys, S. Montangero, and B. Vermersch, A randomized measurement toolbox for an interacting rydberg-atom quantum simulator, New Journal of Physics 25, 103006 (2023).
  * Gärttner  _et al._ [2017] M. Gärttner, J. G. Bohnet, A. Safavi-Naini, M. L. Wall, J. J. Bollinger, and A. M. Rey, Measuring out-of-time-order correlations and multiple quantum spectra in a trapped-ion quantum magnet, Nature Physics 13, 781 (2017).
  * Gärttner  _et al._ [2018] M. Gärttner, P. Hauke, and A. M. Rey, Relating out-of-time-order correlations to entanglement via multiple-quantum coherences, Phys. Rev. Lett. 120, 040402 (2018).
  * Li and Benjamin [2018] Y. Li and S. C. Benjamin, One-dimensional quantum computing with a ‘segmented chain’ is feasible with today’s gate fidelities, npj Quantum Information 4, 25 (2018).
  * Wurtz  _et al._ [2023] J. Wurtz, A. Bylinskii, B. Braverman, J. Amato-Grill, S. H. Cantu, F. Huber, A. Lukin, F. Liu, P. Weinberg, J. Long, S.-T. Wang, N. Gemelke, and A. Keesling, Aquila: Quera’s 256-qubit neutral-atom quantum computer (2023), arXiv:2306.11727 [quant-ph] .
  * Kim  _et al._ [2023] Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. van den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme, and A. Kandala, Evidence for the utility of quantum computing before fault tolerance, Nature 618, 500 (2023).
  * [117] Details of analytical derivation and a numerical code used for obtaining the Bell correlator for concrete values of $N$ and $r$ are available at the following address: https://github.com/mplodzien/Many-body-quantum-correlations-in-spin-chains.

Supplemental Material for
“Generation of scalable many-body Bell correlations in spin chains with short-
range two-body interactions”
Marcin Płodzień,1 Tomasz Wasak,2 Emilia Witkowska,3 Maciej Lewenstein,1,4 and
Jan Chwedeńczuk5

1ICFO - Institut de Ciències Fotòniques, The Barcelona Institute of Science
and Technology, 08860 Castelldefels (Barcelona), Spain
2Institute of Physics, Faculty of Physics, Astronomy and Informatics, Nicolaus
Copernicus University in Toruń, Grudzia̧dzka 5, 87-100 Toruń, Poland
3Institute of Physics PAS, Aleja Lotnikow 32/46, 02-668 Warszawa, Poland
4ICREA, Pg. Lluis Companys 23, ES-08010 Barcelona, Spain
5Faculty of Physics, University of Warsaw, ul. Pasteura 5, PL-02-093 Warsaw,
Poland

##  I Derivation of formula from Eq. (7)

In this section, we derive Eq. (7) from the main text. We begin by invoking a
general expression for the density matrix of $N$ qubits which is

$\displaystyle\hat{\varrho}_{\rm
in}=\sum_{\vec{s},\vec{s}^{\prime}}\varrho_{\vec{s},\vec{s}^{\prime}}|\vec{s}\rangle\\!\langle\vec{s}^{\prime}|.$
(S1)

The double sum over $\vec{s}$ and $\vec{s}^{\prime}$ is a shortened notation
for $2N$ sums over $s_{1},\ldots,s_{N}=\pm 1$ and
$s^{\prime}_{1},\ldots,s^{\prime}_{N}=\pm 1$. The considered Hamiltonian

$\displaystyle\hat{H}=\sum_{k,l=1}^{N}J_{kl}\hat{\sigma}_{z}^{(k)}\hat{\sigma}_{z}^{(l)},$
(S2)

where $J_{kl}=1$ for $0<|k-l|\leqslant r$, and $0$ otherwise, can be expressed
in the following form

$\displaystyle\hat{H}=\sum_{r=1}^{N-1}\sum_{k=1}^{N-r}\hat{\sigma}_{z}^{(k)}\hat{\sigma}_{z}^{(k+r)},$
(S3)

Since the Hamiltonian contains only $\hat{\sigma}_{z}^{(l)}$ operators, it is
natural to use as the basis the product of eigenstates of $N$ $z$-axis Pauli
matrices. Hence the action of the evolution operator on each such ket replaces
all Pauli operators in Eq. (S3) with a corresponding set of $\vec{s}$ numbers,
i.e.,

$\displaystyle e^{-i\hat{H}\tau}|\vec{s}\rangle$
$\displaystyle=e^{-i\tau\sum_{r=1}^{N-1}\sum_{k=1}^{N-r}\hat{\sigma}_{z}^{(k)}\hat{\sigma}_{z}^{(k+r)}}|\vec{s}\rangle_{z}=$
(S4)
$\displaystyle=e^{-i\tau\sum_{r=1}^{N-1}\sum_{k=1}^{N-r}s_{k}s_{k+r}}|\vec{s}\rangle_{z}\equiv
e^{-i\hat{H}_{\vec{s}}\tau}|\vec{s}\rangle_{z}.$

Therefore, the output state becomes

$\displaystyle\hat{\varrho}_{\rm
out}(\tau)=e^{-i\hat{H}\tau}\hat{\varrho}_{\rm
in}e^{i\hat{H}\tau}=\sum_{\vec{s},\vec{s}^{\prime}}e^{-i(\hat{H}_{\vec{s}}-\hat{H}_{\vec{s}^{\prime}})\tau}\varrho_{\vec{s},\vec{s}^{\prime}}|\vec{s}\rangle\\!\langle\vec{s}^{\prime}|_{z}$
(S5)

Since the correlator $\mathcal{E}_{N}$ is calculated with the product of $N$
operators rising the spin projection along the $x$-axis, it is convenient to
change the basis using

$\displaystyle|s_{k}\rangle_{z}=\frac{1}{\sqrt{2}}(|{+1}_{k}\rangle+s_{k}|{-1}_{k}\rangle),$
(S6)

where the absence of the subscript $z$ indicates that we are working in the
eigenbasis of $x$-axis Pauli operators. The above expression can be expressed
in a compact way as follows

$\displaystyle|s_{k}\rangle=\frac{1}{\sqrt{2}}\sum_{m_{k}=\pm
1}s_{k}^{\frac{1-m_{k}}{2}}|m_{k}\rangle.$ (S7)

Substituting this result into Eq. (S5) we obtain

$\displaystyle\hat{\varrho}_{\rm out}(\tau)$
$\displaystyle=\frac{1}{2^{N}}\sum_{\vec{s},\vec{s}^{\prime}}e^{-i(\hat{H}_{\vec{s}}-\hat{H}_{\vec{s}^{\prime}})\tau}\varrho_{\vec{s},\vec{s}^{\prime}}\times$
$\displaystyle\times\sum_{\vec{m},\vec{m}^{\prime}}|\vec{m}\rangle\\!\langle\vec{m}^{\prime}|\prod_{k=1}^{N}s_{k}^{\frac{1-m_{k}}{2}}s_{k}^{{}^{\prime}\frac{1-m^{\prime}_{k}}{2}}.$
(S8)

The correlator couples the two extreme elements of the density matrix: the one
where all spins are up: $|+1\rangle^{\otimes N}$ with that where all are down
$|-1\rangle^{\otimes N}$. Hence it can be expressed by a single element of
$\hat{\varrho}_{\rm out}(\tau)$, namely the coherence term between these two
elements. This corresponds to setting all $m$’s to $-1$ and all
$(m^{\prime})$’s to $+1$, giving

$\displaystyle\mathcal{E}_{N}$ $\displaystyle=\big{|}\varrho_{+1^{\otimes
N},-1^{\otimes N}}\big{|}^{2}=$ (S9)
$\displaystyle=\big{|}\frac{1}{2^{N}}\sum_{\vec{s},\vec{s}^{\prime}}e^{-i(\hat{H}_{\vec{s}}-\hat{H}_{\vec{s}^{\prime}})\tau}\varrho_{\vec{s},\vec{s}^{\prime}}s_{1}\ldots
s_{k}\big{|}^{2}$

as invoked in Eq. (7) of the main text.

##  II Diagrammatic spin-inversion asymptotic expansion theory

To capture the physical mechanism of the correlations’ enhancement beyond the
universal critical range $r=4$, we resort to an asymptotic diagrammatic
expansion of the many-body state coefficients in the number of spin-inversions
throughout the evolution. The proposed method provides information about the
time dependence of $\mathcal{E}_{N}$, enabling the estimation of the critical
time and indicating the origin of enhanced correlation at the critical range
$r=4$.

We consider first the part of the Hamiltonian from Eq. (1) which couples a
single pair, i.e.,
$\hat{H}_{kl}=J_{kl}\hat{\sigma}_{z}^{(k)}\hat{\sigma}_{z}^{(l)}$. Then, under
the action of $\hat{H}_{kl}$, we obtain

$\displaystyle
e^{-i\tau\hat{H}_{kl}}|1_{k}\rangle|1_{l}\rangle=\hat{f}_{kl}|1_{k}\rangle|1_{l}\rangle,$
(S10)

where
$\hat{f}_{kl}=\cos(J_{kl}\tau)-i\sin(J_{kl}\tau)\hat{\sigma}_{z}^{(k)}\hat{\sigma}_{z}^{(l)}$.
The action of the second term in the operator $\hat{f}_{kl}$ inverts the spins
due to the property $\hat{\sigma}_{z}^{(k)}|1_{k}\rangle=|{-1}_{k}\rangle$
unless the distance $d=|k-l|$ exceeds the range $r$. Since the evolution of
the state is given by a product of the terms from Eq. (S10),

$|\psi(\tau)\rangle=\prod_{k,l}{}\hat{f}_{kl}|\psi_{\mathrm{in}}\rangle,$
(S11)

we can represent the final state diagrammatically as a sequence of dots (each
corresponding to the initial state of the individual spin) and lines
connecting pairs of dots. To each line connecting a pair we assign either the
amplitude $\cos(\tau)$ and an unchanged state of the pair, or the amplitude
$-i\sin(\tau)$ to a pair with inverted spins. The final state is the sum over
all possible assignments to all the lines and all possible ways of connecting
the pairs of spins that are compatible with the distance of the interaction.
In such a case, the nontrivial lines can connect only points that are not
further apart than the distance $r$, and there are in total $K\equiv
r(2N-r-1)/2$ of such lines. We now apply this diagrammatic method to calculate
the $\mathcal{E}_{N}$ at short times and identify at which $r$ the correlator
crosses the Bell limit.

####  II.0.1 Diagrammatic approach

Here, we describe the diagrammatic method allowing for physical insights
regarding the origin of the critical interaction range $r=4$. We begin by
analyzing the number of possibilities of forming the diagrams $P_{r}$ for
small values of $r$. We show that the number of diagram classes in the case of
$r=4$ exceeds those of $r=3$. As a result, new additional diagrams appear for
$r=4$, which describes a correlation that spreads over the whole chain and
thus increases the exponent of $P_{4}(N)$.

Figure S1: (a) An example of diagrams for $r=3$ and $N=10$. Three disconnected
cluster diagrams are visible. (b) The classes of clusters including 2, 4 and 6
spins ($A_{2,4,6}$) for $r=3$. $R^{a/b}_{x\to y}$ and $E_{xy}$ with
$x,y=1,2,3$ mark extended and closing cluster diagrams as explained in the
main text. (c,d) The classes of cluster diagrams with one (two) unconnected
spins are marked with blue (red) points. The green lines close the cluster,
the blue lines exchange the classes between one and two unconnected last
spins. The red line extends the diagram within the same class.

To proceed, we first investigate the case $r=3$, see Fig. S1(a). As can be
seen from the diagram which is an example for $N=10$, there are clusters or
blocks of spins, between which there are no lines. The positions of such
different clusters can be exchanged leading to an exponential scaling with $N$
due to many permutations when $N$ grows. These clusters including $k=2$, 4 and
6 spins, which form classes $A_{2,4,6}$, see Fig. S1(b), can appear anywhere
in the chain.

For comparison, in the case $r=1$ only a single class $A_{2}$ from Fig. S1(b)
contributes, which contains just a single element. As a result, there are
$(N/2)!$ such permutations built from the same diagram, and thus
$P_{1}=(N/2)!/(N/2)!=1$. In the case, $r=2$, in addition to the class $A_{2}$,
the first diagram from the class $A_{4}$ from Fig. S1(b) contributes. Here, we
can perform the summation analytically and, after straightforward
calculations, we obtain $P_{2}(N)=F_{N/2+1}$, where $F_{n}$ is the $n$-th
Fibonacci number. For large $N$, we have $P_{2}\propto e^{\ln\phi\,N/2}$,
which confirms the exponential scaling with the exponent $(\ln\phi)/2\approx
0.241$ determined by the golden ratio $\phi$.

Importantly, in the case of $r=3$, there is a single cluster diagram, the
class denoted by $R_{1\to 1}$ in Fig. S1(b), which can be extended up to
arbitrary $k\leqslant N$ spins by adding a single line shown in red starting
from the blue point. The notation $R_{p\to q}$ is meant to denote a recursive
type of diagram that brings $p$ open points to $q$ open points in the cluster.
This extended cluster diagram can be terminated at any point by adding the
green line as shown in the diagram denoted by $E_{11}$. Interestingly, such a
cluster can include all spins. The appearance of such extended diagrams
increases the exponent of $P_{r}$ to $0.427$, which we found by a simple fit,
and so approximately by a factor of 2 compared to the case $r=2$. This
increase, however, is insufficient to reach the Bell limit.

For the case $r=4$, the situation changes qualitatively. Similarly to the
$r=3$ case, small cluster diagrams appear in the expansion, but now the class
of extended diagrams is much larger. In Fig. S1(b), we show examples of
diagrams which end with one unbound spin (blue). Not only the cluster can be
closed by adding one of the green lines from classes $E_{1,j}$ with $j=1,2,3$,
but also it can be further extended by using the red line from class $R_{1\to
1}$. Alternatively, the diagram can be converted to a cluster ending with two
spins, as in class $R_{1\to 2}$. Next, such diagrams can be closed (using
green lines from classes $E_{2,1/2/3}$ in panel (c)), or converted to a
diagram with one unbound spin (as in panel (b)) by using blue lines from class
$R^{a/b}_{2\to 1}$. As a consequence, by fitting an exponential, we find
$0.563$ for the exponent $P_{4}$. For completeness, we mention that we found
the exponent $0.670$ for $r=5$.

The increase of the exponent of $P_{4}(N)$, relative to the $r=3$ case, is
sufficient to surpass the Bell limit. Such long-range cluster diagrams have a
much larger contribution in the case $r=4$ than for $r=3$. For instance, for
$N=10$, there are 15 different cluster diagrams embodying all the $N$ spins
for $r=4$ compared to a single diagram for $r=3$. As a result, the rate of
increasing $P_{r}(N)$ as a function of $N$ is larger with increasing $r$ and
leads to crossing the Bell correlations limit when $r=3$ is increased to
$r=4$.

####  II.0.2 Expansion of correlator $\mathcal{E}_{N}$

The correlator $\mathcal{E}_{N}(\tau)$ from Eq. (3) can be conveniently
rewritten in the form

$\mathcal{E}_{N}(\tau)=|C_{+}(\tau)|^{2}\,|C_{-}(\tau)|^{2},$ (S12)

for the pure state from Eq. (S11), where $C_{+}$ $(C_{-})$ is the amplitude of
finding all spins up (down) along the $x$-axis. Therefore, in order to
calculate the amplitude $C_{-}$ we need to count all the processes that lead
to the spin-inversion of all $N$ spins, and in order to obtain $C_{+}$ we have
to take into account all the processes that resulted in no net spin-inversion
of the initial state. Our approach is to consider only those diagrams in which
each spin inverted the fewest number of repetitions – which is justified for
short times.

We first determine $C_{-}(\tau)$, focusing on the dominating process where all
the spins have inverted only once. Diagrammatically, each pair of points is
connected only once with a line with the amplitude $-i\sin{\tau}$, and there
are $N/2$ such lines, while all the remaining $K-N/2$ lines have amplitudes
$\cos\tau$. Therefore, the coefficient $C_{-}$ is given by

$C_{-}\approx
P_{r}(N)\big{(}\cos\tau\big{)}^{K-N/2}\big{(}-i\sin\tau\big{)}^{N/2},$ (S13)

where $P_{r}(N)$ denotes the number of possibilities of forming the diagrams.

A similar expansion holds for $C_{+}$, but now we count the diagrams with the
least number of lines with the amplitude $-i\sin{\tau}$ that result in no net
spin inversion. In the zeroth order, no spin inversions occur resulting in the
scaling of $\cos^{K}\tau$. In the first order, such diagrams contain three
lines where the three spins $i$, $j$, and $k$ are mutually connected and
subject to the constraint that all the distances between the pairs fall within
the range of the interaction. By the symbol $R_{r}(N)$, we denote the number
of such diagrams for the given range $r$ and the number of spins $N$. We thus
arrive at the following expression

$C_{+}\approx\big{(}\cos\tau\big{)}^{K}+R_{r}(N)\big{(}\cos\tau\big{)}^{K-3}\big{(}-i\sin\tau\big{)}^{3}.$
(S14)

For $r=1$ we only have a single diagram with only nearest-neighbours being
connected, and, thus, $P_{r=1}=1$ and $R_{r=1}=0$. As a result, we recover the
exact formula cited in the main text, namely

$\displaystyle\mathcal{E}_{N}(\tau)=\sin^{N}(\tau)\cos^{3N-4}(\tau).$ (S15)

For small values of $r$ or $N$, the actual numbers for $P_{r}$ and $R_{r}$ can
be found by counting the number of permutations or from analytical solutions
[code_repo]. In the large $N$-limit, $R_{r}$ scales linearly with $N$, and
$P_{r}$ scales exponentially with an $r$-dependent exponent. This exponential
growth of $P_{r}$ with a sufficiently large exponent is responsible for
surpassing the Bell limit and observation of nonlocal correlations.

In Fig. S2(a), we compare our asymptotic expansion of $\mathcal{E}_{N}$ from
Eq. (S12) calculated with Eq. (S13) and Eq. (S14) for $N=16$ and the range
$r=3$ (red), 4 (blue), 5 (black). From these examples, we see that the
proposed approach via the asymptotic spin-inversion expansion works
sufficiently well for the correlator $\mathcal{E}_{N}$ as it captures the
short time behaviour as well as the approximate time of crossing the Bell
limit for $r\geqslant 4$. For $r=3$ the theory quantitatively shows that the
limit is not surpassed.

Figure S2:  (a) The comparison of exact (solid) and approximate (dashed)
numerics for $N=16$ and $r=3$ (red), 4 (blue), 5 (black). The parameters used
are $(P_{r},R_{r})=(491,40)$ for $r=3$, (3116, 76) for $r=4$ and (12483, 120)
for $r=5$. (b) The behaviour of $\tilde{\mathcal{E}}_{\mathrm{max}}$ as a
function of the range $r$ for various values of $N=12$, 14, 16, 18, 20. The
solid line is a guide to the eye. The value of
$\tilde{\mathcal{E}}_{\mathrm{max}}$ is positive for $r=4$ which signals the
existence of $\tau_{\mathrm{crit}}$ for the Bell inequality violation as given
by Eq. (S19).

####  II.0.3 Gaussian approximation

To gain further insight into the expansion method, we approximate the powers
of cosines by Gaussian functions, sines with the first order in $\tau$, and
$K-3\approx K$ in the second term of $|C_{+}|^{2}$, which is valid for large
$N$. In this way, we obtain

$\displaystyle|C_{-}|^{2}$ $\displaystyle\approx
P_{r}^{2}\tau^{N}e^{-(K-N/2)\tau^{2}},$ (S16) $\displaystyle|C_{+}|^{2}$
$\displaystyle\approx e^{-K\tau^{2}}(1+R_{r}^{2}\tau^{6}).$ (S17)

Furthermore, in the region of $\tau$ where the Bell inequality is violated,
and around the maximum of $\mathcal{E}_{N}$ we can neglect the first term in
$|C_{+}|^{2}$ which yields our intermediate-time approximation

$\displaystyle\mathcal{E}_{N}$ $\displaystyle\approx
P_{r}^{2}R_{r}^{2}e^{-\beta_{N}\tau^{2}}\tau^{N+6},$ (S18)

where we introduced the exponent $\beta_{N}=2K-N/2$.

Now, we estimate the time of crossing the Bell limit, and, as a result, we
show that the approximate expression in Eq. (S18) describes the observed
numerical behaviour, i.e., that for $r=3$ the Bell threshold is not surpassed
and for $r\geqslant 4$ we observe the violation of the inequality. However,
our approach, being asymptotic in nature, does not capture accurately the
critical time $\tau_{c}$, where the Bell limit is surpassed, due to its
systematic underestimation which manifests in a small offset. In Fig. S2, this
offset is not visible but it can be larger for larger values of $r$.

To determine $\tau_{c}$, we expand $\ln(2^{N}\mathcal{E}_{N})$ around maximum,
which is reached at the time $\tau_{\mathrm{max}}=\sqrt{(N+6)/(2\beta_{N})}$,
up to quadratic terms in $\tau-\tau_{\mathrm{max}}$. The Bell limit
$\mathcal{E}_{N}=1/2^{N}$ is reached at

$\tau_{c}\approx\tau_{\mathrm{max}}-\sqrt{\frac{\tilde{\mathcal{E}}_{\mathrm{max}}}{2\beta_{N}}},$
(S19)

which is meaningful, i.e., $\tau_{c}\in\mathbb{R}$, only when
$\tilde{\mathcal{E}}_{\mathrm{max}}\geqslant 0$, and where the maximum of the
scaled logarithm is

$\frac{\tilde{\mathcal{E}}_{\mathrm{max}}}{N}=\frac{2\ln(P_{r}R_{r})}{N}-\bigg{(}1+\frac{6}{N}\bigg{)}\ln\frac{1}{\tau_{\mathrm{max}}}+\ln
2-\frac{1+\frac{6}{N}}{2}.$ (S20)

Due to the linear scaling of $R_{r}$ with $N$, the first term on the right-
hand side extracts the exponent of $P_{r}$ in the first approximation. The
second term, which decreases $\tilde{\mathcal{E}}_{\mathrm{max}}$, comes from
the factor $\tau^{N}$ in the expansion of $C_{-}$ originating from the
requirement that all the $N$ spins have to be inverted during the dynamics. In
Fig. S2(b), we present the function $\tilde{\mathcal{E}}_{\mathrm{max}}$ as a
function of the range $r$ for various $N$’s. The time $\tau_{\mathrm{crit}}$
from Eq. (S19), when the Bell threshold is reached by $\mathcal{E}_{N}$,
exists if $\tilde{\mathcal{E}}_{\mathrm{max}}\geqslant 0$. We find that for
$r\leqslant 3$, $\tilde{\mathcal{E}}_{\mathrm{max}}$ is negative and for
$r\geqslant 4$ it is positive resulting in a physically meaningful $\tau_{c}$.

Physically, $\mathcal{E}_{N}(\tau)$ quantifies how fast the initial state,
given by $|C_{+}(\tau)|^{2}$, is depleted, and in order to surpass the Bell
inequality, the state with all the spins inverted, given by
$|C_{-}(\tau)|^{2}$, has to be populated fast enough. The population in
$|C_{-}|^{2}$ is exponentially enhanced as $r$ increases due to a larger
number of diagrams with spin-pairs inverted by interaction. On the other hand,
$\tau_{c}$ decreases with increasing $r$ causing significant reduction
$\propto\tau^{N}$ of $|C_{-}(\tau)|^{2}$. The two processes compete, and in
order to surpass the Bell limit, $P_{r}$ has to be large enough. We note that
$R_{r}$ also plays a role in the above considerations, since the first term in
the expansion of $|C_{+}|^{2}$, which leads to Gaussian decay, is insufficient
to describe the violation of the Bell limit.

###  II.1 Measuring many-body Bell correlator with the Multiple Quantum
Coherences

The multiple quantum coherences (MQC) technique provides extensive information
about the structure of a many-body state. It allows for the relation of the
quantum correlators with other physical quantities like the out-of-time-order
correlations [garttner2017measuring, PhysRevLett.120.040402]. The MQC is
defined upon picking some observable $\hat{A}$ and expressing a density
operator in terms of the eigenstates and eigenvalues of $\hat{A}$ as follows:

$\hat{\varrho}=\sum_{m}\sum_{\lambda_{i}-\lambda_{j}=m}\varrho_{ij}|\psi_{i}\rangle\\!\langle\psi_{j}|\equiv\sum_{m}\hat{\varrho}_{m}$
(S21)

where $\hat{A}|\psi_{i/j}\rangle=\lambda_{i/j}|\psi_{i/j}\rangle$. The MQC is
defined as a norm of the fixed-$m$ part of the density matrix, namely
$I_{m}(\hat{\varrho})={\rm
Tr}\left[\hat{\varrho}^{\dagger}_{m}\hat{\varrho}_{m}\right]$. In the multi-
qubit case considered in the current work, we consider

$\hat{A}=\frac{1}{2}\sum_{k=1}^{N}\hat{\sigma}_{x}^{(k)}$ (S22)

and $m=N/2$ to obtain $I_{N}(\hat{\varrho})=\mathcal{E}_{N}$. The
$I_{N}(\hat{\varrho})$ is directly accessible in the laboratory, as shown in
Ref. [garttner2017measuring] for $N=6$ qubits.

For a more detailed discussion of the physical significance of the MQCs, their
measurements, and the relation to the out-of-time-order correlations, see Ref.
[PhysRevLett.120.040402] and the references therein.

## References