Quantum point spread function for imaging trapped few-body systems with a quantum gas microscope

Quantum gas microscopes, which image the atomic occupations in an optical lattice, have opened a new avenue to the exploration of many-body lattice systems. Imaging trapped systems after freezing the density distribution by ramping up a pinning lattice leads, however, to a distortion of the original density distribution, especially when its structures are on the scale of the pinning lattice spacing. We show that this dynamics can be described by a filter, which we call in analogy to classical optics a quantum point spread function. Using a machine learning approach, we demonstrate via several experimentally relevant setups that a suitable deconvolution allows for the reconstruction of the original density distribution. These findings are both of fundamental interest for the theory of imaging and of immediate importance for current quantum gas experiments.

Introduction. Imaging with high resolution is a cornerstone for understanding the structure, dynamics and functionality of matter [1][2][3]. In the field of ultracold atoms, quantum gas microscopes have opened new avenues for studying lattice systems [4][5][6][7] and led to remarkable progress and insights, such as density correlations and string order [8], long-range anti-ferromagnetic correlations [9] or entanglement growth [10] in Mott insulators. Naturally, it is of equal interest to study trapped, i.e. non-lattice systems, where imaging with single-atom sensitivity is also vital for exploring beyond mean-field physics, i.e. for probing correlation effects [11]. Singleatom resolved imaging in free space has been demonstrated for metastable helium atoms, which can be detected using a multi-channel plate with a typical resolution of 60 µm [12], and recently for lithium atoms using a short fluorescence pulse, where the position spread due to scattering recoils can be reduced to 4 µm [13]. In order to reach sub-micron resolution, the positions of the atoms have to be frozen by ramping up a pinning lattice before the fluorescence imaging and detection of the atoms takes place. Such a capture of atoms in a pinning lattice was demonstrated starting from a larger scale lattice [14] or a larger scale continuous system [15], but freezing and measuring of density structures on the scale of the pinning lattice spacing was so far not considered and achieved.
Alternative schemes to reach sub-lattice resolution of quantum gases, inspired by related imaging techniques in other fields, have been proposed. Stimulated emission depletion microscopy [16], which breaks the diffraction limit set by the imaging wavelength, can be adapted to quantum gases using the position-dependent dark state of a Lambda-system [17]. A scanning tunneling microscope could be realized by coupling to a single ion [18] or by using dispersive couplings to a cavity [19]. Momentum mapping in combination with phase retrieval should allow imaging with 1-2 orders of magnitude better than the lattice spacing [20]. Finally, scanning electron microscopy was successfully applied to quantum gases reaching a resolution below 150 nm [21]. However, the combination of sub-micron resolution and single-atom sensitivity has so far only been achieved by fluorescence imaging in a pinning lattice.
Here, we propose to perform repeated measurements with shifted positions of the pinning lattice relative to The measurement signal is obtained as the average over many realizations, but it contains the deformation due to the ramp dynamics. By repeatedly preparing a realization of the system and freezing with different phases of the pinning lattice ϕ, the density after the ramp can be sampled with a resolution below the lattice spacing a l . (f) By deconvolution with the quantum point spread function (qPSF), the original density can be recovered. All sub-figures show sketches. the trapped physical system, such that a resolution below the lattice spacing becomes possible and we provide an in-depth analysis of this protocol. We show that the density structures on the scale of the lattice spacing will be distorted due to the dynamics taking place during the ramp-up of the pinning lattice. The lattice ramp has to be sufficiently fast to avoid an adiabatic loading of the ground state of the lattice, but sufficiently slow to avoid projections onto very high bands, where the atom positions are not frozen due to large tunneling rates. The proposed scenario is illustrated in Fig. 1. We show that the distortions during the ramp-up can be captured by a quantum point spread function (qPSF). Using deconvolution techniques, these distortions can be removed, which enhances the resolution of the overall measurement sequence. We find that the deconvolution is both relevant and effective for density structures on the scale of the lattice spacing and provides a sub-wavelength resolution. Our approach and technique suggests itself for immediate applications, because a tight confinement and resulting small structures of the original trapped system allow for strongly interacting quantum systems, while the spacing of the pinning lattice is fixed to typically 0.5 µm by the optical wavelength of the interfering laser beams.
Quantum point spread function. We first derive the qPSF for the measurement of a single particle in the premeasurement state |φ [22] and then extend the concept to many-body systems. The measurement is modeled as a two-step process: the ramp-up of the pinning lattice and the read-out of the state occupations. In the following, we keep the phase off-set of the pinning lattice ϕ fixed and thereafter vary it for resolving fine density structures. During the ramp-up, we assume that all external potentials but the pinning lattice are either switched off or negligible such that the quantum dynamics is governed by the Hamiltonian Here, the lattice depth V (t) is ramped up from zero to its maximal value V f within the time-scale T f by using a tanh-like ramping protocol and k l denotes the pinning lattice wavenumber corresponding to the lattice spacing a l . The lattice sets the recoil energy E r = 2 k 2 l /(2m) as typical energy scale. Directly after ramping up the lattice, the system is in the stateÛ ϕ |φ withÛ ϕ = T exp(−i/ T f 0 dτĥ ϕ (τ )) andT denoting the chronological time-ordering operator.
The occupation of the site i is then read out via fluorescence imaging, which we describe within the established framework of measurement operatorsR i;ϕ and positive operator-valued measuresR † i;ϕR i;ϕ [23]. Being only interested in the probability for finding the particle at site i given the phase off-set ϕ we have to specify the operatorR † i;ϕR i;ϕ . For this purpose, we assume that a particle that ends up in the Wannier state |w α i;ϕ of the pinning lattice Hamiltonianĥ ϕ (T f ) after the ramp-up, where α denotes the band index, is measured with the detection efficiency η α ∈ [0, 1], which can be modeled bŷ Here, a high detection efficiency η α is ensured, if the tunneling rate of the band J α is small compared to the imaging time scales. As J α increases very rapidly for higher bands α, we can approximate the efficiencies by a step function, i.e. η α = 1 for a finite number of 'non-tunneling bands' and η α = 0 for all higher bands. Then the opera-torR † i;ϕR i;ϕ becomes a projector. Atoms in higher bands or continuum states [24] lead to loss and the lattice ramp has to be chosen such that this loss remains small. Deep lattices and not-too-fast ramps keep this loss negligible. Finally, the measurement signal s(x) is obtained by averaging over the pinning lattice shifts ϕ.
In order to define the qPSF we consider the analogy to classical optics, where the exact "object" density ρ(x) becomes blurred in the image plane via the point spread function f (x) according to the convolution , where s(x) denotes the signal in the imaging plane. Given this relationship and the precise form of f , there are various deblurring techniques for (approximately) restoring ρ(x). Our aim here is to reformulate the probability p i;ϕ as a convolution to define qPSF for our imaging protocol. By means of the translation symmetryÛ † ϕR † i;ϕR i;ϕÛϕ = T ia l +ϕÛ † 0R † 0;0R 0;0Û0T † ia l +ϕ with the translation operator T z = exp(−izp/ ), we arrive at our central result: where the kernel Q(x, y) = −x|Û † 0R † 0;0R 0;0Û0 | − y describes both the quantum dynamics during the ramp-up of the pinning lattice and the subsequent fluorescence imaging. Thus, we find that the probability of detecting the particle at site i given the pinning lattice offset ϕ is provided by the diagonal of the two-dimensional (2D) convolution of φ * (x)φ(y) with Q(x, y), which we therefore name quantum point spread function (see Fig. 2). Eq. (4) moreover shows that the probability p i;ϕ can be expressed by a continuous function s(z) evaluated at discrete positions, p i;ϕ = s(ia l + ϕ). By repeating the experiment for various offsets ϕ one effectively samples this pseudo-probability s(z), which shall be called signal in the following. For practical calculations, one can spectrally decompose the qPSF and finds that the signal s(z) equals an incoherent superposition of 1D convolutions of the pre-measurement state with the back-propagated x [a l ] x' whereπ denotes the parity operator. As a side remark, the quantum dynamics during the ramp is non-adiabatic such thatπ|χ α does not coincide with the corresponding Wannier state of the shallower lattices of the ramp. In order to extend the qPSF to many-body systems, which is accomplished in the supplementary material [25], we make the following two assumptions: first, we assume that all interactions are switched off, e.g., via a Feshbach resonance before ramping up the pinning lattice and reading off the site occupations. Second, the atomic density should be small enough such that lightassisted collisions can be neglected during the imaging [4,5]. Then, we can factorize the qPSF describing the statistics of single-shot measurements upon an N -body ensemble by N single-particle qPSF [25]. Moreover, the ensemble average over many such single-shots results in the signal s(z) = dxdy ρ 1 (x, y) Q(z − x, z − y) with ρ 1 (x, y) denoting the pre-measurement reduced one-body density matrix, if the population of unobserved pinninglattice bands and continuum states after the ramp-up is negligible [25]. We note that the spatial dimensionality does not play a role and the framework equally applies to higher spatial dimensions.
Deconvolution. Inverting the relationship Eq. (4) is quite a difficult task: deconvolution in general is an illposed problem and, moreover, we have to cope with the intriguing situation that the measurement signal s(z) constitutes only the diagonal of the 2D convolution (Q * ρ 1 )(z, z ′ ) [26]. By scanning over different lattice ramps we find for suitably chosen V f and T f that the real part of the qPSF Q(x, x ′ ) acquires a dominant diagonal pattern with a fast decay of the off-diagonal elements, while the imaginary part is significantly smaller [see Fig. 2(a),(b)]. This motivates us to express the signal s(z) as a 1D convolution of the one-body density ρ(x) = ρ 1 (x, x) with some yet unknown 1D filter q(x): However, it is a priori not clear, whether such a filter, which is independent of the underlying density, exists and if it does, how to obtain it. Yet, if one has found such a filter, Eq. (6) allows for applying established deconvolution algorithms for obtaining the pre-measurement density ρ(x). Making the most obvious choice by taking the diagonal of the qPSF, q(x) = Q(x, x), turns out to be numerically unstable and inaccurate. To compensate for the complexity of the 2D convolution, this diagonal needs to be readjusted. To this end we call upon a machine learning approach. Inspired by the multi-frame deconvolution technique, which is applied, e.g., in astronomy [27], we pursue the following machine learning approach to learn the unknown filter q(x). Our training set consists of a small number n t of one-and many-body states with known (reduced) one-body density matrix ρ (k) 1 (x, y), k = 1, ..., n t [25]. For each training sample, we calculate the corresponding measurement signal s (k) (z) by the full 2D convolution with the exact qPSF Q(x, y) [28]. Then, we determine the best 1D filter q(z) by minimizing the meansquared error on the training set nt k=1 dz [s (k) (z)− (q * ρ (k) )(z)] 2 via batch gradient-descent with line-search [25] [see Fig. 2(c)]. We finally apply a classical deconvolution algorithm [25] to Eq. (6) for several unseen cases s(z) to reconstruct the underlying pre-measurement density ρ(x).
Applications. We showcase the performance of our qPSF approach and deconvolution strategy using three physical example setups (see Fig. 3): excited harmonic oscillator eigenstates featuring a rapidly oscillating density, two identical bosons with infinite repulsion in a harmonic trap [29] and a Fermi polaron. For details on the implementation of these systems see [25]. The examples are chosen to cover a broad range of different situations: single particle, weakly-and highly-correlated few-body physics.
The deconvolution uses a multi-frame filter q [see  7)] of the measurement signal Ds (blue solid line) and the reconstructed density D d (red solid line) from the genuine singe-particle density as a function of the structure size (see [25] for definitions).
ton in a small BEC as additional training example (see [25]). We find that applying it to the unknown signals from the two-Boson and Fermi polaron problem, yields very good results, emphasizing the power of the method. We stress that we learn the multi-frame filter from singleparticle and mean-field cases and then apply it to the unseen situations, which involve correlated many-body states.
The insets in Fig. 3 show the genuine single-particle density ρ(x), the measurement signal s(x) and the deconvolved signalρ(x) for different physical examples. In all cases, the structure of the genuine density is washed out in the measurement signals, but almost completely recovered by the deconvolution. In particular, we recover all of the many oscillations for the harmonic oscillator with their full original contrast [ Fig. 3(a)]. In the two-Boson example, it is striking that the deconvolution successfully reproduces the original density although the two humps have merged into a single one in the measurement signal [ Fig. 3(b)]. In the Fermi polaron example, both the sharp dip in the center of the trap and the Friedel oscillations around it are fully recovered in the deconvolved signal, although they seemed to be lost in the measurement signal [ Fig. 3(c)]. These examples showcase the power of the deconvolution method using the qPSF.
To judge on the quality, we introduce a dissimilarity measure between two normalized functions g and h as It takes a value of zero for coinciding functions and increases up to one as the absolute deviation becomes more pronounced. Further, we define the dissimilarity between the measurement signal s [30] and the genuine density ρ as D s = D(s, ρ) and the dissimilarity between the reconstructed densityρ and the genuine one as D d = D(ρ, ρ). In Fig. 3, we show how this dissimilarity depends on the typical structure size σ of the genuine density. For structures that are large compared to the lattice spacing (σ > 2a l ), the dissimilarity D is negligible both for the measurement and the deconvolved signal. When the structures are on the scale of the lattice spacing, the measurement signal starts to deviate due to the dynamics during the ramp-up of the pinning lattice. The dissimilarity of the deconvolved signal, however, remains negligible due to the successful deconvolution. Only for structures smaller than about half the lattice spacing (σ < 0.5a l ), we observe an increase of D d , indicating the limitations of the method. Using the deconvolution via the qPSF, we can therefore shift the accessible structure sizes from about 2a l to about 0.5a l , which is a significant improvement that is crucial for many physical examples in quantum gas physics.
Outlook. Our work opens the research direction of high-resolution imaging with single-atom sensitivity also for trapped, i.e. non-lattice systems. We propose to apply a pinning lattice for imaging and to sample the reduced one-body density with a resolution below the lattice spacing by performing repeated measurements with shifted positions of the pinning lattice relative to the physically trapped system. We have shown that density distortions resulting from the dynamics during the ramping up of the lattice can be compensated by deconvolution with a quantum point spread function for a wide range of parameters. Our findings are of immediate relevance for ongoing quantum gas microscope experiments. A reliable measurement of small density structures will allow accessing new regimes and imaging of the corre-sponding physical processes such as the shape of a vortex core taking into account beyond mean-field effects [31] or discrete few-body structures in arbitrary traps. For simplicity, we have focused here on one-dimensional systems, but our framework equally applies to higher spatial dimensions. Further extensions of our work would be the fate of correlation measurements [32] and blurring effects in the measurement of the dynamics. Another important aspect is the imaging after release from a driven system, e.g., for producing artificial gauge fields [33], where switching off the drive can induce further effects. Releasing from lowest Landau levels yields a self-similar expansion of the wave function [34], which could be used before freezing the distribution. P. S. and C. In the following, we extend the qPSF theory to many-body systems by first inspecting the case of distinguishable particles, then discussing single-shot measurements of indistinguishable particles, and finally deriving the relationship between the ensemble average of such single-shot measurements and few-body correlations.

