Measuring Non-local Brane Order with Error-corrected Parity Snapshots

Exotic quantum many-body states, such as Haldane and spin liquid phases, can exhibit remarkable features like fractional excitations and non-abelian statistics and offer new understandings of quantum entanglement in many-body quantum systems. These phases are classified by non-local correlators that can be directly measured in atomic analog quantum simulating platforms, such as optical lattices and Rydberg atom arrays. However, characterizing these phases in large systems is experimentally challenging because they are sensitive to local errors like atom loss, which suppress its signals exponentially. Additionally, protocols for systematically identifying and mitigating uncorrelated errors in analog quantum simulators are lacking. Here, we address these challenges by developing an error correction method for large-scale neutral atom quantum simulators using optical lattices. Our error correction method can distinguish correlated particle-hole pairs from uncorrelated holes in the Mott insulator. After removing the uncorrelated errors, we observe a dramatic improvement in the non-local parity correlator and find the perimeter scaling law. Furthermore, the error model provides a statistical estimation of fluctuations in site occupation, from which we measure the generalized brane correlator and confirm that it can be an order parameter for Mott insulators in two dimensions. Our work provides a promising avenue for investigating and characterizing exotic phases of matters in large-scale quantum simulators.


INTRODUCTION
Conventional phases of matter can be characterized by measuring local order parameters, which represent the degree of symmetry breaking.However, it has become clear that the concept of "non-local order parameters (NLOs)" is crucial to distinguish different types of exotic quantum orders [1][2][3][4][5].The Haldane chain [6] is a paradigmatic example, where a string-type non-local correlator [2] can reveal the hidden quantum phase.Due to the direct accessibility of atomic distribution with a single-site resolution [7,8], ultracold atoms in optical lattices have provided an unprecedented opportunity to study NLOs.For example, the unity filling Mott insulating (MI) phase has been a testbed for studying NLOs [9,10] because the phase hosts bounded particle-hole pairs as virtual excitations on top of the constant density distribution due to a finite tunneling strength.The string correlator has been measured in a one-dimensional Mott insulating phase [11], and recent experiments using Fermi gases have applied the stringtype correlator to reveal hidden antiferromagnetic correlations in doped Fermi-Hubbard chains [12] and to probe the Haldane phase in Fermi-Hubbard ladder [13].Efforts have been made to extend the NLO to higher dimensions, with the recent suggestion of generalized brane correlators in two-dimensional Hubbard models [14,15].Moreover, Wilson loops have been exploited to probe the topological Z 2 spin liquid phase in Rydberg atom arrays [16,17].
However, these non-local correlators can be easily destroyed by experimental imperfections, such as detection atom loss, limiting its practical usage.When measuring the brane parity correlator for N sites, as an example, each site may experience a small loss (error) rate η 1 during imaging process.The atom loss can occur in all lattice sites independently, so this can exponentially suppress the fidelity F of the N-site parity measurements F ∼ (1 − η) N .Hence, to reliably evaluate the NLO in a large-scale quantum simulator, it is essential to alleviate the effect of the incoherent error as much as possible.In literature, the systematic method to identify and reduce the effect of various types of errors from the measurement data is collectively called as the error mitigation protocol [18].Unfortunately, there has been little progress [19,20] on such protocols in ultracold atom simulators.
Here, we present a new error correction (EC) method for atomic quantum simulator and demonstrate its efficacy by measuring non-local order parameters in the two-dimensional Bose-Hubbard (BH) model.The EC method is based on mapping the parity snapshot data of the Bose-Hubbard system to the spin configurations of a 2D Ising error model, enabling the identification of the correlated particle-hole pairs and remove uncorrelated holes from the snapshot.The brane parity correlators with error-corrected snapshots then can successfully distinguish the Mott insulator and superfluid phases in a large-scale system containing more than 100 lattice sites.We further find that the brane parity correlator satisfies the expected perimeter scaling laws.Moreover, we are able to infer the fluctuations in the site occupation number in the MI because the EC method can assign the particle-hole pairs from the parity snapshots.This enables us to evaluate the generalized brane correlator, and confirm recent predictions that it can serve as an order parameter for the two-dimensional Mott insulators [15].

THE BOSE-HUBBARD MODEL
As an experimental platform, we employ ultracold 7 Li atoms in a square optical lattice [Fig.1(a)] to realize a twodimensional BH model [23,24]: The bosonic creation (annihilation) operator at lattice site i is b is the site occupation number, J is the tunneling strength, U is the on-site interaction energy, and ε i is the energy offset from harmonic confinement in the lattice beam.Moreover, we employ a high-resolution fluorescence imaging system [25] and detect the number parity P i = (−1) n i at lattice site i.Notably, we can prepare a large-sized unity filling Mott insulator (40 sites diameter with more than 1000 atoms) by tuning the scattering length of the atoms using a Feshbach resonance.In this study, we focus on the central area of the D c = 20×20 lattice sites to minimize trap inhomogeneity.For the brane correlators, we take all the possible L × L domains within D c [21].
However, the detection of the parity and preparation of MI states are not perfect because of the various errors, such as particle loss during the imaging process or free holes generated from the thermal fluctuations.Here, we attempt to simulate the ground state of MI, the finite temperature is also the source of errors.We will collectively call these holes as "uncorrelated errors" because they occur independently at each site.It contrasts with the correlated parity flips due to the virtual particle-hole excitations over the MI.We estimate the rate η of such uncorrelated errors to be ∼ 3% [21].Despite being small, this can dramatically diminish the ability to measure multi-point correlators in MI including the brane parity correlator O D = ∏ i∈D P i , where ... refers to ensemble average.For example, see Fig. with EC (blue) and without (lightblue) the EC method.The experimental data agree well with the QMC simulation.They can compare with the QMC simulation with no uncorrelated error η = 0 at low temperature T = 1/32U (darkblue dotted line).The gray dotted line represents the quantum critical point in 2D, estimated by the QMC simulation (J/U) c = 0.059 [22].Each data point is obtained over 40 independent experimental runs, and the shaded region denotes the maximum standard error (see supplemental material for details [21]).ity correlator signal only for a small system size with L = 6 [Fig.2(b) inset, open circle].Hence, it is imperative to remove the effect of the uncorrelated local errors in experiments to observe the non-local correlators in noisy environments.

