Detecting Hidden Order in Fractional Chern Insulators

Topological phase transitions go beyond Ginzburg and Landau's paradigm of spontaneous symmetry breaking and occur without an associated local order parameter. Instead, such transitions can be characterized by the emergence of non-local order parameters, which require measurements on extensively many particles simultaneously - an impossible venture in real materials. On the other hand, quantum simulators have demonstrated such measurements, making them prime candidates for an experimental confirmation of non-local topological order. Here, building upon the recent advances in preparing few-particle fractional Chern insulators using ultracold atoms and photons, we propose a realistic scheme for detecting the hidden off-diagonal long-range order (HODLRO) characterizing Laughlin states. Furthermore, we demonstrate the existence of this hidden order in fractional Chern insulators, specifically for the $\nu=\frac{1}{2}$-Laughlin state in the isotropic Hofstadter-Bose-Hubbard model. This is achieved by large-scale numerical density matrix renormalization group (DMRG) simulations based on matrix product states, for which we formulate an efficient sampling procedure providing direct access to HODLRO in close analogy to the proposed experimental scheme. We confirm the characteristic power-law scaling of HODLRO, with an exponent $\frac{1}{\nu} = 2$, and show that its detection requires only a few thousand snapshots. This makes our scheme realistically achievable with current technology and paves the way for further analysis of non-local topological orders, e.g. in topological states with non-Abelian anyonic excitations.


I. INTRODUCTION
Over the last decades, interacting topological systems have provided an exciting opportunity to explore exotic states of matter, most prominently states exhibiting intrinsic topological order.A prime example for this paradigm is the fractional quantum Hall (FQH) effect, where strongly interacting particles in two dimensions are subject to a strong magnetic field [1,2].In contrast to symmetry breaking phases of matter, such topologically ordered states cannot be described by a local order parameter.Nevertheless, since the early days of FQH physics researchers have been intrigued by the idea of using non-local extensions of conventional order parameters in order to characterize topological states of matter and enable insights into their remarkable properties [3][4][5][6].Specifically, based on variational trial states like the Laughlin state [7], the emergence of so-called hidden offdiagonal long-range order (HODLRO) has been established theoretically in continuum systems.This multiparticle variant of conventional off-diagonal long-range order (ODLRO), characterizing Bose-Einstein condensation associated with a global U(1)-symmetry breaking, will play a central topic in this article.As we will discuss, in order to detect HODLRO in a FQH system a measurement of long-range one-particle coherence, ∼ ψ † (r 1 +d) ψ(r 1 ), has to be combined with a simultane-ous measurement of the positions r 2 , ..., r N of all remaining N − 1 particles.This makes a direct measurement of HODLRO unfeasible in traditional solid state experiments.
In this article we show that state-of-the-art quantum simulation platforms, and cold atoms in particular, can be used to directly detect HODLRO in FQH states.Through large-scale numerical density-matrix renormalization group (DMRG) simulations [37][38][39] of the isotropic Hofstadter-Bose-Hubbard model at mag- ○ Possible experimental realization of the proposed measurement protocol using a synthetic bilayer Hofstadter-Hubbard system with internal states |σ = ±⟩. 2 ○ On one side of the system, the edges are coupled to obtain an effective single layer system.3 ○ Measurements of off-diagonal elements of the one-particle reduced density matrix, â † x,y,+ âx,y,− are possible using a local pulse protocol involving both layers, as described in the main text.(b) Correlations extracted from Nsnaps = 10 4 numerical snapshots without any post-processing and with experimentally relevant errorbars from shot-noise.The algebraic decay of the two-site correlations of the composite bosons |ρ x,x ′ ;y | ∼ b † x,y bx ′ ,y (red squares) can be clearly distinguished from the exponential decay of the original bosonic correlations |ρ x,x ′ ;y | ∼ ⟨â † x,y âx ′ ,y ⟩ (blue circles).(c) Tensor-network diagram of the extended perfect sampling scheme which allows for numerical snapshots of composite boson correlations |ρ x,x ′ ;y | equivalent to the ones obtained in the proposed experiment.
netic filling ν = 1 /2, we demonstrate that fractional Chern insulators in lattice systems exhibit the same form of universal HODLRO as the corresponding continuum FQH states.To this end, we perform projective measurements on the matrix-product state (MPS) representation of the ground state, building on previous works introducing a perfect sampling algorithm for MPS [40] to emulate snapshots of local observables [41] via a successive collapse of the many-body wave-function.We extend this scheme in such a way that it enables snapshots of non-local quantities (see Fig. 1(c)) to sample correlation functions like those needed to detect the multi-particle HODLRO.Our numerical prediction can be experimentally tested using two hyperfine manifolds in ultracold atoms in a quantum gas microscope (see Fig. 1(a)), for which we propose a specific experimental implementation enabling long-range measurements of the one-particle coherence combined with simultaneous single-site resolved measurements of the density on the whole lattice, granting direct access to HODLRO (see Fig. 1(b)).

The remainder of this paper is structured as follows:
In Section II we present an overview of the main results of this article.We briefly introduce the key concepts of HODLRO and elucidate explicitly how to obtain snapshots necessary for the detection of this putative multiparticle off-diagonal long-range order from measurements on cold atoms in optical lattices.In Section III we extensively discuss the results of our large-scale numerical analysis of the Hofstadter-Bose-Hubbard ground state showing the absence of conventional long-range order and the emergence of HODLRO using the MPS-based sampling method.Finally, in Section IV we present the MPS-based projective sampling scheme in a self-contained review which can be read independently.

II. OVERVIEW OF MAIN RESULTS
This section summarizes the main findings of the paper.We begin with a short review of HODLRO in continuum FQH systems.We then proceed by generalizing this concept to lattice systems hosting fractional Chern insulators.After introducing a scheme to resolve long-range coherence in quantum gas microscopes, exploiting a bilayer geometry, a detailed experimental setup is proposed which yields all information needed to obtain HODLRO.

