Quantum Virtual Cooling

We propose a quantum information based scheme to reduce the temperature of quantum many-body systems, and access regimes beyond the current capability of conventional cooling techniques. We show that collective measurements on multiple copies of a system at finite temperature can simulate measurements of the same system at a lower temperature. This idea is illustrated for the example of ultracold atoms in optical lattices, where controlled tunnel coupling and quantum gas microscopy can be naturally combined to realize the required collective measurements to access a lower, virtual temperature. Our protocol is experimentally implemented for a Bose-Hubbard model on up to 12 sites, and we successfully extract expectation values of observables at half the temperature of the physical system. Additionally, we present related techniques that enable the extraction of zero-temperature states directly.

In this work, we develop a novel approach to address this issue by introducing a measurement scheme that enables to access system properties at fractions of its temperature T . This "virtual" cooling to a temperature T virtual = T /n (n = 2, 3, . . . ) is facilitated by joint measurements on n copies of the system. A schematic is given in Fig. 1(a). Further, we detail implementations tailored = FIG. 1: (a) Schematic representation of the virtual cooling protocol. Collective measurements on two copies of a thermal state ρ β at temperature T = 1/(kBβ) correspond to standard measurements at half the temperature, T /2. (b) Diagrammatic representation. Two copies are evolved with the unitary F2, and a subsequent measurement of X2 is performed. In combination this gives the expectation value tr{ρ 2β X} corresponding to half the original temperature.
to cold-atom systems in optical lattices, and illustrate our protocol in an experimental quantum simulation of the Bose-Hubbard model. Finally, we show how these ideas can be generalized and discuss protocols to distill the many-body ground state from multiple copies of thermal many-body states.
We are interested in quantum many-body systems described by a thermal state ρ(T ) = e −βH /Z, where H is the Hamiltonian of the system and Z(T ) = tr{e −βH } is the partition function at inverse temperature β = 1/(k B T ). The measurement of an observable X in the state ρ gives the expectation value X T = tr {Xρ}. Be-low we will discuss a protocol that allows us to effectively measure X T /n . The central idea is based on the ability to express the thermal density operator at T /n by the n-th power of ρ(T ) ρ(T /n) = ρ(T ) n /tr{ρ(T ) n }. (1) In order to access the higher powers of the thermal state, we require n copies of the state ρ(T ) prepared in parallel as well as the capability to implement operations that exchange the n copies. More specifically, we have tr {Xρ n } = tr {X s S n ρ ⊗n }, where S n cyclically permutes quantum states in the n copies, i.e. S n |ψ 1 ⊗ |ψ 2 ⊗ · · · ⊗ |ψ n = |ψ 2 ⊗ |ψ 3 ⊗ · · · ⊗ |ψ 1 , and X s is the symmetrized embedding of X on the n-fold replicated Hilbert space X s = 1 n n m=1 S m n (X ⊗ 1 ⊗(n−1) )S m n † . Therefore, the virtual measurement of X T /n at temperature T /n via Eqn. (1) can be reduced to determining the expectation values X s S n and S n on the n copies of the state at temperature T . This is illustrated in Fig. 1(b). Measurements of expectation values of S n can be achieved with auxiliary qubits [38,39], or directly via many-body state interferometry [40][41][42], as recently demonstrated with cold atoms [27]. We also note that our protocols apply to subsystems which are locally thermally, even if the global system is not thermal. In our experiments below, we leverage 'eigenstate thermalization' [43][44][45][46][47] to obtain thermal reduced density matrices from globally pure states of finite energy density in a chaotic system. Earlier theoretical work provided numerical evidence that a chaotic eigenstate or a reduced density matrix of a thermal state encodes correlations at all temperatures [48,49].
Below, we discuss protocols to measure X s S n for arbitrary n and detail the procedure for the simplest example n = 2. We first focus on an interferometric measurement scheme and demonstrate that it can be implemented in current experiments with cold atoms. Alternative virtual cooling schemes using ancillary atoms are discussed in the Supplementary Materials. Finally, we show that schemes with ancillary atoms can be generalized to not only virtually cool a many-body system, but directly distill and prepare the many-body ground state from a thermal state. Importantly, all of the discussed protocols are agnostic to the temperature T of the physical system, and thus can be used to obtain additional, virtual cooling even after all available physical cooling methods have been deployed.