ERROR CORRECTION METHOD
Our EC method is designed to circumvent the difficulty efficiently and allow the calculation of non-local brane correlators [14,15], which can accurately distinguish the two phases of the BH model.It is based on mapping the parity snapshots of the BH model to the spin configurations of the 2D Ising error model We map a bit σ = 0, 1, which parameterizes the number parity P = (−1) σ +1 of the BH model, to the spin-Z in the Ising model, i.e., Z i |0 i = |0 i and Z i |1 i = − |1 i .Here i, j is a link in a 2D square lattice and X i is the Pauli-X operator, which flips the parity, e.g., X i |0 = |1 .Z is a normalization factor, and β is a single fitting parameter in our EC.The properties of |Ψ(β ) , including the phases and phase transition, have been well investigated in the literature [26].For example, |Ψ(β ) is known to be in the paramagnet phase for β < β c ≈ 0.22, and in the ferromagnetic phase (spins are aligned along the X direction) for β > β c .
Our key observation is that |Ψ(β ) at small β < β c can be naturally identified with the parity snapshots of MI.For example, |Ψ(β → 0) represents the uniform parity configuration.It naturally corresponds to the MI with homogeneous site occupation, where the boson tunnelling strength J/U is strongly suppressed.When J/U is finite but small, there are exponentially short-ranged, virtual particle-hole excitations above the uniform background [Fig.1(b)].They will appear as a pair of bitflips in the snapshots.Such virtual excitations can be well captured by a non-zero but small β in Eq. ( 2) which allow the pair of parity flips above the uniform spin configuration [Fig.1(c)].The higher-order terms in β capture the longer-distance pair excitations, which are however exponentially short-ranged in the paramagnetic phase.Intuitively, β controls the strength of the virtual particle-hole fluctuations, which approaches β c as J/U approaches the transition toward SF.In practice, β can be efficiently calibrated for the given experimental snapshot data [21].
The error model allows us to extract the desired information hidden in the snapshot data.In particular, it can reveal the origin of the parity flips.The important quantity in this task is the spin-spin correlator, which approximates the probability for the two parity flips at the sites i E and j E to appear from the virtual particle-hole fluctuation.Intuitively, the virtual particle-hole pairs in the BH model [Fig.1(b)] correspond to XX-operators [Fig.1(c)] in the Ising model, and thus p(i E ↔ j E ) is naturally related to the spin-spin correlator.A more rigorous statement can be found in supplementary information [21].This provides a systematic route to identify the uncorrelated holes for a given experimental snapshot [21].For each snapshot [Fig.1(d)], we first calculate the probability p(i E ↔ j E ) for all the possible pairs of the parity flips [Fig.1(e)].We then pair each parity flip at i E with another at j E , which generates the maximum p(i E ↔ j E ) among many others [Fig.1(e)].When a parity flip at i E satisfies max Scaling of the brane parity correlator.Log plot of the brane parity correlator with domain length L at MI phase with J/U = 0.0391.Without the EC method (open symbol), the brane parity correlator drops quadratically with increasing domain length.On the other hand, the brane correlator with the EC method drops nearly linearly with the domain length.Inset: the exponential decay parameter a as a function of J/U.Solid line is a power-law fit curve with a ∝ (J/U) 1.8 , at an interval, 0.02 ≤ J/U ≤ 0.04.The data is obtained over 40 independent experimental runs, and the error bar denotes one standard error.
with the uncorrelated error rate η, then we conclude that the error is not correlated with all the other parity flips.Hence, it must be uncorrelated local errors, and we correct such errors from the snapshot data [Fig.1(f)] achieving our goal.
To demonstrate the effectiveness of our error correction (EC) protocol in removing uncorrelated errors, we utilize numerical images generated from Quantum Monte Carlo (QMC) simulation.After preparing a numerical image, we add random holes with its population rate η and compare the errorcorrected images with the original image.In the Mott insulator (MI) phase, we observe a negligible difference between the two images [21] It is worth noting that all the calculations and data processing involved in our EC procedure can be carried out very efficiently on a classical computer.It is more efficient than the classical simulation of the BH model, such as QMC simulation of the ground state.Additionally, even if one has a method that exactly simulates the ground state, it is not a priori clear how to identify the local errors in experimental parity snapshots.Only when one has a proper error model for the system, the pairing of the parity flips and identification of the incoherent local errors can be correctly done.
It should be also remarked that our work does not pursue achieving the quantum error correction [27] of logical qubits.Although our protocol is based on the pair matching of errors as the quantum ECs in toric code models [28,29], we solely focus on removing the errors from measurement data and identifying correlations between the parity flips.A comparison with previous quantum error correction protocols is provided in the supplementary information [21].To properly simulate the experiments, the QMC simulation contains the uncorrelated holes with the rate η ∼ 3%.The gray dotted line represents the quantum critical point of the Bose-Hubbard model in 2D,(J/U) c = 0.059 [22].The data is averaged over 40 different experimental realizations, and the shaded area denotes the maximum standard error at L = 3.For comparison, the integer parity correlators with L = 6, 12 (Fig. 2) are drawn in gray symbols and lines.(b) Domain length dependence of the generalized brane correlators at Mott insulating phase (J/U = 0.027) with different parameter α ∈ [0, 1].The integer parity correlator (α = 0) decays with L, while the fractional parity correlators are domain size independent when α ≥ 0.5 for MI states.The solid line represents the QMC data.

