Bounds on detection of Bell correlations with entangled ultra-cold atoms in optical lattices under occupation defects
Tanausú Hernández Yanes Youcef Bamaara Alice Sinatra Emilia Witkowska






















§ ABSTRACT 
Bell non-locality stems from quantum correlations effectively identified using inequalities. 
Spin chains, simulated with ultra-cold atoms in optical lattices, Rydberg atoms in tweezer arrays, trapped ions, or molecules, allow single-spin control and measurement. Therefore, they are suitable for studying fundamental aspects of these correlations and non-locality. 
Occupation defects, such as vacancies or multiple atoms occupying a single site due to imperfect system preparation, limit the detection of Bell correlations. We investigate their impact using a simplified toy model parameterized by the probability of a site being singly occupied.
We derive the corresponding Bell inequality and identify the smallest probability that establishes a lower bound for detecting Bell correlations.
We relate the bound to two physical parameters leading to defects in occupations: non-zero temperature and filling factor, focusing on entangled ultra-cold atoms in optical lattices. 
Finally, we numerically validate the predictions of the toy model by full many-body simulations.



Bounds on detection of Bell correlations with entangled ultra-cold atoms in optical lattices under occupation defects

    E. Witkowska
    September 9, 2024
=======================================================================================================================



*These authors contributed equally to this workfootnote





§ INTRODUCTION



Quantum mechanics introduced groundbreaking concepts of non-local correlations and entanglement that challenged well-established principles of classical physics, including realism, causality and locality <cit.>. 
 A classical picture can be recovered at the price of introducing local hidden variables, and the corresponding (LV) theories are shown to satisfy Bell inequalities <cit.>. 
Quantum correlations that violate Bell inequalities are inconsistent with LV theory and are, therefore, referred to as non-local <cit.>.

The most robust method for entanglement certification is provided by violating Bell's inequalities, as it avoids assumptions about the physical nature and degrees of freedom to be measured or the calibration of measurements <cit.>. 
This is known as the device-independent scenario, which is a powerful resource in many quantum information tasks, such as self-testing <cit.>, randomness amplification and expansion <cit.>, quantum key distribution <cit.>, and quantum sensing and metrology <cit.>.
Recent experimental demonstrations have employed various platforms, including entangled photons <cit.>, spins in nitrogen-vacancy centers <cit.>, superconducting circuits <cit.>, pairs of Josephson phase qubits <cit.>, and neutral and ultra-cold atoms <cit.>.
Many studies have concentrated on few-particle scenarios; however, correlations also naturally emerge in quantum many-body systems <cit.>.

Atom assembly in optical lattices and optical tweezer arrays <cit.>, as well as trapped ions and molecules <cit.>, offer excellent approaches for studying many-body entanglement, with capabilities for single-atom preparation, control, and detection.
In the above-mentioned platforms, entanglement can be generated dynamically from the initial coherent state  well approximated by the one-axis twisting (OAT) model <cit.> where spin-squeezed and GHZ states are produced.
Given the control and local measurements possible in such architectures and their scalability, the systems offer the powerful setup to generate and study fundamental aspects of Bell correlations and non-locality, enabling the exploration of quantum information concepts. 
Imperfections, however, e.g. in system preparation, may introduce occupation defects, limiting Bell inequalities' application and entanglement certification.

In this paper, we study the role of occupation defects in detecting Bell correlations focusing on N two-level bosonic atoms in a spin-entangled state distributed among M sites of an optical lattice, such that N≤ M.
In the Bell scenario, we consider the measurement of two local collective spin observables at each lattice site where 0, 1 or 2 atoms are present. For each measurement result there are three possible outcomes: ±1/2 when the site is singly occupied, and 0 whenever a defect such as an empty or doubly occupied site is present, as illustrated in Fig. <ref>.
We introduce a simplified toy model where the atoms' internal and external degrees of freedom are decoupled. The dynamics of the internal spin degrees of freedom, including possible entanglement among the spins, is assumed to be given by the OAT model with N atoms. The external degrees of freedom are parameterized by the probability p of having a single occupied site allowing vacant, double, etc. occupancy sites when p<1. 
We derive the corresponding p-parameterized Bell inequalities based on M- and two-sites Bell correlations <cit.>. We analytically estimate the lowest (critical) value of the probability p for violation of Bell inequality.
We obtain p_c≈ 1/√(2) for the M-site Bell correlations. For the two-site Bell correlations, we obtain  p_c=4/5 and p_c=√(3)/2 when N=pM<M and N=M, respectively. These results are general and can be relevant for any platform where individual addressing of spins is
possible 
<cit.>, for example as in recent experiments using an array
of trapped ions <cit.> and Rydberg atoms <cit.> demonstrating generation of entanglement in terms of scalable spin-squeezing with tens of spins.





We explore the predictions of the toy model by performing full many-body numerical simulations of ultra-cold atoms in an optical lattice. We consider two distinct protocols for Bell correlations generation, which differ by the source of the occupation defects.

In the first one, Bell correlations between individual atoms are generated in the Mott insulating regime via a weak inhomogeneous magnetic field or interaction anisotropy <cit.>, and the considered occupation defects are vacant sites (holes). 
The probability p of single occupancy sites is the filling factor f=N/M. 
Our numerical simulations confirm predictions of the toy model of the lower bound on the filling factor concerning the detection of the M-sites Bell correlations.
However, the detection of the two-site Bell correlations in the lattice system is hindered whenever f<1 for fixed positions of holes. The movement of holes turned out to be crucial to lower the numerical bound on the filling factor, allowing for the detection of two-site Bell correlations up to f≈ p_c, as predicted by the toy model.

In the second protocol, entanglement is generated in the superfluid regime using atom-atom interactions <cit.> and next is transferred to the Mott phase by increasing the lattice height. The main source of defects is imposed by non-zero temperature T of an initial sample when N=M. We relate analytically the probability p with the initial temperature T of the system in the limit of an isoentropic transformation that leaves the system in thermal equilibrium at each moment. In this case, we determine analytically the upper bound for the initial temperature corresponding to p_c below which M- and two-site Bell correlations can be detected. 

The paper is organized as follows: In Secs. <ref> and <ref>, we provide the theoretical framework for the Bell scenario and toy model. This is applied to the system composed of ultra-cold atoms in an optical lattice presented in Sec. <ref> via analytical and numerical approaches. In Secs. <ref> and <ref> we present the results for the first and second protocols, respectively, discussing the role of defects introduced by holes and non-zero initial temperature of the system.
We conclude with further discussion and outlook in Sec. <ref>.








§ BELL SCENARIO