II.1. A brief Review of HODLRO
The ordering in superfluids and superconductors can be famously characterized by the existence of off-diagonal long-range order (ODLRO) associated with the breaking of the particle number preserving U(1) gauge symmetry of the system signaling Bose-Einstein condensation.Consider the one-particle reduced density matrix In general, one expects ρ r,r+d to fall off rapidly and decay to zero for increasing values of d.In systems exhibiting ODLRO however, the one-particle correlations saturate to a finite value in this limit Contrary, one of the hallmarks of topologically ordered states of matter is the absence of any kind of one-particle long-range order.For example, for Laughlin's trial wave function describing the ground state of a continuum FQH system at magnetic filling ν = 1 /m, where z j = x j +iz j is the jth particle position in the complex plane, it can be shown that the one-particle reduced density matrix decays exponentially Nonetheless, there is a characteristic type of order hidden in the density matrix of the Laughlin state.It can be revealed by introducing a singular gauge transformation [42][43][44] of the form The effect of this transformation is a contribution of the form to the wave function which can be understood as endowing each of the particles with m quanta of magnetic flux Φ 0 .Independent of m being odd/even, Eq. ( 5) maps the fermionic/bosonic FQH problem to composites of flux quanta and charge with bosonic statistics which are hence referred to as composite bosons ψ( †) CB (z j ).For a more detailed derivation, see Appendix A. In this composite boson basis the continuum one-particle reduced density matrix decays algebraically, i.e.
This emergence of quasi long-range order, signaling a condensation of the composite bosons in the Laughlin state [4,6] defines hidden off-diagonal long range order (HODLRO).

II.2. HODLRO in Lattice Systems
We can readily expand the concept of HODLRO to discrete systems.Inspired by earlier studies of lattice anyons [45,46], in a lattice system with local bosonic degrees of freedom â( †) k , attaching m flux quanta to each boson can be achieved by the following lattice gauge transformation where the bold letters represent the coordinates of a particular lattice site, e.g.k = (x, y) and nk = â † k âk is the local density operator.Evidently, Φj is the lattice analog of the continuum emergent gauge degrees of freedom of A j .Hence, using Eq. ( 7), we can define composite boson operators b( †) j := e (−)i Φ( †) j â( †) j , describing the desired composites of m flux quanta and a boson.
In turn, we can introduce the one-particle correlation function for the composite bosons, ρj,l = b † j bl , which can be expressed purely in terms of the original bosonic annihilation (creation) operators â( †) j and the boson number operators nk ρj,l = k̸ =j,l (10) Consequently, the multi-particle composite boson correlation function ρj,l picks up an additional non-trivial phase for each occupied lattice site k ̸ = j, l with respect to the usual correlator ρ j,l ∼ ⟨â † j âl ⟩.Therefore, in order to probe this lattice analog of HODLRO one not only has to measure the long-range bosonic one-particle coherence ρ j,l but combine it with a simultaneous measurement of local densities operators {n k } k̸ =j,l .This enables us to examine lattice models described by the Hofstadter-Bose-Hubbard Hamiltonian Ĥ = −t ⟨j,k⟩ e iϕ j→k â † j âk + H.c. (11) which hosts fractional Chern insulators.The phase ϕ j→k is determined by the choice of the magnetic vector potential and the sum runs over adjacent sites.

II.3. Measuring Off-diagonal Long-range Order
Despite being highly adjustable and offering singlelattice-site and single-particle resolution, so far quantumgas microscopes only allow for a direct measurement of physical quantities which are diagonal or near-diagonal in the occupation number basis.To resolve phase coherence one relies on a combination of a time-of-flight expansion followed by absorption imaging granting access to the momentum distribution [47].That means, even though there has been recent progress in resolving coherence and currents between adjacent sites [48,49], measuring longrange single-particle correlation functions of the form from single-site measurements is a non-trivial task.
Here, we propose to use a bilayer geometry consisting of two hyperfine manifolds of ultracold atoms which are coupled along one edge of the layers, see Fig. 1(a).Using such a setup an effective single-layer is achieved in which we can realize a two-dimensional lattice model.Exploiting this double layer structure, non-local inter-lattice-site correlations ∼ â † r âr ′ of the effective model, where r and r ′ denote sites in different layers but are on top of each other -i.e.r = (x, y, +), r ′ = (x, y, −) -, can be accessed by measuring local correlations between the layers.In order to extract this inter-layer correlations, we propose to perform a local π /2-pulse on the probe sites r and r ′ to rotate the measurement basis into a coherent superposition of both layers.In this manner, the inter-layer correlations can be extracted using standard single-site resolved quantum gas microscopy.In the following subsection, we will discuss in detail how such measurements can be performed and how the bilayer structure enables us to naturally treat bosonic atoms in an artificial, homogeneous magnetic field, making it a prime candidate to probe fractional Chern insulators.

