Detecting high-dimensional entanglement in cold-atom quantum simulators

,


I. INTRODUCTION
Since its initial conception inspired by the EPR paradox [1], quantum entanglement has been identified as a key aspect in the understanding of a plethora of physical phenomena, such as the dynamics of disordered spin systems [2], the thermalization of closed quantum systems [3,4], and even in the context of the black-hole information paradox [5].In recent years, much attention has been directed toward the effects of entanglement in condensed matter, where it has been linked to topological properties of quantum states [6,7] and quantum phase transitions [8][9][10], among others [11].Studying entanglement in these macroscopic systems directly is oftentimes too challenging due to limited experimental control and measurement capabilities.
The development of experimental systems offering quantum control on the level of single particles over recent decades has enabled an alternative approach to studying such genuine quantum phenomena.To simulate complex quantum systems, one constructs simpler synthetic systems, called quantum simulators, which mimic, or emulate, the dynamics of the system of interest.In particular, cold atoms trapped in lattice geometries have evolved into the leading platform for quantum simulation of condensed matter systems [12][13][14][15][16][17][18].Through the application of external fields, model parameters can be tuned within a broad regime ranging from strong repulsive to attractive interactions, equipping the system with an ideal framework to simulate highly entangled quantum states with single-atom-resolved readout [19][20][21][22].The capability to detect entanglement in these platforms is crucial for the investigation of the aforementioned phenomena, but still faces challenges [23].Many experimentally available criteria can, in fact, only indicate ("witness") the existence of entanglement in a state qualitatively [24].
In this work, we want to go beyond detecting the mere presence of entanglement and instead make statements about the entanglement structure.The standard measure of entanglement for bipartite pure quantum states ρAB = |ψ⟩⟨ψ| is the entanglement entropy, defined as S(ρ A ) = S(ρ B ) = − d i=1 p i log p i , with the reduced density matrix ρA = Tr B (ρ AB ) (ρ B analogously) and its eigenvalues p i [4].Even though in many cases much can be learned from this quantity, it contains less information than the full eigenvalue spectrum, also known as the entanglement spectrum, from which it is derived.Therefore, more recently, the entanglement spectrum itself has been used extensively to investigate the role of entanglement in various phenomena, including fractional quantum Hall states [25], topological insulators and superconductors [26], one-dimensional (1D) systems in the scaling regime [27], emergent irreversibility [28,29], and many-body localization transitions [30,31], leading to new insights.Furthermore, the ability to prepare and certify states with a broad entanglement spectrum would arXiv:2305.07413v2[quant-ph] 9 Feb 2024 enable the execution of quantum algorithms that exploit this property for enhancing run time and robustness [32][33][34].
The number of nonvanishing terms in the entanglement spectrum is known as the entanglement dimension, or Schmidt rank, of the state.It represents the number of terms needed to faithfully represent the quantum state in the product Hilbert space (with generalizations established for mixed states).Standard methods to obtain the entanglement dimension for cold-atom systems available today are based on full state tomography, or on efficient fidelity-measurement schemes, for which the number of required measurement bases scales quadratically, or linearly, respectively, with the local Hilbertspace dimension L [24].Recently, advanced methods for accessing information about the entanglement spectrum have been proposed, including Hamiltonian learning [35][36][37], random measurement schemes [38,39], and ancillary-system-based readout protocols [40].However, these approaches either make assumptions about the prepared states potentially leading to bias, or pose stringent requirements on experimental capabilities (for a more detailed discussion, see Sec.VII B).
We propose an alternative approach to detecting highdimensional entanglement in systems of lattice-confined ultracold atoms.Our method is inspired by earlier findings for entangled photon pairs in different polarization states [41].In that work the authors construct a measurable lower bound on the state fidelity to a highly entangled reference state.This approach provides a powerful tool as one can define a set of fidelity thresholds with each threshold corresponding to a matching minimum entanglement dimension of the measured state [42].Bounds on the fidelity to the reference state thus naturally translate to bounds on the entanglement dimension of the prepared quantum state.One can construct such a bound by measuring in only two mutually unbiased bases (MUB) |i⟩ m and |j⟩ n , i.e., ∀m, n : ⟨i|j⟩ m n = L −1 , simplifying the experimental procedure significantly.However, implementation of two such MUB measurements for coldatom systems is a challenging problem.
Our main contribution is to derive lower bounds on the fidelity to highly entangled reference states that only require position-and momentum-correlation measurements, generalizing previously reported bounds in several ways.Both the position and momentum bases can be accessed by measuring the atom positions in situ and after TOF expansion [22,43], techniques that are well established experimentally [44].The fidelity bounds directly yield bounds on the entanglement dimension and thus measurable Schmidt-number witnesses.Furthermore, we show that this protocol is applicable to a large class of reference states, to bipartite systems with multiple indistinguishable particles per species (party) for both fermions and hard-core bosons, and even to a multipartite setting.One might expect that a bound based on the fidelity to a reference state gives satisfactory results only for experimental states close to that reference, i.e., for states the reduced density-matrix spectrum of which is similar to that of the reference state.Our findings indicate, however, that our bound detects high-dimensional entanglement for a broad range of quantum states, even in the presence of strong decoherence.The bound turns out to be robust against typical experimental noise sources and its tightness decreases at most linearly with the noise strength, i.e., with the impurity of the prepared state.
In the remainder of this work, we first establish a fidelity bound for a pair of two entangled atoms in an optical lattice in Sec.II and test its robustness regarding typical experimental noise using a Hubbard model in Sec.III.Subsequently, we generalize the method to multiple indistinguishable atoms per species (Sec.IV) and to a multipartite setting, where more than two different atomic species are entangled (Sec.V).In Sec.VI we derive fidelity bounds for extended classes of reference states.Our conclusions and a discussion of our results are provided in Sec.VII.