We consider M spatially separated subsystems labelled i=1,...,M, sharing a quantum M-partite state described by the density operator ρ̂. In each subsystem, k observables labeled α=0,...,k-1, can be measured, each having d possible outcomes. A Bell experiment, as shown in Fig. <ref>-(a), 
involves choosing an arbitrary observable for each of the M subsystems α={α_i}_i=1^M,
and recording the measurement results r={r^(i)}_i=1^M. The goal is to determine the M-partite probability distribution P_M(r|α) for the outcomes r given the settings α for all possible choices of settings α. In the LV theory <cit.>, any probability distribution P_M(r|α) can be written as

    P_M^( LV)(r|α)=∫ dλ q(λ)∏_i=1^MP^(i)(r^(i)|α_i,λ).

The local nature of this equation resides in the fact that the probability distribution of the outcomes for any given subsystem i depends only on the setting within that same subsystem i. The correlations between different subsystems here take their origin from a dependence relation that was established in the past, when the state ρ̂ was generated. This dependence can be fully described by some variable λ, which may be random with a probability distribution q(λ), affecting simultaneously all the M subsystems. Furthermore, if we consider a local realistic model where the measurement result r^(i) is deterministically determined by the setting α_i and the variable λ in each subsystem i, we have

    P^(i)(r^(i)|α_i,λ)=δ_r^(i),r^(i)(λ).


In the Bell experiment described above, the possibility of decomposing the measured probability distribution P_M(r|α) into the form (<ref>) constitutes the locality condition for correlations between the M subsystems under study. Conversely, if the measured probability distribution P_M(r|α) cannot be written in the form (<ref>),  the correlations present in the system are non-local. This decomposition corresponds to the most general locality condition, as defined by Bell. In practice, as we will see later, this condition is often formulated as a Bell inequality. 



In this paper, we consider an optical lattice with M sites, each containing zero, one or two two-level atoms. 
Each atom with two internal states a and b is an effective spin 1/2, and we can express local collective spins in the quantized form as
ŝ_x^(j) = (ŝ_+^(j)+ŝ_-^(j))/2, 
ŝ_y^(j) = (ŝ_+^(j)-ŝ_-^(j))/(2i), 
ŝ^(j)_z=(â_j^†â_j - b̂_j^†b̂_j)/2 
with ŝ_+^(j) = â_j^†b̂_j, 
ŝ_-^(j)=(ŝ_+^(j))^†.

Within the Bell scenario, we consider in each site j, the measurement of two local collective spin observables ŝ^(j)_α on the j-th site,
with α=0,1, each giving d=3 possible results

    r_0^(j)=0,-1/2,+1/2  ;   r_1^(j)=0,-1/2,+1/2,

as illustrated in Fig.<ref>-(a), where the result r=0 is assigned to any measurement result different from ±1/2. The probability distribution P_M(r|α) to obtain the results r given the settings α can be theoretically calculated, in terms of the density operator ρ̂ describing the system's state, as

    P_M(r|α)= tr[ρ̂⊗_j=1^MΠ̂_α_j,r^(j)]
    ,

where Π̂_α_j,r^(j)=±1/2 projects onto the eigensubspace of ŝ^(j)_α_j with eigenvalue r^(j)=±1/2, and Π̂_α_j,r^(j)=0=1̂-∑_r=±1/2Π̂_α_j,r projects onto the subspace perpendicular to both eigensubspaces of ŝ^(j)_α_j corresponding to the eigenvalues r^(j)=+1/2 and r^(j)=-1/2. 

In the next Section, we derive two Bell inequalities using this framework. The first one, relying on M-site correlations, is mainly useful for highly entangled states with a small number of lattice sites (see, e.g., Fig. <ref>-(b)). The second one, relying on two-site correlations, is more applicable to spin-squeezed states with a large number of lattice sites (see, e.g., Fig. <ref>-(c)). 
For both cases, we introduce a simplified toy model accounting for occupation defects that result in measurements outcomes r^(j)± 1/2.





§ TOY MODEL


In the lattice system considered one starts from a 
product state |x⟩^⊗ N where each atom is in a coherent superposition of two internal states, i.e. coherent spin state (CSS) along the x direction.
The entanglement between spins is dynamically obtained using many-body interactions by different protocols (see e.g., Sec. <ref> and Sec. <ref>) that can be effectively described, in some limit, by the OAT model

    Ĥ_ OAT = ħχŜ_z^2,

where χ quantifies the strength of interactions and Ŝ_σ represents the collective spin operator of N two-level atoms with σ = x,y,z. In (<ref>) we have σ =z. The OAT model generates spin-squeezed states as well as non-Gaussian entangled states, including the GHZ state <cit.>.


In order to study the role of occupation defects in the detection of Bell correlations, we introduce an approximate toy model where the internal and external degrees of freedom of the atoms are decoupled 

    ρ̂ = ρ̂_ext⊗ρ̂_SS.
 
Here ρ̂_ SS describes the internal degrees of freedom of the atoms including entanglement among the spins. In particular, we consider the OAT evolution (<ref>)

    ρ̂_SS=|ψ_t⟩⟨ψ_t|  with  |ψ_t⟩=e^-iĤ_ OATt/ħ|x⟩^⊗ N,

and ρ̂_ ext describes the external degrees of freedom, which we assume to be factorized over the different lattice sites 

    ρ̂_extern = ⊗_j=1^Mρ̂_j  with  ρ̂_j=p |1⟩_j _j⟨ 1| +(1-p) ρ̂^⊥_j,

where |1⟩_j is the Fock state with one spin at site j and p is the probability of having a single occupancy state. The probability of having an empty, double, etc.,  occupied site, as described by ρ̂^⊥_j, is given by (1-p). In practice, in ρ̂^⊥_j we only consider empty sites or doubly occupied sites. In this framework, when measuring a local spin observable, we assign the value r^(j)=0 to any measurement result different from ±1/2. 



 §.§ M-site Bell correlations



We now derive a Bell inequality whose violation allows the detection of non-local correlations. For this, from the measurement outcomes, in the two settings α=0,1, in each subsystem we introduce the complex quantity 

    r_+^(j) = r_0^(j)+i r_1^(j)

and we consider the product

    c_M = ∏_j=1^M r_+^(j),
 
In a local realistic theory, where the probability distribution P_M(r|α) has the form of (<ref>) and (<ref>), the average of c_M over many Bell experiment realizations is given by 

    ⟨ c_M⟩=∫ dλ q(λ)∏_j=1^M r_+^(j)(λ).