Single-shot measurements of distinguishable particles
In order to extend the qPSF theory to many-body systems, we make the following assumptions: (i) As in the single-particle case, we assume that all external traps but the pinning lattice are either switched off or negligible during the full measurement protocol. (ii) Moreover, we assume that the inter-atomic interactions are either switched off by means of a Feshbach resonance or negligible during the full measurement protocol. (iii) Finally, we regard the fluorescence imaging of the pinning lattice sites to be a pure one-body process, i.e. neglect few-body effects such as loss via light-induced collisions [4,5]. The latter approximation is valid, if the pre-measurement atomic density is so low that the likelihood of finding more than one atom in a pinning-lattice site after ramp-up is strongly suppressed.
Under the above assumptions, an N -body system in the pure [35] pre-measurement state |Ψ evolves into the stateÛ |Ψ during the ramp-up of the pinning lattice with given phase off-set ϕ. HereÛ ϕ (τ )] denotes the single-particle time-evolution operator acting on the κ-th particle. Since the fluorescence imaging is modeled as a single-particle process, we can directly transfer the positive operatorvalued measure for the single-particle case [see Eq. (3) of the main text] to the many-body realm and obtain for the probability to detect the 1st, 2nd, ..., N -th particle in the pinning-lattice site i 1 , i 2 , ..., i N , respectively: shall indicate that the whole operator acts on the κ-th particle). Note that for a fixed phase the probability i P ϕ (i) ≤ 1 due to the possibility of detection efficiencies being smaller than one. Making use of the translation symmetry ofM (κ) iκ;ϕ as in the single-particle case discussed in detail in the main text, we may express Eq. (S1) as where the spatial positions of the N particles are abbreviated as x ≡ (x 1 , ..., x N ) and the integrals are taken w.r.t. all particle coordinates, i.e., d N x ≡ dx 1 dx 2 ...dx N . Moreover, the pinning-lattice sites, in which the particles are detected, are abbreviated as i ≡ (i 1 , ..., i N ) and Ψ(x) ≡ x 1 , ..., x N |Ψ refers to the position representation of the N -body pre-measurement state |Ψ . Finally, the N -body qPSF turns out to be the N -fold product of the one-body qPSF, which has been derived for the single-particle case in the main text: where Q(x κ , y κ ) = −x κ |M (κ) 0;0 | − y κ . Thereby, we obtain the N -particle post-measurement distribution P ϕ (i 1 , ..., i N ) for a given pinning lattice phase off-set ϕ by evaluating the signal function at the discrete positions z = ia l + ϕ, i.e., P ϕ (i) = S(ia l + ϕ). Repeating the N -body measurement for various phase off-sets ϕ effectively means sampling from the pseudo-probability S(z).