II. BOUND ON ENTANGLEMENT DIMENSION
Any bipartite pure quantum state on a product Hilbert space H = H A ⊗ H B can be represented in Schmidtdecomposed form |ψ⟩ AB = k i=1 λ i |i⟩ A ⊗ |i⟩ B , a basis choice that minimizes the number of contributions k (also known as Schmidt rank or entanglement dimension, D ent ) needed to represent a given quantum state [45].Any separable state can be written through one tensorproduct term alone and therefore k = 1.The defining feature of entangled states is that this no longer holds true and thus k ≥ 2, as can be seen on the example of the singlet state |Ψ⟩ EPR = (|↑↓⟩ − |↓↑⟩)/ √ 2. The composition of these tensor-product contributions and their weights defines the entanglement structure of a quantum state.Gaining information on the full structure is an exceedingly hard problem for multidimensional states, i.e., states entangled in several internal degrees of freedom, or multipartite states, i.e., states made up of three or more entangled parties.Determining the entanglement dimension instead is both insightful and experimentally feasible, as we show in this paper.We start with the case of attractive interactions and later, in Sec.VI, we generalize to the repulsive case.
The entanglement dimension of a bipartite system is bounded by the size of the smaller of the two local Hilbert spaces, k max = min[dim(H A ), dim(H B )] =: L. In the remainder of this work we take the local Hilbert-space dimensions to be equal and finite.One can choose a maximally entangled state (MES) of the system that has equal coefficients for all L terms, give some intuition, the two subsystems will later be the two atoms in the lattice, where m labels the lattice sites.The fidelity of the experimental state ρ to the reference MES, implies a convenient state-distance measure to compare the two states, as it is bounded as a function of the entanglement dimension of ρ.One can explicitly construct a set of bounds B k on the fidelity to the MES, F (ρ, Ψ MES ), given by which hold for any experimental state ρ with entanglement dimension D ent ≤ k [42,46].In the case of mixed states, ρ = i p i |ψ i ⟩⟨ψ i |, the notion of an entanglement dimension has to be extended to the so-called Schmidt number.This is defined as the maximum entanglement dimension of the pure parts |ψ i ⟩ of the state, minimized over all possible pure-state decompositions: . The violation of the relation in Eq. (3) for given k therefore indicates that ρ is entangled with a dimension of k + 1 or higher.This not only gives a robust entanglement witness, since the lowest threshold B 1 = L −1 already indicates entanglement, but also bounds the width of the entanglement spectrum and hence gives insight into the entanglement structure.We exemplify this in Sec.III C for localized dimer states in a Hubbard model where the detected entanglement dimension correlates with the number of macroscopic Schmidt coefficients.Additionally, one can use the fidelity to construct lower bounds on the entanglement of formation, as has been shown in Refs.[41,47], establishing F (ρ, Ψ MES ) as a versatile source of information about the entanglement content of ρ.Nonetheless, fidelity measurements come with significant experimental complexity, in general requiring measuring in L + 1 different bases for an L-dimensional local Hilbert space [41].
In the following, we establish a lower bound F (ρ, Ψ MES ) on the fidelity F (ρ, Ψ MES ) accessible to experiments with cold atoms in optical lattices (or arrays of optical tweezers).It only requires measurements in two bases, independent of the local Hilbert-space dimension given by the lattice size.We start with the case of two distinguishable atoms here and develop generalizations to higher atom numbers and other reference states in later sections.The two atoms constitute the two entangled parties and their local Hilbert spaces are spanned by the discrete position states (sites) of the atoms in the opticallattice potential.Consequently, the MES for this product basis is a superposition state with both atoms located at the same lattice site, in a superposition summing over all L sites with equal probability [Eq.( 1)].
It is insightful to split the fidelity into two sums, dividing the contributions into state populations (lefthand sum) and two-particle coherences F coh (right-hand sum).The state populations of the two distinguishable species, in the following labeled as ↑ and ↓ with their corresponding number operators n↑ and n↓ , can be obtained by spatially discretizing the joint density distribution ⟨n ↑ (x 1 )n ↓ (x 2 )⟩.It can be probed directly through single particle resolved fluorescence imaging, realizing highprecision in situ measurements [22,44,48].A representation of ⟨n ↑ (x 1 )n ↓ (x 2 )⟩ for the ground state of a Hubbard Hamiltonian with L = 6 at U/J = −12 is displayed in Fig. 1 Such direct experimental access is not available for the two-particle coherences F coh , but one can instead bound F coh from below by measuring in a second basis.A natural choice for cold atoms is the momentum basis, as the system comes with a native implementation of the corresponding basis change, the Fourier transformation.It can be applied efficiently by rapidly switching off the lattice potential and interactions and subsequently letting the atoms propagate in a weak harmonic potential for t = T /4 with trap oscillation period T before taking a fluorescence image [22,44,49].By repeatedly preparing and measuring a state with this scheme, one acquires samples from the momentum correlation function ⟨n ↑ (k 1 )n ↓ (k 2 )⟩.
To show how to utilize this to bound state coherences from below, we construct the corresponding measurement operator by stating the effect of the Fourier transform on the localized Wannier basis functions of the lattice potential.The basis function for the nth lattice site can be expressed as ω(x − nd) due to the discrete translational invariance of the lattice, where d is the lattice spacing.Any shift in position space causes a phase factor in momentum space, so one obtains for the single-atom wave function in momentum space with ω(k) being the Fourier transform of the Wannier envelope.[44].Using the field operators, defined via lattice-site creation (annihilation) operators ĉ † j,σ (ĉ j,σ ), one can represent the particle number operator in momentum space in the position-space basis {|j⟩ | j ∈ {1, ... L}} as |m⟩⟨n| e id(m−n)k .
(7) As only one particle per species is present in the lattice, no differentiation between fermions and bosons has to be made here.The full expectation value ⟨n ↑ (k 1 )n ↓ (k 2 )⟩ in the density-matrix picture is given by the trace over the product of the two momentum number operators and the density matrix, ⟨n Finally, by exploiting the cyclic property of the trace, one arrives at the following expression: and a phase factor obtained through the Fourier transformation [48].The above-given description is naturally rewritten in terms of a new set of basis functions by bundling terms with the same complex phase factors and their conjugate counterparts into trigonometric basis functions of the two lattice momenta k 1 and k 2 .The full momentum correlation function then reads The above-mentioned basis weights Re(g αβ ) and Im(g αβ ) in Eq. (10a) are sums over the real and imaginary parts of the coherences of the density matrix ρ [see Eq. ( 10b)].Each coefficient g αβ is defined by the pair of position-space distances for all contributing coherences The set of all coherences contributing to a given coefficient can simply be constructed by shifting all atom positions of one of the coherences along the lattice.Since we have already combined coherences and the corresponding phase factors with their complex conjugates, we have to introduce the index set M in Eq. (10c) to avoid double counting of coherences.For additional information regarding Eq. ( 10), we refer the reader to Ref. [44].Obtaining the coefficients g αβ is not directly straightforward, as the basis {φ R αβ , φ I αβ } is nonorthogonal due to the modulation of the periodic basis functions through the envelope [44].Projecting the measured distribution [cf.Eq. (10a)] onto the basis function set therefore yields smeared-out coefficients c αβ , where each coefficient also contains small contributions coming from the nonvanishing overlap with other basis elements.To overcome this problem, we explicitly construct the linear transformation Q that maps the set of actual basis weights ⃗ G to the measured coefficients where each element of the matrix Q is given by an overlap integral between a pair of basis functions (for details, see Appendix A).These integrals factorize since the Fouriertransformed Wannier envelope factorizes as well; consequently, only a small number of 1D integrals linear in the number of lattice sites must be computed to construct Q.
The actual basis weights g αβ are then extracted by formally inverting Q and rewriting Eq. ( 12) as Numerically, we employ a conjugate-gradient method to determine ⃗ G.The two projection integrals in Eqs.(11a) and (11b) can be evaluated in a simplified way using Monte Carlo importance-sampling techniques.By treating the momentum correlation function as a normalizable multivariate probability density, it can be absorbed in a redefinition of the integration variable.The remaining integrals, are then directly evaluated through the measured or simulated momentum correlation samples.This evaluation method enables scalability to higher atom numbers introduced later, as the Monte Carlo integration error scaling is independent of the integral dimension, while also reducing the variance of the integrand at the same time.
We make some additional comments regarding synthetic data generation and efficient computation of Q in Appendix A.
At this point one has obtained access to basis weights g αβ equal to sums over subsets of coherences of ρ.However, not only the two-particle coherences relevant for F coh in Eq. ( 4) are contained within the basis weights g αβ but also different-site two-particle coherences that do not contribute to the fidelity F (ρ, Ψ MES ).We note that any general density-matrix element is bounded from above by using Cauchy-Schwarz inequality with the right-hand side containing only already measured state populations and thus adding no new experimental complexity.For pure states, the second inequality in Eq. ( 15) is obviously tight but it grows looser with increasing mixedness of the state.In the twoatom case presented here, the subset of (α, β) ∈ M that carries relevant two-particle coherences reduces to α = β =: δ ∈ {1, . . ., L − 1}.The desired sum of relevant coherences can then be lower bounded by subtracting the bounds in Eq. ( 15) for all noncontributing coher-ences from the sum of relevant basis coefficients, where the second sum of the last expression covers all nondimer coherences.Together with the same-site populations displayed in the first sum of Eq. ( 4), we formulate the complete experimentally accessible lower bound on the fidelity of the experimental state ρ to Ψ MES as Inserting this bound in Eq. (3) yields our first main result, where the fidelity bound F constitutes an entanglementdimension witness and is obtainable directly through fluorescence measurements, in situ and after TOF.Thus, if F (ρ, Ψ MES ) exceeds B k (Ψ MES ) for some given k, the state ρ is certified to be entangled in at least k + 1 dimensions.
The proposed experimental protocol can be summarized as follows: One prepares an ensemble of two atoms of different species in a periodic potential in some state of interest.
1.By single-atom-resolved detection, one measures the position-space correlation function ⟨n ↑ (x 1 )n ↓ (x 2 )⟩.The signal is discretized by identifying the atom positions obtained in each shot with a pair of lattice sites, which yields the positionspace populations ⟨mn| ρ |mn⟩ entering in Eqs. ( 16) and (17).
2. The momentum-space distribution is probed through ballistic TOF expansion, resulting in an effective Fourier transformation of the wave function.
The coefficients c αβ are obtained from the measured momentum correlation function, ⟨n ↑ (k 1 )n ↓ (k 2 )⟩, by computing the overlap between the measured distribution and the trigonometric basis functions {φ R αβ , φ I αβ } [Eq.( 9)], i.e., by evaluating the basis functions using the sampled momenta.From these, the corrected expansion coefficients g αβ are obtained via Eq.( 13) and inserted into Eq.( 16), which yields the desired lower bound on the reference-state fidelity in Eq. (17).
The statistical requirements for confident certification are discussed in Sec.III A and a study of the robustness of the protocol with respect to typical experimental noise effects is given in Secs.III B and III C.

III. CERTIFICATION ROBUSTNESS UNDER REALISTIC CONDITIONS
We study the performance of our method under realistic experimental conditions through numerical simulations.Our model system is a 1D open-boundary Hubbard model, realized by cold atoms in a deep optical lattice (V 0 = 8E r [50]) in the tight-binding approximation [51].Due to limited wave-function overlap between sites, all tunneling going beyond adjacent sites is suppressed.Through the application of external magnetic fields, Feshbach resonances can be utilized to implement an effective on-site atom-atom interaction with a highly tuneable interaction strength [52].The dynamics of the system are captured by the Hamiltonian with the tunneling strength J, interaction strength U , creation (annihilation) operator ĉ † i,σ (ĉ i,σ ) for an atom on site i and in spin state σ ∈ {↑, ↓}, and their corresponding atom number operators ni↓ , ni↑ [53].This system is characterized by the ratio U/J (J > 0), where negative values correspond to attractive and positive values to repulsive interactions.In the simple case of two distinguishable particles, both Fermi-Dirac and Bose-Einstein statistics produce the same dynamics.The Hubbard model was chosen due to its simplicity and widespread use in numerical modeling [16,51] but our readout scheme is also applicable to other lattice Hamiltonians.
In the remainder of this section, we consider the ground state of the two atoms in a lattice of size L = 6 with attractive interactions at U/J = −12 and use 2.5×10 4 momentum-space and 1×10 4 position-space samples for certification, unless specified otherwise.Later, in Sec.VI, we will also consider repulsive interactions, where robust entanglement certification is achieved by adapting the employed reference state.In the configuration given above, ρ is entangled in all six lattice degrees of freedom, meaning that D ent = 6, and thus serves as a suitable test state for our entanglement-detection scheme.Figure 1(a) shows a representation of the position-space probability distribution.Through exact diagonalization, we find that the fidelity F (ρ, Ψ MES ) increases with growing attractive interaction strength (blue line in Fig. 2) but the ground state does not converge to Ψ MES (F (ρ, Ψ MES ) < 1).Knowledge of the exact fidelity would enable us to certify five out of the six entanglement dimensions for moderately attractive interactions [ F (ρ, Ψ MES ) > B 4 for U/J ≲ −6].The offset in fidelity with the MES is an effect of Figure 2. The dependence of the fidelity F and the fidelity bound F on the interaction-to-tunneling-strength ratio U/J for pure (r = 0) and dephased (r ∈ {0.05, 0.15}) ground states.The B k thresholds are the horizontal dashed lines such that fidelities above any B k indicate at least k + 1 entanglement dimensions.Both F and F increase with growing attractive interaction strength before saturation.The tightness of F decreases with increasing mixing rate r.The statistical error bars are small and barely visible.
the finite system size, as central sites are energetically favored for open boundary conditions, since more tunneling pathways are available [see Fig. 3(a)], making the distribution of populations nonuniform.Since the Schmidt coefficients are given by the double-occupation probabilities in the strongly attractive limit, this behavior translates to a nonuniform entanglement spectrum.In the case of the pure ground state, we find our fidelity bound to be tight (blue markers in Fig. 2).The use of our protocol therefore yields the highest certifiable entanglement dimension achievable using the fidelity to the MES.
In the following subsections, we discuss the requirements of our detection scheme in terms of its robustness with respect to typical experimental imperfections and noise sources, starting with finite measurement statistics in Sec.III A. In Sec.III B we study the effect of generic dephasing noise arising in experiments, which can be caused by, e.g., fluctuations of trapping light parameters or control fields during state preparation.Finally, in Sec.III C, we consider fluctuations of the depth of individual lattice sites that are characteristic of realizations using arrays of optical tweezers.Randomized potentials lead to localization of atom pairs and thus potentially to a reduction of the entanglement dimension, an effect that becomes observable through our detection scheme.Moreover, in Sec.III D, we discuss the behavior of our approach in the limit of large lattice sizes.Additionally, in Appendix C, we have simulated the performance of the bound on thermal ensembles .Readers more interested in generalizations of the scheme may jump to the last paragraph of Sec.III, where our results on noise robustness are summarized.For the simulation results presented in Fig. 2, the first two effects are already addressed within the simulation.Finite measurement statistics induce fluctuations of the certified fidelity and thus impact entanglement detection.Additionally, experimental quantum state realizations are in general not pure wave functions |ψ 0 ⟩ but face mixing and decoherence.A simple model to account for this is to replace |ψ 0 ⟩ with a dephased density matrix ρ = (1 − r) |ψ 0 ⟩⟨ψ 0 | + rL −2 1, with a mixing parameter r related to the state impurity p = 1 − p, washing out the probability distribution.Both effects have been in-cluded to produce dephased ground states in Fig. 2, each with sampled correlation functions at mixing strengths r = 0.05 (p ≈ 0.095) and r = 0.15 (p ≈ 0.270), respectively (orange and green data sets throughout this work).These model alterations clearly lead to loss of bound tightness and add random noise to the fidelity bound F , as visible in Fig. 2. We discuss these matters in more detail in the following.