By introducing the functions f(λ)=1 and g(λ)=∏_j=1^M r_+^(j)(λ), and the scalar product ⟨ f,g⟩=∫ dλ  q(λ)f^*(λ)g(λ), the application of the Cauchy-Schwarz inequality to the functions f(λ) and g(λ), gives

    |⟨ f,g⟩|^2≤⟨ f,f⟩⟨ g,g⟩.

This leads to the following Bell inequality

    ℰ_M≡|⟨ c_M⟩|^2≤∫ dλ  q(λ)∏_j=1^M|r_+^(j)(λ)|^2≤ 2^-M,

and ℰ_M is the M-site Bell correlator.
We choose, in each site j, the settings α=0 and α=1 corresponding to ŝ_y^(j) and ŝ_z^(j) respectively <cit.>, to form 

    ŝ_+^(j)=ŝ_y^(j)+i ŝ_z^(j).

It is worth stressing here, that the above setting is optimal for the GHZ state and is less effective for other entangled states generated by the OAT dynamics.
The Bell inequality (<ref>) takes the form

    ℰ_M=|⟨ŝ_+^(1)ŝ_+^(2)...ŝ_+^(M)⟩|^2≤ 2^-M.




When considering the toy model (<ref>), the only nonzero contribution to E_M (<ref>) comes from the sites occupied by a single spin. We have

    ℰ_M^(p1) = p^2Mℰ_M^(p=1),

where p^M represents the probability that the M sites are occupied with a single atom. In the non-Gaussian regime of the OAT dynamics (<ref>), a macroscopic superposition of coherent states is created at some particular instants χτ=π/q labelled by an even integer q=2, 4, 6, …, M, with q=2 for the GHZ state <cit.>. The M-site Bell correlator corresponding to these non-Gaussian spin states, for p=1, reads ℰ_M^(p=1) 1/q^2, see Eq.(17) in <cit.>. By replacing in (<ref>), one can derive a critical value of p below which the M-site Bell correlations present in the states generated at χτ=π/q cannot be detected

    p_ c = q^1/M/√(2),

i.e. for which ℰ_M^(p=p_ c) = 2^-M.
We note that in the large M limit we have p_c≈ 1/√(2) for all q.





 §.§ Two-sites Bell correlations




In the case of the two-sites Bell correlations, we introduce the vector M⃗ and the matrix C̃, whose elements, for α,β=0,1, are respectively given by

    M_α   = ∑_j=1^M⟨ r^(j)_α⟩
    
    C_αβ   =∑_i,j≠ i⟨ r^(i)_αr^(j)_β⟩,   C̃_αβ = C_αβ-M_α M_β .

It can be shown <cit.> that for any input data M⃗ and C̃ compatible with a LV theory (<ref>), any 2× 2 positive semi-definite matrix A and any 2× 1 vector h⃗, the following Bell inequality holds

    L(A,h) = tr(AC̃)+h⃗·M⃗ + E_ max≥ 0,
 
where the classical limit is

    E_ max=max_r⃗[r⃗ ^TAr⃗-h⃗·r⃗],

with r⃗=(r_0,r_1)^T being the vector of all possible pairs of outcomes corresponding to the two settings α=0 and α=1 for a single subsystem. The positive semi-definite matrix A and the vector h⃗ that minimize L(A,h) can be found using a data-driven method as in Ref. <cit.>. We find that the Bell inequality
 [Alternatively, the Bell inequality (<ref>) can be represented using the C matrix, namely
L=C_00+C_11-C_01-C_10
-(M_0-M_1)^2-M_0-M_1 + M≥ 0.]

    L=C̃_00+C̃_11-C̃_01-C̃_10-M_0-M_1 + M≥ 0,

first established in Ref. <cit.> for spin-squeezed states in the case of only two possible measurement outcomes, is optimal also in our case with the three possible outcomes (<ref>) and under occupation defects.


For the case of two-sites Bell correlations, we choose α=0 and α=1 corresponding to the measurement, at each site j, of the spin components <cit.>

    ŝ^(j)_0    =  ŝ^(j)_n⃗cos(θ)+ŝ^(j)_m⃗sin(θ),
    ŝ^(j)_1    =  ŝ^(j)_n⃗cos(θ)-ŝ^(j)_m⃗sin(θ),

where the unit vector n⃗ is in the spin direction and m⃗, perpendicular to the spin direction, is in the best squeezing direction. The choice of the two settings (<ref>)-(<ref>) is dictated by the geometry of the spin-squeezed states. Non-zero expectation values come from the averages calculated in the plane spanned by the n⃗ and m⃗ vectors.

The averages (<ref>) and the correlations (<ref>) for α,β∈{0,1}, that form the Bell inequality (<ref>), take respectively the form

    M_α   = ∑_j=1^M∑_r=±1/2r  tr[ρ̂ Π̂_α_j,r] 
    C̃_αβ   =∑_i,j i^M(∑_r,s=±1/2rs  tr[ρ̂ Π̂_α_i,r⊗Π̂_β_j,s])-M_αM_β.
 


To account for the occupation defects, we evaluate M_α and C̃_αβ to determine the form of the Bell inequality using the toy model density matrix (<ref>) under the two scenarios investigated in the physical systems presented in Secs. <ref> and <ref> when N≤ M and N=M, respectively.






  §.§.§ Non-unit filling



The first scenario, relevant to the system described in Sec. <ref>, assumes that N=pM≤ M.
Under this condition, p<1 indicates a non-unit filling of the lattice whose sites are empty or singly occupied. The density matrix ρ̂^⊥_j in (<ref>), in this scenario, becomes

    ρ̂^⊥_j=|0⟩_j _j⟨ 0|.

The averages (<ref>) and the correlations (<ref>) for α,β∈{0,1} are in this case

    M_α   = ⟨Ŝ_n⃗⟩_ SScos(θ),
    C̃_αβ   =N-p/N-1{
    (ΔŜ_n⃗)^2_ SScos^2(θ)-(-1)^δ_αβ(ΔŜ_m⃗)^2_ SSsin^2(θ) . 
       . 
    -N/4cos[ 2 θ (1-δ_αβ) ]
    }
    +1-p/N-1⟨Ŝ_n⃗⟩^2cos^2(θ).
 

The Bell inequality (<ref>) is given, for this scenario, by

    L^ (V)_θ/M   =p(N-p)/N-1sin^2(θ)(4(ΔŜ_m⃗)^2_ SS/N-1)
         -2p⟨Ŝ_n⃗⟩_ SS/Ncos(θ)+1≥0.

