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

Spin chains, the archetypal many-body quantum systems, are used for prototype quantum computing, ultra-precise sensing, or as quantum simulators. Hence, they attract theoretical and experimental attention, focusing on analyzing their dynamic and static properties, their structure, and the “quantum content”, namely how — and in what sense — a particular chain is non-classical. We address these points and show how to generate a potent resource for quantum technologies, namely the many-body Bell correlations in spin chains, with controllable short-range two-body interactions. Subsequently, we classify the depth of produced Bell correlations. We identify a critical range necessary to generate many-body Bell correlations in the system and provide the physical mechanism behind this critical behaviour. Remarkably, the critical range is short and universal for every length of the chain, establishing the viability of this protocol. In contrast, the critical time, at which the many-body Bell correlations emerge in the system depends on its size, as our analysis reveals. Importantly, we show, that these Bell correlations are fully determined by just a single element of the density matrix, and can be measured by the existing state-tomography methods. The fully analytical findings presented here provide novel insight into the methodology of generating strong many-body non-classical correlations with short-range two-body interactions, offering promising prospects for quantum technologies.


I. INTRODUCTION
It is hard to fully comprehend how inventions like lasers, semiconducting transistors, or magnetic-resonance devices changed our world.They are all fruits of the advancing knowledge and control over many microscopic processes and the last century is often referred to as the "first quantum revolution".Nowadays, emerging quantum technologies push the boundaries of the precise manipulation and control of quantum systems even further-to the level of single atoms, photons, and molecules.They are expected to drive the "second quantum revolution", changing medicine, computation, information storage and transfer, or material science in unimaginable way [1][2][3].
Many of these solutions can be implemented using the archetype of a many-body quantum system schematically shown in Fig. 1, namely a chain of two-level objects, such as spin-1/2 particles (qubits).What makes these systems so much more potent than their classical analogues is the richness of quantum resources.These are, above all, the quantum coherence (at the single-particle level) and the triad of many-body non-classical correlations: entanglement [4], Einstein-Podolsky-Rosen steering [5] and the Bell nonlocality [6].To have a running quantum device, we must master the creation and manipulation process ticularly promising method of the dynamical generation of a many-body quantum-correlated state utilizing twobody interactions in an Ising-type spin chain starting from a fully separable state of N spins.We focus on generating the strongest, most potent, and most exotic form of quantum correlation, namely the Bell correlations.Nevertheless, all the methods presented here are suitable for the characterization of other members of the quantum triad.
Most importantly, we demonstrate that the range of the two-body interactions plays a crucial role in generating Bell correlations, and our analysis reveals the existence of a critical range, above which they are created.We show that this critical range is independent of the number of spins, hence it is universal for all sizes of spin chains.We employ diagrammatic techniques to grasp the physical mechanism behind this critical behaviour.We also characterize the depth of Bell correlations, namely how many spins exhibit non-classical correlations depending on the range of interactions.We show that the strength of these correlations can be fine-tuned to the level desired for a particular application.Also, we identify the critical time at which these correlations emerge and relate this instant of time to the interaction range in the form of a scaling law.Last but not least, we show that the many-body Bell correlations in spin chains can be extracted from the measurement of just a single element of the density matrix.By means of the quantum state tomography, we demonstrate that the measurement of this element is within the reach of the modern experimental tomographic methods.
Our findings suggest one possible route towards scalable quantum devices reinforcing experimental effort on generating entangled cluster states, which support manybody correlations of spin chains on a global scale [7][8][9][10].Permutationally invariant Bell inequalities and entanglement criteria for many body systems based on the first and second-order moments of local observable were intensively studied since [11,12].These were generalized to translationally invariant inequalities for up to 10 parties, [13] and to studies of entanglement depth [14,15], and Bell correlations depth [16].Energy (in particular in spin chains) was used as a Bell correlations detector [17].Similar Bell inequalities were used to detect criticality in simple spin chains [18].Finally, this approach has been generalized to arbitrary qdits (spin) systems [19][20][21].Individual addressing and control of single spin are offered by experimental platforms composed of trapped ions [8], ultra-cold molecules [22], atoms in optical lattices [23] or Rydberg tweezers arrays [24].It has been shown that all-to-all interactions, which can be simulated in the one-dimensional Bose-Hubbard models, generate many-body Bell correlations of arbitrary depth [25][26][27].Among the others, the long-range XY spin model also provides interactions that favour spin alignment and stabilize collective behaviour [28][29][30][31].Recent experiments demonstrated scalable two-body correlations (in terms of spin squeezing) using a 2D array of trapped ions [32] and Rydberg atoms [33] when simulating long-range XY spin model.Apart from the experimental realization of the long-lived Bell states [34][35][36][37][38], a generation of scalable entanglement using Ising spin chains is still awaiting proof-of-principle experiments.
We believe that a thorough analysis of the interaction range's role in generating quantum correlations presented in this work could advance the journey towards this direction.