II.4. Experimental Protocol
In order to detect HODLRO in a cold atom quantum simulator we first need to achieve a long-range measurement of one-particle coherence and secondly, perform projective measurements on the site-local Fock spaces of the remaining sites simultaneously.For a particularly simple measurement of long-range two-point correlations ∼ â † r âr ′ , we exploit the bilayer approach described in Section II.3.In particular, we propose to realize a synthetic bilayer Hofstadter-Bose-Hubbard model using internal states |±⟩ of the atoms to realize the different manifolds.As already demonstrated experimentally [13,14], artificial gauge fields can be realized with opposite signs for the two internal states of the atom.Thus, it is possible to realize a Hamiltonian of the form Ĥexp = − t x,y,σ â † x+1,y,σ âx,y,σ +e σ2πiα(x+1/2) â † x,y+1,σ âx,y,σ + H.c.
where â( †) x,y,σ are the bosonic annihilation (creation) operators and nx,y,σ = â † x,y,σ âx,y,σ are the boson number operators.Here, (x, y) are the lattice positions on a square lattice of size L x × L y parametrizing a single layer while σ = ± labels the internal state of the atoms.The Hamiltonian is written in Landau gauge, i.e. the physical vector potential which generates the homogeneous flux per plaquette α only affects the hopping amplitudes in ydirection which pick up a complex phase factor.
While Eq. ( 13) describes a true bilayer system, by locally driving the internal transition and hence introducing a local coupling between the internal states on the x = 0 edge of the system, i.e. introducing a term Ĥexp = Ĥexp − t y â † 0,y,+ â0,y,− + H.c. , we can realize an effective single-layer model.Note that the additional edge term only consists of hopping contributions in x direction and thus, due to the gauge choice, does not contain complex phase factors.In this manner, an effective single-layer Hofstadter-Bose-Hubbard model of dimension 2L x × L y with homogeneous artificial magnetic field of α flux quanta per plaquette can be realized, see Fig. 1(a). 1  The main advantage of such a synthetic bilayer system is that it allows for a direct measurement of â † x,y,+ âx,y,− by using local tunneling pulses between the layers as follows.Switching on an additional coupling between the two internal states at site (x, y, ±) only and for a time corresponding to a π /2-pulse essentially rotates the bosonic operators to a new basis â± x,y = (â x,y,+ ± âx,y,− ) 1 Another option is to realize a true spatial bilayer system and to couple the edge by driving a Raman transition.
Now, applying a local detuning ∆ between the internal states [50] for some time τ allows us to furthermore rotate the measurement basis to where the relative phase φ − π/2 ∝ ∆τ can be controlled by the offset ∆.
Measuring the local occupation numbers with site and internal state resolution gives the density in the transformed basis n± x,y at the site where the pulse was applied, while at all other sites the original densities are recovered.We note that in the original basis the transformed occupation numbers translate to n± x,y = nx,y, where we introduced the generalized current operator ĵx,y (φ) = −ie iφ â † x,y,− âx,y,+ + H.c..
From Eq. ( 17) it follows immediately that That is, from occupation number measurements in the rotated basis, we can recover the occupation numbers in the initial basis, while simultaneously obtaining ĵx,y (φ) for any value of the offset phase φ resolving the absolute value of the coherence between the two sites via In general, exploiting the mapping of the bilayer system to an effective single-layer model, long-range correlations can be extracted using this protocol while reading out Fock basis snapshots on all remaining sites simultaneously.Consequently, the protocol proposed here furthermore allows for measurements of more complex correlators of the form In this case Eq. ( 20) generalizes to which especially contains all the necessary information to compute HODLRO as defined in Eq. ( 10) if we identify the function f (n xj ,yj ,σj ) with the phase contribution due to the flux attachment.Therefore, our protocol makes it feasible to extract this indicator of intrinsic topological order in state-of-the-art cold atom experiments.
The remainder of this paper will numerically explore the possibility to use our approach for the bosonic Laughlin state at magnetic filling ν = 1 /2.We find that already N snaps = 10 3 numerical snapshots (see Fig. 1(c)) of a system of 10 × 10 sites are sufficient to obtain not only a qualitative agreement with theoretical predictions from the continuum, but also an accurate quantitative estimate for the exponent of the power-law scaling of the HODLRO parameter (see Fig. 1(b)), giving direct insights about the topological order of the probed quantum state.

III. NUMERICAL ANALYSIS
In this section, we present an overview of the key results of the extensive numerical analysis of HODLRO we performed for the Hofstadter-Bose-Hubbard model.First, we introduce the explicit model and the parameters considered in our simulations.After pointing out key features of the ground states like incompressibility and the presence of a homogeneous density droplet in the bulk, we continue by confirming the absence of ordinary ODLRO.Finally, we present strong indications for the emergence of HODLRO in the composite boson one-particle correlations of the obtained ground states.

III.1. Model
We study interacting bosons on an L × L square lattice with open boundary conditions subject to a perpendicular magnetic field, modeled by the Hofstadter-Bose-Hubbard Hamiltonian Here, â( †) x,y are the bosonic annihilation (creation) operators and nx,y = â † x,y âx,y are the boson number operators.Analogously to the proposed experimental protocol, we perform our simulations using the Landau gauge.However all results discussed below, including the correlations indicating HODLRO, also apply to other gauge choices, see Appendix B. We fix the flux per plaquette to 2πα = 2π /6 in all our simulations.We study the hard-core bosonic limit, U /t → ∞, and consider dilute systems of few bosons N , with N /L 2 ≪ 1.Similar models have been explored in the literature for a while and it was established that the ground state of Eq. ( 23) considering the set of parameters listed above close to magnetic filling factor ν = N /N ϕ = 1 /2 is a lattice analog of the Laughlin state [7, 16-18, 24, 26, 32-34] where N ϕ is the number of flux quanta piercing the sample.Furthermore, it has been shown that a finite but strong Hubbard repulsion U /t is already sufficient to stabilize the Laughlin state.Therefore, we expect the results obtained in the hard-core limit U /t → ∞ to carry over to the experimentally relevant case of strong but finite interactions.Local particle number density ⟨nx,y⟩ in units of the flux density α for different system sizes and particle numbers.We find an extended bulk region for sufficiently large systems (L = 10, 12) exhibiting a bulk density consistent with the prediction n/α ≈ 1 /2 for the Laughlin state (indicated by the dotted line in the upper row).In the upper panel it is also demonstrated that these findings are robust to small changes in the particle number N , indicating the (bulk) incompressibility of the studied ground states.
We use the single-site DMRG method [37][38][39] implemented in the SyTen-toolkit [51] to find an MPS representation of the ground state of Eq. ( 23) with maximum bond dimension χ = 1000, where we exploit the U(1) symmetry associated with particle number conservation of the system.We consider different system sizes (L = 7, 10, 12) and different particle numbers N , however remaining in the dilute limit, N /L 2 ≪ 1.The magnetic length, and thus the correlation length of the Laughlin state, is given by ℓ B = a / √ 2πα ≈ a, where a is the lattice constant.
We calculate the local density ⟨n x,y ⟩ for the DMRG ground states and, for sufficiently large systems (L = 10, 12), identify a bulk region of density n ≈ να = α /2 as expected for a ν = 1 /2-Laughlin state, see Fig. 2 (lower panel).In particular, this behavior is robust for a range of particle numbers satisfying N /L 2 ≈ α /2, which can be interpreted as a signature of the incompressibility of the Laughlin state, see Fig. 2 (upper panel).For the smallest systems studied here (L = 7), the situation is less clear because of significant finite-size effects, and no extended bulk region is formed.However, earlier studies found evidence of an approximate Laughlin state even in systems of 4 × 4 sites [16], so that it might still be possible to find signatures of the topological nature in very small systems.