MEASUREMENT OF BRANE CORRELATORS
We turn our attention to the main result of this work, where we implement the EC to the experimental data and compute the brane parity correlator O D across the MI-SF phase transition (Fig. 2).Upon applying the EC, we observe a dramatic increase of the O D in the MI phase, and the parity correlator can well distinguish two phases even for the L = 12.The experimental results also show excellent agreement with QMC results after the error corrections.The error-corrected values of the brane parity correlator are also very close to that of the correlator at a low temperature (obtained from QMC) as ex-pected.The free holes from thermal fluctuations are part of the errors in simulating the ground state of MI, which our EC method targets to correct.Thus, our EC effectively lowers the temperature of the snapshots.
The effect of EC is also reflected in the scaling behaviors of O D with respect to the domain length L. First, we consider the case without the EC.In this case, O D is expected to scale as follows in which the first term ∝ L is the expected perimeter law for O D in MI [30].Whenever the domain boundary intersects the short-ranged particle-hole pairs, the brane operator returns a value of O D = −1 that exponentially suppress its expectation value with the perimeter length.The second term ∝ L 2 is due to the uncorrelated holes.Such holes are uncorrelated and randomly occur in each site, and hence its overall effect on the O D is proportional to the area of the domain D. Indeed, our experimental data (Fig. 3) is well fitted with a ≈ 0.21( 5) and b ≈ 0.034 (9).Note that the fitted value b is comparable to a.This explains why the bare measurement of O D (Fig. 2) does not show any signals regardless of the underlying many-body states.On the other hand, when EC is implemented, we immediately find that the expected perimeter law is extremely well followed (Fig. 3).Here, the EC-assisted data is fitted with a ≈ 0.058(1) b ≈ 0.0012 (1).We also observe that the decay coefficient a scales with tunneling energy a ∝ (J/U) 2 .This observation indicates the correlated particle-hole pairs in the MI phase from the quantum fluctuations, which is also represented in the parity correlation function C(i, j) = P i P j − P i P j [21].
Although the scaling behaviors of O D can distinguish both MI and SF phases in a finite system size, the value of O D itself cannot be used as an order parameter in the thermodynamic limit because of the perimeter law.Recently the generalized brane correlator has been proposed [14,15] to resolve this problem by considering the fractional version of the brane parity correlator where the angle θ = π/L −α depends on the domain length with a power law exponent α ∈ [0, 1].When α ≥ 0.5, the generalized brane correlator can have a finite value for MI in the thermodynamic limit and become zero in SF, serving as an order parameter for the MI [15].
The fractional parity correlator O θ D , however, cannot be evaluated from the standard parity projected fluorescence imaging system since it requires the information of the occupations number in each lattice site n i .By using our EC method, we can evaluate the O θ D even without additional experimental techniques to resolve the site occupation [31][32][33][34].Since our error model can identify the correlated particle-hole pairs, we can statistically infer the site occupation n i in the domain D, and thus estimate the generalized brane correlator in the MI regime.It is based on the fact that the number fluctuations is small inside the MI, such that the n i does not deviate much from its average n = 1.To measure the generalized brane parameter, we first assign the correlated particle-hole pairs for a given parity snapshot with probability computed within our EC model Eq.(2).Then, to a given pair, the site occupation number is randomly specified by either a doubly occupied site n i = 2 or an empty site n i = 0 and evaluate the O θ D .Although the number fluctuations can become large in the SF phase, the generalized brane correlator is zero in the SF so that we can statistically measure the O θ D across the SF-MI phase transition.
Figure 4 shows the EC-assisted measurements of generalized brane correlator.The fractional parity correlator with δ N 2 within a Gaussian approximation.Since the number fluctuations linearly increase with the domain length L in the MI (not shown), the angle θ = π/ √ L can normalize the brane correlator and remove the system size dependence.Indeed, varying the exponent α, we find that O θ D show a negligible dependence on L for MI only for α ≥ 0.5 [Fig.4(b)], which is consistent with the previous theory prediction [15].While in the superfluid phase, the number fluctuations are much stronger, and the fractional parity correlator rapidly vanishes over the domain length.Consequently, the phase transition curve becomes sharper with increasing the domain length, and the quantum critical point can be also fairly well identified.The result confirms the recent prediction that O θ D can be an order parameter of the MI phase in 2D [15].