II. SPIN CHAINS AND BELL CORRELATIONS
A spin chain is a collection of N two-level quantum systems, typically equally distributed along a line.Each entity is characterized by a triad of Pauli matrices σ(k) i , with i = x, y, z and k = 1 . . .N .When spins are connected by two-body interactions, the dynamics of the chain can be realized in many different ways, depending on the architecture of the experimental setup.Here, we consider the natural generalization of the Ising model to the case when spins are allowed to interact not only with their nearest neighbors but also with more distant members of the chain.In the absence of external magnetic fields, the evolution of such a system is described by the following Hamiltonian where J kl is the coupling strength (in units of ℏ) of the k-th spin interacting with its l-th partner.Recent experiment [9] showed a high degree of control over the distance selective interactions J kl .Here, we take the simplest finite-range interaction potential with a single length-scale given by the range of the interaction, i.e., the rectangular function The interactions range r takes the values from r = 1 (nearest-neighbor couplings) to r = N − 1 (all-to-all couplings), see the scheme in Fig. 1.
A natural way of using the Hamiltonian from Eq. (1) to generate quantum correlations is to take N uncorrelated spins, polarized in a direction orthogonal to the z-axis distinguished by the Hamiltonian, for instance and let the interactions correlate the particles.Here, |1⟩ denotes an eigenstate of the x-axis Pauli operator with positive eigenvalue.Regarding this procedure, there are a number of questions that arise.Is there some minimal value of the interaction range r, at which quantum features emerge or become more pronounced?What is the relation between the range and the pace at which these correlations are generated?Or, how should J kl depend on the distance |k − l| to maximize the effect of the dynamics?All these points can be exhaustively addressed using the correlation function where the rising operators are taken along the x-axis, i.e., σ(k) . The E N or, equivalently, Q N (we use these quantities interchangeably, depending on the interpretation convenience) is called the Bell correlator, because E N ⩽ 2 −N (or Q ⩽ 0) is the many-body Bell inequality, see Refs.[39][40][41][42][43][44] and Appendix A for details.Consequently, the observation of E N > 2 −N (or Q N > 0), signals the presence of many-body Bell correlations in the system.Furthermore, the larger the inequality violation, the stronger the correlations as quantified by the nonlocality depth (see discussion below).
For the initial state as in Eq. ( 3), evolving under the Hamiltonian from Eq. ( 1), the Bell correlator is where τ is the evolution time and H ⃗ s is the c-number counterpart of Eq. ( 1) with operators σ(k) z replaced by s k (see Appendix B for the derivation of this formula).The sums over ⃗ s and ⃗ s ′ run over 2 N combinations of s 1 , . . ., s N and s ′ 1 , . . ., s ′ N which take values ±1.At first glance, Eq. ( 5) seems uncomplicated because sums over ⃗ s and ⃗ s ′ separate.Nevertheless, apart from the case when the Hamiltonian (1) couples only the nearest neighbours, no closed formula for Eq. ( 5) is known.Still, extended information about the many-body Bell correlations can be analytically extracted from the above expression, as we show in the following.