III.2. Absence of ODLRO
As discussed above, topologically ordered states of matter are prominently characterized by the absence of conventional long-range order, and in particular ODLRO.We confirm this to be true for the states under study here by analyzing the behavior of the normalized two-point correlation function along the x-direction ρ x,x ′ ;y = â † x,y âx ′ ,y / ⟨n x,y ⟩ ⟨n x ′ ,y ⟩. contracted MPS FIG. 3. Normalized two-point correlation function for the system sizes and particle numbers of Fig. 2 with the reference site on the edge, i.e. x ′ = 0.While we observe significant finite size effects, we find the characteristic exponential decay in all systems using a snapshot based scheme (circles, Nsnaps = 2 × 10 3 , 10 4 ) with an exponential fit (dotted line) in good agreement.
The errorbars of the snapshots correspond to a single standard deviation ± σρ x,0;y/ √ N .Furthermore, a direct contraction of the ground state MPS (diamonds) confirms these findings, in particular at short and intermediate length scales.The correlation length of the fitted decay is ξ ≃ 3a /2 which is significantly lager than the magnetic length lB ≃ a < ξ.
Analytical continuum calculations for the Laughlin states at general filling factor ν = 1 /m for a fixed value of y found a characteristic exponential decay at long distances, which is in particular independent of the filling factor [4].
In our numerical analysis we evaluate the expectation value for the ground state wave function, both using an exact MPS contraction, and analyzing Fock basis snapshots taken by the projective sampling method introduced below in Section IV.
To study long-range correlations we fix the reference site to be on the edge of the system, i.e. x ′ = 0. Furthermore, we additionally fix y = y ′ = y 0 close to the center of the system (y 0 ≈ L /2) to ensure that we eventually probe the bulk of the system.
The quasi-exact MPS expectation value for the oneparticle reduced density matrix ρ x,0;y0 of the ground state reproduces the continuum prediction of Eq. (25) for the system sizes L = 10 and L = 12 up to intermediate-range correlations, see Fig. 3 (green diamonds).As ℓ B ≈ a, we expect only short-range correlations to be affected by the competition of the characteristic length scales set by the magnetic and the lattice length.For long distances of the order |x − x ′ | > L /2, we observe evident signatures of finite-size effects weakening the decay.We believe the presence of an edge mode to be the reason for the deviation from the continuum prediction in this limit.In support of this, a probe of the correlator along the systems edge (see Appendix D) yields conclusive evidence for an algebraic decay characteristic for such an edge mode.For the smallest system, L = 7, finite size effects dominate and no clear signature of an exponential decay can be found.
We additionally analyze ρ x,0;y0 using the projective two-site sampling from MPS.We find that already N snaps = 2 × 10 3 snapshots are sufficient to qualitatively reproduce the exponential decay of ρ x,0;y0 on intermediate scales obtained via the MPS contraction, see Fig. 3 (light-blue circles with errorbars).The plotted values are the means of the sampling given with an errorbar corresponding to a single standard deviation ± σρ x,0;y/ √ N .In Appendix C we justify this choice and analyze the convergence behavior of the sampling to the expectation value in the limit N snaps → ∞ and evaluate its statistical accuracy.Using the data points of the sampling, we perform an exponential fit for all system sizes yielding a correlation length ξ ≈ 3a /2 which is slightly larger than the magnetic length ℓ B ≈ a.