CONCLUSIONS AND OUTLOOKS
We present an error correction method for a neutral atom quantum simulator using an optical lattice, which enables the measurement of non-local multi-point correlators in twi dimensions.Our method is based on mapping the parity snapshot data of the Bose-Hubbard system to the spin configurations of a 2D Ising error model, whose properties can be calculated efficiently.The error model allows us to systematically compute the probability for the two holes to be paired, and thus to identify and remove uncorrelated holes from the experimental image.Using the error-corrected images, we can confirm that the brane parity correlator shows a clear perimeter scaling law to the relatively large scales.Moreover, we successfully measure the generalized brane correlators with negligible dependence on domain size L, which can serve as the order parameter for MI.Our work opens up a number of promising directions for further study.One possible extension is to explore the non-local order parameters of the topological phases, such as a Haldane insulator in an extended Hubbard model [35,36], gapped spin SU(N) chain [37], and Kitaev-Heisenberg ladder [38].In particular, it will be interesting to generalize our EC method to the ultracold atom simula-tion of the Fermi Hubbard model [12,34,39].In these systems, the correlations in snapshot data between different spin species are essential in understanding the physics of manybody states.One may attempt to build an analogous Z 2 × Z 2 error model to capture the correlated virtual excitations and uncorrelated errors of the Fermi Hubbard model.Another immediate extension of our work is to apply our EC model to enhance the visibility of Wilson loop operators in quantum simulations of toric codes [4,16,40] as in [41,42].In the toric code, the Wilson loop operators are non-local operators that can detect and diagnose the topological order [4].Given the duality between the Z 2 gauge theory and the Ising model [43], it seems natural to use our approach to improve the diagnosis of topological orders [16,40].Note added: While preparing the manuscript, we became aware of experiment that measures brane parity correlator in 2D MI phase [44].Error correction estimation using the QMC data.Uncorrelated holes can be successfully removed using the EC method.From the QMC simulation, we estimate the incorrect decision of the error correction scheme.(b) Type-1 error of the error correction scheme.The type-1 error is the probability that we correct the wrong hole that does not come from the incoherent error source.We found that this probability is at most 1.2%, less than the experimental error rate of η ≈ 0.03 from the thermal excitation and the imaging loss.It implies that the measured brane correlator in the MI is not from the over-correction of the holes from the incorrect error correction method.(c) Type-2 error of the error correction scheme.Contrary, the type-2 error in MI is much larger than the type-1 error.The error correction method prefers not to correct the incoherent error when it is not certain.In the SF regime, the type-2 error is larger than 0.1 in the whole J/U values and it became 1.

Estimation of Error rate η
To measure the loss during the imaging process, we took two consecutive images with 1s exposure and 0.1s interval.By comparing two images, we calculated the image fidelity by the Raman sideband imaging.The error during the imaging process mostly appears in the form of the loss with probability η f id = 0.01.Another source for the uncorrelated error is the thermal excitation of the atoms.We measured the average filling of the Mott insulator to be n = 0.97.We compared it with the QMC simulation (see below) and computed that η thermal = 0.02 which makes total uncorrelated error η = η thermal + η f id = 0.03.

QMC simulations
We use the directed loop QMC simulation proposed in Ref. [46,47].The system size is set to be 50 × 50.The maximum boson occupation is set to be n max = 3.The harmonic curvature is set to be ω = 3.156 × 10 −2 U with the potential V = 1 2 ωr 2 .For each J/U, we calibrate the chemical potential so that the mean total boson number in QMC matches that of the experiment.
The total particle number is ∼ 1220.The temperature is calibrated to T = 1/12U so that we have η thermal ≈ 0.018.We sample 3 × 10 4 snapshots to compute the expectation values of brane correlators.The error bar represents the standard error of mean.The QMC simulation with random noise is drawn by using the error correction (solid line) and without the error correction (dotted line).The density correlation is maximum near the critical point.

Error investigation of the error correction method with QMC data
We investigated the Type-1 and Type-2 errors of the error correction method by correcting the added uncorrelated errors with the QMC simulation in Fig. S1.The probability of Type-1 error is at most 1.2%, less than the experimental error rate of η ≈ 0.03 from the thermal excitation and the imaging loss.It implies that the measured brane correlator in the MI is not from the overcorrection of the holes from the incorrect error correction method.Contrarily, the type-2 error in MI is much larger than the type-1 error.The error correction method prefers not to correct the incoherent error when it is not certain.

SECTION B: TWO-SITE PARITY CORRELATOR
We measured the two-site correlator and the density fluctuations with different values of J/U and demonstrate the presence of particle-hole pairs in the Mott insulator.The two-site correlation is maximized near the critical point at J/U = 0.059.The experimental result is well matched with the numerical simulation result considering the harmonic curvature.Using the error correction method, the two-site parity correlator is increased in the Mott insulator phase where the error correction is successfully done.In the SF phase, the error correction method cannot distinguish the uncorrelated holes from the correlated holes, and the parity correlator remains the same.