A. Sampling statistics
In experiments, both the joint position-space distribution ⟨n ↑ (x 1 )n ↓ (x 2 )⟩ and the momentum-space distribution ⟨n ↑ (k 1 )n ↓ (k 2 )⟩ are probed by repeated state preparation and measurement, each experimental run providing one sample point drawn from the respective distribution.The finite sample numbers are the cause of statistical errors in our fidelity-bound estimation.In this section, we systematically explore the scaling of the standard error (SE) of our bound F (ρ, Ψ MES ) with regard to the sample size to determine how many samples are required for acceptable error margins.The position-space distribution can be obtained directly in discretized form with L 2 different outcomes, whereas the momentum-space distribution is continuous in k 1 and k 2 and needs to be processed via Monte Carlo integration, demanding more samples.We therefore put special emphasis on the momentum distribution in the following and fix the number of positionspace samples to N pos = 1×10 4 .
To analyze scaling properties with regard to available measurement statistics, we compute the fidelity bound F for a wide range of synthetic momentum-space sample sizes N s .The results for different values of r ∈ {0, 0.05, 0.15} are presented in Fig. 4(a).For p = r = 0 (blue data set), the average of the distribution (dashdotted line) and the true state fidelity coincide, indicating that we can reconstruct the right fidelities without bias.The SE σ F of the distribution for different impurities and sample numbers is shown in Fig. 4(b).We report no significant dependence of σ F on the impurity and find a power-law behavior with exponent b = (−0.48± 0.02) [Fig.4(c), computed with the r = 0 data set], consistent with the expectation of Monte Carlo error scaling The complete set of all fitting parameters for this and all following numerical fits can be found in Appendix B. At high momentum-space sample numbers, we observe a saturation of the error, as the number of position-space samples has been kept constant and the corresponding statistical fluctuations start to dominate.We conclude that 1×10 4 position-space samples and 1.2×10 4 momentum-space samples are sufficient to reduce the SE to σ F < 0.01, independent of the state impurity.
When the fidelity lower bound is used to certify the entanglement dimension of the experimental state, the statistical requirements for faithful certification solely depend on the distance to the next threshold value B k .Fidelity-bound values directly in the middle of two B k lines maximize this distance and have the highest admissible margin of error, whereas fidelities close to thresholds call for ever-increasing sample sets to provide the needed accuracy.The measured bound value can be monitored on the fly to adapt the number of samples taken in order to fulfill the statistical requirements.
Our data indicate that surprisingly low sample numbers can be sufficient for robust entanglement detection.For example, the distance to the next relevant threshold B 5 (red dashed line in Fig. 4(a) for the pure ground state at U/J = −12 is given by B 5 − ⟨ F ⟩ > 2σ F , even at only N s = 2000 samples, a statistically significant statement.The fidelity to the MES drops with decreasing attractive interaction strength and, with it, the distance to the next lower fidelity threshold.The sample set sizes should thus be adapted for less attractive interaction strengths.

B. Dephasing noise
Both the fidelity F and the fidelity bound F decrease linearly with the mixing parameter r, as shown in Fig. 5.The bound declines faster than the actual state fidelity; the linear fit slopes are a = −0.76(fidelity) in contrast to ã = −1.15± 0.03 (fidelity bound) [54].Consequently, the tightness gap widens linearly with a slope of a Gap = 0.39±0.03with r.Our bound certifies the same entanglement dimension as could be certified with the actual fidelity for most of the investigated impurity regime r ≤ 0.25 and with only one dimension less in the regime 0.1 ≲ r ≲ 0.16 (Fig. 5).Certification of high-dimensional entanglement therefore remains possible even for significantly mixed states.One concrete decoherence effect potentially arising during state preparation is the presence of a thermal bath, resulting in a Gibbs thermal state with finite temperature, i.e., a mixture of ground and excited states.We discuss this case in Appendix C, finding that our certification scheme is quite robust in the sense that the fidelity bound remains tight up to rather high temperatures.Thus, the generic white noise considered here may well overestimate the typical impact of decoherence on the proposed method.

C. Lattice-potential disorder
Next, we investigate the tightness of our fidelity bound in the presence of lattice-potential disorder, which typically arises in experiments with arrays of optical tweezers where the relative intensities, and thus the depths, of the individual tweezer traps are difficult to stabilize.We introduce a modified Hamiltonian based on Eq. ( 19) including a normally distributed potential depth fluctuation for each lattice site, with the tunneling strength J as the energy scale.It should be noted that the fluctuations are modeled to be uncorrelated, a realistic assumption in the case of optical tweezer arrays but not necessarily for optical lattices.Imperfections in the potential landscape cause localization in the ground-state wave function and decreased fidelity to the reference state Ψ MES , as seen in Fig. 3(c).The resulting composition of the localized state is quite different compared to that of Ψ MES , with strongly peaked double-occupation probabilities around some localization center.Even though the two distributions differ significantly, our method still enables one to certify an entanglement dimension of D ent = 4, demonstrating the wide applicability of the protocol, as we can detect major components of the entanglement spectrum of a state not close to the reference.In particular, it allows us to track the reduction in entanglement due to disorder-induced pair localization, as we discuss in the following.
In the strongly attractive regime of U/J = −12, the states with both atoms at the same lattice site make a contribution of 95.7% to the pure undisturbed groundstate populations.It is therefore a reasonable simplification to treat the atom pair as a bound dimer moving through the lattice.The ground-state localization is then in agreement with the predictions of Anderson localization for disordered potentials, where the occupation probability is suppressed exponentially when going away from the localization center [55].1D systems are expected to localize for any nonzero potential disorder, with the localization length depending on the disorder strength [56].
To investigate the effect of shot-to-shot latticepotential fluctuations on the (detected) MES fidelity, we simulate single experiments on ground-state mixtures of 1×10 3 individual disorder realizations configured according to Eq. ( 20) and compute the disorder ensemble average F over 1×10 3 experimental runs.Both F and F decrease according to a stretched exponential law ∝ exp{−b(Jσ V ) c } with increasing potential depth fluctuation Jσ V and approach the B 1 entanglement threshold, shown in Fig. 6.The bound tightness does not decrease significantly compared to the disorder-free lattice, even for the strongest simulated fluctuation strengths.This is quite remarkable, as each investigated state is a mixture of thousands of individually localized disorder realizations.Consequently, our bound certifies the same entanglement dimension or Schmidt number as the true fidelity for a large regime of disorder strengths.This behavior is markedly different from dephasing noise, where we have found a linearly widening gap between fidelity and bound.The breakdown of the fit at small disorder strengths is caused by the finite size of the lattice.For very weak disorder, the localization length exceeds the lattice scale.In this regime, the fidelity therefore only decreases slowly with increasing disorder strengths, up until single disorder centers can be fully resolved in the lattice.The Jσ v = 0 data point additionally marks the critical point of the localization phase transition in 1D, so anomalous behavior is expected here.Consequently, small disorder strengths do not significantly decrease the fidelity F (ρ, Ψ MES ), and thus our bound also decreases at a reduced rate.

D. Lattice-size dependence of the state fidelity
The scalability of entanglement certification with respect to the lattice size L is of significant concern for experimental implementations.To systematically investigate this, we repeat our previous simulation for a range of different lattice sizes.Our data for a system with finite sampling statistics and dephasing noise show an algebraic asymptotic decline of the fidelity with growing lattice size, as shown in Fig. 7(a) (the fit model and parameters are given in Appendix B).The fidelities and fidelity bounds approach constant nonzero values for L → ∞, depending on the state-mixing strength r.Consequently, certified entanglement dimensions continue to grow as L → ∞.We find that the scaling behavior of our bound depends on the level of dephasing noise.
The addition of lattice disorder changes the situation.From previous data (see Sec. III C), we expect localiza-tion into dimers but the dependence on the lattice size is not immediately evident.Our investigation of a lattice with fixed disorder strength Jσ v = 0.05 yields an exponential fidelity drop-off ∼ exp(−bL), as shown in Fig. 7(b).As our investigated states are mixtures of ground states of different disorder realizations, we do not expect the bound to be tight in the limit of large L. To account for this, we include an offset c in the exponential fit to F .All fits match the data very well at large lattice sizes but significantly underestimate the fidelity in double-and triple-well configurations.Again, finite-size effects offer a plausible explanation for this behavior: the localization length of the system can exceed the lattice size, making it impossible to resolve a localization center fully in small systems.
Since the entanglement-dimension thresholds scale only linearly, B k ∼ L −1 , as compared to the exponentially decaying fidelity, the certified entanglement dimension decreases to D ent = 1 for L → ∞.Consequently, after an initial increase of certifiable entanglement, the entanglement dimension accessible through the bound starts to decline.Based on the reported fit, we extrapolate a maximum certifiable entanglement dimension of D ent = 7 for pure states and D ent = 5 for r = 0.15 for a fixed disorder standard deviation of Jσ v = 0.05.
The above simulations demonstrate the strong impact of site-to-site lattice-potential fluctuations on the scaling behavior of fidelity and thus on the certifiable entanglement dimension for large lattice sizes.While in the case of a disorder-free lattice the certifiable entanglement dimension increases indefinitely with lattice size, disorderinduced localization effects lead to a maximal certifiable dimension reached at some finite lattice size, depending on the disorder strength.We note that site-to-site potential fluctuations may be present for arrays of optical tweezers, while in the case of optical lattices, realized by a single retroreflected laser beam, intensity fluctuations will lead to correlated potential fluctuations not affecting the ground-state properties.Also, the precise properties of the prepared state may depend on the experimental preparation scheme, not discussed in this work.Furthermore, viewing disorder as a feature and tuning its strength deliberately, our method allows the study of dimer localization through the lens of the entanglement spectrum.
Finally, we consider the dependence of the statistical errors on the lattice size, again with fixed sample numbers.For this purpose, we analyze 1×10 4 bootstrap resampling realizations to estimate the SE σ F for lattice chains with lengths up to L = 20, as displayed in Fig. 8.All three simulated mixing rates give qualitatively and quantitatively similar errors.To extract the general trend, we average over the three mixing-rate data sets for improved statistics and fit with an algebraic asymptotic growth model, which shows good agreement in the investigated regime (the model is also listed in Appendix B).Therefore, we find σ F to be largely independent of the lattice size.Scaling up to extended lattice chains is there- The scaling of the fidelity F of the ground state of a flat optical lattice as a function of the number of lattice sites L and the mixing parameter r.The numerical data, including statistical noise, fit well to the asymptotic algebraic behavior.The certifiable entanglement dimension continues to grow with increasing lattice size, as the fidelity asymptotically approaches its infinite-system-size value.(b) A log-linear plot of the fidelity for a disordered lattice of size L with fixed disorder strength JσV = 0.05, showing an exponential fidelity decay.The number of entanglement dimensions accessible to certification has a maximum of Dent = 7 before decreasing again with growing system size.The numerical fits have been computed using data points with L ≥ 6 (nondotted lines) and contain an offset for the bounds.The data for r = 0.05 are not included in the figure for better visual clarity.The statistical error bars, especially on the true fidelity F (L, r), are barely visible.Both configurations are evaluated with an increased 2.5×10 4 position-space and 5×10 4 momentum-space samples.Figure 8.The dependence of the SE of the fidelity bound σ F on the lattice size L at a fixed number of samples.All mixing rates r ∈ {0, 0.05, 0.15} show similar initial growth of σ F before saturation.The combined data set is well described using an algebraic asymptotic growth model, showing little variation in σ F for lattices with L ≳ 10 fore not statistically prohibitive, opening up the possibility of feasibly preparing and certifying states with very high-dimensional entanglement.
In summary, state dephasing and lattice fluctuations have different impact signatures on both the true fidelity F and on our fidelity bound F .While the bound tightness is loosened by growing dephasing, with a linearly widening gap between F and F , it remains mostly tight in the presence of lattice-potential fluctuations.The statistical errors follow the expected Monte Carlo scaling ∝ 1/ √ N s ; very moderate sample numbers of ≈ 1×10 4 both in momentum and position space are sufficient to reduce the SEs to the subpercent range for all investigated lattice sizes.The bound is therefore robust with respect to typical noise sources and the entanglementcertification capability comes close to that of the true fidelity.

IV. MULTIPLE PARTICLES PER SPECIES
In the context of quantum simulation of condensed matter physics problems, the two-atom configuration discussed so far presents a somewhat unphysical low-density limit.Eventually, one would like to access the entanglement structure near half-filling, meaning atom number N = L/2, where true many-body effects emerge.
However, our method still relies on the measurement of coefficients of trigonometric basis functions, the number of which scales with the local Hilbert space size and thus exponentially in the particle number.Hence the true many-body regime stays out of reach for the scheme presented in this work.Nonetheless, studies of few-body cold-atom systems in the past decade have revealed that the few-body dynamics approache the many-body limit even at very moderate particle numbers [57,58].Fewbody systems are thus interesting candidates for quantum simulation and, in particular, entanglement certification, and give experimentalists capabilities beyond that of simpler two-particle systems such as entangled photon pairs.Recent success in the preparation and control of indistinguishable atom systems motivate this ansatz [59,60].Here, we want to extend our method to multiple particles in each of the two subsystems and present numerical simulations for up to N = 4 particles per species.

A. Theoretical considerations
Systems in which the number of atoms per species is increased to N > 1 can conveniently be described in a second quantization picture with different Fock modes.These modes are labeled by the spin of the particles and their lattice positions and are occupied by a given number of particles.For the fermionic atoms in the Fermi-Hubbard model, each lattice site can at most be populated by one atom per species due to Pauli exclusion.Hard-core bosons have the same exclusion rule, here enforced by strong repulsive on-site interactions.The resulting dimension of the local Hilbert space, i.e., the Hilbert space available to each species, which determines the maximal entanglement dimension, thus becomes A half-filling configuration gives the highest-possible entanglement dimension for a given lattice size, with with the normalization adapted to reflect the changed Hilbert-space size.In this notation |m 1 . . .m N ⟩ A/B designates the Fock state of species A or B, where sites m 1 . . .m N are occupied by one atom each.Here, the m i are in ascending order and are mutually different due to the aforementioned exclusion rules.
The general approach of bounding the fidelity to the MES to bound the Schmidt number remains unaltered.All experimental tools for single-atom and spin-resolved detection are still applicable but one has to take care to properly address the indistinguishability within the subspecies.While state populations can be extracted in a straightforward extension to the two-atom case, some changes have to be applied to access the coherences in Eq. ( 4) in the second quantization picture.Here, due to different commutation relations of fermions and bosons, our bound behaves differently for the two cases.In this work, we focus only on hard-core bosons and fermions, as their Hilbert spaces are identical and can thus be treated analogously.The detailed construction of the fidelity lower bound from multiparticle real-space and momentum correlation functions is presented in Appendix D. The crucial difference between fermions and bosons is the appearance of signs in the coherence terms stemming from the fermionic anticommutation relations.This diminishes the tightness of the estimate used in Eq. ( 16), even for pure states, and makes high-dimensional entanglement certification more challenging for fermionic systems than for hard-core bosons, as we show in Sec.IV B.
Lastly, we briefly address the scalability of the method toward larger particle numbers.Increasing the system size in terms of the number of atoms in the system requires significant computational resources, both for synthetic data generation and data processing.Furthermore, the necessary measurement statistics also increase for systems with higher atom counts.Our data processing is based on Monte Carlo techniques, which do not inherently scale with the dimension of the momentum space, i.e., the number of atoms in the system, but scaling can be introduced through the variance of the joint momentum distribution.We investigate these statistical scaling properties in Sec.IV C.

B. Numerical results
We simulate the ground state of N = 3 particles of both species in a lattice with L = 6 for both fermions and hard-core bosons.This setup enables a maximum entanglement dimension of D max ent = 6 3 = 20 [see Eq. ( 21)].To compare the behavior of this few-body system with that of two atoms, we repeat the interaction-strength sweep shown in Fig. 2. Our numerical data show that the fidelity to the MES in the strongly attractive regime is lower than in the two-atom case [see Figs.9(a) and 9(b)] but the behavior is otherwise qualitatively the same.The fidelity reduction is caused by a combination of same-site exclusion, which increases the distance between atoms of the same species, and finite-size effects, penalizing occupation of sites close to the edges.The combination of both effects leads to a very nonuniform distribution of dimer populations.However, as anticipated, the fidelity bound F shows a strong dependence on the underlying quantum statistics; for bosons, much higher and thus tighter fidelity bounds were achieved compared to fermions.This also leads to a large difference in terms of the certified entanglement dimension; in the pure case at U/J = −15, we certify D ent = 7 for fermions and D ent = 13 for hard-core bosons.The trend also carries We find good agreement with exponential decay for both the true fidelity F and our bound F for Jσv ≥ 0.01 disorder strengths.The simulation has been conducted at U/J = −12.All measurements have been simulated using 1×10 5 momentum-space and position-space samples each.
over to dephased states with r > 0, where we observe a stronger impact of dephasing than in the two-atom case.
In the fermionic case, at the strongest investigated dephasing of r = 0.15, no entanglement is witnessed in the ground state and only D ent = 2 is found for weaker dephasing of r = 0.05.The impact is less severe for bosons, where we find D ent = 3 for r = 0.15 and D ent = 13 for r = 0.05, respectively.These findings contrast with our data for a half-filling configuration for N = 2 atoms per species with L = 4, shown in Fig. 16 in Appendix E. Here, bosons and fermions show very similar results, with minor deviations only visible for pure states in the weakly attractive regime.Additionally, the effect of dephasing is much more comparable to our initial findings for N = 1 atom per species in Fig. 2. With higher numbers of particles and lattice sites present in the system, an increasing amount of coherences have to be subtracted using the bound in Eq. 16, explaining the difference in performance for different system sizes.Finally, we also investigate the case of N = 4 hardcore bosons per species on a lattice with L = 8.Using 3×10 5 samples in both position and momentum space, we are able to estimate F = 0.50 ± 0.06 at U/J = −15.Given D max ent = 8 4 = 70, this fidelity translates into a certified entanglement dimension with respect to a 1σ confidence interval of D ent = 31 (3σ confidence: D ent = 23).In accordance withe the above-described dephasing characteristics at r = 0.05, we see a strongly reduced fidelity bound of F = 0.10 ± 0.07, which gives D ent = 3 at the 1σ level.
Interestingly, the addition of disorder on the lattice reveals some key differences compared to the 1 + 1-atom case.Instead of the stretched exponential approach toward B 1 we find standard exponential decay of F with the bound decreasing significantly below the entanglementdetection threshold [cf.Fig. 9(c)].The matter wave function cannot converge to one localization center but is distributed across the entire lattice due to the abovediscussed exclusion rules.A nonzero number of states with unpaired atoms retain nonvanishing populations and connected coherences, which in turn induce tightness loss.Nonetheless, the initial transitional phase is comparatively short, with a good numerical fit agreement already for Jσ V ≥ 0.01.To summarize, the onset of localization effects is found for weaker disorder in few-body systems and quickly converges toward the sensible infinite-disorder ensemble, i.e., the perfect mixture of localized dimer states.However, state dephasing is the dominant effect, as even the strongest disorder strength Jσ v = 0.15 results in a deviation ∆ F := |F − F | ≈ 0.08, half of the difference caused by minor dephasing at r = 0.05.

C. Scaling of sampling complexity with N
To gain a better understanding of the complexity in terms of the required experimental runs N s , we con-Figure 10.A Log-linear plot of the fidelity bound SE σ F for N + N atoms in a lattice at half-filling (N = L/2).For coherence reconstruction, 3×10 5 position-and momentumspace samples each were used.The SE is well described using an exponential numerical fit.Because of the higher computational cost for 4 + 4 atoms, a smaller number of resamples was taken, leading to higher uncertainty in σ F .
ducted a scaling analysis of N + N atoms at half-filling for N ∈ {1, 2, 3, 4}.We have found that the SE of the fidelity bound σ F increases exponentially with N (see Fig. 10).From the fit coefficients one can extract an expected increase of σ F by a factor s N = 3.31 ± 0.13 for every additional atom pair introduced into the system.Using the earlier-confirmed (cf.Fig 4) σ F (N s ) ∼ 1/ √ N s relation, an increase of the sample size by a factor of s 2 N = 11.0±0.9 is necessary to keep statistical errors constant while increasing the particle number per species, N , by one.It shows that scaling deep into the many-body regime remains infeasible but configurations with a few atoms per species are realistically achievable.

V. MULTIPARTITE ENTANGLEMENT
In the previous section, we showed that when extended to few-body systems of multiple atoms per atomic species, our method still succeeds in the certification of high-dimensional entanglement, even for mixed states.Since only two spin states are populated, the system is fully described by the atom number N , the interspecies interaction strength U , and the tunneling strength J, so the same experimental toolbox can be used as for the case of one atom per species.When, instead, a higher number of spin states and thus entanglement parties is involved, a plethora of experimental and theoretical complications arise for entanglement detection but we can adapt the bound F to be able to certify genuine multipartite entanglement.In the following, we first formulate the theoretical framework needed for the classification of multipartite entanglement in this system and then describe a possible setup certifying high-dimensional tripartite en-tanglement.Finally, we present simulation results of entanglement certification for three atomic species in an optical lattice of L = 6 sites.

A. Multipartite-entanglement certification
While bipartite entanglement of pure states is fully developed theoretically, many questions are still open concerning the characterization and certification of multipartite entanglement.For states consisting of three entangled qubits, two sets of nonequivalent states sharing genuine tripartite entanglement have been identified, those equivalent under local operations and classical communication (LOCC) to the Greenberger-Horne-Zeilinger (GHZ) state and those LOCC equivalent to the W state [4,61,62].Because of this nonequivalence of entanglement, the Schmidt decomposition is no longer defined for general multipartite states.Different methods are therefore needed to obtain and describe the entanglement structure of multipartite quantum states.Numerous different approaches have been proposed to define canonical forms of tripartite and multipartite states with a minimal number of nonzero coefficients.However, to uniquely define any given quantum state through these methods, a number of parameters significantly higher than the local Hilbert-space dimension is required [47,[63][64][65].For some states, most notably also for generalizations of the GHZ state to higher local dimensions, it is still possible to define a generalized Schmidt decomposition, as every contribution to |ψ⟩ ABC combines orthogonal basis vectors |i⟩ for all three subsystems, assumed to have the same local Hilbert-space dimensions.No basis transformation can therefore reduce the number of terms used for the representation of Eq. ( 23) any further [4,66].One can now define a multipartiteentanglement dimension with the properties of an entanglement monotone in analogy to that of bipartite states [67], also with a maximum value of D max ent = dim H A .A generalized GHZ state with equal contributions on all sites given by therefore represents a suitable generalization of the twoatom MES [Eq.( 1)] as a reference state.It should be noted that this state is not maximally entangled in the sense that it has the maximum number of terms needed to be faithfully represented among all states with genuine tripartite entanglement but, rather, has the highest number of terms possible for it to also have a generalized Schmidt decomposition of the given form.(a) The dependence of the three different interaction-to-tunneling-strength ratios U12/J, U13/J, and U23/J on the external magnetic field (in gauss) based on Ref. [68] and gauged to fit experimental data published in Ref. [44].(b) Magnification of the narrow s-wave Feshbach resonance at 543 G [69,70].
The entanglement dimension of an experimental state ρ can be bounded by a set of fidelity thresholds B k to that reference state analogous to Eq. ( 3), opening up in principle the same certification route taken for bipartite entanglement.We prove these bounds in Appendix F. The algorithm given by Eqs. ( 8)-( 16) can be adapted straightforwardly to include three or more atomic species (for details, see Appendix G).

B. Experimental model
There are several different possible cold-atom implementations in which multipartite entanglement can be realized.Here, we choose a generalized Hubbard model with three distinct atomic species.The interaction strength between different spin states is usually regulated through the use of a magnetic Feshbach resonance [22,52].When a third spin state is added to the system, each of the three possible atom pairs is now governed by their individual interaction strengths U ij .To experimentally realize control over a mixture of three different spin states, an isotope with three overlapping Feshbach resonances connecting three low-energy eigenstates can be used.One possible choice is 6 Li, for which Feshbach resonances for the three lowest energy states at 690 G, 811 G, and 834 G are experimentally accessible and have been realized before [71][72][73].Since all three Feshbach resonances are magnetic, it is no longer possible to control the individual interaction strengths independently, but, rather, all three values U 12 , U 23 , and U 13 are tuned at the same time through shifts of the external magnetic field.For field strengths up to 527 G, all three scattering lengths are negative, delivering a broad regime of attrac-tive interactions between all atom species, as shown in Fig. 11.States close to the MES [Eq.( 24)] can thus be realized by preparing the Hubbard-model ground state in this regime.A three-particle extension to the Hubbardmodel Hamiltonian can be constructed as with σ, σ 1 , σ 2 ∈ {1, 2, 3} labeling the different hyperfine states [73].We base our numerical simulation of tripartite entangled systems on precise scattering lengths for 6 Li published in Ref. [68].From these measurements, we derive U/J values for different magnetic field strengths gauged to fit the interaction-strength data for U 13 reported in Ref. [44] [Fig.11(a)] to establish experimental comparability.This provides access to the interactionstrength triplet for a wide field-strength regime and thus enables one to study high-dimensional tripartite entanglement in Hubbard-model ground states.An alternative approach could be based on ultracold fermionic atoms in optical lattices with SU(N )-symmetric interactions [74,75].They have recently been shown to feature strong effective multibody interactions, making them a promising atomic platform for the preparation of multipartite entanglement in the future [76].

C. Numerical results
To assess the effect of the new intricate triplet structure of interaction strengths, we perform a sweep across the accessible range of magnetic field strengths B for three distinguishable atoms in the ground state of Eq. ( 25).The result is presented as a function of U 13 in Fig. 12(a).All presented true fidelities F (ρ, Ψ MES ) are again computed through exact diagonalization of the Hamiltonian.The signal found at U 13 ≈ −3.8J relates to a narrow s-wave Feshbach resonance at B = 523 G [magnified in Fig. 11(b)], which was earlier reported in Refs.[69,70].The observed fidelities are of similar magnitude as values reported for the 3 + 3 atom configuration in Figs.9(a) and 9(b) but with significantly higher fidelity bounds for the tripartite configuration.The impact of dephasing is of similar strength, as observed for simple two-atom configurations in Fig. 2.
The analysis of lattice disorder reveals differences compared to our results for two-species settings.At U 13 /J = −12, we see a stretched exponential decay in both F and F with lasting bound tightness, matching our results for two-atom configurations.However, we find a much steeper fidelity reduction and a smaller initial plateau.The fit yields stretch powers of c = 0.300±0.013for the true fidelity and c = 0.489 ± 0.011 for our fidelity bound.The three strongly attractive interaction strengths drive the atoms into triple-occupation states, which dominate the pure ground state at these values of U/J with 96.62% triple-occupation (trimers) and 3.35% The dependence of the disorder-averaged fidelity F on the normalized strength of the lattice-potential fluctuations, JσV , for the pure ground state at U13/J = −12 (U12/J = −3.67,U23/J = −6.66).The fidelities are adequately described through a stretched exponential fit approaching the B1 boundary in the strong-disorder limit.The two simulations used to create the displayed data sets each utilized 5×10 4 momentum-space and 1×10 4 position-space samples.The statistical errors caused by disorder averaging and limited statistics are relatively small.Data points in the dotted part of the fit have been excluded from fitting.
double-occupation (dimers) contributions.In a lattice with nonvanishing disorder, the wave function therefore localizes solely around a small number of lattice sites.Such bunching is prohibited for bipartite settings with indistinguishable particles, where Pauli exclusion enforces a maximum of two atoms per site (see Sec. IV), explaining the greater disorder susceptibility in the tripartite system.Decreasing the attraction strength diminishes the triplet bunching effect, leading to a less strong impact of disorder.This property comes at the cost of lower fidelity at very low disorder strengths, since the singleoccupation and double-occupation probabilities rise accordingly.
At vanishing disorder, robust certification of fourdimensional tripartite entanglement is possible and for very strong disorder, Jσ v ∼ 0.25, two-dimensional tripartite entanglement can still be confidently certified.We thus witness multipartite entanglement for an extended disorder regime.In the regime shown in Fig. 12(b), the contribution from lattice disorder to bound tightness is less than, or of the same order of magnitude as, the contribution from state dephasing.However, the reduction in true fidelity through disorder dominates all other error sources considered.

VI. GENERALIZATION TO OTHER REFERENCE STATES
Up to this point, we have shown the application of our method to ground states of Hubbard models with attractive interspecies interactions.We now develop generalizations to repulsive models of two or more atoms.As the entire process of measuring the g αβ coefficients is agnostic with regard to the measured state, the steps outlined in Eqs. ( 10)-( 14) can be applied in an identical manner, leaving the experimental procedure unchanged.However, the reference state Ψ ref must change and therefore one must extract different coherences.The employed scheme for deriving fidelity lower bounds can, in principle, be applied to any reference state.However, the bound tightness, especially in the presence of dephasing noise, will generally depend on the properties of the chosen reference state, leaving room for optimization in a given experimental scenario.

A. Two repulsively interacting atoms
A suitable reference state for the ground state of a repulsive Hubbard model of two atoms in a lattice of L sites may be given by an equal superposition of all nondimer states, which in turn means that the fidelity to that reference state is given by The procedure to extract coherences by subtracting bounds on all other noncontributing coherences presented in Eq. ( 16) can then be adapted to remove coherences not of the type of Eq. ( 27).This can be done without added complexity, as all coherences have the same weight in the fidelity and can be homogeneously extracted from contributing g αβ coefficients.This delivers a valid and accessible lower bound on F (ρ, Ψ ref ).However, since Ψ ref is not natively given in a Schmidt-decomposed form, one first has to compute the Schmidt decomposition i .In the case of L = 6, we find λ 1 = 5/6 and λ 2 = . . .= λ 6 = 1/30.This results in a high barrier of B 1 = 5/6 to detect entanglement at all, while the higher thresholds are equally spaced between B 1 and 1.This is, in essence, caused by the fact that our initial guess for a reference state is simply not that highly entangled, as can be seen through the entanglement entropy S(Ψ ref ) = 5/6 log(6/5) + 1/6 log(30) ≈ 0.719 compared to the MES used in the attractive case with S(Ψ MES ) = log(6) ≈ 1.792.
We can compensate this shortcoming by exploiting the additional structure of reference states of the form (26).As we show in Appendix H, uniform nondimer reference states always have an associated Schmidt basis vector ( The bound can then be derived as follows.First, one bounds the nondimer-nondimer contributions from below by taking the sum of all coefficients and subtracting the bounds of all other terms as originally shown in Eq. ( 16), where the second sum on the right-hand side includes dimer populations as well as all dimer coherences.Second, we have to bound the dimer-dimer terms.In principle, this can be done analogously to Eq. ( 29) and would amount to subtracting bounds for all nondimer terms, which for a repulsive model are much larger than dimerdimer terms.Even slight dephasing would cause a significant underestimation of these coherences and a corresponding loss of bound tightness.A more controlled approach is to bound the sum as, where we have replaced the dimer-dimer coherences with the negative of their upper bound.This introduces a small bias in the case of w d d ≥ 0, meaning that all states with nonvanishing dimer populations can no longer deliver a tight bound, but the loss of bound tightness from this term is largely independent of the level of dephasing and thus is much more stable.
While fidelity bounds can be derived analogously for arbitrary reference states, this susceptibility to dephasing renders reference states with no structure in their entanglement spectrum less suitable in practice.Coherences would appear in the fidelity with widely varying weights and would require a large number of bounds in the style of Eq. ( 29).Dephasing leads to underestimation of many of those bounds as described above and ultimately to the loss of any usable bound.
A similar argument can be made for the the dimernondimer coherences as shown below: Consequently, one can first perform the measurement scheme as originally introduced in Sec.II and then optimize the certified entanglement dimension by varying the value of λ 1 of |Ψ ′ ref ⟩.The results of this optimization procedure for the two-atom ground state with L = 6 are displayed in Fig. 13.The highly peaked entanglement spectrum of the "naive" initial guess |Ψ ref ⟩ offers the best fidelity but insurmountable entanglement thresholds.On the other hand, a uniform spectrum, i.e., a maximally entangled state, delivers ideal bounds but at the cost of loss in fidelity.We find the highest certified entanglement dimension of D ent = 4 at λ 1 ≈ 0.706.

B. Multiple atoms per species with repulsive interactions
Next, we also investigate a repulsive system of N atoms per spin state at half-filling (L = 2N ).Here, the adaptation of our method is much more straightforward than we saw before.In the attractive case, the Hubbard-model ground state was close to a superposition of states in Figure 14.The numerical results for the entanglementdimension certification of 3 + 3 hard-core bosons in a lattice of size L = 6 with repulsive interactions.The represented data mirror the above results for the same setting but with attractive interactions, in Fig. 9(b).The dotted lines again represent the infinite-measurement-statistics limit F∞ computed using exact coherences of ρ.All measurements were simulated using 1×10 5 momentum-space and position-space samples each.
which all atoms were bound in dimers, across all latticesite combinations [cf.Eq. ( 22)].Therefore, we used this MES as a reference state.In the repulsive regime, the ground state is close to a superposition of states with no dimers.At half-filling, for a given configuration of sites being occupied by species-A atoms, there is a unique configuration of species-B atoms realizing a dimer-free state, namely all B atoms occupying the sites not occupied by A atoms.This means that the standard choice of perfectly anticorrelated atom positions, is also a MES, equivalent to Eq. ( 22) up to a permutation of the species-B basis states and thus ideally suited for our entanglement-detection scheme.One therefore only has to extract coherences of perfectly anticorrelated states instead of correlated ones, while leaving the rest of the scheme unaltered.
We have simulated the extraction procedure for the ground state of a system of varying repulsive U/J with N = 3 hard-core bosons per species in a lattice with L = 6.The results are shown in Fig. 14.Depending on the level of dephasing, certification of up to D ent = 13 is feasible.Our data mirror previous data from our investigation for attractive systems in Fig. 9(b), where we have observed matching fidelities and bounds for exchanging U ↔ −U and replacing the attractive MES in Eq. (22) with the repulsive MES from Eq. ( 33), reminiscent of a particle-hole symmetry.This means that the method can be applied to Hubbard-model ground states over the entire range of interaction strengths, with particularly favorable properties in the half-filling case.But, also, multiatom scenarios away from half-filling can be treated.There, we do not have a unique particle-hole matching, but several contributions with nonvanishing configurations have to be considered in a reference state.For states close to half-filling the entanglement spectrum remains mostly flat but the lower the density in the lattice, the less informative is the knowledge of the positions of species A about species B. This causes a more strongly peaked entanglement spectrum, as we have observed in the case of two atoms on L = 6 sites.To demonstrate that certification is still possible we performed, as an example, a simulated application of the method for the ground state of N = 2 hard-core bosons per species in L = 6 lattice sites at U/J = −12, using 5×10 4 samples in momentum and position space each.The naive uniform superposition of all nondimer states as 6, so D ent = 4 can be certified at 3σ confidence.Also, here one could conceive of a scheme to design better-suited reference states in the spirit of Fig. 13, which we leave for future investigations.

A. Summary
We have presented a new method to bound the fidelity of few-body states of ultracold-atom systems to a highly entangled state.High fidelity indicates the presence of high-dimensional entanglement in the experimental state and can be used to bound entanglement quantifiers such as the entanglement dimension or the entanglement of formation.We have constructed lower bounds on the fidelity that are measurable in systems of ultracold atoms in optical lattices utilizing only position-and momentumspace measurements.A detailed study of the statistical significance and tightness of these bounds under realistic assumptions about experimental measurement conditions and noise sources indicates manageable experimental and statistical requirements.Interestingly, states that are highly mixed due to lattice-potential fluctuations retain their bound tightness to a high degree, allowing one to observe the disorder-induced reduction of ground-state entanglement.Generic white noise has been identified to cause linear decline of the tightness of our fidelity bound, while finite temperature has a comparably mild impact on bound tightness.We have generalized this method to certify entanglement in multipartite systems and con-figurations with several atoms per spin state, requiring alterations to the coherence-extraction framework to account for partial indistinguishability.In these settings, we have demonstrated the feasibility of certifying up to D ent = 31 entanglement dimensions for 4 + 4 hard-core bosons and up to D ent = 4 of genuine tripartite entanglement.Furthermore, by using reference states beyond the canonical MES, we have demonstrated the wide applicability of our method to quantum simulation experiments with itinerant particles in lattice geometries.

B. Literature context
Our work should be considered in the context of research lines focusing on efficient state tomography schemes or, leaving tomography out as an intermediate step, direct entanglement-detection methods.Here, we briefly review these research lines, commenting on their strengths and weaknesses compared to our method, without any claim of completeness.
Full quantum state tomography in a bipartite system generally requires a number of different measurement bases that scale quadratically in the local Hilbert-space dimension [77].This limits its applicability to very small system sizes [78], despite significant advances in the efficiency of maximum-likelihood estimation [79] and Bayesian tomography methods [80].A more economic scaling of experimental cost can be reached by restricting the state space in which the reference state is being searched.Examples for such approaches are compressed sensing tomography [81][82][83], assuming that the prepared state has reduced rank, and methods using variational ansatz functions, such as neural-network quantum state tomography [84][85][86][87] or matrix-product state tomography [88][89][90], which restricts its search space to weakly entangled states-operating exactly in the opposite regime to the one targeted in this work.The drawback of this class of methods is that restricting the state space necessarily leads to bias, as it is generally not known whether the experimentally prepared state lies in the class of states representable by the ansatz.
For extracting properties of the entanglement spectrum, it is often not necessary to fully reconstruct the quantum state.For example, if the global state of the system can be assumed to be pure, the entanglement spectrum can already be extracted from the state ρA of a subsystem.A variational approach to determine the entanglement Hamiltonian, i.e., the logarithm of the reduced density matrix, has recently been demonstrated experimentally [35][36][37].Another notable approach is the use of random measurements to detect entanglement [91] within the framework of shadow tomography [92].This framework can be applied to extract Schmidt-number witnesses by probing correlation matrices [38,39].However, the method requires large sample sizes and the implementation of Haar-random unitary operators, an open challenge for systems of itinerant particles.If one is only interested in Rényi entanglement entropies, methods using multiple copies of the quantum state can be employed [93].A related protocol uses ancillary particles to measure the entanglement spectrum directly for cold latticeconfined bosons [40].While this method is very elegant, it poses stringent requirements on experimental capabilities.
The approach pursued in our work relies on measurable lower bounds on the fidelity to a highly entangled reference state for probing the entanglement dimension.A number of works have studied efficient methods for estimating fidelity, or at least bounding it, ranging from correlation-measurement-based approaches [41,94,95] to variational methods [96] and random Pauli-string measurements [97,98].Often, these schemes are tailored to a specific experimental system, such as entangled photon pairs in the case of Refs.[41,95], where the capability to measure in a pair of MUBs is exploited.This makes it difficult to apply these methods to other platforms, where these capabilities are not given.The strength of our proposal lies in the development of an entanglementcertification scheme that relies on techniques readily available to cold-atom experiments and is generally applicable to bipartite and multipartite scenarios realizable with this versatile quantum simulation platform.

C. Scalability
The term quantum simulation often entails the notion of scalability to system sizes that are beyond the reach of classical simulation methods, i.e., reaching the regime of quantum advantage.Here, we summarize our findings on the scaling of both experimental and computational cost of our method as a function of lattice size L and particle number per species N .
We have found an algebraic saturation of the SE of our fidelity bound for growing lattice sizes for N = 1.Statistical requirements for faithful entanglement certification thus remain approximately constant for an extended regime of lattice sizes, which opens up one pathway to prepare highly entangled states in large lattices.However, increasing the number of atoms (keeping the density constant) in bipartite configurations results in an exponential increase of the SE, which consequently necessitates an increase of samples taken by nearly one order of magnitude to add an additional pair of atoms to the system.We note, however, that experimental sampling of momenta is achieved through fluorescence imaging, where all momenta of one atomic species are captured in a single image.Consequently, there is no inherent connection between the sampling rate and the system size.This allows comparably fast sample production in systems with several atoms compared to the creation of such samples by numerical simulations, where the computational complexity is linked to the number of atoms.
The data-processing routine used in this work consists of several steps: projection of the sampled momentum distribution onto modes of the momentum-basis expansion, basis change via formal matrix inversion of the matrix Q to correct for nonorthogonality, and coherence extraction.All these steps have exponential computational complexity scaling in N , which also prohibits application to genuine many-body systems.By contrast, the local Hilbert-space size, and thus the processing complexity, is only polynomial in the lattice size L.

D. Outlook
The detection scheme proposed here can be generalized and extended in various ways.First, the momentumspace measurement, achieved by completely switching off the lattice potential, operates in the continuous domain.One could also envision only tuning the interparticle interaction strength to zero and allowing the particles of each species to undergo a tunneling evolution in the lattice before they are imaged.This would correspond to a measurement in a discrete basis complementary to the in situ measurement, similarly giving access to coherences as the current scheme but avoiding the step of projecting measured data onto a function basis in continuous space.Furthermore, it would also relax the resolution requirements in momentum space, making the method accessible to a even broader range of contemporary experimental setups.Second, entanglement-dimension witnesses based on measurements in two complementary bases may be developed analogously for other quantum simulation platforms.Examples are trapped ions, superconducting qubits, or Rydberg atoms, realizing spin systems, where the entangled subsystems consist of multiple spins, or qubits, with native local unitary transformations available to each specific platform.Furthermore, one could include a small number of additional measurement bases to give further constraints on state coherences, improving bound tightness, especially for higher atom numbers.Finally, while investigating the impact of experimental imperfections on our ability to certify entanglement in realistic settings, we have found distinctly different signatures for pure state dephasing and lattice disorder.This implies that this bound could be used as a probe of disorder and localization in the prepared state.More generally, we would like to apply the developed method to more quantum states of interest, beyond Hubbard All numerical results presented in this paper require a number of processing steps, ranging from synthetic sample generation to nonorthogonality corrections and coherence extraction, that come with computational complexity scaling exponentially with the size of the system.In the following we will briefly describe the numerical methods we used for optimizing the performance of classical processing steps, which have been crucial for reaching the largest reported system sizes.
The key technique for synthetic data generation is the sampling process for high-dimensional probability distributions.Computing the full distribution on a grid becomes prohibitively expensive, so a Monte Carlo type algorithm must be employed instead.Realizing that the momentum integral over Eq. (10a) separates for each term, one can utilize a so called ancestral sampling procedure [99]: Since the integrals can be split up, one can integrate out all but one of the momenta to obtain the marginal p(k 1 ).After sampling k 1 from that distribution, one can fix k 1 and integrate out the rest, now leaving k 2 open to obtain the conditional probability distribution p(k 2 |k 1 ).This scheme can be repeated until all momenta are fixed and a complete sample is generated.Replacing d-dimensional integrals with the product of d 1D integrals, that can be computed beforehand, greatly increases accessible system sizes.The remaining 1-dimensional integrals I(δ) are of form We find that | w(k)| can be well approximated through a Gaussian g(µ = 0, σ), which also agrees with experimental findings [44].Replacing the Wannier envelope yields an expression of the form I(δ) ∼ exp −(dσδ/2) 2 .With this, the expression for p(k 1 ) becomes a sum over terms weighted by the integral over the remaining l − 1 momenta, I(δ 2 , . . ., δ l ) ∼ exp −(dσ/2) 2 l i=2 δ 2 i .As these are exponentially small in l i=2 δ 2 i , we can define a cutoff δ c and neglect all terms for which l i=2 δ 2 i > δ c .A similar technique can be used to simplify the computation of the matrix elements of the basis overlap matrix Q [see Eq. ( 12)].The matrix is needed to correct the measured coefficients c αδ for overlap with different nonorthogonal basis elements.In the two-atom case, each element is given by By similar manipulations of the trigonometric functions under the integral, one can obtain the following factorized form Using the fact that f (γ) = f (−γ), one only needs to evaluate f (γ), where γ ∈ {0, . . ., 2(L − 1)}.
The formal inversion ⃗ G = Q −1 ⃗ C can be efficiently approached by exploiting that Q is hermitian and positivedefinite and using Cholesky decomposition, Q = LL † [100], which gives the lower-diagonal matrix L acting as a preconditioner for Q −1 .For large Q, saving the dense L can become too costly, so that the iterative conjugate-gradient algorithm [101] becomes the more practical solution.Both algorithms are available within the numpy and scipy scientific computing libraries in Python [102,103], which we have used for our numerical simulations presented in this work.

Appendix B: Fitting parameters
In Table I we list the numerical fit models and fit parameters for all conducted fits appearing in this paper.L dependence order  The bound rapidly starts to loose tightness for βJ ≲ 0.5.As before, the data has been taken using 1×10 4 position-space and 2.5×10 4 momentum-space samples.

Appendix C: Bound performance for thermal states
To ascertain the susceptibility of the derived fidelity bounds to thermal excitation, we have conducted additional simulations for thermal states, ρT ∼ exp(−βH), of two attractively interacting distinguishable atoms in a lattice of size L = 6.The resulting fidelities in dependence of the normalized inverse temperature, βJ, are displayed in Fig. 15.In contrast to our results using white noise as a generic decoherence model, presented in Fig. 5, we observe no significant loss of bound tightness for a broad temperature range of βJ ≳ 0.5.This value of βJ translates to a ground-state fraction of ≈ 19.3% and an ensemble purity of ≈ 0.164.For even higher temperatures (smaller βJ) the bound starts to deviate from F , as shown in the inset, but this is inconsequential for entanglement-dimension certification, as both F and F are lower than B 1 and entanglement can no longer be witnessed.
These results indicate that the proposed bound is resilient against dissipation through coupling to a finite temperature bath.Thus, the method can be applied in an experimental setup where one cools directly into the Hubbard-model ground state, ending up at some finite temperature.We note that ground states may also be prepared by adiabatic deformations of the optical potential, starting with a localized dimer.Realistic modeling of the preparation process will depend on the concrete experimental setup and is left for future investigation.

Appendix D: Details on indistinguishable atom bipartite entanglement certification
Here, we want to give a full derivation of our fidelity bound for systems of multiple indistinguishable particles per species.We begin by giving the field operators Ψ↑ (k), Ψ↓ (k), in terms of their lattice-site creation and annihilation operators, where we have introduced shortened notation âj := ĉj,↑ and bj := ĉj,↓ .Operators within one bosonic (fermionic) species are subject to commutation (anticommutation) relations.Both species are distinguishable via their spin degree of freedom, meaning that operators from different species always commute.Formally, these statements then read where we have used [•] ± to denote anticommutator and commutator, respectively, and ω(k) is the Fourier transform of the Wannier envelope.We can then rewrite the 2N -atom momentum correlation function using the field operator form given in Eq. (D.1a) and obtain in analogy to the two-atom result in Eq. (8).Note that we have used the normal ordered correlation function ⟨: n↑ (k 1 )n ↑ (k 2 )n ↓ (k 3 )n ↓ (k 4 ) :⟩, as it correctly represents single-atom-resolved measurements in momentum space [105,106].This subtle distinction was not necessary for fully distinguishable particles, as the field operators Ψσ commute with each other.By inserting the coefficient expansion of ρ in the Fock basis and exploiting the relations in Eq. (D.1d), the expectation values in Eq. (D.2) become where we split the expectation value into a product between the two subsystems, or species.They can be evaluated independently using the commutation (anticommutation) relations from Eqs. (D.1b)&(D.1c),resulting in Inserting the above results from Eqs. (D.3) and (D.4) into Eq.(D.2) gives a rather long expression.Therefore, for the moment, we leave out all terms connected solely to one of the two subsystems as shown in Eq. (D.5a) and reintegrate them later once all terms relating to the remaining subsystem have been sufficiently simplified.First, we resolve all open δ terms to eliminate {p, q, p ′ , q ′ } from the summation, which produces a sum of phase factors, each with a different permutation of site indices appearing the in bra and ket in Eq. (D.5a).The sign of the different terms is determined by the underlying quantum statistics.Treating this sum of complex phase terms within the brackets in Eq. (D.5b) (in combination with their complex conjugate counterparts) as the new basis functions would give the desired combinations of coherences as weights.However, we find that these functions are linearly dependent, so that Q is generally rank deficient and thus cannot be inverted.Unambiguous reconstruction is therefore not possible.We overcome this issue by reorganizing terms.In a first step, we relabel all site indices such that all phases are of the same form, noted below the complex phases in Eq. (D.6a).One has to take care to properly transport the conditions k < l and k ′ < l ′ , which implement a second quantization picture, by introducing the corresponding Heaviside step functions θ(x) in Eq. (D.6b).Using the symmetry (antisymmetry) of the density matrix ρ under particle exchange, one can rejoin all four terms into one term but without any restrictions regarding an ordering of {k, l} or {k ′ , l ′ }, as seen in Eq. (D.6c).Finally, we replace two of the summation variables by the differences of the index pairs ∆k = k − k ′ and ∆l = l − l ′ to group together phases with the same factors appearing in the exponent Repeating above-described procedure for the second subsystem yields Eq. (D.7).This grouping of complex phases in combination with the Wannier envelope make up a complete basis and can thus be used to unambiguously reconstruct the corresponding basis weights, (D.7) The summation over the differences ∆k and ∆l also allows for configuration where k > l, terms that were previously excluded in Eq. (D.6a).This swapping of indices is equivalent to a particle exchange, which is accompanied by an additional minus sign for fermions in the corresponding matrix elements in ρ.The alternating signs in Eq. (D.5b) have effectively been shifted into the definition of the density matrix ρ.This means that unlike the simple two-atom case, some coherences inherently acquire a sign here, which can lead to "destructive interference" between coherences.The direct consequence is a loss in bound tightness for fermions, as observed in Sec.IV B.
Appendix E: 2+2 Atoms Figure 16 shows data for the case of N = 2 atoms per species in analogy to Fig. 9.We observe that the loss of tightness for fermions is far less pronounced than in the case of N = 3 atoms per species.The subtle differences between Fermi-Dirac and Bose-Einstein statistics can be most notably observed for very weak interactions strengths and pure states, where our bound underestimates the true state fidelity only for fermions.F on the interaction-to-tunneling-strength ratio U/J for pure and dephased states with (a) fermions and (b) hard-core bosons.The dotted line represents the infinite-measurementstatistics limit F∞ computed using exact coherences of ρ.Differences between the two plots are most discernible for pure states at U/J ∼ 0. All measurements were simulated using 5×10 4 momentum-space and position-space samples each.

Appendix F: Tripartite entanglement-dimension bounds
Here, we extend the concept of entanglementdimension bounds B k to tripartite reference states with generalized Schmidt decomposition, as given in Eq. ( 23), in close analogy to original work for bipartite states given in Refs.[42,107].The general idea is again to give bounds on the maximal fidelity between some generalized reference state |Ψ⟩ L = L i=1 λ i |iii⟩ and some state ρk with generalized Schmidt number k.This comparison can be made in a sensible way, as |Ψ⟩ L is already given in a form similar to the Schmidt decomposition of bipartite systems.All contributions are combinations of orthogonal basis states |i⟩ on the three subsystems and only appear once each.No unitary basis transformation therefore can reduce the number of states appearing in |Ψ⟩ L , giving it the same role as the bipartite entanglement dimension.It is not necessary to consider general mixed states, as convexity of the fidelity guaranties that fidelity is maximized through a pure state, so we restrict the proof to pure states only [42].The highest possible fidelity between the reference and a pure state as used in Sec.V.This resemblance is directly related to the restriction to Schmidt-decomposable states as reference states.Multipartite states in general cannot be brought into a form where each subsystem basis vector appears only once through some basis transformation.
Therefore, this technique can never be expected to be able to detect all terms for a generic multipartite quantum state but, at most, the minimum of all local Hilbertspace dimensions.
Appendix G: Details on multipartite-entanglement certification Extending the original scheme for entanglement certification to multipartite entanglement is straightforward but tedious.Here, we briefly want to give a starting point of how this extension is derived and present the final bound Fcoh .Like before, we decompose the momentum correlation function of three atoms ⟨n 1 (k 1 )n 2 (k 2 )n 3 (k 3 )⟩, in terms of coherences and consider phases picked up due to the Fourier transformation.This results in This is necessary to be able to do the full reconstruction of the momentum correlation function, since the true coefficients g αβγ have to be obtained from the full distribution of measured coefficients c αβγ first.The redefined set M for three atomic species is given in Eq. (G.3c).
All remaining steps outlined in Eqs. ( 11) to ( 15) can be adapted analogously, such that one arrives at the final result for the bound of the coherence contributions as follows: Fcoh (ρ, Ψ MES ) = |ij⟩ . (H.3) Consequently, we can split up our reference state as where we have included all remaining contributions in a normalized state, |ii⟩ , (H.5) and λ ′ = 1/ √ L. Note that both |λ 1 ⟩ and |λ ′ ⟩ are symmetric under lattice-site exchange, such that any superposition of these states will have the same symmetry.This means that one can define the one-parameter family of reference states by varying the relative weight between them.Finally, we compute the weights w of coherences in the fidelity, introduced in Eq. ( 28), as a function of the parameter λ.Using the above-discussed symmetry, we know that all dimer-dimer terms ⟨ii| ρ |jj⟩ must contribute equally, giving them a shared weight w . (H.8c) It is clear from Eq. (H.8a) that, for all L ≥ 2, one has w nd nd ≥ 0, since λ ∈ [0, 1].

Figure 1 .
Figure 1.(a) The position-space correlation function ⟨n ↑ (x1)n ↓ (x2)⟩ of the two-particle attractive Hubbard-model ground state for L = 6 lattice sites with lattice spacing d at U/J = −12.(b) A graphical representation of both particles occupying the same lattice site (top) or adjacent lattice sites (bottom) with the respective signals in (a).(c) The momentum correlation function ⟨n ↑ (k1)n ↓ (k2)⟩ corresponding to the position correlation function of (a).All values smaller than 1×10 −5 in both (a) and (c) have been masked.
(a) (for details on the model and numerical implementation, see Sec.III).Each grid point represents a twoparticle state contributing to ρ.The signals on the diagonal represent dimer population probabilities, whereas off-diagonal elements correspond to configurations with atoms on different sites [Fig.1(b)].The wave-function envelope is determined by the on-site Wannier basis of the lattice and depends on the lattice depth V 0 and site spacing d.Since L m=1 ⟨mm| ρ |mm⟩ /L ≤ 1/L, it is clear that the populations contribute at most ∝ 1/L to F .Their impact therefore becomes negligible compared to coherences for large systems.

Figure 3 .
Figure 3.The lattice-potential configurations in units of the recoil energy ER (left) and the respective dimer-occupation probability distributions for the ground state (right).The dashed lines indicate the potential baseline depth V0 (orange) and the uniform probability distribution of the MES ΨMES (red).In all three cases, the ground state has the maximum entanglement dimension of Dent = 6.(a) An even unaltered lattice potential.The dimer population is heavily centered on the central lattice sites.One can certify up to Dent = 5.(b) A lattice with increased potential depth at the outlying sites, resulting in a uniform distribution among all lattice sites and full certification of Dent = 6.(c) A lattice with potential fluctuations ∆E ∼ N (0, 0.08J) on each lattice site.The dimer population shows strong localization and is far away from ΨMES.Nevertheless, Dent = 4 can be certified.

Figure 4 .
Figure 4. (a) The fidelity bound F as a function of the number of momentum-space samples for mixing strengths r ∈ {0, 0.05, 0.15}.The gap between the true fidelity F (solid lines) and the fidelity-bound average ⟨ F ⟩ (dash-dotted lines) increases with growing impurity.Each marker represents one sampling realization.(b) The dependence of the SE of the fidelity bound σ F on the number of momentum-space samples Ns for mixing strengths from (a).(c) The linear regression of log-log-represented σ F data for r = 0 from (b), demonstrating power-law scaling.

Figure 5 .
Figure 5.The systematic linear scaling of both F and F with the state-mixing parameter r.Since the true fidelity F is computed exactly, linear regression errors [O(1×10 −16 )] are solely caused by machine precision and are omitted here.The error of F is of a statistical nature and is barely visible.

Figure 6 .
Figure 6.The disorder ensemble average of the fidelity F as a function of the strength (standard deviation JσV ) of lattice depth fluctuations.After an initial transitional phase, the fidelity and the bound follow a stretched exponential decay for potential fluctuations with JσV ≳ 0.07.The data points in the dotted part of the line are excluded from the fit, indicating deviating behavior for very weak fluctuations due to finite-size effects.

Figure 7 .
Figure7.The lattice-size dependence of the fidelity.(a) The scaling of the fidelity F of the ground state of a flat optical lattice as a function of the number of lattice sites L and the mixing parameter r.The numerical data, including statistical noise, fit well to the asymptotic algebraic behavior.The certifiable entanglement dimension continues to grow with increasing lattice size, as the fidelity asymptotically approaches its infinite-system-size value.(b) A log-linear plot of the fidelity for a disordered lattice of size L with fixed disorder strength JσV = 0.05, showing an exponential fidelity decay.The number of entanglement dimensions accessible to certification has a maximum of Dent = 7 before decreasing again with growing system size.The numerical fits have been computed using data points with L ≥ 6 (nondotted lines) and contain an offset for the bounds.The data for r = 0.05 are not included in the figure for better visual clarity.The statistical error bars, especially on the true fidelity F (L, r), are barely visible.Both configurations are evaluated with an increased 2.5×10 4 position-space and 5×10 4 momentum-space samples.

Figure 9 .
Figure 9.The numerical results for the entanglementdimension certification of 3 + 3 indistinguishable atoms in a lattice with L = 6.Scaling of the fidelity F and the fidelity bound F for different interaction-to-tunneling-strength ratios U/J for pure (r = 0) and dephased (r ∈ {0.05, 0.15}) states with (a) fermions and (b) hard-core bosons.The dotted line represents the infinite-measurement-statistics limit F∞ computed using exact coherences of ρ.(c) The scaling of the disorder ensemble average F as a function of the normalized optical-lattice depth fluctuation Jσv for the pure ground state of hard-core bosons.We find good agreement with exponential decay for both the true fidelity F and our bound F for Jσv ≥ 0.01 disorder strengths.The simulation has been conducted at U/J = −12.All measurements have been simulated using 1×10 5 momentum-space and position-space samples each.

Figure 11 .
Figure 11.(a)The dependence of the three different interaction-to-tunneling-strength ratios U12/J, U13/J, and U23/J on the external magnetic field (in gauss) based on Ref.[68] and gauged to fit experimental data published in Ref.[44].(b) Magnification of the narrow s-wave Feshbach resonance at 543 G[69,70].

Figure 12 .
Figure 12.The numerical results for entanglement-dimension certification of a tripartite-state configuration in a lattice with L = 6.(a) The scaling of the fidelity F and the fidelity bound F for different interaction-to-tunneling-strength ratios U/J for pure (r = 0) and dephased (r ∈ {0.05, 0.15}) states.(b) The dependence of the disorder-averaged fidelity F on the normalized strength of the lattice-potential fluctuations, JσV , for the pure ground state at U13/J = −12 (U12/J = −3.67,U23/J = −6.66).The fidelities are adequately described through a stretched exponential fit approaching the B1 boundary in the strong-disorder limit.The two simulations used to create the displayed data sets each utilized 5×10 4 momentum-space and 1×10 4 position-space samples.The statistical errors caused by disorder averaging and limited statistics are relatively small.Data points in the dotted part of the fit have been excluded from fitting.
j=1 |ij⟩, an equal superposition of all states in the entire Hilbert space.If one now defines a new reference state |Ψ ′ ref ⟩ in terms of the same Schmidt basis but varies the value of λ 1 and uniformly adapts the remaining Schmidt coefficients to preserve normalization, one obtains a family of highly entangled states symmetric under lattice-site exchange.This is important, as the weight w of coherences now only depends on whether they are of type dimer-dimer (⟨ii| ρ |jj⟩, w d d ), dimer-nondimer (⟨ii| ρ |jk⟩, w nd d ), or nondimer-nondimer (⟨ij| ρ |kl⟩, w nd nd ≥ 0), independent of the specific lattice sites.This makes their extraction much simpler and the process more resilient against dephasing effects, as we show below.Explicit expressions for the weights w are also given in Appendix H.The fidelity for a generic reference state from that family then reads F (ρ, Ψ ′ ref ) ρ |jk⟩ + ⟨jk| ρ |ii⟩ .

Figure 13 .
Figure 13.entanglement spectra of the adapted MES with equally spaced B k thresholds in (a) and of the referencestate choice with maximal fidelity in (b).(c) By varying the first coefficient λ1 of Ψ ′ref in postprocessing, we find the optimal Dent = 4 at λ1 ≈ 0.706 for the ground state for the two atoms at U/J = 30.The measurement was simulated using 5×10 4 shots in position and momentum space.The blue-shaded area represents the 1σ confidence interval of the bound.

1 LF
Fig. Description Fitting Model Fitting Parameters a b c 4(c) SE Ns dependence σ F (Ns) = aN b s

1 LFFigure
Figure Numerical study of a thermal of two atoms at U/J = −12 in a lattice with L = 6 sites.(a) Fidelity F (ρT , ΨMES) in dependence of the normalized inverse ensemble temperature βJ on the left axis.On the right, ensemble ground-state participation amplitude.(b) Zoom in on values of βJ ≤ 1.The bound rapidly starts to loose tightness for βJ ≲ 0.5.As before, the data has been taken using 1×10 4 position-space and 2.5×10 4 momentum-space samples.

Figure 16 .
Figure 16.Numerical results for entanglement-dimension certification of 2 + 2 indistinguishable atoms in a lattice with L = 4. Dependence of the fidelity F and the fidelity bound F on the interaction-to-tunneling-strength ratio U/J for pure and dephased states with (a) fermions and (b) hard-core bosons.The dotted line represents the infinite-measurementstatistics limit F∞ computed using exact coherences of ρ.Differences between the two plots are most discernible for pure states at U/J ∼ 0. All measurements were simulated using 5×10 4 momentum-space and position-space samples each.