INTERFEROMETRIC MEASUREMENT
To simplify the presentation we first discuss a virtual cooling scheme for bosonic atoms in optical lattices. The key idea is to represent the permutation operator S n in the bosonic Hilbert space as S n = F n R n F † n , where the unitary F n denotes the discrete Fourier transformation and R n = j e −i2π/n n p=1 p np,j [56]. Here a p,j denotes the bosonic annihilation operator on site j in copy p, and n p,j = a † p,j a p,j is the corresponding number operator. Note that F n can be realized by simply introducing tunnel coupling between neighboring copies [41], and R n can be directly measured with a number-resolving quantum gas microscope. This representation of the permutation operator suggests that we introduce an operator X n = F n X s F † n , which is the discrete Fourier transform of the observable X that we want to measure. With this definition we can express A measurement of X at the virtually reduced temperature T /n thus consists of a measurement of X n R n and R n after application of the discrete Fourier transform across the copies. For many interesting observables one finds [X n , R n ] = 0 so that R n and X n can be measured independently.
As a specific example, we consider the experimentally simplest case n = 2 and the measurement of the on-site density by choosing X ≡ n j . The corresponding protocol consists of three steps. (i) We prepare n = 2 identical instances of the thermal many-body state ρ(T ). This can be achieved, for example, by preparing two identical states in neighboring 1D tubes, or 2D planes. It is essential that the copies are decoupled at this stage, which can be achieved by using a large optical potential between the tubes or planes to suppress any inter-copy tunneling. (ii) We then freeze the dynamics within each copy, and lower the potential between the two copies, e.g. using an optical superlattice. This induces tunneling between the two copies via the Hamiltonian H BS = −J BS j (a † 1,j a 2,j + h.c.), which allows us to realize the so-called beamsplitter operation F 2 that maps ρ ⊗2 → F 2 ρ ⊗2 F † 2 . Interactions between the atoms need to be turned off (e.g. via a Feshbach resonance) or made negligible as compared to J BS during this step. (iii) Finally, we measure the on-site occupation number on all sites in both copies using a number-resolving quantum gas microscope. This gives direct access to R 2 = (−1) j n1,j and X 2 = F 2 1 2 (n 1,j +n 2,j )F † 2 = 1 2 (n 1,j +n 2,j ). Averaging the results over multiple experiments gives the expectation value of the local density at T /2 via Eqn. (3) (for a schematic of a single measurement trial, see Fig.  2(a) below). Remarkably, this experimental procedure parallels the one employed to determine the second order Rényi entropy of cold atoms, with atom numberresolved measurements being the only additional requirement. Such measurements were first demonstrated for one-dimensional systems using full-atom-number-resoved imaging in quantum gas microscope [28].