III.3. Emergence of HODLRO
The central physical question addressed in this paper is the existence of HODLRO in fractional Chern insulators.As introduced in Section II, HODLRO can be understood as the algebraic long-range order of emergent objects.To achieve this, the original bosons are transformed to composite bosons via a singular gauge transfor- x,y bx ′ ,y , squares and triangles) for the system sizes and particle numbers of Fig. 2 with reference at the edge, x ′ = 0.That means, for ρ, f (nx j ,y j ) ≡ 1 while for ρ, f (nx j ,y j ) is to be identified with the phase contribution due to the flux attachment in Eq. (10).Already relatively few snapshots (Nsnaps = 2 × 10 3 , orange triangles) are sufficient to distinguish the algebraic decay of the correlations ρx,x ′ ;y of the composite bosons from the exponential decay of the ordinary correlations ρ x,x ′ ;y .Furthermore, for sufficiently large systems and Nsnaps = ×10 4 (red squares) a fit of the algebraic decay gives an exponent − 2 /m consistent with m = 2.The fit is obtained using only intermediate-range points within the bulk of the system.The errorbars for the sampling of both the bosonic (blue) and the composite bosonic (orange, red) two-point correlations correspond again to a single standard deviation ± σρ x,0;y/ √ N and ± σ ρx,0;y/ √ N respectively.mation [42,43,52], which corresponds to the attachment of flux quanta to the fundamental bosons.In the case of the Laughlin state at ν = 1 /m, each particle is endowed with exactly m flux quanta.The attached flux quanta cause the composite object to acquire an additional statistical angle of mπ.Hence, independent of the filling and the underlying particles the composite particles obtained by this procedure are always bosonic.
For varying x, x ′ and fixed y = y ′ = y 0 the composite boson correlations in the continuum Laughlin state follow a power-law scaling of the form indicating HODLRO in the Laughlin state [4,6].In contrast to the ordinary bosonic correlation, the transformed composite bosonic expression cannot readily be evaluated using direct contractions of the MPS representation of the ground state, because one explicitly needs to know the occupation numbers n k on all lattice points which do not coincide with j = (x, y 0 ) and l = (x ′ , y 0 ), while simultaneously determining the bosonic one-particle coherence on the two remaining sites.However, the projective two-site sampling algorithm we introduce in this paper grants access to exactly this information.Applying this algorithm, we can generate snapshots to evaluate the expectation value with in principle arbitrary accuracy.
Moreover, this scheme precisely emulates the measurement protocol for a quantum gas microscopy experiment as proposed in Section II above.
Analogously to the analysis of the bare correlator, we compute the normalized transformed, multi-particle twopoint correlation function ρx,x ′ ;y = b † x,y bx ′ ,y / ⟨n x,y ⟩ ⟨n x ′ ,y ⟩, for all system sizes with the reference site fixed at the center of one edge, x ′ = 0, y 0 ≈ L /2.Already for N snaps = 2 × 10 3 snapshots we find a clear difference between the exponential decay for the bare correlator ρ and the algebraic decay of the transformed expression ρ for all system sizes, see Fig. 4.This clearly indicates the presence of HODLRO in the ground state of our lattice model for a number of snapshots realistically achievable in cold atom experiments.Moreover, the emergence of the power-law scaling is robust to the choice of the reference site x ′ , which is shown in more detail in Appendix E.
Building on those qualitative observations, we perform a fit of the analytically predicted decay in Eq. ( 26) to our sampled data and find an excellent agreement of the exponent with the analytically predicted value for N snaps = 10 4 snapshots in the largest systems, L = 12.We find a power-law scaling with an exponent m L=12 = 2.01 ± 0.05 which agrees perfectly with the predicted value m = 2.
For the smaller systems, L = 7, 10, we find that we can still distinguish algebraic and exponential decay, however the fit to the analytical prediction becomes less accurate.
In the case of the smallest system studied here (L = 7), the situation gets even more complicated as an obvious bulk region is hardly identifiable and edge effects are believed to significantly influence the results.Interestingly, while the exponential decay is destroyed by finite size effects, the algebraic decay in the composite bosonic correlations is robust against them, signaling an underlying universal property of the state responsible for its emergence.We conclude that such experiments constitute a prime platform giving insights about the nature of the intrinsic topological order of the state by pure analysis of the ground state.

IV. METHODS
We now turn to the discussion of the sampling protocol to simultaneously generate one-and two-site snapshots from a MPS which was instrumental for our preceding numerical analysis of HODLRO.First, we provide the reader with a brief review of some features of MPS which are essential to our snapshot set-up and the idea to numerically emulate experimental cold atom quantum gas microscope measurements.Afterwards, as the singlesite/two-site algorithm is a straightforward generalization of the single-site sampling, we continue by shortly reviewing the idea of the single-site method.Building up on this, in the last subsection, we discuss the full singlesite/two-site variant of the algorithm.

IV.1. Matrix Product State Set-up
We consider a tensor product Hilbert space H of L sites with local basis states {|σ j,1 ⟩ , |σ j,2 ⟩ , ..., |σ j,p ⟩}, where p denotes the local physical dimension, so that the full Hilbert space is given by The wavefunction of an arbitrary pure state |ψ⟩ ∈ H can be represented by a product of matrices such that which yields the MPS representation of the state.Here, the site-local matrix M σj j is unique only up to an invertible linear map.We can exploit this gauge freedom to rewrite the site tensor in its left-canonical /right- Left-/right-canonical MPS are defined by requiring that they consist of left-/right-normalized tensors only.We can also construct a mixed-canonical MPS by fixing a site j for which the tensor M σj j remains unchanged while we demand all site tensors to the left/right to be left-/right-normalized tensors.We then call site j the active site.The canonical form of the MPS tensor-network is visualized in Fig. 5.The canonical form is especially advantageous if we want to compute the expectation value of an arbitrary site-local operator ôj (31) as the one-site reduced density matrix ρj becomes trivial and grants access to the full probability distribution for the eigenvalue spectrum of the operator given by the diagonal elements ρ oj oj .For a detailed discussion of MPS and its technicalities we remit to [38].

IV.2. Quantum Simulator Set-up
As we are interested in Fock basis snapshots, from now on we work in the occupation number basis and choose ôj = nj , but in principle the following scheme works for arbitrary operators.
Taking Fock basis snapshots of the pure state |ψ⟩ corresponds to simultaneously measuring the site-local occupation number operator nj on each lattice site.We can decompose each nj using local projection operators Pnj n 1 P (n1) = 1 via the diagonal elements of the one-site reduced density matrix ρ n 1 n 1 = ⟨n1|ρ1|n1⟩ = P (n1).Upon drawing a certain eigenvalue n1 from the probability distribution the state is projected in the corresponding site-local product state.By moving the active site of the projected state to the second site one can read out the conditioned probability distribution n 2 P (n2|n1) = 1 via ρ n 2 n 2 = ⟨n1, n2|ρ2|n1, n2⟩.Finally, one obtains a single-full lattice snapshot of the form (n1, n2, n3) from onto the eigenspaces of the corresponding measurement outcomes n j , where The global particle number operator N = L i=1 nj can be rewritten as and consequently, the measurement outcome of a single snapshot is given by a tuple (n 1 , . . ., n L ), which accordingly corresponds to a pattern of projectors of the initial state into site-local eigenspaces of the density operator.