SECTION C: RELATION WITH PREVIOUS QEC METHODS
Here we present a brief comment on the relation of our EC method to the previous works.We will mainly compare our EC method to quantum EC protocols on the toric codes [27][28][29].We first briefly review the quantum error correction (QEC) protocols in the toric code.We then compare them with our method.In the toric code, the errors always occur in pairs, and they may propagate far from each other.The QEC then attempts to pair these errors from the measurement data and remove them from the wavefunctions by performing proper quantum operations.In this process, it is important to define properly how the pair of errors are correlated.Dannis and Kitaev [27] assigned an error probability on each link of the toric code and mapped the QEC problem to the problem of finding the minimum energy configuration of the random bond Ising model.The (free) energy minimization problem is then mapped to other mathematical problems, such as the minimal weight perfect matching (MWPM) problem in graph theory, which can be efficiently solved [28,29].Note that these methods can be used to detect correlated errors in our setup.One can also identify uncorrelated measurement errors by repeated measurements.In summary, QEC protocols target to identify and correct both the uncorrelated and correlated errors on quantum states.In this respect, one may find certain similarity between our EC method and the QEC protocols.However, there are several important differences.
First of all, we note that pairs of the parity flips in the parity snapshots of the MI are not errors that should be removed.Instead, we use this information to isolate and remove the uncorrelated holes, which will have no partner in the snapshot.The correlation between the parity flips also provides an estimate for the number fluctuations, from which we can evaluate the generalized brane correlators.In this regard, the goal of our method and that of the QEC protocols are different.
In addition, the QEC protocols for the toric code cannot perform the tasks that we need to do in our setup.To better appreciate this, we remind that the QEC protocols concern how to find the single most probable pairing pattern between the errors.This is legitimate for the QECs, which have to determine a pairing configuration and correct them.However, this does not take other probable patterns into account, which may appear in our problem.On the other hand, our EC method takes other probable patterns into account by finding and keeping the probability of each pair of holes in a parity snapshot to be correlated.By doing so, our method can pair long-ranged parity flips, which are hardly found by QEC protocols but relevant to our purpose.Holes in a parity snapshot could be correlated even when they are not nearby.In particular, such a long-ranged correlation gives a scaling of the brane correlator super-exponential in the perimeter of the membrane in the SF phase.This scaling would not be well captured by QEC protocols, for example the MWPM, which pairs parity flips locally.On the other hand, our method can capture the scaling well as demonstrated in the main text.

APPENDIX D. DETAILS OF ERROR CORRECTION
This section contains details of our EC methods, which were schematically discussed in the main text.Let us start with a brief note.If one wishes to develop a new EC method, there are a few (natural) requirements.First, the new EC scheme should be designed to reconstruct correlations appropriately between holes in each parity snapshot.In addition, it should be efficient than the full simulation of the BH model.In this Appendix, we will see that our error correction scheme satisfies the two requirements nicely.
The rest of this section is organized as follows.First, we define the correlation between holes in a snapshot in Sec.I, which will be related with the spin-spin correlation function of transverse field Ising model (TFIM) in Sec.III.We then show in Sec.II that the correlation of TFIM can be classically efficiently estimated.We next present our detailed arguments why the correlation of BH model between two holes can be approximated by that of TFIM in Sec.III.To approximate the correlations of the BH model as that of the TFIM, a single parameter of TFIM, which is the "temperature" β of TFIM, should be calibrated.The calibration is done by matching the expectation values of the brane parity correlators of BH model and TFIM, which is discussed in Sec.V.In Sec.IV, we introduce how to efficiently compute the expectation values of the brane correlators in TFIM.In Sec.V, we finally discuss how to fix the parameter β of TFIM and complete our EC protocol.

I. DEFINITION OF CORRELATION BETWEEN HOLES
Here, we define what the correlation between holes in a snapshot means.Let us assume that we have a parity snapshot σ of a wavefunction |ψ .σ is a vector of 0 and 1 whose elements represent the parity of each site, i.e., P i = (−1) σ i (0 and 1 represent odd and even parities of the site occupation of the BH model, respectively).On top of this parity snapshot pattern σ , we want to introduce two parity flips (two holes) at i and j with σ i = σ j = 1.See Fig. S3.We will write the relation between σ and σ as Here, ⊕ is the modulo two summation, and 1 i is the unit vector having the unit i-th element.Fig. S3 illustrates the definition of σ ⊕ 1 i ⊕ 1 j .The two parity flips 1 i ⊕ 1 j in σ can be introduced either from a paired virtual excitation of |ψ or uncorrelated errors, e.g.free holes.The two holes are correlated if and only if they are introduced by a paired virtual excitation.The likelihoods of the two cases are given by with the probability p(σ ) = | σ |ψ | 2 , the parity flipped snapshot σ = σ ⊕ 1 i ⊕ 1 j , and the measurement error probability η.
For a given parity configuration σ , we now ask the origin of the two parity flip.We may answer this question by comparing the two likelihoods.More precisely, we may conduct the maximum likelihood test: the two holes originate from measurement errors if and only if L corr < L err .Equivalently, we may compare p(σ )/p(σ ) with η.To simplify the notation, we will define f (i, j|σ ) ≡ p(σ )/p(σ ).We interpret that f (i, j|σ ) measures how easy to introduce a parity flip 1 i or 1 j on σ via a virtual excitation, which will be compared to measurement errors.
In general, a snapshot σ has more than two holes, so each hole can pair with one of the other holes or can be left uncorrelated.
In this case as well, we will conduct the maximum likelihood test per each parity flip.More precisely, for a given hole of σ at the site i, we first compute f (i, j|σ ⊕ 1 i ⊕ 1 j ) for all the other holes at the site j.We then compare its maximum value with the error probability η.The hole at the site i is uncorrelated if and only if Note that since the number of holes in a snapshot is ∼ O(L x L y ), the number of computations of f (i, j|σ ) for each snapshot scales O(L 2 x L 2 y ).Thus, if f (i, j|σ ) is hard to compute, then we may not be able to conduct the maximum likelihood test.
In fact, the computation of f (i, j|σ ) requires the ratio between p(σ ) and p(σ ⊕ 1 i ⊕ 1 j ), which is almost equivalent to the full computation of the ground state wavefunction.Thus, the direct computation of f (i, j|σ ) for the error correction is not meaningful because it is already as expansive as the full simulation of the BH model.Hence, we instead approximate f (i, j|σ ) by its average Before explaining the factor 1/2 in the above, we note that the above quantity measures how easy, on average, to introduce a parity flip at 0 i or 0 j on a typical snapshot of |ψ via a virtual particle-hole pair excitation.This can be compared to the uncorrelated error rate η to perform the maximum likelihood test.We will also regard this quantity as the square root of the likelihood p(i ↔ j) of the pair of the parity flips at the two sites i and j being correlated.In the next section, we will see that of the TFIM is equivalent to the spin-spin correlator, which can be easily computed.Hence, the maximum likelihood test can be done classically efficiently.
The factor of 1/2 in Eq. (S3) comes from the fact that we always choose a pair of parity flips in a snapshot σ and try to reconstruct σ via undoing the virtual particle-hole excitation between the two sites.This introduces an ordering between snapshots, and the ordering can be approximately encoded by the condition p(σ ) > p(σ ).This condition comes from the fact that the introduction of a virtual particle-hole excitation in MI requires finite energy cost, so the probability assigned on σ is typically smaller than that on σ .Then, the average of f (i, j|σ ) over configurations satisfying p(σ ) > p(σ ) is given by Here, the sum over σ satisfying p(σ ) > p(σ ) is equivalent to the sum over all possible pairs of σ with the 1/2 factor.
The validity of our approximation above can be confirmed explicitly in the perturbative regime.Let us assume that we have a snapshot σ of a MI state having a single pair of particle-hole excitations in a local region (see Fig. S4).Then, the probability of getting σ is much smaller than that of the configuration σ .In addition, the probability of getting the snapshot σ from the ground state wavefunction is almost equal to one.Also probabilities for getting other configurations are negligibly smaller than p(σ ) in the same reasoning.Thus, we have In other words, we have