The minimization of L^ (V)_θ with respect to θ gives

    L^ (V)/M   =1-p(N-1)/N-p⟨Ŝ_n⃗⟩^2_ SS/N[N-4(ΔŜ_m⃗)^2_ SS]
         -p(N-p)/N-1[1-4(ΔŜ_m⃗)^2_ SS/N],

for an optimal θ such that cosθ_ opt=N-1/N-p⟨Ŝ_n⃗⟩_ SS/N-4(ΔŜ_m⃗)^2_ SS[Similar to the previous case, this solution should verify |cos(θ_ opt)|≤1, which introduces the restriction p≤ N-(N-1)|⟨Ŝ_n⃗⟩/[N-4(ΔŜ_m⃗)^2]|. Here, L^ (V) can be written in terms of the squeezing parameter ξ^2=N (ΔŜ_m⃗)^2_ SS/⟨Ŝ_n⃗⟩_ SS^2  and the normalized mean spin v=⟨Ŝ_n⃗⟩_ SS/(N/2) as

    L^ (V)/M= 1-p(N-1)/N-pv^2/4/1-ξ^2 v^2-p(N-p)/N-1(1-ξ^2v^2).

]. 
The equation (<ref>) is represented in Fig. <ref> as a function of time of the OAT dynamics for a given N and for different values of the occupation probability p. This reveals that for a given N and at each moment of time τ, there is a critical value of p below which, the two-site Bell correlations present in the system cannot be detected. This critical probability p=p_ c is the value for which, the following 3^ rd order equation holds

    p^3-2N(1-4(ΔŜ_m⃗)^2_ SS/N)p^2
    +[N^2 + N-1 ] p 
       +[4N (ΔŜ_m⃗)^2_ SS+(N-1)^2⟨Ŝ_n⃗⟩^2_ SS/N[N-4(ΔŜ_m⃗)^2_ SS]]p 
       -N(N-1)=0.

In the limit of large atom number N, the minimal value of (<ref>), over χτ, converges to the p-dependent constant 

    L^ (V)_ min/M≈ 1-5/4p,

and the corresponding critical value of the occupation probability tends to p_ c=4/5.









  §.§.§ Doubly occupied sites


 

In this scenario, relevant to the system described in Sec. <ref>, we assume that N=M. Under this condition, p<1 indicates a redistribution of the N atoms among the M lattice sites resulting in the emergence of both empty and doubly occupied sites. The filling factor is one, f=1. Thus, the density matrix ρ̂^⊥_j in (<ref>) can be written as

    ρ̂^⊥_j=|0⟩_j _j⟨ 0|+|2⟩_j _j⟨ 2|/2.

The calculation of the averages (<ref>) and the correlations (<ref>) for α,β∈{0,1} gives

    M_α   = p ⟨Ŝ_n⃗⟩_ SScos(θ),
    C̃_αβ   =p^2 (ΔŜ_n⃗)^2_ SScos^2(θ)-(-1)^δ_αβp^2 (ΔŜ_m⃗)^2_ SSsin^2(θ)
         -p^2N/4cos[2θ(1-δ_αβ)],
 
where the subscript SS refers to an expectation calculated in the ρ̂_ SS state (<ref>), and the Bell inequality (<ref>) reads

    L_θ^ (VI)/M   =p^2sin^2(θ)(4(ΔŜ_m⃗)^2_ SS/N-1) 
       -2p⟨Ŝ_n⃗⟩_ SS/Ncos(θ)+1≥0.

By minimizing L^ (VI)_θ with respect to θ we obtain

    L^ (VI)/M=1-⟨Ŝ_n⃗⟩^2_ SS/N[N-4(ΔŜ_m⃗)^2_ SS]-p^2[1-4(ΔŜ_m⃗)^2_ SS/N],

for an optimal θ such that cosθ_ opt=1/p⟨Ŝ_n⃗⟩_ SS/N-4(ΔŜ_m⃗)^2_ SS[This solution should verify |cos(θ_ opt)|≤1 which introduces the restriction p≥ |⟨Ŝ_n⃗⟩/[N-4(ΔŜ_m⃗)^2]|.],[An alternative formula of L^ (VI) can be obtained by introducing the spin squeezing parameter ξ^2 and the normalized mean spin v as

    L^ (VI)/M= 1-v^2/4/1-ξ^2 v^2 - p^2(1-ξ^2v^2).

]. 
The equation (<ref>) is represented in Fig. <ref>-(a) as a function of time τ of the OAT dynamics for a given N and for different values of p. We note that for a given N and at each moment of time, there is a critical value of p below which, the two-site Bell correlations present in the system cannot be detected, from (<ref>) we obtain

    p_ c = √(N[N-4(ΔŜ_m⃗)^2_ SS]-⟨Ŝ_n⃗⟩^2_ SS)/N-4(ΔŜ_m⃗)^2_ SS.

In the limit of large atom number N, the minimal value of (<ref>), over χ t, approaches a p-dependent constant value

    L^ (VI)_ min/M≈3/4-p^2,

and the corresponding critical probability p_ c approaches p_ c=√(3)/2.  A lower value of p can be obtained with a Bell inequality including onsite and two-site correlations up to the fourth-order. In the limit of large N, this leads to a critical occupation probability of p_ c=1/2 <cit.>.











§ APPLICATION WITH ULTRA-COLD ATOMS IN OPTICAL LATTICES
 



We test the predictions of our toy model using a system consisting of N ultra-cold bosonic atoms confined in an optical lattice. We focus on rubidium-87 atoms occupying two internal states, labelled a and b, and loaded into an optical lattice potential, akin to recent experiments employing quantum gas microscopes <cit.>. The optical lattice comprises M lattice sites and is considered to be in one dimension.


The system, in the lowest energy band, is
conveniently considered in the Wannier functions basis <cit.>. In the tight-binding limit, where the lattice potential exceeds the recoil energy E_R=ħ^2k^2/(2m), the system can be conveniently described by the Bose-Hubbard model,

    ℋ̂_ BH   = - t ∑_i, j=i± 1(â_i^†â_j + b̂_i^†b̂_j) + U_aa/2∑_i=1^M n̂^a_i (n̂^a_i -1) 
        + U_bb/2∑_i=1^M n̂^b_i (n̂^b_i -1) 
    + U_ab∑_i=1^M n̂^a_in̂^b_i ,

where t and U_σσ' are the tunneling and interaction parameters.
â_i (b̂_i) is the annihilation operator of an atom in internal state a (b) in the i-th site of the lattice, and n̂^a_i=â_i^†â_i, n̂^b_i=b̂_i^†b̂_i are the corresponding number operators. 