IV.3. Sampling Single-Site Operators
Given the MPS representation of an arbitrary quantum state, we can emulate state-of-the-art quantum simulator measurements by drawing independent snapshots employing the perfect sampling scheme first introduced by Ferris and Vidal [40].Contrary to other sampling procedures such as Markov chain Monte Carlo sampling, this algorithm produces perfectly uncorrelated samples (and hence the name perfect sampling).This is computationally advantageous, as we do not need to account for additional equilibration and autocorrelation times.The perfect sampling scheme for a global lattice operator, which is decomposable in site-local observables is a straightforward successive application of the canonical form and tensor contractions which has been discussed in the framework of MPS in [41].
For a treatment of the ground state of the Hofstadter-Bose-Hubbard model, we also need to account for the U(1)-symmetry of the system associated with the conservation of the particle number N constraining tensor manipulations.To clarify this, we briefly review the key steps of the single-site algorithm here.
The idea of projectively sampling is best understood, if we first consider the expectation value of the global observable N on our lattice Hilbert space H, which we can express as Here N denotes the set of all p L possible outcome configuration of the outcome tuple ⃗ n ≡ (n 1 , ..., n L ) and Evidently, ⃗ n∈N P (⃗ n) = 1 and |N | = p L .In order to approximate ⟨ψ| N |ψ⟩, consider only a subset of configurations Ñ ⊂ N with | Ñ | < p L .Then upon normalizing with Z = ⃗ n∈ Ñ P (n) we can use the following approximation of the exact expectation value If we now suppose that the | Ñ | configurations of the subset Ñ were drawn randomly from all possible configurations |N | in accordance to the probability distribution P (⃗ n), we can rewrite this approximation as That is, we estimated ⟨ψ| N |ψ⟩ by means of | Ñ | independent samples from the random variable (P (⃗ n), E(⃗ n)).FIG. 7. Schematic sketch of the modification of the one-site algorithm towards the mixed one-site/ two-site sampling to evaluate expectation values of the form ⟨n1...nj−1 qj,j+1 nj+2...nL⟩.After j − 1 single-site steps we bring the projected MPS state |ψj−1⟩ in a canonical gauge where both site j and site j + 1 are the active sites (round blue nodes).We are now interested in the probability distribution q j,j+1 P (qj,j+1|nj−1, ..., n1) = 1 which is now accessible via the diagonal elements of the two-site reduced density matrix ρ q j,j+1 q j,j+1 j,j+1 = ⟨n1, ..., nj−1, qj,j+1|ρj,j+1|n1, ..., nj−1, qj,j+1⟩ = P (qj,j+1|nj−1, ..., n1).After drawing qj,j+1 in accordance to P (qj,j+1|nj−1, ..., n1) the scheme is continued by performing single-site measurements starting at site j + 2 until we obtain a full-lattice snapshot of the type (n1, ..., nj−1, qj,j+1, nj+2, ..., nL).
Instead of performing the full contraction in order to obtain Eq. ( 34), we can sample the observable N by drawing a finite number of samples from the probability distribution.Based on this result, constructing the probability distribution P is the key to the projective sampling algorithm.As elicited above, MPS are an ideal starting point to construct P , due to the isometric character of the site tensors which enable us to express the MPS in the canonical/unitary form giving access to the full probabilities through the one-site density matrices.
Hence, the first step of the algorithm is bringing |ψ⟩ in the right-canonical form Subsequently, we select the first site-local density operator n1 and break it down in its p projector constituents On the tensor-network level, this can be achieved following Ref.[53] by manipulating the local rank 4 site tensor of the density operator in accordance with the constraints given by particle number conservation until we have a compactified rank 2 tensor.The rank 2 tensor is a Hermitian matrix and thus, we can apply an eigenspace decomposition.
We obtain the probability P (n 1 ) to draw n 1 via where |n 1 ⟩ ≡ Pn1 1 |ψ⟩.Doing this for all p possible measurement outcomes yields the full probability distribution satisfying n1 P (n 1 ) = 1. ( We randomly choose one of the p eigenvalues representing the measurement outcome in accordance to their probability distribution P (n 1 ).
Finally, we project the MPS site-tensor in the drawn site-local product-state and normalize it by means of a singular value decomposition in order to shift the canonical center towards the second site This one-site sampling step is now repeated for each of the L lattice sites.The initial state |ψ k ⟩ of the kth sampling step is given by The sampling probabilities at each step are conditioned by all already projected site tensors and therefore the final probability distribution P (n 1 , ..., n L ) ≡ P (⃗ n) is given by ) Consequently, by sweeping over the full lattice we obtain a single snapshot-tuple (n 1 , ..., n L ) drawn according to P (⃗ n).The whole perfect sampling scheme for MPS is illustrated in Fig. 6.