EXPERIMENTAL DEMONSTRATION
In order to demonstrate our protocol, we experimentally realize it in a one-dimensional Bose-Hubbard model. In the experiment, a Bose-Einstein condensate of 87 Rb atoms is loaded into a two-dimensional optical lattice positioned at the focus of a high-resolution imaging system. The dynamics of the atoms is well-described by a Bose-Hubbard Hamiltonian parametrized by tunneling strength J and on-site interaction energy U (see Ref. [28] for details).
The experimental protocol consists of four steps: initialization, quenched thermalization dynamics, beamsplitter operations, and measurements. During initialization, optical potentials are sequentially manipulated in order to isolate an initial product state, |ψ 0 , with a single atom on the central 2 × 6 sites of a 2 × L plaquette in the deep 45E r lattice where the tunneling between the sites is negligible [28]. Each 1 × L tube represents an identical copy of the system. Next, the lattice potential along the chains is suddenly lowered, allowing particles to tunnel and interact within each chain. It has been previously shown [28] that this quenched dynamics drives the thermalization of small subsystems within the chain. Hence, after sufficiently long time evolution, the state of the subsystem can be described by an effective temperature T and chemical potential µ, which are determined by the total energy and particle number density of |ψ 0 . (See also [49].) After the desired time evolution the dynamics of the system is frozen by suddenly increasing the lattice depth along the chains, and a beamsplitter operation F 2 is implemented by lowering the potential barrier between the two chains, such that particles can tunnel (in the transverse direction) for a prescribed time. Finally, the number of particles on each individual lattice site is measured. This procedure is repeated multiple times in order to obtain sufficient statistics.
We apply our virtual cooling protocol in three regimes (A, B, and C), with differing initial states |ψ 0 , system size L, and Hamiltonian parameters U/J. For the data sets A and B, each of L = 6 sites is initially occupied by one particle, whereas for the data set C, only the middle six out of the total L = 12 sites are occupied by one particle per site. The tunneling rates are set such that U/J ≈ 1.56 (data set A) or 0.33 (data sets B and C). These combinations lead to the effective temperatures and chemical potentials (T /J, µ/J) ≈ (3.5, −1.0), (11.5, −6.3), and (18.3, −17.7) of subsystems for data sets A, B, and C, respectively. Based on our protocol, we extract the average particle number density n i for thermal ensembles at reduced temperature. Fig. 2(b) shows the resulting single-site particle density Schematic for a single measurement trial of X2 = 1 2 (n1,i + n2,i) and R2 restricted to the ith site, after F2 has been applied to the two copies. (b) Measured singlesite density, averaged over all but the edge sites of the chain, after virtual cooling has been applied to the system (blue circles with vertical error bars). Red discs show the single-site density of the state before our protocols are utilized (the actual density of particles in each experiment), whereas light blue discs correspond to the prediction of the effective thermal ensemble (see text) at half the temperature. Agreement of the data with the reduced-temperature ensemble validates the applicability of our method in the experimental system. Error bars denote the standard error of the mean.
after virtual cooling for all three cases. We compare these results with the initial single-site density at the original temperatures as well as theoretical predictions from an ideal thermal ensemble ρ 2β at half of the original temperatures. All data points are in good agreement with the reduced temperature ensemble indicating that our virtual cooling scheme works in the experimental system.