We explore two protocols for the dynamical generation of Bell correlations within this system as detailed in Sections <ref> and <ref>, both resulting in Mott entangled states where atoms exhibit spin entanglement and are evenly distributed across the lattice. 
In both cases, the initial state ρ̂_a with all atoms in the internal state |a⟩ is turned into a
coherent superposition of |a⟩ and |b⟩ by applying a π/2-pulse, which is equivalent to 

    ρ̂_ ini = e^-i Ŝ_y π/2 ρ̂_a e^i Ŝ_y π/2  .

The two protocols differ in their specific mechanisms for inducing entanglement in the system and hence in the relevant source of occupation defects.

In the first protocol, Bell correlations are generated directly in the strongly interacting Mott regime U_σσ'≫ t where individual spins interact via spin-exchange interactions between neighbouring spins (<ref>). 
These interactions allow for the generation of entanglement within the system either through weak anisotropy or by coupling with an inhomogeneous field.
The mechanism for generating entanglement was explained and demonstrated in prior works such as <cit.>.
In the second protocol, Bell correlations are induced by atom-atom interactions in the superfluid regime t≫ U_σσ' with all-to-all individual spin connections. Subsequently, these correlations are transferred to the Mott phase by an adiabatic increase of the optical lattice depth, a process described in detail in previous studies <cit.>.

The final state of both protocols allows measuring spin components at a specific lattice site. In the Bell scenario, the outcomes of local measurements across all lattice sites can be collected, and used for the detection of Bell correlations using inequalities (<ref>) and (<ref>).
In the ideal realization of the Mott phase protocol, when ρ̂_a is the ground state where each lattice site hosts precisely one atom, two distinct measurement outcomes are possible.
However, imperfections introduce non-unit filling throughout the lattice, leading to additional measurement outcomes, such as when a site is either empty or double occupied. These failures arise from imperfect preparation of initial state ρ̂_a or to non-zero temperature.

In the following sections, we consider the role of imperfections and we establish a connection between the factor p of the toy model (<ref>) and the filling factor f=N/M or the initial temperature T for these two protocols.




§ ENTANGLEMENT GENERATION IN THE MOTT PHASE




In the strongly interacting limit U_σ, σ'≫ t and in the ground state manifold with at most one atom per lattice site, the system (<ref>) is approximately described by the t–J model 

    Ĥ_t-J=- t ∑_i, j=i± 1P̂_0
     (â_i^†â_j + b̂_i^†b̂_j)P̂_0
        + Ĥ_ XXZ ,

where P̂_0 is a projector operator over the manifold of at most single occupancy states, and where

    Ĥ_ XXZ=
        - J ∑_j=1^M-1( ŝ_x^(j)ŝ_x^(j+1) + ŝ_y^(j)ŝ_y^(j+1) + Δŝ_z^(j)ŝ_z^(j+1) - 1/4),

is the Heisenberg XXZ model with the spin-exchange interaction parameter J=4 t^2/U_ab and the anisotropy parameter Δ= U_ab/U_aa + U_ab/U_bb - 1 <cit.>. 
When Δ=1 the Hamiltonian takes the form of the isotropic Heisenberg XXX model and it is a natural case for rubidium-87 where U_aa≈ U_bb≈ U_ab. The anisotropy parameter Δ can be tuned by changing the values of interaction strengths using either Feshbach resonances or by shifting optical lattice potentials for states a and b.
The collective spin operators are just a summation over the individual ones,
Ŝ_σ = ∑_j=1^M ŝ^(j)_σ for σ = x,y,z, ±.
The tunnelling term in (<ref>) is relevant whenever the filling factor f=N/M is not one, meaning there are holes (empty sites) in the system and N≤ M. 

Let us consider the case of zero temperature when the initial state ρ̂_a is in the Mott regime. 
In the ideal case, we have ρ̂_a = ⊗_j=1^M |a⟩⟨ a|_j, indicating that at each lattice site there is an atom in the internal state |a⟩. However, in the presence of holes, the on-site state can be |0⟩⟨ 0|_j if it is not occupied. In the given experimental realisation, the number of holes can be arbitrary as well as their positions.
In our simulations, we evolve single realizations with initially each site being empty or singly occupied, and then average over many realizations to evaluate expectation values. Specifically,
 we generate a random number x_j∈ (0,1] for each lattice site, obtaining

    |Ψ_x⟩_j = θ̅(f-x_j) |a⟩_j 
        + ( 1-θ̅(f-x_j) ) |0⟩_j,

where θ̅(x) is the Heaviside step function.
The density matrix ρ̂_a describing N_r realizations is obtained through a simple averaging,

    ρ̂_a=
        1/N_r∑_r =1^N_r[ ⊗_j=1^M 
        |Ψ_x⟩⟨Ψ_x |_j
        ],

over the set of on-site random numbers r={x_j^1, x_j^2,  ⋯, x_j^N_r} for all j ∈{1,...M}.
One can show, that in the continuous limit, we obtain

    ρ̂_a = ⊗_j=1^M ( f |a⟩⟨ a|_j + (1-f) |0⟩⟨ 0|_j ),

given that all the sites are equivalent.
Therefore, the filling factor f can be identified with the parameter p of the external state of the toy model (<ref>) provided that ρ̂_j^⊥=|0⟩⟨ 0|_j. 
Owing to the presence of holes, the average number of atoms N differs from the number of lattice sites M when f 1, namely N=f M.

It is worth noting here, that the density matrix ρ_a in (<ref>) can also be cast in the following way 

    ρ̂_a    = 
        f^M 
        ⊗_j=1^M| a ⟩⟨ a |_j 
       +f^M ( 1-f/f)^1( 
        |0 ⟩⟨ 0 |_1 ⊗_j=1^M-1 | a ⟩⟨ a |_j + 𝒫_1) 
       +f^M ( 1-f/f)^2( 
        |0 ⟩⟨ 0 |_1 ⊗
        |0 ⟩⟨ 0 |_2 ⊗_j=1^M-2 | a ⟩⟨ a |_j + 𝒫_2) 
       + ⋯,

where 𝒫_σ are all other configurations with σ holes.
The state (<ref>) can be written using the shorter notation

    ρ̂_a = 
        f^M [
        ⊗_j=1^M | a ⟩⟨ a |_j + ∑_σ=1^M( 1-f/f)^σΣ_σ[  ρ̂_σ] ],

where ρ̂_σ is the product state of a given configuration having σ holes and 1-σ occupied sites while Σ_σ represents summation over permutations of all possible configurations. 