IV.4. Sampling Two-Site Operators
In the context of revealing HODLRO we are interested in sampling an operator expectation value of the type â † j âl k̸ =j,l nk (45) where we adopted the notation of Sec.II, i.e. a bold letter corresponds to a discrete coordinate tuple (x, y).
Eq. ( 45) underlines that we want to simultaneously sample Fock basis snapshots on all but two sites, and on the two remaining sites we are interested in the bosonic oneparticle correlation function.Such operators are not accessible through the introduced one-site perfect sampling scheme, however, by generalizing the method in fact any kind of multi-site observable can be probed.
Here we will present the modification needed in order to sample operators of the type of Eq. ( 45), i.e. we need to additionally extract probabilities for eigenvalues of non-local two-site operators.Analogously to the single-site case, these probabilities are contained within the two-site reduced density matrix ρj,l = Tr k̸ =j,l ρ. (46) We can simplify sampling Eq. ( 46) by exploiting the invariance of the trace under permutations P , which are involutory, i.e. fulfill P P = 1.Hence, in order to sample an arbitrary non-local two-site observable qj,l with |j − l| > 1, we can use the permutation P to change the ordering of the sites {j, l} → {j ′ , l ′ }, such that |j ′ − l ′ | = 1.From a MPS perspective this is beneficial, as we do not need to account for non-local operations.That means, we can safely assume that we are interested in sampling the following set of local operators (n 1 , ..., nj−1 , qj,j+1 , nj+2 , ..., nL ) where the coordinate arithmetic j ± 1 corresponds to the two adjacent lattice points of site j in an arbitrary mapping of the two-dimensional lattice to a one-dimensional chain.In this way, up until to the jth sampling step we can follow the single-site procedure, i.e. we begin by sampling a pattern of site-local occupation number eigenvalues of the form (n 1 , ..., n j−1 ) in accordance to the probability distribution P (n 1 , ..., n j−1 ).
Next, we take the projected state |ψ j−1 ⟩ and bring it in a mixed canonical form, but this time both site j and j + 1 form the canonical center.This enables us to trivially compute the two-site reduced density matrix granting access to the conditioned sampling probabilities.
The next step is to spectrally decompose qj,j+1 .On the tensor network level qj,j+1 is now a rank 6 tensor which increases the complexity in fusing it down to its compactified rank 2 form yielding qj,j+1 = q j,j+1 q j,j+1 Pq j,j+1 j,j+1 .
The two-site step is completed by the projection on the corresponding two-site product state: and ultimately, by moving the canonical center towards the site j + 2. The extended perfect sampling algorithm in its tensor network notation is shown in Fig. 7.
From here onward, the procedure is just the plain one-site algorithm again until we finally obtain (n 1 , ..., n j−1 , q j,j+1 , n j+2 , ..., n L ).How a snapshot of this form can be used to probe HODLRO is discussed in detail in Appendix F and G.

V. CONCLUSION AND OUTLOOK
The main result of our work is the demonstration that state-of-the-art quantum simulation platforms, and cold atoms in particular, can be used to directly detect HODLRO in fractional Chern insulators.This finding is based on strong numerical evidence for the existence of HODLRO in a lattice analog of the ν = 1 /2-Laughlin state, which has far-reaching consequences for topologically ordered lattice systems.We explicitly generalized existing continuum results to lattice systems and described how correlations of the emergent composite bosons give rise to HODLRO in this context.Such correlators are of interest to probe the condensation of emergent composite bosons and to deepen our understanding of the intrinsic topological order in the quantum state.
In view of the current abilities of cold atom quantum simulators, we proposed an experimentally accessible scheme to extract long-range correlations of the form ⟨â † j âl k̸ =j,l nk ⟩.Our proposal is based on realizing a bilayer system which is coupled along one edge to realize an effective extended single-layer Hofstadter-Bose-Hubbard model.This allows to measure non-local coherences ∼ â † j âl in current setups.Our numerical results are based on the extension of the existing perfect sampling scheme [40] to allow to projectively sample more complex, multi-particle correlation functions like HODLRO.The application of this sampling algorithm allowed us to probe HODLRO, thus demonstrating that this concept persists in lattice systems accessible to near-term quantum simulators.In particular, we showed that already relatively few snapshots (N snaps = O(10 3 )) are sufficient to distinguish the exponentially decaying correlations of the underlying bosons from the quasi-long-ranged power-law decay of the composite boson correlations.Furthermore, we were able to resolve the exponent of the algebraic decay which is directly related to the filling factor of the Laughlin state and gives direct qualitative insights about the intrinsic topological order of the state, specifically its K-matrix.
The recent progress in studying FQH states in cold atom experiments [16] calls for additional ways to probe these states directly.The snapshot-based protocol discussed here provides a way to explore the unconventional correlators needed, for example, to reveal the intrinsic topological order of such states.Having revealed the exotic correlations present in such systems, a next step would be further investigations of their origin.Mi-croscopically, computing the entire one-particle reduced density matrix for the composite bosons employing the snapshot protocol could give further insights into the speculated condensation of the composite bosons in the Laughlin state.Another intriguing possibility is to explore the fate of HODLRO in systems out-of equilibrium.While we restricted our analysis to a particularly simple FQH state, we believe that the concepts introduced here are also applicable to more exotic FQH states, like for example the Pfaffian state [54] or parafermion states [55].Furthermore, applying similar ideas to other systems exhibiting topological order, like chiral spin liquids [56], might deepen our understanding of such exotic states of matter.These more complex states should be in principle accessible to quantum simulator platforms like ultracold atoms or superconducting qubits where our approach can be readily applied.
We introduced the Hofstadter-Bose-Hubbard model in Landau gauge and hence, all calculations have been carried out using this particular gauge choice.It is straightforward to show that the emergence of HODLRO found in Section III is independent of the gauge.Consider the original two-point correlation function ρ x,x ′ ;y,y ′ = â † x,y âx ′ ,y ′ .(B1) Under a gauge transformation on the physical vector potential A of the form the bosonic annihilation (creation) operators transform as â( †) x,y −→ e (−)iχx,y â( †) x,y .(B3) For instance going from Landau gauge A x,y;L = 2πα (−y, 0) to the symmetric gauge A x,y;S = πα (−y, x) corresponds to a gauge transformation χ x,y = −παxy.
Consequently, the correlation function of Eq.B1 simply picks up an additional local phase under general gauge transformations χ x,y maintaining a homogeneous magnetic field ρ x,x ′ ;y,y ′ −→ e −i(χx,y−χ x ′ ,y ′ ) ρ x,x ′ ;y,y ′ , (B4) which in particular implies that the absolute value of ρ is preserved.The singular lattice gauge transformation which promotes the bosonic operators to composite bosonic operators on the other hand had the form â( †) x,y −→ e (−)i Φ( †) x,y â( †) x,y .
Since Φx,y directly depends on the site-local density operator nm,n where (m, n) ̸ = (x, y) we cannot treat it as a pure phase which we can pull out of the ensemble average However, as e (−)iχx,y is just a complex number it commutes with any Φx,y and it follows that the composite bosonic correlation function ρx,x ′ ;y,y ′ = −i Φ † x,y â † x,y âx ′ ,y ′ e i Φx ′ ,y ′ , (B6) also simply acquires a trivial phase under gauge transformations on the physical vector potential ρx,x ′ ;y,y ′ −→ e −i(χx,y−χ x ′ ,y ′ ) ρx,x ′ ;y,y ′ .studied here is supposed to exhibit a chiral Luttinger liquid edge mode with a phase field operator that is expected to have considerable overlap with the creation and annihilation operators on the lattice.The existence of a chiral edge mode can be probed by considering two-point correlation functions along the edge of the system.While the two-point correlations are expected to decay exponentially deep in the bulk of the system, the presence of an edge mode weakens the falloff moving closer to the edge.Measuring correlations explicitly along the edge one should find an algebraic decay signaling the presence of the gapless edge mode.Additionally, in systems with finite size (which we are considering in our studies) one would also expect the revival of the correlation function once we approach the opposite edge of the system [25].
In Fig. 9 we study the behavior of the two-point correlation functions along the x-direction for varying values of y while keeping the reference site fixed on the edge, i.e. x ′ = 0. We find that for values of y ≲ 2 close to the edge the decay is indeed suppressed resembling signatures of a power-law decay.For values of y ≃ L/2 deep in the bulk the correlations are vanishing exponentially.Furthermore, the revival of the correlation function can be observed for all values of y.Since the edge mode has a shorter distance to cover between x = 0 and x = L we find more pronounced revivals for smaller y.Hence, we find conclusive evidence for the coexistence of a gapless edge mode living on the boundary of the system and a gapped bulk.considerations within the errorbars for x ′ = 0.As discussed in the results in Section III, points x close to x ′ and points x close to the edge might be effected by length scale, finite size and edge effects and are not considered for the fit.Hence, the bigger x ′ the fewer data points can be taken into account and the fit quality becomes less reliable.