Single-shot measurements of indistinguishable particles
Next, we concentrate on the special case of N indistinguishable particles. In this case, the outcome of a single-shot measurement is a pinning lattice occupation-number histogram (n 1 , ..., n L ) ≡ n, where n i denotes the number of particles found in the i-th lattice site, i = 1, ..., L and L refers to the number of lattice sites. Here, we have in particular few-body situations in mind, where one can easily probe the full distribution of the N particles in the pinning lattice.
Obviously, the N -body qPSF Q N (x, y) remains invariant under simultaneous permutation of the particle labels in both x and y. Given a system of indistinguishable particles, the probability Eq. (S2) does not depend on the concrete particle labeling but only on how many particles are found in a certain site. Taking this combinatorically into account, one finds for the probability of the histogram n for a given phase off-set φ where i n denotes some N -dimensional lattice-site index vector, which features n r -times the entry r with r = 1, ..., L.
We remark that while Eq. (S5) describes the (within the considered measurement model) correct probability of detecting the histogram n given the number of particles N and the phase off-set ϕ, these probabilities do not sum up to unity in general when considering all conceivable histograms n with N particles. In fact, the probability for not detecting all N particles in the pinning lattice due to the occupation of higher bands with detection efficiencies smaller than unity or continuum states after ramp-up [see Eq. Ensemble averages over single-shot measurements and few-particle correlations Having taken many single-shot measurements of identical copies of the many-body system, one may evaluate the corresponding ensemble average of certain n-particle observables. In classical absorption imaging of atomic samples for instance, one obtains the reduced one-body density by averaging the spatial particle number distributions of many single-shot measurements. Density-density correlations can be inferred from absorption images by averaging the product of occupation-number fluctuations at two spatial positions over many single-shot measurements. Here, we stress that in classical absorption imaging the average of an n-particle quantity over many single-shot measurements is directly connected to the pre-measurement reduced n-body density matrix, whereas in the case of the quantum gas microscope measurement protocol pursued in this work, this relationship is more complicated in general and shall be derived here.
First, let us derive the probability p (1) ϕ (r) to find an atom in the pinning-lattice site r when averaging over many single-shot measurements with the same phase off-set ϕ, i.e., many different histograms n distributed according tō P ϕ (n). Using Eq. (S5), we find Apparently, only histograms n with n r > 0 contribute to p (1) ϕ (r). Substituting n = m + e r , where m denotes an arbitrary (N − 1)-particle histogram and e r an occupation number vector with all components zero except for the r-th one being set to unity, one obtains Next, we rewrite the summation over (N − 1)-particle histograms as a summation over N − 1 lattice site indices i;ϕ , we finally obtain Since particles in higher bands or continuum states after the pinning lattice ramp-up are not detected,K (κ) ϕ = 1 1 and thus p (1) ϕ (r) is not given as the expectation value of a one-body observable. So the one-particle quantity p (1) ϕ (r) may depend on up to N -particle corrections and cannot be represented as the trace of a one-particle observable times the pre-measurement reduced one-body density operator in general, which is in contrast to the case of absorption imaging.
Our simulations in the main text, however, show that the probability to populate higher bands or continuum states by the pinning lattice ramp-up is negligibly small for suitably chosen experimental settings (see the discussion on the impact of higher bands on the qPSF in Section B as well as Table S1). Under these circumstances,K (κ) ϕ effectively acts as the identity operator on the κ-th particle in |Ψ and we obtain the relation where ρ 1 (x, y) = x|ρ 1 |y denotes the position representation of the pre-measurement reduced one-body density operatorρ 1 , which one obtains from the pre-measurement many-body state by a partial trace over all but one particle,ρ 1 = tr 1 |Ψ Ψ| . Similarly, one can derive the corresponding expressions for the ensemble average of an n-particle quantity with n > 1 over many single-shot measurements. Here, we only explicate this relationship for the case n = 2, i.e., the probability p (2) ϕ (r 1 , r 2 ) to detect a particle at site r 1 and another particle in site r 2 in the ensemble average: (S11) Similar manipulations as above can be applied and one arrives at If the population of higher bands and continuum states after the pinning-lattice ramp-up may be neglected, we end up with where r = (r 1 , r 2 ) and ρ 2 (x, y) denotes the position representation of the pre-measurement reduced two-body density operatorρ 2 , which one obtains by a partial trace over all but two particles,ρ 2 = tr 2 |Ψ Ψ| .
Section B: Numerical procedure to obtain the quantum point spread function According to Eq. (4) of the main text the qPSF is an operatorQ = α η α |χ α χ α | with |χ α =πÛ † 0 |w α 0;0 . Therefore, we have to calculate the Wannier state |w α 0;0 of the band α at site i = 0 for the pinning lattice with the final potential depth V f and phase off-set ϕ = 0 and propagate it with the time-evolution operatorÛ † 0 =T exp(−i/ 0 T f dτĥ 0 (τ )), which describes the lowering of the pinning lattice from V f to zero depth. Finally, the parity operatorπ is applied to reformulate the measurement signal in terms of a convolution.
The Wannier states are obtained by representing the position operatorx in the basis ofĥ 0 (T f ) and then diagonalizing it. Afterwards, we set a band limit α max . Modeling the detection efficiencies η α for energetically high lying bands, however, is more involved as these depend on both the tunneling and fluorescence imaging time scale. The tunneling rates grow exponentially with the band index, such that they can be divided into tunneling and non-tunneling bands within the fluorescence imaging time to a good approximation. For the sake of simplicity, we therefore assume that all (bound) bands lying energetically below V f are detected with unit detection efficiency, meaning η α = 1 ∀ α ≤ α max , and η α = 0 for the continuum states, since atoms in these states are not pinned during the fluorescence imaging. As a consequence, this model does only give a lower bound on the loss in the measurement signal [Eq. (S25)] due to unobserved channels. We use a lattice containing L = 99 sites with 33 grid points to resolve each site, unless stated otherwise, while the potential depth is varied in the range V f ∈ [50, ..., 300] E r .
The back-time propagation of relevant Wannier states is performed with the Multi-Layer Multi-Configuration Time-Dependent Hartree for bosons (ML-MCTDHB) approach [36,37] to obtain |χ α . The ramping times cover T f ∈ [1, ..., 9] /E r and the ramping protocol V (t) is a logistic function of sigmoid form: parameter ensures that V (T f ) = V f does not deviate much from the saturated value V max . With η = 10 −3 fixed, T f alone determines the adiabaticity of the ramping protocol (see Fig. S1). In the case of an adiabatic preparation of many-body ground states in optical lattices, such as the bosonic Mott insulator, the optimal shape of the ramp function has been extensively discussed [38]. In contrast, for pinning the distribution on the lattice in quantum gas microscopes, simple s-shaped ramps have proven sufficient [4,5]. We note that in our setting, the dynamics during the ramp will be strongly non-adiabatic in order to avoid a loading of the ground state of the lattice, but freeze the atoms in their original position. Therefore, we expect that the precise shape of the ramp should not be important.
We show the real and imaginary part of the spatial representation of the qPSF Q(x, x ′ ) = x|Q|x ′ for a quick ramp T f = /E r with a deep lattice V f = 200E r and for a slow ramp T f = 9 /E r with a comparatively shallow lattice V f = 50E r (Fig. S2). In the first case we observe a diagonal pattern in the real part with a fast decay of the offdiagonal, while in the second case the real part displays a Gaussian profile with the imaginary part being suppressed by an order of magnitude. Both cases are rather localized around a small region of approximately 5a l . The diagonal pattern can be induced and enhanced by choosing deeper lattices, meaning that higher bands are responsible for this effect, although by successively adding bands for the qPSF calculation we found that approximately only the first half of the bands α max is responsible for the pattern formation. Going to ramp times beyond T f = 9 /E r requires large lattices with more than L = 100 lattice sites, because the Wannier states, propagated back in time, almost reach the boundaries of the grid. x [a l ] x' The qPSF has no direct relation to the classical PSF of the imaging system with finite numerical aperture NA, which is used for the fluorescence imaging after the pinning of the atoms. As long as the NA is large enough to allow for a reconstruction of the lattice occupation (typically NA=0.6-0.8), it drops out of the problem. If one repeats the measurement with varying positions of the pinning lattice with respect to the initial system via the displacement by ϕ, even the lattice constant a l does not pose a fundamental limit to the resolution. In the numerical examples a sampling with resolution 0.03a l was used and similar relative positioning of 0.1a l between the pinning lattice and further traps were reported experimentally [39]. The distortion from the dynamics during ramp-up, which is captured by the qPSF and is relevant for structures on the order of a l , is therefore the fundamental limitation on the resolution. The deconvolution with the qPSF can then lead to density measurements with a resolution significantly higher than a l .

Section C: Examples of application
For the numerical implementation we make use of recoil units x r = 1/k l , E r = 2 k 2 l /(2m), T r = /E r with the wavenumber k l = 2π/λ l of the laser beam of wavelength λ l = 1064 nm to create the lattice potential and m being the mass of the trapped particles, here 87 Rb. The lattice constant is a l = λ l /2. i) Harmonic oscillator (HO) eigenstates where a ho = mω the harmonic oscillator length, ω the frequency of the trap and n ∈ N 0 the excitation level. To characterize the structure size of the HO modes with respect to the lattice we consider the variance of the position operator divided by the number of peaks in the density profile σ n /a l with σ n = 1 n+1 ( ψ n |x 2 |ψ n − ψ n |x|ψ n 2 ) 1/2 . ii) Dark soliton We prepare a dark soliton within the mean-field approximation placed in a reflection-symmetric box with an extension L b smaller than that of the pinning lattice L l = La l . We position the soliton in the center of the box and ensure that it is sufficiently separated from the walls: where the prefactors c i are chosen such as to ensure the continuity and the normalization of the wave function, ξ = 1/ √ 8πρa sc is the healing length of the condensate, a sc the scattering length andρ the constant background density. The structure size is chosen as σ ξ /a l = 2ξ/a l , which is approximately the full-width-at-half-maximum (FWHM) of the soliton profile. The soliton example is used for training of the filter only.
iii) Impurity in a Fermi sea We put N spin-polarized fermions in a reflection-symmetric box of length L b < L l . A stationary impurity positioned in the middle of the potential acts as a repulsive delta-potential of infinite strength, inducing a density profile of fermions similar to that of a soliton, but with an oscillatory background. The eigenstates have a defined parity: The density operator for an even number of fermions is then given by a mixed statê Here, the structure size is assigned by an average extension of a peak in the one-body density σ N /a l = (L b /N )/a l . iv) Two bosons with infinite repulsion in HO The highly correlated problem of two bosons trapped in a harmonic trap and interacting with each other via a delta-potential of infinite strength can be solved analytically in the relative frame [29]. By transforming the solution back into the laboratory frame and tracing out one of the coordinates one obtains the following one-body density matrix: g(x, y) = √ π xy + 1 2b 2 (erf(bx) − erf(by) + 1) + and b = 1/a ho . The correlated two-body system requires the full 2D convolution to create the signal, which is cumbersome to achieve on a large grid with fine resolution. So we consider very large trapping frequencies and reduce the grid to L = 33 lattice sites. The structure size is defined similar to the HO case: σ/a l with Section D: Simulation of the measurement signal In the most general formulation the distorted signal s(z) can be obtained directly via a 2D convolution of the one-body density matrix ρ 1 (x, x ′ ) of the initially prepared system (Section C) with the kernel Q(x, x ′ ) (Section B): However, the (L · 33) × (L · 33) matrices lead to approximately (L · 33) 4 numerical operations, which renders the direct calculation inefficient for larger grids. One way to circumvent this issue would be to make a smaller support for the density by confining it more tightly and for the filter by defining a cutoff, when the amplitudes drop below a certain value. Here, we just verified that a spacing ∆x = (1/33)a l provides converged signals by doubling the site resolution.
Independently of the above statements we can reduce the numerical effort to ∝ (L · 33) 2 , namely for weakly correlated systems the spectrally-decomposed one-body density operator has a finite number of natural populations λ γ with considerable weight:ρ with |φ γ natural orbitals. Inserting this relation and additionally the expansion of the qPSF into Eq. (S22) we obtain the signal as a sum of 1D convolutions of the natural orbitals φ γ with the 'band' filters χ * α : There is another important point worth mentioning, namely the padding. Since we are working with finite systems, a convolved function spans a larger region than the input functions. Thus, we need to provide values for chosen densities outside the grid and padding with zeroes is the most natural choice for trapped systems, while periodic padding would be suitable for ring geometries. Also, we ensure that the distortion of the signal does not reach the boundaries of the grid.
For a given lattice realization, meaning fixed V f , T f and phase ϕ, the signal s(z = i · a l + ϕ) with spatial sampling period a l sums up to unity only when η α = 1 ∀ α, becauseR † i;ϕR i;ϕ then forms a positive operator-valued measure. In our case, neglecting continuum states lying energetically above V f results in a particle loss Ω. In other words, Ω = 1 − L i=1 s(i · a l + ϕ) is the probability of finding a particle in none of the sites, but in the unobserved channels. Averaging over multiple lattice realizations ϕ ∈ {0, ..., π} we can estimate the mean particle lossΩ expected for the given pre-measurement reduced one-body density ρ 1 : In Table S1 we show the average lossΩ for densities and ramp-up parameters discussed in the main text, which is indeed very small and has a tendency to decrease for larger structures.
structure size σ/a l 0.5 1.0 2.0 HO n = 10 0.030 0.014 0.013 ρ 1 two-Bosons 0.021 0.015 0.014 Fermi polaron 0.033 0.031 0.017 TABLE S1. Average particle lossΩ(V f , T f , ρ1) for the densities, discussed in the main text, and different structure sizes σ (see Section C) relative to the lattice spacing. The ramp-up parameters are V f = 200Er and T f = /Er.