III. RESULTS
For the aforementioned r = 1 case, the Bell correlator is given by the analytical formula see Appendix C for the derivation.This expression has a maximum max τ (E N ) ≈ 2 −1.6N +1 at τ ≃ π/6, which is exponentially smaller than the Bell limit 2 −N .Therefore, the nearest-neighbors interactions are not capable of generating many-body Bell correlations in the scheme considered here.The other extreme case is the all-to-all type of interactions with r = N − 1, which was considered in detail in [26].The method of generating quantum correlations by means of Eq. ( 1) is then called the one-axis twisting and the Bell correlations are present starting from τ ≃ 1.77/N , and the correlator reaches its maximal value E N = 1/4 at τ = π/4, where the N -body Greenberg-Horne-Zeilinger (GHZ) state is formed [25,27,45,46].
In Fig. 2(a), we display the time evolution of the correlator Q N calculated using Eq. ( 5) for N = 64 spins and various ranges of interaction r.We observe that the correlator breaks the Bell limit, i.e., Q N > 0, only when r ⩾ 4.This is the first indication for 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 i.e., maximized Q N with respect to τ calculated for different N 's.In Fig. 2(b) we show Q max N , as a function of the number of spins N and for various values of r.These maximal values, for a particular r, lie on a straight line which translates to the exponential scaling of E N , namely E N ∝ 2 γN .In panel Fig. 2(c), we display the value of this exponent as a function of the interaction range.
The observed behaviour of Q max N and γ(r) > 0 for r ⩾ 4 demonstrates that r = 4 is indeed a critical range, at which the Bell correlations are detected.Crucially, Q max N increases with growing N , indicating scalable Bell correlations in this system.Remarkably, although r = 4 in the limit N ≫ 1 is a local interaction encompassing an intensive number of neighboring spins, the system still exhibits an increasing degree of violation of the Bell inequality.Finally, the invariance of the critical range with respect to the number of spins proves that this range is universal.In Section IV, we discuss the origin of the universal range using the diagrammatic methods and the asymptotic expansion theory [47].In Appendix D, we discuss the many-body Bell correlations generation for power law interactions decaying as J kl = |k − l| −α with α > 0. We merely mention here that in this typical setup, the potential cannot fall off too rapidly with the distance to observe the Bell correlations, namely, the critical exponent of the potential leads to the bound α ≲ 1.5 [28,32,33].

A. Bell correlations depth as a function of range r
The value of the correlation function Q N provides information about the structure of the quantum state of the spin chain.While Q N > 0 signals the presence of Bell correlations among the spins, the extent to which the bound is surpassed informs about the "Bell correlation depth" -namely, what is the fraction of spins that, for this particular value of Q N can still be modeled with local hidden variables.Specifically, if where ν ∈ N, then up to N − ν − 3 spins are not Bellcorrelated.Hence the Bell correlations depth is ν + 3.
Consequently, for the maximally entangled GHZ state, for which, according to Eq. ( 4) we obtain Q N = N −2, the depth is N , i.e., all the spins exhibit the Bell correlations.In Fig. 3(a) we display the dependence of the maximal fraction of uncorrelated spins, i.e., n = (N − ν − 3)/N as a function of the normalized range r/N , deduced from Q max N .The result approximately follows the line n = 1 − r/N , shown with a dashed black line.Hence, the Bell correlations depth is naturally related to the range of the spin-spin interactions -the higher the r, the more particles form a Bell-correlated cluster.

B. Critical time
The interaction range r not only determines the Bell correlations depth in a spin-chain but also the time at which these correlations emerge in the system.This is hinted by the Fig. 2(a) -the time τ at which the correlator Q N crosses the Bell limit depends on r, and the generation of correlations seems to accelerate for higher ranges.To identify this dependence, consider the expression given by Eq. ( 5).Take some internal spin, which couples with 2r neighbours.The sum over s k = ±1 will vanish unless the corresponding phase-term e −iτ s k (s k−r +...+s k+r ) oscillates quickly enough.The sum in the parenthesis is at best equal to 2r, hence the phase factor will vary significantly between s k = −1 and Thus the Bell correlator becomes significantly non-zero if τ ≳ β/r, where β is some constant.This is confirmed by the exact solution of the dynamics generated by Eq. ( 1).We take different N 's and identify the critical time when Q N overpasses the Bell limit.By changing the range from the critical value r crit = 4 to r = 10 we recover the τ crit ∝ 1/N dependence, see the right panel in Fig. 3(b).