FIG. 1 .
FIG. 1.(a) 1 ○ Possible experimental realization of the proposed measurement protocol using a synthetic bilayer Hofstadter-Hubbard system with internal states |σ = ±⟩.2○ On one side of the system, the edges are coupled to obtain an effective single FIG.2.Local particle number density ⟨nx,y⟩ in units of the flux density α for different system sizes and particle numbers.We find an extended bulk region for sufficiently large systems (L = 10, 12) exhibiting a bulk density consistent with the prediction n/α ≈ 1 /2 for the Laughlin state (indicated by the dotted line in the upper row).In the upper panel it is also demonstrated that these findings are robust to small changes in the particle number N , indicating the (bulk) incompressibility of the studied ground states.

L = 12 , N = 10 ρx,0; 5 Nsnaps = 10 4 ρx,0; 5 Nsnaps = 2 × 10 3 ρx,0; 5
FIG. 4.Decay of correlations for the original bosons ( |ρ x,x ′ ;y | ∝ â †x,y âx ′ ,y , circles) and the composite bosons(|ρ x,x ′ ;y | ∝ b †x,y bx ′ ,y , squares and triangles) for the system sizes and particle numbers of Fig.2with reference at the edge, x ′ = 0.That means, for ρ, f (nx j ,y j ) ≡ 1 while for ρ, f (nx j ,y j ) is to be identified with the phase contribution due to the flux attachment in Eq.(10).Already relatively few snapshots (Nsnaps = 2 × 10 3 , orange triangles) are sufficient to distinguish the algebraic decay of the correlations ρx,x ′ ;y of the composite bosons from the exponential decay of the ordinary correlations ρ x,x ′ ;y .Furthermore, for sufficiently large systems and Nsnaps = ×10 4 (red squares) a fit of the algebraic decay gives an exponent − 2 /m consistent with m = 2.The fit is obtained using only intermediate-range points within the bulk of the system.The errorbars for the sampling of both the bosonic (blue) and the composite bosonic (orange, red) two-point correlations correspond again to a single standard deviation ± σρ x,0;y/

jFIG. 5 .FIG. 6 .
FIG. 5. (a) Schematic sketch of the tensor-network of the right-canonical MPS representation of an arbitrary state |ψ⟩ ∈ H.(b) Now the state is shown in its mixed-canonical form where all sites to the left/ right of the active site are in the left-/ right-canonical gauge (blue left-/ right-orientated triangle nodes).(c) The one-site reduced density matrix ρj in its tensor network graph representation.Due to the mixedcanonical form of the state the tensor contraction left/ right to the active site (round blue node) will give the identity (lightorange boxes).ρj gives access to the full site-local probability distribution.
gauge transformations is only unique up to a translation on the connected component covered by the exponential map of the group U(1), we can choose the extension of the singular gauge transformation to be Φx,y − χ x,y means, not only the norm is preserved but generally the emergence of HODLRO.

FIG. 9 .
FIG. 9. One-particle correlation function along the x-direction.Clear signatures of a slower decay close to the edge, signaling the presence of an edge current.