OBSERVABLES
In our experiment, we measured the single-site density of bosons in an optical lattice at half the physical temperature. If we measure the single-site density on the jth site, we have X ≡ n j , X s = 1 2 (n 1,j + n 2,j ), and X 2 = F 2 1 2 (n 1,j + n 2,j )F † 2 = 1 2 (n 1,j + n 2,j ). Notice that X s = X 2 , since each operator is invariant under the many-body Fourier transform F 2 between the two system copies. Here, X 2 is easily measured by averaging the number of atoms on the jth site in the two copies. Furthermore, X 2 commutes with R 2 , and so we can measure the observables in either order. In fact, X 2 and R 2 commute with all single-site densities n 1,j , n 2,k , and so we can simply measure the individual particle numbers and combine them to compute the expectation values of X 2 and R 2 .
For more complicated observables such as densitydensity correlators X ≡ n j n , the situation is more subtle. A direct application of the procedure outlined above requires a measurement of While the first term in Eqn. (4) (i.e., the final equality) is easily measurable with standard quantum gas microscopy, the second term requires additional interferometric apparatus. However, the first term of Eqn. (4) by itself contains interesting information about the system at half of its temperature. This first term of Eqn. (4) is easy to measure, since it commutes with R 2 and all of the number operators. Doing so would output the unconventional correlator The term on the left here is the desired equal-time density-density correlator at half the system temperature, whereas the term on the right is peculiar. In fact, this peculiar term is equal to the unequal imaginarytime correlator 1 2 tr{n j (1/T ) n ρ(T /2)} where n j (τ ) = e Hτ n j e −Hτ is the number density evolved in imaginary time. If our system is translation invariant and at sufficiently low temperature, we expect tr{n j n ρ(T /2)} to depend on |j − |, whereas the peculiar term should not strongly depend on |j − |. This is because at low temperatures, the large imaginary time evolution of the operator n j scrambles it strongly, destroying the memory of its initial position j. Indeed, in the limit of T → 0, the peculiar term is just ψ 0 |n j |ψ 0 ψ 0 |n |ψ 0 which is clearly independent of |j − |. At high temperature and small |j − |, both terms in Eqn. (5) have a nontrivial dependence on |j − | and so we are unable to extract each term separately. Nevertheless, it is interesting to note that our protocol yields some information about the unequal imaginary-time correlator in this regime. We note also that when |j − | is much larger than the thermal correlation length, both terms in Eqn. (5) approach 1 2 tr{n j ρ(T /2)} tr{n ρ(T /2)}. We explore the dependence of tr{n j ρ(T ) n ρ(T )} on |j − | as a function of T in the Supplementary Materials, and confirm that there is essentially no dependence at sufficiently low temperatures.
Another approach to measuring density-density correlators is to use a different protocol involving an ancillary qubit to implement a controlled swap operation (see Supplementary Materials). In this protocol, one does not need to measure R 2 by observing the particle number distribution with a quantum gas microscope, and so we have more flexibility in our ability to measure X 2 . Another interesting observable which is easy to measure in the ancillary qubit setting is the hopping operator X ≡ a j a † +h.c., and accordingly X s = 1 2 (a 1,j a † 1, + a 2,j a † 2, ) + h.c. which satisfies X 2 = X s .
We are often interested in local observables X, which in turn correspond to the local observables X s . Suppose that X is supported on a subregion R. Then X s is supported on the joint region R 1 ∪ R 2 of the two corresponding system copies. For concreteness, suppose our system is one-dimensional. We desire to measure where ρ R (T /2) = tr R {ρ(T /2)} is the reduced density matrix of ρ(T /2) on R. Naïvely, it seems that we only need to perform our procedure on the subsystem R 1 ∪ R 2 of the two copies. However, this is not correct, since Nonetheless, suppose we extend R by buffering each of its boundaries by a number of sites corresponding to the correlation length of the system at temperature T /2. Let us denote this extended region by B. Here, R ⊂ B, but B is smaller than the whole system. The corresponding joint region of the two system copies is B 1 ∪ B 2 . If we perform our procedure on the subsystem B 1 ∪ B 2 of the two copies, we can access the state σ B ≡ ρ 2 B /tr B {ρ 2 B }, which satisfies σ B ≈ tr R {ρ(T /2)}, and therefore So if we choose B large enough (but in most cases, smaller than the size of the entire system), we can still approximately measure our desired observable.
In an experiment, the performance of a measurement protocol is limited by the number of repetitions required to achieve sufficiently high precision. In our setting, the measurement statistics required to precisely measure the denominator Z n = tr{ρ(T ) n } in Eqn. (1) may be a limiting factor. In a many-body system, Z n is directly related to the Rényi-n entropy S n = 1 1−n log(Z n ), which scales with volume for local systems. Z n is therefore often exponentially small in the system size. Hence, one would generally need a large number of measurements N m ∼ 1/Z 2 n ∼ exp{2 s(T )|R|}, where s(T ) is the entropy density at temperature T and |R| is the size of the subregion on which ρ(T ) is supported. In the limit of low temperature, this scaling becomes favorable since s(T ) generally decreases. However, the thermal correlation length ξ(T /n) can increase as T is lowered, requiring a larger subregion size |R| ≥ ξ(T /n). Together, the number of measurements required to achieve some fixed precision scales as N m ∼ exp{2 s(T ) ξ(T /n)}.
Of course, if ρ is only approximately thermal, then expectation values of ρ n /tr{ρ n } for larger values of n can have amplified deviations from thermality. However, if we are interested in the physics of the ground state |ψ 0 , then ρ n /tr{ρ n } ∼ |ψ 0 ψ 0 | for larger values of n so long as |ψ 0 is the dominant eigenstate of ρ. We discuss a related idea for ground state distillation in the following section.