II. VIRTUAL EXCITATION RECONSTRUCTION OF TFIM
Here, we show how the virtual excitation reconstruction of TFIM can be done classically efficiently.The ground state of TFIM that we will consider is given by the following imaginary-time evolved wave function: where Z is given by We want to emphasize that this wavefunction has a single free parameter β .The free parameter will be used to match the statistical nature of the parity snapshots of BH model.The parent Hamiltonian of this wave function is given by [48] For small β , it becomes the conventional transverse field Ising model.This model also has the paramagnet and ferromagnet phases at roughly 2β c = 0.441, which is the well-known critical temperature of the 2D Ising model.
One particularly important property of |ψ in the above is that it is a "classical" wavefunction, i.e., all the coefficients of |ψ can be set to be positive and are given by the square root of the probability.Even more, a straightforward calculation shows that E[ f (i, j|σ )] becomes the two-point X i X j correlator: Here, p(σ ) p(σ ) with σ = σ ⊕ 1 i ⊕ 1 j is the product of the coefficients of |σ and |σ , and σ is a configuration that can be obtained by applying spin flips on iand j-th sites.Thus, the sum of p(σ ) p(σ ) over all configuration σ is equivalent to 1 2 X i X j .Using this, we can conduct the error correction for TFIM efficiently, if we can compute X i X j easily.Now we show that X i X j can be indeed computed classically efficiently.In fact, it is given by the two-point correlator of the classical 2D Ising model: 0| e 2β ∑ i, j X i X j |0 = ∑ x x i x j e 2β ∑ i, j x i x j ∑ x e 2β ∑ i, j x i x j = x i x j Ising,2β , (S10) which is the thermal expectation value of the two-point correlator x i x j of the classical 2D Ising model.Note that the two-point correlator x i x j of the classical 2D Ising model for all i, j can be classically efficiently computed from a single Monte Carlo simulation [49].

III. REASONING BEHIND APPROXIMATION BY TFIM
Here, we provide some intuitive explanaions why dynamics of virtual particles of BH model can be well mimicked by that of TFIM.We summarize our mapping between the two models in Table S1.To support the mapping, we will first show that the probability distribution of parity configurations of BH model (with the mean filling of n = 1) and TFIM can be matched perturbatively.This implies that E[ f (i, j|σ )] can also be matched for the two models.In addition, we also argue that the dynamics of virtual particles of BH model and TFIM are qualitatively similar.Finally, we show that expectation values of diagonal observables of the both models, i.e., observables in Z or parity basis, can be matched simultaneously.
First, we note that the wavefunctions of both models at the stable fixed points give the same probability distribution for the spin σ or number parity P at each site: (S11) The spin and the number parity is related by P = (−1) σ +1 .Second, the energy cost of introducing a pair of two excitations X i X j or b † i b j to the fixed-point configuration σ = 0 or P = 1 of MI is given by ∆E = 2 or ∆E = U, respectively.Note that each Boson hopping operator has a direction, while each Ising operator has no direction.Thus, one may expect that some perturbative configurations made by Boson hopping operators cannot be made by Ising operators.However, in the MI phase, the energy cost of making such a configuration is large.For example, |030 number state can be made by applying b † 2 b 1 and b † 2 b 3 operators to |111 number state.Obviously, |030 number state does not have the corresponding state in TFIM.However, its energy cost in MI is 3U, so it is unlikely to be seen for large U.It means that their parity distribution can be matched at least within the perturbation theory in MI.Fig. S4 illustrates this in more detail.
In addition, the two-point correlators b † i b j (in BH model) and X i X j (in TFIM), which measure how easy to introduce a pair of virtual particles, behave similarly.First, they both exponentially decay to their length |i − j| in the MI and paramagnetic phases, i.e., b † i b j and X i X j follows ∼ e −|i− j|/ξ with the correlation length ξ .Second, they both converge to finite numbers in the SF and ferromagnetic phases at the zero temperature.
Based on these facts, we expect that the values of the observables of both the models, which are diagonal in the parity basis or spin-Z basis, can be well matched for MI.A particularly interesting diagonal observable in both sides is the brane parity correlator ∏ i∈D P i in BH model or equivalently the domain wall operator ∏ i∈D Z i in TFIM.The brane parity correlator decays exponentially in |∂ D| in MI (which corresponds to the paramagnetic phase in TFIM) and super-exponentially in |∂ D| in SF (which corresponds to the ferromagnetic phase in TFIM).As the diagonal observables can be well matched in both sides, we may calibrate the parameter β of TFIM by matching the expectation value of the brane parity correlator in both sides.In the next Here, l i 1 and l i 2 are sites connected by the link l i .Each spin σ i at the i-th site is given by the parity of the total number of links connected with the site, i.e., σ i = ∑ l∈i l i (mod 2).These link configurations are redundant, and multiple link configurations can give a single spin configuration.We will refer to the set of all links giving a spin configuration σ as link(σ ).The probability associated with a configuration |σ is then given by We can think of it as products of two classical systems with link variables l and s whose statistical weights are β n n! where n is the degeneracy of each link.We note that since l and s give the same spin configuration σ , they should be related by a closed loop of links o, which does not change the spin configuration.Our next step is to develop a cFigS3lassical Monte Carlo scheme simulating Eq.S14.The classical MC protocol comprised of two updates: the parity update and the loop update.The parity update changes σ by adding or removing one link term from a link.The statistical weight associated with adding (removing) a link term to l i and s i is β 2 /(n l i + 1)(n s i + 1) (n l i n s i /β 2 ).The loop update adds a closed loop of links on either l or s.To do so, we first choose arbitrary position and inject a worm which propagates along links.It adds (removes) a link term on each link based on the statistical weight β /(n + 1).The worm moves until its head and rear meet.We find that two update schemes are sufficient to simulate the probability distribution p(σ ) of the wavefunction |ψ .