Having the above expression, the initial spin coherent state (<ref>) can be expressed as follows:

    ρ̂_ ini =
        f^M [ 
        ρ̂_ ini^(0)
        +
        ∑_σ=1^M( 1-f/f)^σΣ_σ[  ρ̂_ ini^(σ)] ],

where the rotations act on the density matrix corresponding to specific configurations of σ holes

    ρ̂_ ini^(σ) =    
         e^-i Ŝ_y π/2 ρ̂_σ e^i Ŝ_y π/2 ,

and ρ̂_ ini^(0)=e^-i Ŝ_y π/2 ( ⊗_j=1^M | a ⟩⟨ a |_j ) e^i Ŝ_y π/2.
Therefore, any expectation values of the operator Ô can be evaluated as 

    ⟨Ô⟩ = f^M [ 
         Tr[ρ̂_ ini^(0)Ô] + 
        ∑_σ=1^M 
        ( 1-f/f)^σ⟨Ô⟩_σ],

where ⟨Ô⟩_σ = Σ_σ Tr[ρ̂_ ini^(σ)Ô],
and where summation runs over all configurations of σ holes on M lattice sites. 

In Subsections <ref> and <ref>, we employ the above-discussed approach for the numerical evaluation of Bell correlations and the critical value of the filling factor allowing for their detection. 
Before, however, let us discuss the mechanism responsible for generating entanglement within the system in the Mott phase.




 §.§ Effective microscopic description



To generate entanglement from the initial spin coherent states (<ref>) driven by the t–J Hamiltonian (<ref>) we consider two methods.
Both of them were thoroughly investigated in <cit.>.

The first (A) uses a weak anisotropy Δ 1 in (<ref>), such that Δ≪ 3 - 2 cos(π/M) and Δ≫ 2cos(π/M) - 1, to generate entanglement in the system.
Under an ideal scenario with N=M,the system can be represented by a single spin chain whose dynamics is described by the OAT model <cit.>.
The presence of holes divides the spin chain into partial chains separated by holes whose configurations are included in (<ref>). For instance, if a single hole is located somewhere in the middle of the chain (not at the borders), the system can be effectively viewed as two partial chains. The system dynamics composed of the partial chains separate only when the positions of holes are fixed. 
The microscopic model describing the dynamics of each specific configuration involved in (<ref>) is effectively approximated by the OAT model when the positions of holes are pinned to the lattice. 
The model describing each of n partial chains for a given configuration of σ holes reads

    Ĥ^(0)_ eff, n = -χ_n^(0)Ŝ_z,n^2 ,     with   χ_n^(0)=J (Δ-1)/(L_n-1),

where L_n is the number of atoms consisting of the partial chain.


The second method (B) uses a weak inhomogeneous magnetic field 

    Ĥ_B=∑_j β_j ŝ^(j)_z,

with β_j ≪ J [cos (π/N ) - 1],
for the isotropic case Δ=1 to generate entanglement in the system. 
In this case, when σ holes are pinned into the lattice, the effective model describing each of the n partial chains is described effectively by the following Hamiltonian

    Ĥ_ eff, n = χ_ nŜ_z,n^2 + v_nŜ_z,n,

where we omitted constant energy terms, and where

    χ_n   = 1/L_n-1∑_q=1^L_n-1|c^(q)_n|^2/E^(q)_n,
        v_n    = 1/L_n∑_l=l_n^l_n+L_n-1β_l,
    
        c^(q)_n   = √(2)/L_n∑_l=l_n^l_n+L_n-1 p_l^(q,n) (β_l-v_n),

with l_n being the location of the first spin in the partial chain and E^(q)_n = J (1-cos(π q / L_n)). 

The effective models for the (A) and (B) methods are derived and explained in detail in <cit.>. Their validity for describing spin-squeezing generation was also demonstrated.

Therefore, one can employ the effective models based on the OAT dynamics for partial chains when evaluating expectation values as given by (<ref>) when considering fixed positions of holes (effectively t=0 in the t–J model (<ref>)). We employ them to calculate the upper bound on the filling factor required for the Bell correlations detection as demonstrated in <ref>.
However, the effective models become inaccurate when the tunnelling of holes occurs. In such a case, we perform full numerical many-body calculations by the exact diagonalization method of the t–J model (<ref>) when demonstrating Bell correlations and validity of the Toy model.

In the next subsections, we demonstrate the generation of Bell correlations in the system using scenarios A and B.
In our numerical simulations, we consider open boundary conditions <cit.> and use the parameters as in the recent experiment of A. Rubio-Abadal et al <cit.> with ^87Rb atoms, lattice spacing d=532 nm, and inter- and intraspecies interactions U_aa∼ U_bb∼ U_ab=U <cit.> with U=24.4t.





 §.§ M-site Bell correlations



We numerically evaluate the many-body Bell correlator 
ℰ_M = |⟨ŝ_1,+…ŝ_M,+⟩ |^2 as defined in (<ref>) for the initial spin coherent state given by (<ref>). 


In the upper panel of Fig. <ref> we show the evolution of the M-site Bell correlator ℰ_M with M=12 and for the two proposed methods, A and B. 
The small shift in time scale and lower magnitude of the correlator maxima in the anisotropic case result from its effective model being approximated to the first order <cit.>.
The only relevant contribution comes from the part of the system state describing all sites occupied, as already discussed in Section <ref>. The tunnelling of holes cannot change the value of this correlator.

In the lower panel of Fig. <ref> we check the scaling of the correlator with the filling factor f for relevant instances in time using a logarithmic scale. 
We see good agreement with the scaling of the toy model result in Eq. (<ref>) when identifying p with f.  
Likewise, we see the GHZ state result crossing the classical limit at the value obtained in Eq. (<ref>) for q=2, identified by the edge of the grey area. 








 §.§ Two-sites Bell correlations



In the short-time dynamics, where spin squeezing is generated, the two-site Bell correlator L given by (<ref>) is more relevant.

In Fig. <ref> we demonstrate the variation of the two-sites Bell correlations generated using the method A (interactions anisotropy, left column) and B (inhomogeneous magnetic field, right column) for different filling factors and effective tunnelling rates t/J.
The color bar indicates the magnitude of effective tunnelling when holes are present in the system relative to the relevant energy scale J. 
While in all instances t  0 to have a non-zero value of J, we assume for the effective tunnelling t/J=0 that positions of spins are fixed for those particular numerical simulations. 
Shaded areas indicate the upper and lower boundaries set by the results of holes fixed in place and holes moving infinitely fast, respectively. 
The upper bound is calculated numerically for the case of the fixed position of holes by using effective models as described in <cit.>.
The lower bound is given by the toy model and the corresponding Bell inequality (<ref>) where expectation values ⟨·⟩_ SS are replaced with the one given by the corresponding OAT model.
We also observe the generation of Bell correlations requires a larger filling factor in the inhomogeneous magnetic field scenario; both when the movement of holes is allowed and when not.
Nevertheless, the numerical results demonstrate the effectiveness of the proposed toy model  (<ref>) in the estimation of the lower bound in the detection of two-sites Bell correlations.