GROUND STATE DISTILLATION
The above ideas can be generalized to schemes that not only allow us to measure a system at reduced temperatures, but further enable the distillation of the ground state from multiple copies of a thermal ensemble. This is akin to entanglement purification proposals for quantum communication over noisy channels [50]. Consider a non-destructive measurement of the swap operator, S 2 , on two systems which are each prepared in the state ρ. Since S 2 is unitary and hermitian, the two possible measurement outcomes are ±1, corresponding to projections into the symmetric or anti-symmetric subspace with respect to the exchange of the two copies. The state after such a measurement is thus given by P ± (ρ⊗ρ)P ± /tr {P ± (ρ⊗ρ)}, with P ± = (1 ± S 2 )/2. If the measurement outcome is −1, we discard both systems. But for those instances that yield a measurement +1 we retain one of the systems, and discard only the other one. The resulting state of this first system ρ 1 is obtained by tracing out the degrees of freedom of the second system, For an initial thermal state ρ(T ), the new state ρ 1 corresponds to a mixture of ρ(T ) and ρ(T /2). Clearly, ρ 1 has the same eigenvectors as ρ, but with different eigenvalues. In particular, ρ 1 is purer than ρ, and the eigenvalue of the largest eigenvector (i.e. the ground state for thermal ρ) is larger. This purification is of course probabilistic, as its success is conditioned on the proper measurement outcome for S 2 . Remarkably, the success probability p + = (1 + tr ρ 2 )/2 is always larger than 1/2 and approaches 1 as the system is purified. Starting with multiple copies one can iterate the above process, which will ultimately converge to a system in the largest eigenstate of ρ. For thermal states, the procedure distills the many-body ground state, i.e. the zero-temperature state. In a cold atom setup, non-destructive measurements of the swap operator are typically more challenging to implement than the direct measurements which form the basis of the virtual cooling scheme discussed earlier. One way to realize such measurements is to use ancillary qubits [51]. For example, one can envision encoding a qubit in two internal states of an ancillary atom. Employing the Rydberg blockade mechanism, one can then use this qubit, e.g. to control the tunneling amplitude between two copies in optical lattices and so realize a controlled exchange operation (see Fig. 3 and [51]).

DISCUSSION
Reaching low temperatures is paramount for studying interesting quantum many-body phases with quantum simulators. In particular, the small energy scales in cold atom systems pose a major challenge for accessing the required temperature regimes. In this work, we proposed and demonstrated novel techniques that enable access to properties of a system at a fraction of its actual temperature. This virtual cooling is enabled by collective measurements on multiple copies of the system.
More generally, our schemes illustrate a connection between thermal physics and entanglement. In particular, the temperature of a system is intimately connected to its entanglement with its surroundings [28,44,45,49,52,53]. Accordingly, measuring correlations of a thermal system at virtually lower temperatures involves manipulating and probing entanglement. This is why the tools for measuring a system at virtually lower temperatures resemble those that allow access to entanglement entropies [27,40,41].
A natural future direction is to experimentally perform quantum virtual cooling for more complicated observables. A particularly interesting application would be to experimentally study a quantum many-body system with a finite-temperature phase transition at some temperature T c . One could prepare the system at some temperature T > T c , and use virtual cooling to probe features at or below the phase transition. (For related theoretical work, see [55].) Understanding the range of applicability of quantum virtual cooling is an exciting theoretical and experimental program, which will require new insights in subsystem ETH and thermalization. We thank Alex Avdoshkin for helpful conversations. JC is supported by the Fannie and John Hertz Foundation and  Here we present detailed quantum virtual cooling schemes, including alternative ones that do not appear in the main text. We analyze the case of two system copies, so that quantum virtual cooling allows us to probe observables at half of the physical temperature. We start with a more general method involving an ancilla qubit, and then specialize to bosons and fermions in optical lattices.
Using an ancilla qubit Suppose we have two identical systems, each in the same state ρ and with Hilbert spaces H, and an ancilla qubit with Hilbert space C 2 . The joint Hilbert space of the entire system is C 2 ⊗ H 1 ⊗ H 2 with H H 1 H 2 . We define the ancilla states By jointly manipulating the ancillas and two identical systems, we will show how to measure tr{X ρ 2 }/tr{ρ 2 } for any observable X.
2. Apply to the joint system a controlled-(U ⊗ 1) gate for some unitary U on the first system copy, followed by a controlled-S 2 gate where the control qubit is the ancilla.
3. Measure the ancilla qubit in the {|+ , |− } and {|L , |R } bases, respectively, to obtain the measurement probabilities: These probabilities can be used to reconstruct tr{U ρ 2 } for any chosen unitary U .
4. Let {U i } be an orthogonal basis (i.e., tr{U † i U j } = dim(H 2 ) · δ ij ) for B(H 2 ). For instance, if the two system copies are comprised of qubits or qudits, one could use strings of Pauli operators or generalized Pauli operators. Then using the first three steps above, for any operator X ∈ B(H 2 ), one can compute tr{X ρ 2 }/tr{ρ 2 } using the relation The last step may seems daunting, since the sum in Eqn. (11) contains dim B(H 2 ) = (dim H 2 ) 2 terms. However, a good choice of {U i } renders only a few terms in the sum non-zero. For instance, suppose that the two identical systems are each spin chains of qubits, and that X is a product two Pauli operators, each on a separate site. If {U i } is chosen to be products of Pauli operators, then only one term in the sum would be non-zero (namely the term for which U † i = X).

