Detecting two-site spin-entanglement in many-body systems with local particle-number fluctuations

We derive experimentally measurable lower bounds for the two-site entanglement of the spin-degrees of freedom of many-body systems with local particle-number fluctuations. Our method aims at enabling the spatially resolved detection of spin-entanglement in Hubbard systems using high-resolution imaging in optical lattices. A possible application is the observation of entanglement generation and spreading during spin impurity dynamics, for which we provide numerical simulations. More generally, the scheme can simplify the entanglement detection in ion chains, Rydberg atoms, or similar atomic systems.

We derive experimentally measurable lower bounds for the two-site entanglement of the spindegrees of freedom of many-body systems with local particle-number fluctuations. Our method aims at enabling the spatially resolved detection of spin-entanglement in Hubbard systems using highresolution imaging in optical lattices. A possible application is the observation of entanglement generation and spreading during spin impurity dynamics, for which we provide numerical simulations. More generally, the scheme can simplify the entanglement detection in ion chains, Rydberg atoms, or similar atomic systems.

I. INTRODUCTION
The role of entanglement for the quantitative understanding of quantum many-body systems has been the topic of a large number of theoretical studies [1][2][3]. In contrast, the experimental detection of entanglement in quantum many-body systems is less developed, which currently hinders the establishment of more direct links between experiments and theory. So far, entanglement witnesses have been extracted from macroscopic properties or diffractive probes, such as magnetic susceptibilities [1,4,5], spin-or atom-number squeezing parameters [6][7][8][9], or time-of-flight imaging [10,11]. Further, experiments using controlled collisions in optical lattices indicated the generation of entangled cluster states [12,13]. However, these experiments did not access the spatial dependence of entanglement measures which is crucial for observing some of the elementary properties of entanglement in many-body systems, such as area laws [3] or the dynamical generation and spreading of entanglement [1,14,15]. A candidate for establishing a direct experiment-theory connection are quantum spin systems [1][2][3]16]. Such spin Hamiltonians can effectively describe the low-energy physics of certain materials [17][18][19][20], for which a local detection of entanglement seems challenging. However, recent atomic physics realizations of quantum spin systems, such as neutral atoms in optical lattices [21][22][23][24][25][26][27][28] and trapped ions [29][30][31][32][33][34][35], offer the possibility of a local read-out of spin correlations. In ion traps, such local detection of spin correlations and entanglement has been the standard for many years but was mostly used in the context of quantum computing [36]. Only recently, these techniques were employed to detect entanglement in a simulation of a spin system showing the first spatially-resolved detection of entanglement spreading after a local quantum quench [34]. For quantum many-body systems in optical lattices, lo- * Electronic address: manuel.endres@mpq.mpg.de cal detection of individual particles and their correlations has only been demonstrated in the past few years using high-resolution microscopy [24][25][26][37][38][39][40]. Proposals have been made to detect the Rényi entropy with this technique [15,[41][42][43] but no experiment has shown the spatially resolved detection of entanglement in such systems to date.
A key difference between quantum magnetism experiments in ion traps and optical lattices is that in the latter, on-site number fluctuations coexist with spin fluctuations. The reason is that spin interactions in optical lattices are typically generated via superexchange as a second order process in the large interaction limit of Hubbard models [21][22][23], where number fluctuations are suppressed but not absent. Particularly in low dimensions, local number fluctuations can be sizable even at zero temperature [38,40], and additionally the currently achievable temperatures lead to thermal activation of defects [39]. In solids, number fluctuations naturally arise through hole-doping of Mott insulators, leading to an effective description in terms of t-J models [18]. For such systems, a detection of spin-entanglement must take the presence of occupation number fluctuations into account. On the theoretical level, the distinction between entanglement in internal and number degrees of freedom has been clarified in Refs. [44][45][46]. However, concrete experimental proposals for detecting the entanglement between spins in quantum many-body systems of atoms with local particle-number fluctuations are lacking. Here we propose an experimentally feasible scheme to detect spin-entanglement between two sites in the presence of number fluctuations in Hubbard systems using singleatom-and single-site-resolved imaging of atoms in optical lattices [37][38][39]. To this end, the key challenges are current limitations in these setups, namely the lack of arbitrary local spin rotations [47], the lack of full spin resolution [25], and the parity-projection problem [38,39]. Fully accounting for these restrictions, we derive detectable lower bounds for the concurrence [48], an entanglement measure, of the spin-degree of subsystems consisting of two lattice sites. Our method can be readily implemented in current high-resolution imaging setups for optical lattices without technical modifications [24,25,[37][38][39][40]. The scheme is immediately applicable to studying the entanglement generation and spreading during single spinimpurity dynamics in one-dimensional Bose-Hubbard chains [25,49,50]. For this scenario, we provide numerical simulations identifying a parameter range where such experiments could be performed. While our focus is on spin-impurity dynamics, the method can be used in a broader context. For example, it could be an important diagnostic tool in the current experimental search for antiferromagnetic order in the fermionic Hubbard model realized with cold gases [51]. For ion trap implementations of quantum magnetism, the bounds derived in Sec. III could lead to a simplified detection of entanglement in impurity dynamics [34] or global quantum quenches [35] without the need for a full state reconstruction. Further, our results also apply to experiments with Rydberg atoms in optical tweezers [52][53][54][55][56], where atom number fluctuations can result from trap loss. Finally, our method could be used to detect the entanglement in spatially ordered structures of Rydberg excitations in optical lattices [57]. The outline of the paper is as follows. In Sec. II, we give an introduction to entanglement generation and spreading during single impurity dynamics in the nearest neighbor spin-1/2 XX-chain [49,50]. The derivation of lower bounds for the concurrence then follows in several steps taking into account the known experimental limitations for high-resolution imaging of quantum gases in optical lattices. In Sec. III, we derive a lower bound neglecting number fluctuations based only on global pulses in order to circumvent the lack of arbitrary local spin rotations [47]. We then give a conceptual introduction to the detection of spin-entanglement in the presence of number fluctuations in Sec. IV, followed by a case study of spin impurity dynamics in the one-dimensional Bose-Hubbard model in Sec. V. We extend the detection scheme to include number fluctuations in Sec. VI assuming fully spin-resolved detection. In Sec. VII, we account for the current inability to detect two different spin states at once [25] and also treat the restriction to local parity imaging [38,39]. We finish with a conclusion and outlook section.

II. ENTANGLEMENT DURING IMPURITY DYNAMICS IN THE XX-CHAIN
To provide a concrete example and target application, we review the entanglement generation and spreading during spin impurity dynamics in a spin-1/2 XXchain [49,50] with Hamiltonian where J ex is the exchange coupling andŜ ± j = 1 2 (σ x j ±iσ y j ) are spin-1/2 raising (lowering) operators. Withσ α j (α = x, y, z) we denote the Pauli operators applied to site j. Hamiltonians of this type are important for describing recent experiments realizing spin-impurity dynamics in one-dimensional Bose-Hubbard systems [25] and ion chains [34]. In the case of Hubbard systems, the spindescription is precise only for a single spin impurity in the deep Mott insulating limit at zero temperature, where on-site number fluctuations are strongly suppressed. We will come back to this point in more detail in Sec. V and first neglect on-site number fluctuations. For the ion trap implementation, the correct description would be a longrange XX-model instead of the nearest-neighbor Hamiltonian (1). Nonetheless, the following discussion still applies to this case with a simple substitution as detailed below. In the following, we will write a state with a single upspin impurity on site j as where L is the total number of sites, and | ↑ (| ↓ ) refers to up-spin (down-spin) states in the z-basis. As an initial state, we choose a single up-spin impurity at the center of the chain |ψ 0 = |j = 0 . For an infinite chain (L → ∞), the time-evolution under Hamiltonian (1) leads to a spreading of this impurity according to with φ j = i j J j (J ex t/ ), where J j (x) is the Bessel function of the first kind, t is the evolution time, and is the reduced Planck constant. For the long-range XX-model, which is relevant for ion chains, φ j must be substituted by a different function that can be calculated numerically [34]. For the experimental observation in a Hubbard system [25], the probability of finding the spin impurity on site j after various evolution times was observed in quantitative agreement with Eq. (2). However, this experiment did not quantify the correlations and entanglement between spins on different sites A and B. This information is encoded in the two-site reduced density operator , where the trace runs over all sites but A and B. The superscript s stands for single spin-impurity. We find writing the two-site density matrix using basis states | ↑, ↑ , | ↑, ↓ , | ↓, ↑ , | ↓, ↓ for the A and B sites. For any state with a single impurity in an otherwise polarized background, the reduced two-site density matrix has the structural form ofρ s A,B (t). The entanglement between sites A and B can be quantified with the concurrence C [48], a commonly used bipartite entanglement measure [58,59]. The concurrence for a general bipartite pure state |ψ 1,2 in a tensor product H 1 ⊗ H 2 of two finite-dimensional Hilbert spaces H 1 , H 2 can be defined as [60,61] C(|ψ 1,2 ) = 2( ψ 1,2 |ψ 1,2 − Tr(ρ 2 1 )), whereρ 1 = Tr 2 (|ψ 1,2 ψ 1,2 |) is the reduced density operator of subsystem 1. The concurrence defined in this way can also be applied to subnormalized states. The concurrence C(ρ 1,2 ) of a bipartite mixed stateρ 1,2 is defined via a convex roof construction [61] using the infimum over all decompositions ofρ 1,2 into pure states |φ i : Even if the global state |ψ 0 (t) is pure, the reduced density operatorρ s A,B (t) is mixed. We are therefore dealing with a mixed bipartite two spin-1/2 system. Due to the X-matrix form ofρ s A,B (t), the concurrence can be easily calculated [62] (see Eq. (7)): a result obtained earlier in Refs. [49,50].
To get a better intuition for this outcome, we can restrict ourselves to sites with A = −B. In this case, the two-site density matrix can be written as a mixture of a Bell-state |Ψ + = 1 √ 2 (| ↑, ↓ + | ↓, ↑ ) and | ↓, ↓ : Therefore, the concurrence amounts to the probability of finding the system in the Bell state. We show C(ρ s A,B (t)) for various times and sites A and B in Fig. 1a and b, which illustrates how entanglement is generated and spreads in a wave-like fashion during the impurity dynamics.
III. SCHEME FOR SPIN-1/2 SYSTEMS Experimentally, we are facing the problem of detecting the concurrence of an unknown two-site density matrix ρ A,B that might be close to but not necessarily equal tô ρ s A,B due to experimental imperfections. Detecting the concurrence of an unknown state is possible using a full state tomography. For two spin-1/2 systems, a full state tomography can be achieved by measuring all nine combinations of Pauli operators σ α Aσ β B with α, β = x, y and z [63]. We assume that the final measurement is always Figure 1: a Density plot of the concurrence C(ρ s A,B (t)) for the single spin impurity dynamics as a function of lattice sites A and B for tJex/ = 0.2, 1, 3 (left, middle, right). b Concurrence C(ρ s A,−A (t)) for the single spin impurity dynamics evaluated at sites A, −A for tJex/ = 3, 4, 5 (open circles, filled circles, open rectangles). Lines are shown as a guide for the eye. c Density plot of the lower bound for the concurrence K(ρ s A,B (t)) (see Eq. (9)) for the single spin impurity dynamics as a function of lattice sites A and B for tJex/ = 3. Note that K(ρ s A,B (t)) = 0 for odd distances A − B, which results in a checkerboard pattern.
performed in the z-basis, for example, by reading out the populations of two atomic energy levels that encode the spin-1/2 system. A measurement in a different basis is possible by applying pulses that rotate the individual spins before the measurement. A pulse on a single spin on site j can be represented with a unitary operator written in the | ↑ , | ↓ basis. For example, a measurement in the x basis can be realized by a θ = π/2, φ = −π/2 pulse becauseσ x j =R (π/2, −π/2) † jσ z jR (π/2, −π/2) j . For the following discussion, it is important to distinguish pulses on individual spins, which allow for the measurement of σ α Aσ β B for all combinations α, β = x, y and z, and global pulses on both spins, which restrict the measurements to equal axes α = β. Pulses on individual spins arranged in a chain are commonly employed in ion trap implementations. For example, in Ref. [34], the authors show the detection of the concurrence generated during spin impurity dynamics in a long-range XX-model using a full state tomography. However, for optical lattice implementations of spinsystems using Hubbard models, only θ = π pulses on individual atoms have been demonstrated using a rapid adiabatic passage [47]. Pulses on individual atoms with arbitrary θ, φ require improved experimental control and are yet to be implemented. This currently restricts the detection to elements σ α Aσ β B with α = β. Therefore, we now present a simplified scheme for the detection of a lower bound for the concurrence using only global pulses.

A. Bound for global pulses with controlled φ
The first step in deriving the bound is to split the unknown two-site density matrix into an X-and O-part according toρ Knowledge of only the X-part is sufficient to detect a lower bound for the concurrence because [64,65]: The concurrence for density matrices in X-form is given by [62] C Because Thus, we find the lower bound which only requires global pulses for the detection of σ x Aσ x B and σ y Aσ y B . The probabilities for having both spins up, P ↑,↑ , and both spins down, P ↓,↓ , can be detected in the z-basis without pulse before the measurement.

B. Bound for global pulses with undetermined phase φ
The phase φ of the applied pulse is difficult to control experimentally. For the case of the impurity dynamics detailed above, controlling the phase φ would require having a defined phase of the applied field for the pulse relative to the starting time of the dynamics. This is difficult to reach for the implementation in Ref. [25] because the spin dynamics occurs in the tens of hertz regime, while the applied pulses are in the gigahertz regime. We will assume that the pulses are not phase-locked to the starting point of the dynamics. In this case, φ is essentially random. All observables after a global pulse with θ are then effectively described by an equal statistical mixture over all angles φ described by a density matrix is the two-site density matrix after a global pulse with angles θ and φ. Let us denote the average value ofσ z Aσ z B after a global pulse with θ = π/2 and random φ by σ z Aσ z B π/2 . Then, we have where we used the invariance of the trace under cyclic permutation in the second line. Therefore, a measurement after a global π/2 pulse with random φ corresponds to a measurement of the mean of σ x Aσ x B and σ y Aσ y B . The bound (8) can then be rewritten as Importantly, a detection of K A,B only requires measurements with and without a global π/2 pulse (with random phase φ), which simplifies the experimental effort dramatically as compared to a full state reconstruction.

C. Quality of the bound for the case of single-spin dynamics
An important question is how tight the bound (9) is for the case of single-impurity dynamics detailed in Sec. II. Using Eq. (3), we find that The reason for this behavior is that spins at even distances have a parallel alignment in the x − y plane in the sense that | σ x ) is plotted as a function of A and B (Fig. 1c). While the fact that K(ρ s A,B ) = 0 for odd distances is a disadvantage on first glance, this checkerboard pattern can serve as an experimental signature on top of noisy experimental data. Without going into details, we note that by applying a magnetic field gradient before the detection, the offdiagonal element ρ ↑↓ acquires a time-dependent complex phase-factor. Tuning this phase to π/2 changes the parallel alignment of the spins into perpendicular alignment and vice versa. As a result, the measured bound would be tight for odd distances and zero for even distances. Using this technique, a tight bound can be achieved for all pairs of spins.

IV. SPIN-ENTANGLEMENT IN THE PRESENCE OF ATOM NUMBER FLUCTUATIONS
Quantum magnetism experiments in optical lattices are typically performed using mixtures of atoms in two different hyperfine states [21][22][23][25][26][27][28]. The local onsite states can be written as |n + l , n − l , where n + l and n − l is the number of atoms in the two hyperfine states on site l. The state of the whole system can be expanded in basis states l |n + l , n − l . We also introduce a notation for the total atom number on site l as n l = n + l + n − l . The connection to spin systems is obtained using the Schwinger representation (see, e.g., Ref. [66]), which maps the onsite states to a total spin j l system with spin-projection m l defined as We will also use the notation |j l , m l = |n + l , n − l . In the large-interaction limit of Hubbard models, the dynamics in subsectors with fixed j l = 1/2 is governed by XXZ-models [21,22]. However, due to the finite temperature of the samples [38,40] and quantum fluctuations [39], number fluctuations are introduced into the system. This results in contributions of on-site states that map to different j l = 1/2. We are facing a situation where both spin fluctuations (i.e., fluctuations of m l for a fixed j l ) and number fluctuations (i.e., fluctuations of j l ) are present in the system. It is both experimentally and conceptually interesting to ask whether entanglement between the spin-projection degree of freedom is detectable in this scenario. Again, we consider a subsystem consisting of two sites A and B, for which the reduced density operator now also includes contributions from different occupation numbers: where we used the Schwinger and occupation number notation in the second and last line respectively.

A. Entanglement of particles
First, we are dealing with the question of how to conceptually differentiate the entanglement in the spin degree of freedom from entanglement that stems from different total local occupation numbers [44][45][46]. For example, superpositions of states with different local atom numbers, such as 1 √ 2 (|1, 0, 0, 0 + |0, 0, 1, 0 ) (corresponding to a single plus atom in a superposition between site A and B), should not appear entangled. An appropriate procedure to achieve this goal is to first project onto states with fixed local atom numbers. To this end, we define projected two site operatorŝ whereΠ n l l = m l |j l = n l /2, m l j l = n l /2, m l | is the projection operator at site l onto local total atom number n l = n + l + n − l , or, in the Schwinger notation, onto local total spin j l = n l /2. The entanglement in the spin degree of freedom can then be captured by the so-called entanglement of particles [44,45] where we used the concurrence C as an entanglement measure, and is the probability of finding the system with n A atoms on A and n B atoms on B.
In the first line of (11), the concurrence is evaluated with the normalized stateρ n A ,n B A,B /p n A ,n B . The second line follows from the definition (5) of the concurrence applied to the subnormalized operatorρ n A ,n B A,B . A trivial lower bound for E p (ρ A,B ) is The projected operatorρ 1,1 A,B describes the subsector with unity filling on both sites, that is, with local total spin j l = 1/2 on both sites. We will refer to this as the spin-1/2 sector in the following. Our goal is to formulate a detectable lower bound for the entanglement contained in the spin-1/2 sector quantified by the concurrence C(ρ 1,1 A,B ), which can eventually be used to bound the entanglement of particles via the previous inequality.
B. Simplified spin-1/2 notation For the following sections, we will introduce a shorthand notation for the density matrix elements in the spin-1/2 sector based on the Schwinger notation: Instead of the cumbersome notation with ± 1 2 , we will use ↑ and ↓ for spin up and spin down. For example, To illustrate these concepts and to investigate the influence of number fluctuations, we carried out a case study by numerically simulating the dynamics of a mobile spin impurity in the one-dimensional two-species Bose-Hubbard model: Hereb ( †) σ,j is the operator that annihilates (creates) a boson of species σ = {+, −} at site j, J is the hopping amplitude and U the interaction strength. Note that for simplicity the inter-and intra-species interaction parameters are taken to be equal, although in usual alkaline gases they assume slightly different values. In the limit U J at N + + N − = L (N ± is the total atom number of the respective species), the system is in a Mott phase with one particle per site, where charge degrees of freedom are frozen, but internal ones are not. They can be described with an XXZ Hamiltonian via second-order perturbation theory [21,22]: where J ex = 4J 2 /U . The local states | ↑ and | ↓ upon which the spin-1/2 operatorsŜ ± j act are identified with the states |n + = 0, n − = 1 and |n + = 1, n − = 0 , respectively, using the Schwinger representation (see Sec. IV). For the case of a single spin impurity in an otherwise polarized chain, the last term ofĤ XXZ is only a constant offset, and the dynamics are described by the XX Hamil-tonianĤ XX as discussed in Sec. II. Due to on-site number fluctuations, this mapping can break down in experimentally relevant parameter ranges. We consider two possibilities in the following. First, for stronger hopping J, significant quantum fluctuations of the on-site particle number are introduced in the form of correlated particle-hole pairs [40] even at zero temperature. One of the open questions here is up to which dimensionless hopping strength J/U the spin description holds. Second, at finite temperature, thermally excited defects can lead to a break down of the spin-description even for values of J/U where the XXZ model would be a very good approximation at zero temperature. In this case, a crucial question concerns the temperature range in which an observation of spin-entanglement is experimentally feasible.

A. Influence of Quantum Fluctuations
To investigate the influence of quantum fluctuations, we studied the situation U J, where the system is in a Mott insulating phase but particle fluctuations are not negligible [40] using algorithms based on Matrix-Product-States [67]. A system of size L = 30 is initialized in the ground state of Hamiltonian (14) in the sector where N − = L and N + = 0. Since we consider the regime 20 ≥ U/J ≥ 3.5, this correspond to a Mott insulating phase of the σ = − bosons in the thermodynamic limit [68]. We subsequently perform a spin flip for the central spin of the chain using the pro- tocol: |n + = 0, n − = 0 L/2 → |0, 0 L/2 and |0, n − L/2 → |1, n − − 1 L/2 . With this protocol, we need to consider a local Hilbert space that has to accommodate at most one σ = + boson per site, simplifying the numerical simulation. For the σ = − bosons, we truncate their local Hilbert space to four occupancies, with the further constraint that there can be at most four particles per site (the state |1, 4 L/2 is thus discarded). The system is then evolved in time with Hamiltonian (14) using a time-evolving block decimation algorithm (TEBD) [69]. During the time-evolution the maximal allowed bond link is D = 3000. In Fig. 2, we show the concurrence C(ρ 1,1 A,B ) for the subsector with a single particle per site for the sites A = L/2 + 1 and B = L/2 − 1 as a function of time for several values of U/J. Oscillations have a clear U/J 2 period, which is the time-scale associated with the typical energy scale of spin dynamics J ex . A clear decrease of the maximum concurrence for lower U is visible. One reason for this decrease is that for lower U , the probability p 1,1 to find a single particle per site is reduced, which corresponds to a reduced trace ofρ 1,1 A,B . To check for this effect, we compare the dynamics to the prediction from Eq. 2 for the XX model rescaled by p 1,1 (t = 0). The curves for the rescaled XX dynamics are shown in Fig. 2 as solid lines. For U/J 8, the dynamics appears to be well described by the rescaled XX predictions, indicating that effective spin dynamics in the sector with one particle per site are undisturbed by the presence of number fluctuations. For lower U/J, the concurrence C(ρ 1,1 A,B ) is smaller than predicted by the rescaled solution. We now inspect the reduced density operatorρ 1,1

A,B
more closely by comparing the quantities P ↑,↓ P ↓,↑ and |ρ ↑↓ |, which are equal in the spin case [see Eq. (3)]. The equality of both quantities signals fully coherent dynamics. In Fig. 3, we show the time evolution of both quantities for several values of U/J. Interestingly, the two quantities take similar values down to U/J ∼ 6. The fact that we observe |ρ ↑↓ | < P ↑,↓ P ↓,↑ for lower values of U/J can be interpreted as effective decoherence dynamics.
Finally, we consider a global measure of entanglement in the system by investigating the sum of the squared concurrences: The motivation for summing over the square of the concurrences stems from the monogamy inequality [70], which holds for spin-1/2 systems. For the ideal spin dynamics in the XX-Hamiltonian, C 2 (t) = 4(1 − A |φ A (t)| 2 ) → 4 for long times. In Fig. 4, we show C 2 (t) for several values of U/J. For U/J 8, C 2 (t) increases with time, and the prediction of the XX chain weighted with the probability p 1,1 (t = 0) (solid lines) captures the behavior. For smaller U/J, stronger deviations are visible, which indicates decoherence.

B. Influence of Thermal Fluctuations
We now turn to the effects of number fluctuations introduced by a finite temperature and a chemical potential. In order to single out these effects, the hopping strength is set to J = 20U , for which the dynamics at zero temperature is well captured by the XX prediction as shown in the previous section. To this end, we performed exact diagonalization of a system of size L = 6 and included only the lowest-energy part of the Hilbert space. As the hopping of particles is only a small perturbation, we consider only Fock states with an interaction energy smaller than a given energy cutoff: (U/2) n σ,j (n σ ,j − δ σ,σ ) < E c . The system is initialized in the grand canonical ensemble of σ = − bosons with temperatures in a range T = 0 − 0.1U/k B (k B , Boltzmann constant) and chemical potentials in a range µ = 0.25U − 0.75U . We perform the same flip protocol as in the previous section and let the system evolve in time with HamiltonianĤ BH (14). The simulation includes a total number of particles N tot = (N + + N − ) ∈ [4,8] and the cutoff energy is E c = 3U + µ(N tot − L). Convergence of the simulations upon inclusion of more particle sectors and more states has been verified and an error on the order of a few percents is estimated, which is better than the expected experimental precision. In order to test the influence of the relatively small size of L = 6 on the time-evolution, we compared the zero-temperature concurrence spreading at L = 6 with the spreading at L = 30 with the TEBD (see previous section). We find that until time t ∼ 0.5 U/J 2 , the two predictions agree within a few percent. Even if finite-size corrections are expected to be more significant at higher temperatures (we compared the data with those at L=5, not shown), the data in Fig. 5 should be sufficiently accurate to predict the behavior of typical experimental systems with L ≈ 15 − 20 within a few percent. In Fig. 5, we show the concurrence C(ρ 1,1 L/2−1,L/2+1 ) for chemical potentials µ = 0.25U, 0.5U, 0.75U and several values of the temperature T . For increasing temperatures the signal drops. This reduction with temperature is relatively small at µ = 0.5U compared with the other two chemical potential values. This can be attributed to the fact that the gap to excited states is largest, and thermal excitations are thus suppressed, at µ = 0.5U in the limit J/U = 0 [39]. This statement holds in approximate form also for J/U = 1/20. Additionally, the effect of an increase of the chemical potential from the optimal value µ ≈ 0.5U to µ = 0.75U is more damaging to the entanglement than a decrease to µ = 0.25U (compare left and right plot in Fig. 5). This dependence on the chemical potential highlights the importance of tuning the chemical potential at the center of a trapped system to µ ≈ 0.5U . Similar to the case of quantum fluctuations, we check whether decreased concurrence can be ascribed to the reduced population of the single-occupancy sector. Solid lines in Fig. 5 represent the prediction of the XX model rescaled by p 1,1 (t = 0). Whereas the XX model captures the features of the entanglement dynamics for low temperatures and for µ = 0.5U , it fails at the highest temperatures considered. Concluding, we provided evidence that the entanglement propagation scheme previously described can be carried out in a realistic parameter range for experiments. For current temperatures of T ≈ 0.1U/k B [25,39], a drop of the concurrence signal by maximally a factor of two compared to the zero temperature situation is to be expected due to number fluctuations introduced by finite temperature in the grand canonical ensemble. Therefore, the signal should be strong enough to be experimentally detectable.

VI. SCHEME IN THE PRESENCE OF NUMBER FLUCTUATIONS ASSUMING FULL SPIN-RESOLUTION
We now turn to the description of an entanglement detection scheme for a lower bound of the concurrence C(ρ 1,1 A,B ). In this section, we assume that the measurement can be done with full spin-resolution, that is, the individual populations n ± of both species can be detected in a single experimental run.

A. Observable with full spin-resolution
Restricting ourselves to global pulses, the most general observable, in this case, is the joint probability of finding n ± A atoms of species ± on site A, and n ± B atoms of species ± on site B after a global pulse with angles θ and φ. This probability can be written in terms of the diagonals of the reduced density operator after the global pulse with angles θ and φ. The rotation operatorR(θ, φ) l is a generalization of the spin-1/2 rotation (6) to arbitrary local total spins j l . It can be obtained using the transformation of the creation operatorsâ † l,± for species ± on site l a(θ, φ) † l,+ = cos(θ/2)â † l,+ + ie −iφ sin(θ/2)â † l,− a(θ, φ) † l,− = ie iφ sin(θ/2)â † l,+ + cos(θ/2)â † l,− , which yields the mappinĝ for the basis states |j l , m l [66].
B. Observable for random phase φ As discussed in Sec.III B, the phase φ of the global pulse is assumed to be random. Therefore, we consider an averaged density operator Additionally, we are interested in the probabilities p n A ,n B for observing the total atom numbers n A and n B . They can be detected by summing over P n + A ,n − A ,n + B ,n − B (0) with the constraint that n + A + n − A = n A and n + B + n − B = n B :

C. Lower bound
We will now derive a lower bound for the concurrence ofρ 1,1 using the probability of finding a spin-up atom on each of the sites A and B.
The reason for focusing on P ↑,↑ (θ) will become apparent in Sec. VII B. Using the rotation formula, we find the important result P ↑,↑ (π/2) = p 1,1 4 where denotes the real part. The key point is that one can still detect [ρ ↑↓ ] in the presence of number fluctuations using With the same reasoning as in Sec. III, we find the lower bound for the concurrence in the spin-1/2 sector Current implementations of single-site resolved imaging in optical lattices do not resolve the individual atom numbers of both species [25]. Instead, the procedure is to push out one of the species using a resonant pulse and to detect the remaining atoms. For concreteness, we assume that the minus component is pushed out. The observed probability for the atom numbers n + A , n + B of the remaining plus-atoms is then The detected signal therefore mixes contributions from different minus-atom numbers.
In addition to the probabilities after push-out, one can also simply image without push-out pulse. The observed probability thus corresponds to measuring the probability p n A ,n B for the total atom numbers n A , n B according to Eq. (17).

B. Analysis of the problem
A key obstacle for formulating a lower bound without spin-resolution is to extract [ρ ↑↓ ] from the detected sig-nalP n + A ,n + B (θ). For deriving bounds for the concurrence, we will useP 1,1 (θ), which does not contain a signal from empty lattice sites. Writing out Eq. (21), we find Using the rotation formula (16), we can express the unwanted contributions in the second line in terms of probabilities before the pulse, P n + A ,n − A ,n + B ,n − B (0), for states with at least one of the sites occupied by two or more atoms of the same species. Consequently, these terms vanish for fermionic atoms in a single band Hubbard model [51]. A suppression of doubly occuped sites for bosons is possible if the local chemical potential µ in optical lattice experiments is tuned to lower values (0 µ 0.5U ) at the expense of increasing the probability for holes [39]. Further, for experiments with Rydberg atoms in optical tweezers [52][53][54][55][56], the filling of the traps is typically only zero or one. In these situations, the terms in the second line of Eq. (22) vanish and the bound (20) can still be used without full spin resolution. For situations when doubly occupied sites of the same species cannot be neglected, modified bounds can be found by making certain assumptions onρ A,B . We outline two methods in the following sections.

C. Lower bound based on subtraction of 1/4
A modified version of the bound (20) can be derived, using the following assumptions: • A1 The probability of finding sites occupied by three or more atoms before applying the pulse is negligible: There is no coherence between a state with two minus-atoms on A and two plus-atoms on B and a state with two plus-atoms on A and two minusatoms on B: ρn+ Assumption A1 is well fulfilled in the deep Mottinsulating regime of the Bose-Hubbard model at unity average filling for realistic experimental temperatures T ≈ 0.1U [25,39]. A breakdown of the second assumption A2 would require a non-negligible probability for having one of the sites occupied by two minus-atoms and the other site occupied by two plus-atoms. Again, for the Bose-Hubbard model at unity average filling for realistic experimental temperatures [25,39], the joint probability of having both sites doubly occupied (independent of the spin) is much lower than one percent. Hence, assumption A2 is typically valid. Using A1 and A2, one can show that Using the fact thatP 1,1 (0) ≥ P ↑,↑ (0) andP 1,1 (π) ≥ P ↓,↓ (0), we arrive at a corresponding bound The bound (23) works with rather weak assumptions but is not particularly tight. The reason is that the subtraction of 1 4 instead of p 1,1 4 leads to a reduction of the bound when p 1,1 is significantly smaller than one. Additionally, there is no absolute value around the first term, which can lead to a negative signal if [ρ ↑↓ ] < 0.

D. Lower bound based on correlations
Therefore, we derive an improved bound compared to Eq. (23) based on evaluating the quantitȳ The subscript c for P c 1,1 (θ) stands for 'connected' because P c 1,1 (θ) resembles the form of a connected correlation function. Using Eq. (19), we find where, for j = A or B, p 1,j is the single-site probability of finding the total atom number n j = 1 and P 1,n − j (θ) is the probability of finding a single plus-atom and n − j minus-atoms after a θ-pulse. The connection to the corresponding joint probabilities is analogous to Eq. (24). The signal P c if the following two assumptions hold: • B1 The probability of finding a single atom on site A is independent of the probability of finding a single atom on site B: p 1,1 ≈ p 1,A p 1,B .
• B2 There are no correlations in the sectors with local occupation number higher than one: With ρ nj j we refer to the single-site reduced density operator for site j projected onto local atom number n j (see Eq. (10)).
A mechanism that would violate these assumptions is the introduction of density-density correlations via quantum fluctuations in the form of particle-hole pairs [40]. However, these correlations are extremely small beyond nearest-neighbor distances. A potential danger arises if the system is brought out of equilibrium, for example, via a fast quench, which can induce longer-range densitydensity correlations [71]. This can be avoided with a careful adjustment of lattice ramps. Additionally, there are experimental checks for the validity of Eq. (26), such as an observation of the checkerboard pattern described in Sec. III C on a zero background signal, that is,P c 1,1 (π/2) ≈ 0 for even distances. Further, the absence of density-density correlations can be checked using imaging without a push-out pulse. Based on Eq. (26) we can formulate a lower bound G c (ρ A,B ) := 4 P c 1,1 (π/2) − 2 P 1,1 (0)P 1,1 (π), (27) which holds if the assumptions B1 and B2 are fulfilled.

E. Influence of parity-projection
In the current experiments with single-site resolution, only the parity of the on-site occupation number can be observed due to a pair-wise loss from light-assisted collisions [38,39]. The parity-projection only occurs during the actual detection of the remaining species but not during the push-out [25]. The observed probabilities after push-out and subsequent parity projection are thus where n + A , n + B < 2. The additional terms that enterP 1,1 (θ) all stem from triply or higher occupied sites. These terms vanish if A1 is fulfilled. Consequently, the bound (23) can still be used with parity-projection. Similarly, parity-projection adds several terms to Eq. (25) which all vanish if B2 holds. Therefore, bound (27) also remains unaffected.

VIII. CONCLUSION AND OUTLOOK
In conclusion, we proposed a scheme for detecting lower bounds for the concurrence of two sites of a lattice many-body system, which could be used for measuring spin-entanglement in quantum magnetism experiments with coexisting spin and number fluctuations. Our analysis showed that a detection of the lower bounds should be possible in current high-resolution imaging setups for quantum gases in optical lattices [24][25][26][37][38][39][40] despite several technical limitations. However, the scheme would simplify if full spin-resolution was achieved experimentally, and the bound (20) could be used. A possible solution for one-dimensional systems is to prepare a single chain of atoms and let the atoms tunnel orthogonally to the chain before the detection. If a magnetic field gradient is applied during the orthogonal dynamics, atoms with positive and negative magnetic moment would spatially separate. The spatial separation could allow a detection of the local occupation numbers of both spin states in a single experimental run. In this sense, an in-situ Stern-Gerlach experiment could be realized with full spatial resolution along the onedimensional chain. Such a scheme could also be useful to detect the correlations induced by impurities in strongly interacting superfluids, enabling the direct imaging of a polaron cloud [25]. Concerning the actual influence of on-site number fluctuations on spin-entanglement, we performed numerical simulations of spin impurity dynamics in the one-dimensional Bose-Hubbard model. The effect of quantum fluctuations within large parts of the Mott insulating phase could be captured by a renormalized XX-spin dynamics normally only valid in the very deep Mott insulating regime. A similar behavior results from thermally activated number fluctuations. Importantly, our simulations showed that the entanglement generation and spreading survives for the temperatures and parameters of current experiments [25,26]. Thus, the application of the proposed detection technique for this type of experiment should be immediately possible. More generally, the experimental detection of spin-entanglement in Hubbard models realized with optical lattices [21][22][23][24][25][26][27][28] is now within reach.