Note, the microscopic effective models (<ref>) and (<ref>) differ by the presence of a linear term, such that their results will be bounded by different limits as illustrated in Fig. <ref>. 
The critical value for the lower bound in this figure is obtained by solving the corresponding cubic equation by setting (<ref>) to zero, for which only one real root exists.
The anisotropy scenario A, while less accurately described by the corresponding effective OAT model (<ref>) 
 as explained in <cit.>, achieves lower bound results predicted by the toy model (<ref>) when tunnelling processes happen.
This is much harder to achieve in the inhomogeneous magnetic field scenario B since the linear term in (<ref>) changes at each configuration of holes and makes the dynamics largely incoherent.
This is a stark contrast with the M-site Bell correlator presented in Fig. <ref> where we find the lower bound set by the toy model (<ref>) to be equally accurate. 














§ ENTANGLEMENT GENERATION IN THE SUPERFLUID PHASE



We consider here the protocol where the quantum correlations are generated using atom-atom interactions within the superfluid phase and then stored in the Mott phase through adiabatic increasing of the lattice height <cit.>. In this protocol, initially, all atoms occupy the internal state a, and at zero temperature, they are in the ground state of the Bose-Hubbard Hamiltonian |ψ^(a)_0⟩ in the superfluid regime. Subsequently, a π/2 pulse is applied to put the atoms in a coherent spin state where the collective spin is along the x-axis. During the dynamics, the lattice height V_0 linearly increases as

    V_0(τ)=V_ ini+(V_ fin-V_ ini)τ/τ_ramp,

such that, at τ=τ_ramp, for an adiabatic evolution and initially zero temperature, the system reaches the Mott-squeezed state where the atoms are spatially distributed with one atom per lattice site and show entanglement in their internal degree of freedom (from spin squeezed state up to a GHZ state according to the chosen ramp duration τ_ramp). This corresponds to a state of the form (<ref>) with ρ̂_ ext=⊗_i=1^M|1⟩_j _j⟨ 1|.  At finite temperature T, the system exhibits initial thermal fluctuations described by the density operator

    ρ̂_T=1/Z∑_n e^-E_n^(0)/k_ BT|ψ_n^(a)⟩⟨ψ_n^(a)|,
 
where k_ B is the Boltzmann constant, E_n^(0) are the eigenenergies of the initial Hamiltonian (in the superfluid phase), all atoms are in the internal mode a so |ψ_n^(a)⟩ are the corresponding eigenstates, and Z=∑_ne^-E_n^(0)/k_ BT is the normalization constant. These initial thermal fluctuations lead, at the end of the ramp, to a non-zero probability of having holes and double occupations of the lattice sites. In this section, we explore the effect of these initial thermal fluctuations on the detection of the M-site and two-sites non-local Bell correlations present in the final state. We first note that numerical observations performed in small 1D lattices reveal that when U_aa=U_bb≳ U_ab> 0, which is relevant to our purposes, the external dynamics is weakly affected by the internal dynamics up to the Mott transition. Consequently, we estimate the occupation statistics of different lattice sites in the final state by restricting the Bose-Hubbard model to a single internal state a. In this case, the spectrum of the Hamiltonian in the deep Mott phase with t→ 0, is simple: a non-degenerate ground state |ψ_0^ MI⟩, with one atom per lattice site, and gapped excited states |ψ_n^ MI⟩ showing holes and doubly occupied sites. The ground state |ψ_0^ MI⟩ is obtained, from the initial ground superfluid state |ψ_0^(a)⟩, by a unitary evolution with the time dependent Bose-Hubbard Hamiltonian

    |ψ_0^ MI⟩=Û|ψ_0^(a)⟩.




 §.§ M-site Bell correlations



Since the excited states |ψ_n^ MI⟩ have at least one hole and one doubly occupied lattice site, they give no contribution to the M-site Bell correlator (<ref>) whose value at non-zero temperature is given by

    ℰ_M^(T0)=P_0^2ℰ_M^(T=0),

where P_0=e^-E_0^(0)/k_ BT/Z represents the probability of occupying the ground state in (<ref>)[It is important to note that in general, the final external configuration of the system does not consist of a factorized state over lattice sites, as in (<ref>), in particular the probability P_0 is not of the form P_0=p^M.]. 
Fig. <ref>-(a) illustrates this relationship for N=M=4 over varying temperatures, revealing the existence of a critical temperature T_ c below which M-site Bell correlations are detectable. In the weakly interacting regime, the probability P_0 can be analytically determined using the Bogoliuobov theory, where, in the initial superfluid regime, when all the atoms occupy the internal state a, the system can be approximately described by the Hamiltonian

    Ĥ_ Bog=E_0+∑_j≠ 0ħω_jd̂_j^†d̂_j.

In equation (<ref>), E_0 is the ground state energy of the system, d̂_j is the annihilation operator of a Bogoliubov quasi-particle associated with the quasi-momentum q_j, and ħω_j is given in terms of the tunneling parameter t and the interaction parameter U=U_aa, by

    ħω_j=4t√(sin^2(π/Nj)(sin^2(π/Nj)+U/2tN/M)).

At non-zero temperature T, the probability P_0 of being in the ground state is given by 

    P_0=∏_j≠0p_j(n_j=0),

where, p_j(n_j=n)=e^-nħω_j/k_ BT/Z_j with Z_j=(1-e^-ħω_j/k_ BT)^-1. By replacing in (<ref>), one obtains

    P_0=∏_j≠ 0(1-e^-ħω_j/k_ BT).

The critical temperature T_ c, below which Bell correlations in states generated at times χτ=π/q of OAT dynamics are detectable, can be determined by setting P_0 (<ref>) equal to the critical probability p_ c^M, where p_ c is given by (<ref>). Fig. <ref>-(b), illustrates this critical temperature for a GHZ state (i.e. q=2) as a function of the atom number N.

 



 §.§ Two-sites Bell correlations