V. CALIBRATION OF β
In this section, we elaborate on how to calibrate the parameter β of the wavefunction |ψ .We calibrate β by matching the expectation value of the domain wall operator of |ψ with the expectation value of the brane parity correlator computed from error-corrected parity snapshots.We scan β to find the appropriate β that optimally matches the two observables in the two models.For each β , we compute the domain wall operator of |ψ as discussed in Sec.VI.We then perform error correction on the parity snapshot measurements using the correlator X i X j .Estimation of X i X j is illustrated in Sec.VII.After scanning β , we find the optimal β .More details on the fixing β can be found in Sec.VIII.The overall process of the calibration is illustrated on Fig. S5.

VI. CALCULATION OF DOMAIN WALL OPERATOR IN TFIM
Here, we discuss how to set β of |ψ .To set β , we match the expectation value of the brane parity correlator of BH model with the expectation value of the domain wall operator of the TFIM.For each β , given that the uncorrelated error rate is non-vanishing, we can first identify uncorrelated errors in snapshots of BH model and compute the brane parity correlator.In addition, we can also compute the expectation value of the domain wall operator of the TFIM.We then try to find a self-consistant β that matches the expectation value of the brane parity correlator (obtained from the error-corrected parity snapshots) with the expectation value of the domain wall operator.FIG.S5.The overall process of the calibration.We try to find β that matches ÔD and ÔD,TFIM .For each β , we compute ÔD from error corrected experimental snapshots and ÔD,TFIM from the transverse field Ising model.We scan β and find the value that matches the two optimally.
The detailed fitting process is as follows.We first compute the domain wall operator as a function of β .We then interpolate β dependence of the domain wall operator using a (cubic) polynomial f (β ).Note that this is equivalent to introduce a cutoff on β , i.e., it cannot be greater than a certain value β * , which is a zero of f .This is particularly important in SF phase.
The reason why we introduce such a temperature cutoff is two-fold.First, let us assume that we do not introduce any cutoff and O D is small, i.e., SF states.Then, β becomes extremely sensitive to the exact numeric values of O D .It means that if O D is not sufficiently accurate compared to the experimental error bars, then β cannot be determined.To fix β in this case, we introduce a cutoff on J/U and corresponding β .This make β not exceed the cutoff even when O D is uncertain.Second, since we compute the brane parity correlator to distinguish the two phases of BH model, practically we do not need to increase β above β c ≈ 0.22 if the brane parity correlator is already vanishing.

VII. CALCULATION OF TWO-POINT CORRELATORS X i X j
We compute the two-point correlator X i X j using the classical Ising correlator x i x j Ising,2β as discussed in Eq. (S10).To compute the Ising correlator, we perform the classical Monte Carlo simulation proposed in Ref. [49] on 64 × 64 square lattice, which is sufficiently larger than 20 × 20 experimental parity snapshots.We note that a single Monte Carlo simulation is sufficient to compute x i x j Ising,2β for all i, j.

VIII. SEARCH FOR OPTIMAL β
We find the optimal β by comparing expectation values of the brane parity correlator in BH model and the domain wall operator in TFIM.For typical J/U, we can find a unique β = 0 that perfectly matches the two observables.In this case, we can unambiguously fix β and error correction is self-consistent.However, it is also possible that there are multiple self-consistent solutions β = 0 that perfectly match the observables in the two sides.In this case, we choose the minimum β , which maximally (but not entirely) removes holes in a self-consistent way.In addition, for small J/U, there could be no consistent solution except the trivial one with β = 0, which eliminates all holes in snapshots.Although the trivial solution is also a self-consistent solution,