C. Many-body Bell correlator measurement
The expectation value of the correlator E N in Eq. ( 4) can be expressed in terms of a single element of the density matrix ρ, namely where the correlator describes the coherence between the states when all the spins are up and all the spins are down, with respect to the x-axis Pauli matrix of the spins.
As such, the problem of the many-body Bell correlator measurement is shifted to performing the quantum state tomography (QST).For a given system size N , interaction range r and time τ , we prepared L = 10 reconstructions of the target density matrix ρ(τ ) via classical shadows tomography technique [48].In each l-th reconstruction, l = 1, . . ., L, we prepared M = 10 circles) with standard deviation (shaded areas) for N = 4, 6, 8, 10 spins with the interaction range r = 4, and all-to-all couplings when r = N − 1.The linear scaling of a number of classical shadows M with system size N allows for estimating the many-body Bell correlator with a very good quantitative agreement with analytical predictions for all-to-all interactions and a good qualitative agreement for r = 4.The observed loss of precision for r = 4 is due to the estimation of very small values of E N , which require more reconstructions.The manybody Bell correlations can be measured in experimental quantum simulator platforms [49] like Rydberg atoms array [9,10].

IV. 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 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., Ĥkl = J kl σ(k) z σ(l) z .Then, under the action of Ĥkl , we obtain where fkl = cos(J kl τ ) − i sin(J kl τ )σ Since the evolution of the state is given by a product of the terms from Eq. (10), 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(τ ) and an unchanged state of the pair, or the amplitude −i sin(τ ) to a pair with inverted spins.
For more information about the diagrams see the discussion in Appendix F. 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 ≡ r(2N − r − 1)/2 of such lines.We now apply this diagrammatic method to calculate the E N at short times and identify at which r the correlator crosses the Bell limit.

A. Expansion of correlator EN
The correlator E N (τ ) from Eq. ( 4) can be conveniently rewritten in the form for the pure state from Eq. (11), where C + (C − ) is the amplitude of finding all spins up (down) along the xaxis.Therefore, in order to calculate the amplitude C − we need to count all the processes that lead to the spininversion 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 − (τ ), 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 τ , and there are N/2 such lines, while all the remaining K − N/2 lines have amplitudes cos τ .Therefore, the coefficient C − is given by 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 τ that result in no net spin inversion.In the zeroth order, no spin inversions occur resulting in the scaling of cos K τ .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 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 obtain the exact formula from Eq. ( 6).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 [50].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 Appendix F, we present a detailed discussion of the diagrammatic approach.In Fig. 5(a), we our asymptotic expansion of E N from Eq. ( 12) calculated with Eq. ( 13) and Eq. ( 14) 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 E N as it captures the short time behaviour as well as the approximate time of crossing the Bell limit for r ⩾ 4. For r = 3 the theory quantitatively shows that the limit is not surpassed.

B. 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 τ , and K − 3 ≈ K in the second term of |C + | 2 , which is valid for large N .In this way, we obtain Furthermore, in the region of τ where the Bell inequality is violated, and around the maximum of E N we can neglect the first term in |C + | 2 which yields our intermediate-time approximation where we introduced the exponent β 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. (17) describes the observed numerical behaviour, i.e., that for r = 3 the Bell threshold is not surpassed and for r ⩾ 4 we observe the violation of the inequality.However, our approach, being asymptotic in nature, does not capture accurately the critical time τ crit due to its systematic underestimation which manifests in a small offset.In Fig. 5, this offset is not visible though but it can be larger for larger values of r.
To determine τ crit , we expand ln(2 N E N ) around maximum, which is reached at the time τ max = (N + 6)/(2β N ), up to quadratic terms in τ −τ max .The Bell limit which is meaningful only when Ẽmax ⩾ 0, and where the maximum of the scaled logarithm is given by (19) 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 Ẽmax , comes from the factor τ N in the expansion of C − originating from the requirement that all the N spins have to be inverted during the dynamics.In Fig. 5(b), we present the function Ẽmax as a function of the range r for various N 's.The time τ crit from Eq. ( 18), when the Bell threshold is reached by E N , exists if Ẽmax ⩾ 0. We find that for r ⩽ 3, Ẽmax is negative and for r ⩾ 4 it is positive resulting in a meaningful τ crit .Physically, E N (τ ) quantifies how fast the initial state, given by |C + (τ )| 2 , is depleted, and in order to surpass the Bell inequality, the state with all the spins inverted, given by |C − (τ )| 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, τ crit decreases with increasing r causing significant reduction ∝ τ N of |C − (τ )| 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.