In large systems, studying M-site Bell correlations can be challenging both theoretically and experimentally. Therefore, for such systems, we focus on two-site Bell correlations. In this context, we assume that the system is brought to the Mott phase through a transformation that leaves the system in thermal equilibrium at each moment and conserves the entropy. Under such conditions, one can explicitly write an approximate external density matrix ρ̂_ ext of the system in the final Mott phase and thus analytically determine the probability p(T_ i) of having exactly one atom in a given lattice site as a function of the initial temperature T_ i. We first calculate the initial entropy S_ SF as a function of T_ i. We then determine the final temperature T_ f of the system in Mott using entropy conservation <cit.>

    S_ SF(T_ i)|_V_0=V_ ini=S_ Mott(T_ f)|_V_0=V_ fin.
 
In the weakly interacting regime, in the superfluid phase, the entropy can be calculated using the Bogoliubov theory. Indeed, By using the Bogoliubov spectrum (<ref>), the system partition function at temperature T can be written as (non interacting bosons)

    Z_ Bog(T)=∏_j(1-e^-ħω_j/k_ BT)^-1.

The free energy defined as F_ Bog=-k_ BTln Z_ Bog is given by

    F_ Bog(T)=k_ BT∑_jln(1-e^-ħω_j/k_ BT).

We thus deduce the entropy of the system as

    S_ SF^( Bog)(T)/k_ B≡-∂_T F/k_ B
       =-∑_j[ln(1-e^-ħω_j/k_ BT)+ħω_j/k_ BTe^-ħω_j/k_ BT/1-e^-ħω_j/k_ BT].




We now calculate the entropy of the system in the Mott phase by considering the limit t=0 and small temperatures. In this case, the partition function, when N=M, can be estimated using a two particles-holes excitation approximation, where, we take into account only the states with at most two holes and two doubly occupied lattice sites, with U=U_ aa

    Z_ Mott(T)   =1+M(M-1)e^-U/k_ BT
         +M(M-1)(M-2)(M-3)e^-2U/k_ BT.
 
By using the free energy, in the Mott phase, one can obtain the system entropy

    S_ Mott(T)/k_ B   =ln Z_ Mott+1/Z_ MottU/k_ BT[M(M-1)e^-U/k_ BT.
         .+M(M-1)(M-2)(M-3)e^-2U/k_ BT].
 
After determining T_ f using (<ref>), we approximate the external density matrix of the system at the end of the ramp as 

    ρ̂_ ext   ≈1/Z_ Mott(T_ f)[|ψ_0^0 h⟩⟨ψ_0^0 h|
    . 
       .
    +e^-U/k_ BT_ f∑_j=1^M∑_k≠ j^M|ψ_jk^1 h⟩⟨ψ_jk^1 h|.
         .+e^-2U/k_ BT_ f∑_j=1^M∑_k≠ j^M∑_l≠ j,k^M∑_m≠ j,k,l^M|ψ_jklm^2 h⟩⟨ψ_jklm^2 h|],
 
where |ψ_0^0 h⟩ is the state with exactly one atom per lattice site (zero holes), |ψ_jk^1 h⟩ is the state with only one hole in the j^ th site and a single double occupancy in the k^ th site and |ψ_jklm^2 h⟩ is the state with two holes in the j^ th and k^ th sites respectively and two double occupancies in the l^ th and m^ th sites respectively. This approximation enables us to analytically calculate, as a function of the initial temperature, the probability p of having a Fock state with one atom in a given lattice site

    p   ≡ tr[|1⟩_i_i⟨1| ρ̂_ ext]
       =1/Z_ Mott(T_ f)[1+(M-1)(M-2)e^-U/k_ BT_ f.
         .+(M-1)(M-2)(M-3)(M-4)e^-2U/k_ BT_ f]
 

 

In numerical simulations of the exact dynamics, p can be calculated as the trace of the projection of the density matrix of the system on the eigenspace of the observable n̂_i associated with the eigenvalue n_i=1. Equation (<ref>) is represented in Fig. <ref>-(a) as a function of the initial temperature, compared to exact results using both one- and two-component Bose-Hubbard Hamiltonian. This reveals that at sufficiently low temperatures in small systems, the probability p(T_ i) can be accurately approximated using (<ref>) and the entropy conservation condition (<ref>). By equating p(T_ i) with the critical probability (<ref>) from the toy model, one can determine an upper bound, dependent on the ramp duration τ_ramp, for the initial temperature T_ c above which two-sites non-local Bell correlations cannot be detected. Fig. <ref>-(b) shows the critical temperature associated to the minimal (critical) p, over the ramp duration τ_ramp, as a function of the atom number N 
[For larger atoms number, beyond the range of N presented in the figure, the two-hole approximation becomes invalid.].









§ SUMMARY AND CONCLUSIONS




In this paper, we considered the detection of Bell correlations using two-level ultra-cold bosonic atoms loaded into optical lattices. 
We focus on identifying Bell violations for dynamically generated entangled states when imperfections related to non-unit filling appear.

Our proposed toy model accounts for imperfections that provide non-unit filling per site in the preparation and measurement stages through the probability p of a site being single-occupied. The M-site Bell correlator shows a violation for this model in the limit of a large system, when p>p_c≈ 1/√(2). 
On the other hand, the two-site Bell correlator relying on collective spin measurements shows the bounds p_c=4/5 and p_c=√(3)/2 when N=pM<M and N=M, respectively. 
We illustrate these general results under two practical methods of entanglement generation varied by the source of imperfections: vacant sites due to initial non unit filling, and vacant and multiply occupied sites due to non-zero temperature of the initial state. 
In the presence of holes, we study the generation of entanglement in the Mott insulating via interaction anisotropy and inhomogeneous field demonstrating the validity of the toy model predictions for the critical value of the filling factor allowing generation of Bell correlations.
For non-zero temperatures, we study the generation of entanglement in the superfluid regime produced via atom-atom interactions transferred to the Mott regime via adiabatic rising of the lattice height. 
We find a connection between the probability p of the toy model and the effective temperature T of the initial state that drives the system from the unit filling regime under study. 
We identified the critical value of the initial temperature allowing violation of Bell inequalities.

Our results reveal the fundamental limits on detecting Bell correlations in the lattice system due to occupation defects. The bounds identified in specific protocols fall within the range of typical experimental realizations, suggesting the feasibility of detecting Bell correlations.





§ ACKNOWLEDGMENTS


We gratefully acknowledge discussions with Irénée Frérot. This work was supported by the Polish National Science Center DEC-2019/35/O/ST2/01873 and DEC-2023/48/Q/ST2/00087. 
T.H.Y acknowledges support from the Polish National Agency for Academic Exchange through the Foreign Doctoral Internship Grant NAWA Preludium BIS 1 No. PPN/STA/2021/1/00080. 
A part of the computations was done at the Centre of Informatics Tricity Academic Supercomputer & Network.