Boson interferometry
If our two identical systems are bosonic, then we can perform quantum virtual cooling along the lines of the main text. In particular, we do not need an ancilla qubit to facilitate the application of the swap operator. Consider the bosonic Hilbert space Sym(H 1 ⊗ H 2 ), comprising of two systems with N sites each. A basis for Sym (H 1 ⊗ H 2 ) is for (2), F 2 is a unitary which maps Furthermore, R 2 = (−1) j n1,j is the total parity operator for the first of the two identical systems. It is easy to check that and so Then if we have an operator X that we wish to measure, the idea is to instead measure X 2 = F 2 1 2 (X ⊗ 1 + 1 ⊗ X)F † 2 so that, in essence, Of course, there is a detailed measurement procedure which realizes the above equations.

Apply
3. Measure the operator X 2 , given by Here, X({a 1,i , a † 1,i }) denotes that the operator is written in terms of sums of products of creation and annihilation operators in the set {a 1,i , a † 1,i } i∈sites , and similarly for X({a 2,i , a † 2,i }). The operator X 2 has the property [X 2 , R 2 ] = 0, which will be utilized shortly. Suppose X 2 = i λ i P i where the {P i } are orthogonal projectors. Then after measurement one is left with 4. Measure R 2 = Π + − Π − (where Π ± is the projector onto the ± eigenspace) to obtain 5. The probability that one measures R 2 as +1, after having measured ρ ⊗2 to be in the subspace corresponding to P i , is denoted by Prob(+ | i). Similarly, the probability that one measures R 2 as −1, after having measured ρ ⊗2 to be in the subspace corresponding to P i , is denoted by Prob(− | i). After obtaining Prob(+ | i) and Prob(− | i), one can compute where we have used [X 2 , R 2 ] = 0 to go from the second line to the third line, and Eqn. (17) to go from the third line to the last line. A similar procedure can be used to determine tr{ρ 2 }, and then one can compute the quotient tr{X ρ 2 }/tr{ρ 2 }.
In an actual experiment, one does not directly measure the parity operator R 2 , but instead measures the number operator on every site. Since the common refinement of the eigenspaces of all of the number operators is a refinement of the eigenspaces of R 2 , one can measure R 2 via the number operators and obtain the same result as above.