V. CONCLUSIONS
In this work, we presented a systematic analysis of the dynamical process triggered by the Ising-type Hamiltonian from the point of view of its suitability as a source of many-body Bell correlations.Taking an uncorrelated state as a starting point, we have derived the analytical expression for the Bell correlator that allows for the extended understanding of the process and for the numerical simulations not restricted by the dimensionality of the Hilbert space.Though this formula seems to be uncomplicated, it requires advanced diagrammatic methods to grasp the mechanism behind the most important result of this work-namely, that the range of the twobody interactions crucially influences the rate of creation of Bell correlations and its depth.Most of all, we have identified a critical range where a spin interacts with four neighbours at each side, which is the minimal distance of non-vanishing interactions that ensure Bell correlations in the system.Additionally, the feasibility of experimental measurement of the many-body Bell correlator was positively examined by us using QST in the current physical platforms realizing spin-1/2 chains.
It is important to note that our analysis reveals the universality of this critical range-it is size-independent.As the number of spins N grows, the critical range becomes negligible concerning the system size.Yet, the degree of the violation of the many-body Bell inequality grows with the number of spins.Hence, the necessary interaction range, though more extended than the standard nearest-neighbor case, is in fact, local for large N .
We hope these results will indicate the modus operandi in some modern quantum-technological setups.The possibility to control and manipulate the archetypal manybody quantum system, as the spin chain, allows for the creation of strong correlations necessary to advance the technology toward the "second quantum revolution".

AUTHOR CONTRIBUTIONS
M.P. indicated the existence of the critical interaction range, performed many-body simulations, analytical calculations in Ref. [50], 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.
Since the Hamiltonian contains only σ(l) z 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. (B3) with a corresponding set of ⃗ s numbers, i.e., Therefore, the output state becomes Since the correlator 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 where the abscence 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 Substituting this result into Eq.(B5) we obtain The correlator couples the two extreme elements of the density matrix: the one where all spins are up: | + 1⟩ ⊗N with that where all are down | − 1⟩ ⊗N .Hence it can be expressed by a single element of ρout (τ ), namely the coherence term between these two elements.This corresponds to setting all m's to −1 and all (m ′ )'s to +1, giving as invoked in Eq. ( 5) of the main text.