Section E: Multi-frame filter
Our starting point is Eq. (6) of the main text. First, we generate a batch of signals s (k) (z) with k = 1, ..., n t by the full 2D convolution of chosen one-body density matrices ρ (k) 1 (z, z ′ ) with the qPSF Q(z, z ′ ). Then we estimate each signal s (k) (z) as a 1D convolution of the corresponding densities ρ (k) (z) = ρ (k) 1 (z, z) with the same density-independent filter q(z): Further, we assume q to be space invariant, but otherwise no priors will be imposed, because it has no physical interpretation and is rather a mathematical tool. Thus, while the densities of the training set vary from signal to signal, the same filter q is common to all signals. Each of them provides additional information on q, thereby restricting the space of possible solutions.
As to the choice of training samples we create a random selection of five soliton samples with different extension ξ ∈ [10, ..., 40] × 10 −8 m in a box L b ∈ [3/4, ..., 10/12]L l as well as a random selection of five HO samples with trapping frequency ω ∈ [200, ..., 800]×2π Hz and excitation level n ∈ {1, ..., 10}. We explicitly don't include correlated examples from Section C to later test the filter q on unknown signals.
In the next step we define a total loss function L, which describes a deviation between the true signals s (k) and their approximations ρ (k) * q, a least squares problem: where in the last step we switch to a numerical grid with s (k) and q being (L · 33)-dimensional vectors and A (k) denoting a (L · 33 × L · 33) Toeplitz matrix, which represents a 1d discrete convolution with zero padding and limited support.
To find the filter q, that is more likely to have created the observed distortions in the signals, we perform the gradient descent algorithm in batch mode, meaning that we take into account all the frames simultaneously. We find that a small amount of samples is sufficient to obtain a well-performing filter. Thus, we do not need to resort to more memory-efficient optimization algorithms such as stochastic or mini-batch gradient descent.
In each iteration step m ∈ N 0 the filter is updated such that we follow a path towards the minimum of L by taking a direction of negative gradient ∇L. As initial guess we take the diagonal of the qPSF q 0 = Q(x, x) and then iterate where the gradient of the loss function reads: The step size or learning rate β can be optimally calculated (accurate line search) for each iteration step as β = argmin β L(q m − β∇L(q m )) = 1 2 .