Fermion interferometry
It is straightforward to adapt the boson interferometry techniques to fermions, although a few modifications to the protocol are required. Our protocol is inspired by the work of [42]. Suppose we have two systems of fermions, and require that states of different fermion number lie in different superselection sectors. Technically, the superselection rule means that for all observables X, we have ψ 1 |X|ψ 2 = 0 if |ψ 1 and |ψ 2 are states of definite, but distinct fermion number.
For fermions, it is not true that tr R 2 F 2 ρ ⊗2 F † 2 = tr{ρ 2 }. Instead, we have where V has eigenvalues ±1 which depend on the total number of fermions N tot , the floor of half of the total number of fermions N tot /2 , and the number of fermions N 2 in the second copy of the subsystem. (There are, in fact, many choices of V which satisfy Eqn. (23), and so we choose a convenient one for our purposes.) The measurement outcomes for V are given in the  The procedure for measuring tr{X 2 ρ 2 }/tr{ρ 2 } is the same as in the bosonic case above, except that now we need (where here the f, f † operators are fermionic) to additionally satisfy So first let us find which operators, in general, commute with V. Suppose we have an operator of the form n1 of these f 1,2 · · · f n 2 ,2 n2 of these (25) where {i 1 , ..., i m1 }, {j 1 , ..., j m2 }, {k 1 , ..., k n1 }, { 1 , ..., n2 } are all sets with non-repeating elements. All of these operators transform multiplicatively by either +1 or −1 after conjugation by V. Letting m = |m 1 − m 2 | and n = |n 1 − n 2 |, the possibilities are tabulated below: For example, letting X 2 = 1 2 (n i,1 + n i,2 ), we have [X 2 , V] = 0. If instead X 2 = F 2 ( 1 2 (n i,1 n j,1 + n i,2 n j,2 )F † 2 , we likewise have [X 2 , V] = 0.

II.
Extracting two point correlations in the low temperature limit In this section, we numerically study the effect of the second term in Eqn. (5) in the main text. More specifically, we have argued that one can extract a density-density correlation from a more experimentally accessible quantity: While the first term is the desired density-density correlation, the second term arises as a consequence of the Fourier transform of local operators n j and n . As described in the main text, however, we expect that at sufficiently low temperatures the second term does not exhibit any systematic dependence on the distance between two points d = |j − |, allowing us to extract physically meaningful quantities such as correlations lengths from fitting C(j, ) as a function of d.
In order to confirm this expectation, we consider a 1D Bose-Hubbard Hamiltonian with nearest-neighbor hopping rate J and on-site repulsive interaction U = 3J. We numerically compute thermal density matrices for N = 4 particles on L = 16 lattice sites with periodic boundary condition at various temperature T /J ∈ { 1 10 , 1 5 , 1 4 , 1 2 , 1}. For each temperature T , we compute each term in C(j, ) as well as their sum as a function of the distance d ∈ {1, . . . , 8}. Fig. 4 below summarizes our numerical results, from which it can be checked that the density-density correlation (the first term in C(j, )) displays strong anti-bunching (Fig. 4a) at low temperature. By contrast, the second term exhibits diminishing distance-dependence as the temperature decreases (Fig. 4b). We find that the distance dependence of the total value C(j, ) is indeed dominated by the density-density correlation (Fig. 4c) at sufficiently low temperatures. Additional contribution to C(j, ) arising from the second term. Crucially, this contribution exhibits decreasing distance-dependence in the low temperature limit. (c) The position dependence of the total value C(j, ) is dominated by the first term in low temperature limit.

III. Experimental methods
Our experiments start from a high fidelity Mott insulator with a single particle per lattice site. Using high-precision, site-resolved optical potentials, created by a digital micro-mirror device (DMD), we isolate two neighboring six-site long chains of atoms with exactly one atom on each site. In order to ensure the high fidelity of the initial state we hold it in the 45E r deep optical lattice in both directions. To obtain a locally thermal state we suddenly drop the lattice depth along the chains, allowing atoms to tunnel, while keeping the lattice high between the chains. We use a pair of DMD beams to offset the sites right outside the region of interest, thereby defining the overall length of the system. After variable evolution time, we freeze the dynamics along the chains by suddenly ramping up the lattice back to 45E r . In order to make sure that the state has thermalized, we pick evolution times for which the entanglement entropy of the region of interest has reached its saturation value. Table III shows the times used in Fig. 2   In order to implement the beamsplitter operation, we drop the lattice depth between the chains and let the atoms evolve for a certain time duration. During this process the lattice depth along the chains stays high, preventing wavefunction evolution in that direction. At the end of this sequence, we read out the state of the system in the particle number basis with single-site and full atom-number resolution. For more details see [28].