Appendix C: Nearest-neighbours interaction
For the nearest-neighbours case and taking a pure state as an input the correlator from Eq. (B9) can be written as a product of two coefficients, namely where We focus on the line (C3a) and perform the sum over s = 2 cos τ (e −iτ s3 + e iτ s3 ) (C5) We now see that the last sum has the same form as that of the r.h.s. of Eq. (C4), however now starting from s 3 rather than s 2 .Hence the summing procedure can be iterated up until s N −1 giving 2 cos τ N − 2 times.The last term is Hence we obtain that Next, we calculate the second coefficient contributing to the Bell correlator, see Eq. (C3b).Again, we focus on the sum and first perform the summation over s 1 to get The minus sign between the two exponents, which is due to the presence of the s 1 term, is a crucial difference with respect to previous case.We now perform the summation over s 2 to get Hence we recover the structure of Eq. (C8).Therefore the Bell coefficient will be an alternating pattern of sines and cosines.If N is even, then the last term will read If it is odd, then the last term is Hence the result is non-zero only if N is even, giving We conclude that the Bell coefficient for nearestneighbors interactions is zero for odd N , and for even N it is as reported in the main text.
Appendix D: Algebraic decay of J kl with distance We now relate our findings for finite-range interactions to the case of long-range, algebraically decaying interaction potentials that were vastly studied theoretically and experimentally in the literature.These types of interactions appear naturally in many contexts, for instance, if J kl comes from the dipole-dipole interactions, then J kl ∝ |k − l| −α , with α ∈ [3,6], depending on the type of dipoles [8,[22][23][24].It can also be engineered for any other values with trapped ions [52].It was already understood for XY model that the generation of at least two-body correlations requires α < 5d/3 in d dimensions [28].The long-range Ising model follows the same laws and the same theoretical prediction is valid here.It is confirmed by a study of propagation of two-body correlations triggered by long-range potentials, and Lieb-Robinson bound in Ising model for the well-defined casual region outside of which correlations are exponentially suppressed [53][54][55].
The experiments with one-dimensional chains of trapped ions confirmed an acceleration of the correlation transfer through the chain [52] for long-range interactions with α ⩽ 1.
In Fig. 6(a) we plot the maximum of the Bell correlator Q max N for N ranging from 8 to 16 and with α = 0, 1, 2, 3, 4 and 5. Indeed, we observe that when α ⩾ 2, the maximum is below zero, signalling the absence of Bell correlations.Nevertheless a more detailed study of the correlator exponent γ shows that even for α ≲ 5/3 the Bell correlations persist, see Fig. 6(b).Hence, the conclusion is that to create many-body nonlocal states in spin chains, the interactions must be tailored to satisfy the condition α ≲ 5/3 for these long-ranged potentials.
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. 7(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 ∝ e ln ϕ N/2 , which confirms the exponential scaling with the exponent (ln ϕ)/2 ≈ 0.241 determined by the golden ratio ϕ.Importantly, in the case of r = 3, there is a single cluster diagram, the class denoted by R 1→1 in Fig. 7(b), which can be extended up to arbitrary k ⩽ N spins by adding a single line shown in red starting from the blue point.The notation R p→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 ap-pear in the expansion, but now the class of extended diagrams is much larger.In Fig. 7(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→1 .Alternatively, the diagram can be converted to a cluster ending with two spins, as in class R 1→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→1 .As a consequence, by an exponential, we find 0.563 for the exponent P 4 .For completeness, we find 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.

FIG. 2 .
FIG. 2. (a) The Bell correlator QN (τ ) for N = 64 as a function of time τ for the interaction range: r = 1 (blue), r = 2 (light-green) r = 3 (dark-green), r = 4 (red), r = 5 (orange) and r = 63 (violet).Thin lines are for r = 61 and 62.The first maximum of QN with respect to τ is marked by a point (dot, triangle, diamond, star and cross for r ∈ [1, 5], respectively).The values QN > 0 mark the region where the many-body Bell inequality is violated.(b) The first maximum of QN with respect to time, Q max N , as a function of the total number of spins N for various r ∈ [1, 5], as marked by the points on the right corresponding to the panel (a).(c) The exponent of the Bell correlator in the scaling with N , approximated as Q max N ≈ γN + const, or, equivalently, 2 N maxτ EN (τ ) ∝ 2 γN .The growth of scalable many-body Bell correlations with N is manifested by positive γ when r ⩾ 4.

FIG. 3 .
FIG. 3. (a)The maximal fraction of uncorrelated spins, n = (N − ν − 3)/N as a function of the normalized range r/N , see Eq.(8).Different curves correspond to r = 4 . . .N − 1 for N = 8 . . .28 with four-spin increment (the longer the curve, the higher the number of spins).For reference, the dashed black line shows the n = 1−r/N dependence.(b) The critical time at which the Bell correlator passes the Bell limit as a function of r = 4 . . . 10 and for different N : 16 (red), 24 (orange), 32 (blue) and 40 (green), shown in log-log scale.For reference, the 1/r curve is shown with a black dashed line.

4 FIG. 4 .
FIG. 4.The dynamics of the estimated many-body Bell correlator QN 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 stand for the value of the reconstructed many-body Bell correlator from the classical shadows tomography.Standard deviation is marked as a shaded area.
(k) z σ(l) z .The action of the second term in the operator fkl inverts the spins due to the property σ(k) z |1 k ⟩ = |−1 k ⟩ unless the distance d = |k − l| exceeds the range r.

s2=±1 s 2
FIG. 6.(a) The maximum of the normalized Bell correlator Q max N in the range τ ∈ [0, 0.1] versus N , with J kl = |k − l| −α for different values of α = 0, 1, 2, 3, 4 and 5 (from top to bottom lines).(b) The rate of change in the maximum of the normalized Bell correlator, approximated as Q max N ≈ γN + const, for given α.The growth in the many-body Bell correlations with N is indicated by the change in the sign of γ.