FIG. 1 .
FIG. 1. Schematic diagram of the error correction protocol.(a) Mott insulating (MI) phase of neutral atoms in a two-dimensional square optical lattice.The entangled particle-hole pairs (yellow circles) can be generated because of the finite tunneling amplitude.An uncorrelated hole (purple dashed circle) can be found in the MI phase because of experimental imperfections.(b),(c) Correspondence between the Ising error model and the Bose-Hubbard (BH) model.In the BH model, the particle-hole excitations can be created by the tunneling operator b † i b j on the background of a unity filling.After the fluorescence imaging, atoms in the doubly occupied sites are lost, resulting in an empty lattice site (dashed circles).These sites have the same even parity.In the Ising model, the Pauli-X correlator X i X j creates a paired spin-up excitation along a link (red line).The particle-hole pairs after the parity projection can be considered as the excitations in the Ising model.(d) Experimental snapshot image of MI phase.The correlated particle-hole pairs are not distinguishable from the uncorrelated hole.(e) The error identification protocol: firstly, for a given single experimental image, we evaluate the probability p(i E ↔ j E ) for the pairs of the parity flips at the sites i E and j E .The pairs with higher (lower) probabilities are marked by dark (light) red arrows.The uncorrelated hole at the site i E can be identified when the probabilities for all possible pairs p(i E ↔ j E ) is smaller than the (pair of) error rate η 2 .(f) Error corrected image.Correlated particle-hole pairs (tied with yellow lines) and the uncorrelated hole (blue-shaded circle) are identified.To determine the fitting parameter β of the error model, we compare the brane correlator in the BH model (gray box) and the domain wall operator in the Ising model (yellow box) [21].

FIG. 2 .
FIG. 2. Brane parity correlator.(a) The error corrected image in the experiment.The error correction method detects and removes the uncorrelated holes in the snapshots.(b) Brane parity correlator in two-dimensional Bose Hubbard model with the domain size D = 12 × 12. Without the EC (open circle), the brane correlator hardly distinguishes MI and superfluid phase.The EC-assisted measurement (closed circle) shows a dramatic signal increase and can identify the superfluid-to-Mott insulating phase transition.Inset shows the brane parity correlator with L = 6.Solid lines are the QMC simulation at T = 0.083U and uncorrelated hole rate η = 0.03 with EC (blue) and without (lightblue) the EC method.The experimental data agree well with the QMC simulation.They can compare with the QMC simulation with no uncorrelated error η = 0 at low temperature T = 1/32U (darkblue dotted line).The gray dotted line represents the quantum critical point in 2D, estimated by the QMC simulation (J/U) c = 0.059[22].Each data point is obtained over 40 independent experimental runs, and the shaded region denotes the maximum standard error (see supplemental material for details[21]).
FIG. 3.Scaling of the brane parity correlator.Log plot of the brane parity correlator with domain length L at MI phase with J/U = 0.0391.Without the EC method (open symbol), the brane parity correlator drops quadratically with increasing domain length.On the other hand, the brane correlator with the EC method drops nearly linearly with the domain length.Inset: the exponential decay parameter a as a function of J/U.Solid line is a power-law fit curve with a ∝ (J/U)1.8, at an interval, 0.02 ≤ J/U ≤ 0.04.The data is obtained over 40 independent experimental runs, and the error bar denotes one standard error.

FIG. 4 .
FIG. 4. Generalized brane correlator.(a)The generalized brane correlator with α = 0.5 is drawn for different domain length L = 3 − 12.The generalized brane correlator for MI stays almost constant even at large system size L = 12.The solid line represents the generalized brane correlator, calculated within the QMC data.To properly simulate the experiments, the QMC simulation contains the uncorrelated holes with the rate η ∼ 3%.The gray dotted line represents the quantum critical point of the Bose-Hubbard model in 2D,(J/U) c = 0.059[22].The data is averaged over 40 different experimental realizations, and the shaded area denotes the maximum standard error at L = 3.For comparison, the integer parity correlators with L = 6, 12 (Fig.2) are drawn in gray symbols and lines.(b) Domain length dependence of the generalized brane correlators at Mott insulating phase (J/U = 0.027) with different parameter α ∈ [0, 1].The integer parity correlator (α = 0) decays with L, while the fractional parity correlators are domain size independent when α ≥ 0.5 for MI states.The solid line represents the QMC data.
can well distinguish MI from SF and has a negligible dependence on the domain length in MI [Fig.4(a)].These results can be understood by noticing the generalized brane correlator in the MI becomes O θ D e − π 2 θ 2 2 FIG. S1. (a)Error correction estimation using the QMC data.Uncorrelated holes can be successfully removed using the EC method.From the QMC simulation, we estimate the incorrect decision of the error correction scheme.(b) Type-1 error of the error correction scheme.The type-1 error is the probability that we correct the wrong hole that does not come from the incoherent error source.We found that this probability is at most 1.2%, less than the experimental error rate of η ≈ 0.03 from the thermal excitation and the imaging loss.It implies that the measured brane correlator in the MI is not from the over-correction of the holes from the incorrect error correction method.(c) Type-2 error of the error correction scheme.Contrary, the type-2 error in MI is much larger than the type-1 error.The error correction method prefers not to correct the incoherent error when it is not certain.In the SF regime, the type-2 error is larger than 0.1 in the whole J/U values and it became 1.
FIG.S2.Parity-projected density correlation.Experimental data is drawn with the error correction (closed symbol) and the without the error correction (open symbol).The error bar represents the standard error of mean.The QMC simulation with random noise is drawn by using the error correction (solid line) and without the error correction (dotted line).The density correlation is maximum near the critical point.

TABLE S1 .
Mappings between Bose-Hubbard model and Transverse Field Ising model.The mappings are exact deep inside the MI phase.