(S30)
Finally, we iterate until the relative change in the total loss function reaches some threshold.

Section F: Deconvolution
During the image acquisition by microscopes in molecular biology [40] or telescopes in astronomy [41] multiple degradation sources can distort the true form of the object: noise, scatter, glare and blur. The blur, caused by the passage of light through the imaging system, leads to a non-random light redistribution and poses a fundamental limitation to the imaging system. The recorded image is usually modeled as a convolution of the object with a filter, also known as point spread function (PSF). There exists a variety of methods to reverse this process and retain the original object, called deblurring or deconvolution algorithms [42]. They can be classified as inverse ) or iterative (Van-Cittert [44], Lucy-Richardson [45, 46], Steepest Descent [47]); with prior knowledge of the filter (non-blind deconvolution) or completely unknown (blind deconvolution [48]); imposing priors such as non-negativity and smoothness or without them; modeling potential noise sources or neglecting them; using a single frame or a batch of sampled frames (multi-frame deconvolution [49]).
For our case we require a package, which implements an iterative algorithm, as they are more stable and provide a better restoration of degraded resolution, although at the cost of longer computation times. Since we obtained the filter in Section E it should be non-blind. The density we are trying to reconstruct is positive, so a non-negativity constraint is a must, but otherwise no pre-filtering is necessary. We also neglect all sources of noise and the measurement signal is considered as a single frame. wolfram mathematica 10.4 provides two suitable algorithms for this purpose: Lucy-Richardson and Steepest Descent. The output is a positive function. We iterate until converged, disable the pre-filtering and don't include any noise. Steepest Descent proves to be more reliable than Lucy-Richardson and the reconstruction is usually of a better quality.