Quantum read-out for cold atomic quantum simulators

Quantum simulators allow to explore static and dynamical properties of otherwise intractable quantum many-body systems. In many instances, however, the read-out limits such quantum simulations. In this work, we introduce an innovative experimental read-out exploiting coherent non-interacting dynamics. Specifically, we present a tomographic recovery method allowing to indirectly measure the second moments of the relative density fluctuations between two one-dimensional superfluids, which until now eluded direct measurements. Applying methods from signal processing, we show that we can reconstruct the relative density fluctuations from non-equilibrium data of the relative phase fluctuations. We employ the method to investigate equilibrium states, the dynamics of phonon occupation numbers and even to predict recurrences. The method opens a new window for quantum simulations with one-dimensional superfluids, enabling a deeper analysis of their equilibration and thermalization dynamics. Quantum simulators give access to quantum dynamics but for large systems it is difficult to read-out informative observables. By using precise modelling the authors demonstrate how to read-out expectation values of non-commuting sets of observables of phonons in cold atomic gases which gives information about their dynamics and thermodynamics.

Q uantum simulators offer entirely new perspectives of assessing the intriguing physics of quantum many-body systems in and out of equilibrium. They are experimental setups allowing to probe properties of complex quantum systems under unprecedented levels of control [1][2][3] , beyond the possibilities of classical simulations. Among other platforms, experiments with ultra-cold atoms involving large particle numbers or even continuous quantum fields have been particularly insightful [4][5][6][7][8][9][10][11][12][13][14] .
And yet, key questions remain open for a highly unexpected reason: the read-out of state-of-the-art quantum simulators is limited. In one-dimensional (1D) superfluids, for example, one can probe the dynamics of equilibration 4 occurring in the presence of an effective light-cone 5 and leading to generalized Gibbs ensembles 6 . The excellent experimental control over that system allowed to observe coherent recurrences in the dynamics of a system of thousands of atoms 8 . However, in that particular setup, further quantifying the recent observations is currently obstructed because only phase quadratures but not canonically conjugate density fluctuations can be measured. On the contrary, if both quadratures could be measured, and hence if genuine quantum read-out was possible, then studies of intricate questions on the role of interactions, or entanglement dynamics after a quench, could become possible.
This situation is by no means an exception: In fact, in any quantum simulation platform, read-out prescriptions are always restricted in one way or another, which constitutes a crucial bottleneck towards studying intricate physical questions. For cold atoms in optical lattices, akin to the development which will be laid out in this work, innovations such as the quantum gas microscope [15][16][17] directly opened up the path towards studying exciting physical phenomena [9][10][11][12][13][14] . Sophisticated read-out methods are therefore highly desirable and key to a platform.
In this work, we show that quenching a single global parameter in the system can enable a genuine quantum read-out and even allow for state reconstructions. We hence open a new 'window' into a quantum simulator for which so far-as is common in quantum simulation-only incoherent, or 'classical', read-out was natively possible.
Tomography of many-body systems is typically limited due to the number of necessary observables and complexity of control. However, often for large systems, the observed dynamics can be described by an effective free field theory capturing the dynamics by means of modes which are long-lived, i.e., not over-damped. We shall demonstrate that observing their dynamics can already suffice to perform reconstructions of the relevant correlation functions for many-body systems. The basic principle is that nonequilibrium evolution can mix 'quadratures' (non-commuting operators describing the dynamics of a mode) in such a way that consistency of observed correlations implies constraints regarding the unobserved ones. We will show that they can be quantitatively reconstructed.
Specifically, in this work, we set out to introduce a novel method of tomographic read-out for quantum simulators, by combining known quantum dynamics and available measurements to obtain more information. Related ideas of exploiting known random or deterministic unitary dynamics to get access to otherwise inaccessible types of measurements have been theoretically explored [18][19][20][21][22][23] . However, closest in spirit are tomographic methods in quantum optics where a harmonic rotation in phase space allows to measure two canonically conjugated quadratures by a detector sensitive to only one of them, and hence perform a 'quantum measurement'.
Here, we consider this basic idea in a genuine multi-mode setting and demonstrate a practical application of the method to 1D superfluids. We acquire data at different times for many modes at the same time and make use of semi-definite programming techniques for achieving reconstructions insusceptible to noise. After introducing our new recovery method, we explore the physics of 1D superfluids studying the properties of quench dynamics and its initial conditions. We use the quantum read-out information concerning both density and phase fluctuations in momentum space to fit the temperature and the global tunnel coupling parameter of the initial state preparation. Concerning out-of-equilibrium dynamics, we are able to predict recurrences by relying solely on data taken at times away from the recurrence occurrence, demonstrating that the system is coherent throughout the evolution. Finally, we monitor dynamics of phonon occupation numbers constraining their growth over an extensive observation time giving further experimental evidence for the validity of the effective model. The successful functioning of our method demonstrates an excellent agreement of the experiment with the theory of elementary excitations of 1D superfluid 24,25 . Our approach is based on very general and ubiquitous ingredients, hence the framework that we establish can be expected to be in a natural way applicable to various quantum simulators.

Results
The system considered. In order to apply our read-out method in practice, we will consider the setting of two adjacent 1D Bose gases realized with ultra-cold atoms. Their low-energy relative fluctuations in phase and density,φ and δρ, are described by the effective Hamiltonian 8,25 with m being the atomic mass and n GP the average density profile defined by the ground state of a 1D Gross-Pitaevskii (GP) equation. The Hamiltonian describes phonons which are the elementary density-phase excitations, satisfying bosonic commutation relations ½δρðzÞ;φðz 0 Þ ¼ iδðz À z 0 Þ. The corresponding operators are defined within the atomic cloud whose spatial extension R 1D is given by the support of n GP . Experimentally, one can engineer the density profile n GP through the trapping potential. This determines the density-density interaction strength, which is functionally dependent on the density profile, gðzÞ ¼ g½n GP ðzÞ 8,26 (see Supplementary Note 1). For typical experimental parameters, the last term in Eq. (1) has less importance than the first two 27 , which together make up the Luttinger model. In the Supplementary Note 1, we describe a general numerical scheme for obtaining approximate eigenfunctions ofĤ for any n GP of interest. For a constant density profile, this model gives a linear spectrum and oscillatory eigenfunctions -which is qualitatively also the case for Eq. (1) even with small GP profile inhomogeneities. As the excitations are confined within the finite atomic cloud, their spectrum fω k ; k ¼ 1; 2; g is discrete 28 . We denote eigenmode operators of phase and density fluctuations byφ k and δρ k , respectively, and use their corresponding wave functions f ρ k ; f ϕ k 2 C 2 ð½ÀR 1D ; R 1D Þ to decompose the real-space fields aŝ Note that here we have ½δρ k ;φ k 0 ¼ iδ k;k 0 where δ k;k 0 is the Kronecker symbol in contrast to the Dirac delta in the commutation relations of the real-space fields above. Written in terms of the eigenmode degrees of freedomφ k and δρ k , the Hamiltonian becomes diagonal Hereρ 0 does not contribute to the visible dynamics because it is conjugate to the global phase which carries no energy and is removed from the data. We shall refer to the eigenmode operators as quadratures as they constitute a discrete set of bosonic observables, which harmonically rotate and do not mix between different modes as can be directly seen from their time evolution within the Heisenberg-picturê In the method presented below, we will make use of approximate eigenmodes obtained from the eigenfunctions of a spatial discretization of the considered model which, for simplicity, we will continue to denote byφ k ; δρ k and f ρ k ; f ϕ k -see the Supplementary Note 1 for a more detailed discussion. Let us stress that by the time evolution in Eq. (4) density fluctuations, which are not accessible by direct measurements, are dynamically mixed into the observed phase sector which is the foundation to our reconstruction approach.
Quadrature tomography. In this section, we turn to describing the reconstruction procedure, exploiting the known and efficiently tractable Hamiltonian dynamics on the one hand and ideas of reconstruction and signal processing on the other. This gives rise to a practical and versatile method of reconstructing correlation functions of a type inaccessible to direct measurement. The atom chip experiment 29 which we are considering here, measures referenced correlation functions of the relative phase through matter-wave interferometry 7,8,22,[30][31][32] Φðz; z 0 ; tÞ ¼ hðφðz; tÞ Àφðz 0 ; tÞÞðφðz 0 ; tÞ Àφðz 0 ; tÞÞi : ð5Þ Here we chose to reference the phase with respect to the middle of the system z 0 ¼ 0 μm, which removes from the data the global phase canonically conjugate toρ 0 . We aim to reconstruct the second moments of the initial state of the quadraturesr ¼ ðφ 1 ; ;φ N ; δρ 1 ; ; δρ N Þ T of the N lowest lying eigenmodes ofĤ satisfying bosonic commutation relations ½r k ; For this, we define the covariance matrix of the initial state as the collection of second moments It is important to note that a matrix V constitutes a collection of physically admissible second moments if and only if the matrix , has nonnegative eigenvalues. This condition reflects the Heisenberg uncertainty principle of canonically conjugated observables 33 . It will be convenient to use the notation Note that here we include the interest in both type of correlations, not only those related to phase fluctuations. Altogether, V allows to provide a Gaussian description of the full unknown state of the system whose validity can be verified by measuring vanishing higher-order connected correlation functions 7 . In this work, we will consider initial states that are approximately Gaussian and under this assumption determining V yields a full state reconstruction. For non-Gaussian initial states, the method presented in the following can still reconstruct the second moments V, but will not provide specific information on the higher moments.
Using the decomposition into eigenmodes in Eq. (2), we obtain for the observable second moments defined in Eq. (5) Note that we have introduced the cut-off N in the summation over the eigenmodes, anticipating that higher energy modes will have a negligible contribution in the measured signal either because they carry too much energy to or due to finite real-space resolution in the experiment.
Next, we exploit that the time evolution of the Hamiltonian in Eq. (3) does not mix quadrature operators of different eigenmodes as stated in Eq. (4) which gives where As summarized in Fig. 1, we now have all ingredients needed for quantum read-out. Specifically, based on the relations (8) and (9), we can recover the density correlations through a leastsquares recovery problem. For this, we collect all measured values of Φðz; z 0 ; t i Þ at different points z and z 0 and times t i in a vector b. Furthermore, we define a linear map AðṼÞ which, given some trial covariance matrixṼ, outputs the values of Φðz; z 0 ; t i Þ sorted as in b via Eq. (8). The time-evolution is implemented using Eq. (9) such that only the covariance matrix of the initial state is used to fit the observed data. If W denotes a weighting matrix, then WðAðṼÞ À bÞ is the vector of the weighted least-squares residues. The covariance matrix optimally fitting the data is then given by the solution to the following optimization problem The first line implements the minimization of the length kÁk 2 of the vector of weighted least-squares residues and the condition in the second line ensures that V is a physical covariance matrix. The  (1) is decoupled on the level of momentum space operators fφ k ; δρ k g involving multiple modes rotating at different frequencies (3). Measurements of real-space continuum fields fφðz; tÞg yield the referenced two-point correlation functions Φðz; z 0 ; tÞ defined in Eq. (5). In this work, by means of sophisticated post-processing using tools from signal processing, we are able to recover from real-space data taken at equidistant measurement times the full covariance matrix V (6) of the nonlocal eigenmodes of the Hamiltonian. This approach is general and can be applied to any other system in which quenches to non-interacting multimode Hamiltonians are available. optimal solution to this convex quadratic problem with a semidefinite constraint yields the covariance matrix V of the initial state with a minimal value of Θ. The optimization can be performed efficiently and reliably numerically with standard methods for semidefinite programming. We use the package cvx 34 , see the Supplementary Note 2 for more details on the implementation. The most basic idea of an algorithm solving (10) is to repeatedly take a steepest-descent step towards minimizing the least-squares residue and impose QðVÞk0. Standard convex optimization packages like cvx solve such a problem in a more sophisticated way ensuring numerical accuracy and converge in a matter of seconds. In the implementation, we chose a diagonal weighting matrix W with entries σ½Φðz; z 0 ; t i Þ À1 , where σ½Φðz; z 0 ; t i Þ denotes the standard deviation of each measured value. This weighting allows to put more emphasis on more precise values and yields with this a more reliable scheme as we find that Φðz; z 0 ; tÞ grows typically for increasing spatial separations jz À z 0 j but Φðz; z 0 ; tÞ=σ½Φðz; z 0 ; tÞ % const.
Note that the above idea and in fact the whole framework of quantum read-out formulated here is independent of the dimensionality of the Hamiltonian and can be applied, e.g., in two dimensions. There is also no restriction to continuum systems so lattice models can be treated similarly.
Experimental data analysis. Let us consider the state preparation procedure used in the recent experiment 8 where recurrent dynamics has been observed. Given an estimate of the average number of atoms per gas N Avg ' 3400 and the shape of the experimental boxlike potential, we can numerically obtain the average density profile n GP ðzÞ from the GP equation. This specifies the HamiltonianĤ (1) with R 1D ' 25 μm. Hence, we have the full information necessary to compute the eigenmode wave functions needed for the reconstruction procedure in Eq. (8). The number of relevant wave functions N % 10 can be upper bounded by considering the finite resolution of the interference images. In our case, the phase fluctuations Φðz; z 0 ; tÞ defined in Eq. (5) can be measured at points z; z 0 spaced by the pixel size of the camera δ % 2 μm 32 . In addition, other effects including diffraction limit the resolution. The measured values can be related to theoretical continuum predictions by implementing a real-space cut-off via a Gaussian convolution with standard-deviation σ % 3:5 μm (see Supplementary Note 2).
Initially the two adjacent gases whose relative phase fluctuations we are studying are strongly coupled. Hence, the state preparation is to a good approximation governed bŷ where the tunnel coupling term of strength J is pulling the relative phase field to zero (i.e., hcosðφÞi ' 1). The initial state can be expected to be a low-temperature thermal state ofĤ Ini . Following that preparation, the system is quenched by suddenly turning off the tunnel coupling. Experimentally, this is realized by separating the two gases over a time of 2 ms until J drops to zero. The middle of this ramp defines the initial time t ¼ 0 ms and the subsequent evolution under the HamiltonianĤ (1) is measured in steps Δt ¼ 2:5 ms. Based on the data from this initial dynamics, we can reconstruct the initial state at t ¼ 0 ms. In Fig. 2a, we plot the reconstructed a b c Fig. 2 Initial state reconstruction. Reconstruction of the full initial state right after decoupling at t ¼ 0 ms based on the phase correlations measured during the dephasing dynamics immediately after the quench (t = 1, 3.5, 6, 8.5, 11 and 13.5 ms). a Comparison of the measured phase correlations Φ Data ðz; z 0 ; tÞ to the reconstructed ones Φ Rec . The time slice presented in the foreground corresponds to t ¼ 1 ms. We find that the reconstruction yields good agreement with the data, evidenced by the respective weighted difference to the data ΔΦ ¼ jΦ Data À Φ Rec j=σ½Φ Data (right). b Reconstructed covariance matrix V of the initial state at t ¼ 0 for the eigenmodes j; k ¼ 1; ; 10. From left to right, the reconstructed phase-phase V ϕϕ , density-density V ρρ and phase-density correlations V ϕρ are plotted. The correlations V ϕϕ and V ρρ are close to diagonal in the numerically obtained wave functions f ϕ k , indicating that the eigenmodes of the system are well captured. For an initial thermal state of the Hamiltonian (11), the cross-correlations should vanish V ϕρ 0, here we find a small contribution. Note that the influence of higher-energy modes is suppressed by the limited spatial resolution in the experiment (see Supplementary Note 2). c Comparison of the diagonal elements of V ϕϕ (blue bullets) and V ρρ (red bullets) at t ¼ 0 with the predictions for a thermal state of the pre-quench Hamiltonian given in (11) (solid lines). The error bars correspond to the 80% confidence intervals obtained from a bootstrap analysis 35 . The thermal predictions are corrected for the suppression due to the finite imaging resolution, see Supplementary Note 2 and ref. 32 . We find the correlations of the first five modes to agree well with a thermal state at T ¼ 52 nK and J ¼ 2π 1:1 Hz obtained by a combined least-squares fit of hφ 2 k i and hδρ 2 k i. The strong suppression of the higher mode signals by the imaging renders a meaningful comparison impossible. ARTICLE COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-019-0273-y phase correlations, showing good agreement with the measured values signifying the consistency of our method. The corresponding covariance matrix of the full initial state is shown in Fig. 2b. Most importantly, note that we are indeed able to infer density fluctuations of the form hδρ j δρ k i. Using the mode transformation from Eq. (2), this information from the eigenmode-space can be also translated to real-space. However, many physical properties of the initial state can be directly extracted from the eigenmode correlations. We firstly observe that the blocks V ϕϕ and V ρρ are close to being diagonal. Hence, we find that the collective modes of the system are well captured by the numerically obtained wave functions f ϕ k . This supports the expectation that the initial state is thermal with respect to the pre-quench Hamiltonian (11). As the eigenmode wave functions are not strongly affected by the quench for the chosen trap geometry (see Supplementary Note 3), if the system was thermal with respect to the initial Hamiltonian, the reconstructed state should remain diagonal even when expressed in terms of the wave functions of the quench Hamiltonian.
On the other hand, we remark that allowing for off-diagonal correlations and cross correlations in V is necessary for an accurate reconstruction as otherwise the comparison in Fig. 2a is significantly worse. One reason for their presence can be small deviations between the assumed eigenmodes and the true eigenmodes of the system, e.g., due to pre-or post-quench trapping potential imperfections not included in the GP profile. Another reason could be that the initial state is genuinely out of thermal equilibrium, which would be interesting from the quantum information perspective in the context of the resource theory of coherence 36,37 .
Independent of these subtleties, our read-out method allows us to study how the energy is distributed in the system based on the measured out-of-equilibrium phase fluctuations. We now have access to the phonon occupation numbers given by n k ¼ 1 2 hφ 2 k þ δρ 2 k i À 1 2 and can study energy expectation values hĤi ¼ P N k¼1 ω k n k withĤ given in Eq. (1). More specifically, we can check from the observed data if the energy is distributed among the modes in a thermal way. The gas is prepared in a double well trap with large tunnel coupling and we expect the prepared state to be thermal with respect to the Hamiltonian (11). Based on Fig. 2b, we find that the reconstructed initial state shows significant suppression of the phase fluctuations which are an order of magnitude smaller than the density fluctuations. This is consistent with the initial energetic penalty on phase fluctuations. In Fig. 2c, we show the quantitative comparison of the reconstructed second moments hφ 2 k i and hδρ 2 k i compared to the thermal theory with the Hamiltonian (11). Due to the finite imaging resolution, we are only able to resolve the lowestlying eigenmodes. We find that their second moments agree with a thermal distribution of the coupled Hamiltonian (11). On the other hand, we also observed that the reconstructed initial state does not agree with the experimental state before the start of the decoupling ramp, as it leads to weaker phase locking than what was measured. This hints at the finite decoupling ramp having a significant influence on the correlations observed in the quench dynamics despite the tunnel coupling decreasing exponentially in the height of the barrier that is being ramped-up. The reconstruction hence extracts an effective initial state of the dynamics. If the physics of the initial Hamiltonian is of particular interest, then this effect might be diminished by performing a faster quench to the free system. On the other hand, let us remark that the physics of quenches of this type is key to the observation of generalized Gibbs ensembles in the considered setup 6 . In general, it is difficult to model the complete process of state preparation theoretically as it involves a strongly correlated phase of the sine-Gordon model out of equilibrium 7 and our method could offer new experimental insights.
Recurrent dynamics. With the reconstruction of the full state of the system, also its evolution beyond the interval of input times can be calculated. Propagating the covariance matrix V forward or backward in time via (9) allows us to pre-and redict the system's dynamics. In ref. 8 , this dynamics was visualized and quantified through the correlator Cðjz À z 0 j; tÞ ¼ cosφðz; tÞ Àφðz 0 ; tÞ The phase-locked initial state corresponds to C % 1 independent of the longitudinal separation z ¼ jz À z 0 j. In Fig. 3a, we show Cð z; tÞ obtained from the experimental data. Due to a linear dispersion relation and an equally spaced spectrum, the involved modes start to rephase after the inital dephasing dynamics leading to partial recurrences of the initial state 8 . Figure 3b-d shows how this rephasing dynamics can be predicted from the reconstructed states. For the reconstruction based on the initial dephasing dynamics, for example, we obtain a good qualitative prediction of the recurrences (red interval). Quantitative agreement is lost over time due to interaction effects between the modes. These interactions are mediated by higher-order terms beyond the effective model assumed in (1) and can therefore not be captured 25 . Nevertheless, the reconstruction method is robust enough such that we can obtain an accurate short-time prediction even using data that are seemingly fully dephased, i.e., where Cð z; tÞ is nearly indistinguishable from the correlations of a thermal state of the quench HamiltonianĤ (blue interval). However, in some cases (green interval), we find that statistical fluctuations can lead to large error bars, an effect reproduced by numerical simulations (see Supplementary Note 3). Note also that the last two intervals were intentionally chosen to be short, including only seemingly dephased data between the recurrences. They cover I ¼ 5 input times, during which the slowest eigenmode performs only about a quarter of a rotation. Therefore, the influence of the finite statistical sample size is more severe in these reconstructions.
Phonon occupation dynamics. Besides providing new insights into the state preparation, access to the full covariance matrix can enable entirely new ways of exploring the effects of interactions. The effective model given in (1) is obtained in a perturbative expansion of the Lieb-Liniger Hamiltonian up to second order 25 . For long evolution times, however, the dynamics can also be affected by the neglected terms that, e.g., can give rise to effects such as Beliaev-Landau damping 25 . It is challenging to obtain the rates of such processes by numerical calculations as interacting bosonic dynamics are notoriously difficult to treat and various approximations are necessary [38][39][40][41] . Therefore, it would be interesting to use the atom chip experiments to measure the damping rates and compare with theoretical predictions to validate different methods.
Here we show how the recovery method described above can be used to investigate these higher-order processes. To that end, we perform the recovery procedure for different input intervals of length I ¼ 8, with varying starting points. For each interval, we obtain an estimate of the central moments of phase and density fluctuations, hφ 2 k ðtÞi and hδρ 2 k ðtÞi, and calculate the phonon occupation numbers n k ðtÞ ¼ 1 2 hφ 2 k ðtÞ þ δρ 2 k ðtÞi À 1 2 . Scanning the starting point of the input interval through the measurement times allows us to investigate the dynamics of these observables, as shown in Fig. 4. The interval length is chosen long enough such that the slowest eigenmode picks up enough dynamical phase to ensure a stable reconstruction. At the same time, it is chosen short enough such that interactions between the modes do not influence the reconstruction.
The occupation numbers n k are constants of motion ofĤ. In Fig. 4a, we show their reconstructed dynamics for the five lowest eigenmodes. We find that overall the occupation numbers do not vary strongly and stay almost constant. This is expected as perturbations to the quench Hamiltonian should be negligible, and in any case they are irrelevant in the sense of the renormalization group. Note, however, that for different measurements with other system sizes, we find indications of a trend of slowly increasing mode occupations (see Supplementary Note 4). While the dynamics of occupation numbers is constant, Fig. 4b shows how at the same time the individual modes rotate between phase and density fluctuations. We find that this dynamics is clearly damped. This hints that the source of the recurrence damping observed in Fig. 3a and ref. 8 is a loss of the initial quadrature 'squeezing' hφ 2 k ð0Þi=hδρ 2 k ð0Þi ( 1 within each mode k rather than changes in their occupations. Our method makes it possible to extract mode-resolved damping rates of the collective excitations: In the future, using a b  k ðtÞi and hδρ 2 k ðtÞi for the first three modes k ¼ 1; 2; 3. We observe a gradual decay of the oscillation amplitude reflecting the apparent irreversible equilibration present in the system. The error bars indicate the 80% confidence intervals obtained from a bootstrap analysis 35 . The lines connecting the data points are a guide to the eye. Using this experimental dataset, we have performed reconstructions based on input data from a given time window (indicated by red, green and blue intervals) and we calculated the phase correlation functions Cð z; tÞ capturing the full spatio-temporal dynamics by evolving the reconstructed states to times beyond the input intervals. b The full spatio-temporal dynamics based on input data indicated by the red interval together with the cut at z c , which shows the quantitative comparison to the measured data. The red-coloured curve corresponds to the propagation of the reconstructed state based on initial times. c Similarly, we show the full spatio-temporal dynamics and a cut at z c obtained from the reconstructed state based on seemingly dephased data after the quench (green interval). d Likewise, based on the times between the first and second recurrences (blue interval), we are able to obtain a reliable reconstruction. While the dynamical prediction based on the initial dephasing dynamics works well, the reconstructions based on data taken in between the recurrences is limited by the finite sample size (green interval)-an effect reproduced by numerical simulations (see Supplementary Note 3). Note that the quantitative discrepancy at times far away from input intervals is due to terms of higher order not captured by the effective theory of Eq. (1). The shaded area around the curves, as well as the error bars of the data, indicate 80% confidence intervals obtained from a bootstrap analysis 35 . smaller time steps Δt and possibly non-equidistant measurement times should allow to study also higher modes and test theoretical predictions concerning the dynamics under perturbations to the non-interacting effective model.

Discussion
We have formulated and demonstrated the functioning of a new quantum read-out method for quantum simulators where we reconstruct the second moments of pairs of conjugated observables by measuring at different times only one of them. The developed scheme allows us to reliably reconstruct the covariance matrix of non-local low-energy excitations of a 1D superfluid based on experimental data from an atom chip experiment which makes phase measurements but does not directly access density fluctuations.
We found several interesting insights into the physics of the system. Firstly, the strong energetic penalty on phase fluctuations present during the state preparation is reflected in the reconstructed correlations as there are significantly less phase than density fluctuations. The reconstructed state is almost diagonal, which underlines that the eigenmodes before and after the quench are closely related and on a higher level demonstrates that the effective theory captures the relevant degrees of freedom of the system. A fit to a thermal model for the initial state allowed us to estimate the temparature and the effective tunnel coupling in the state preparation. In the considered setting, recurrences of the system have been recently observed 8 , which are due to an approximately linear spectrum of the phonons. We have demonstrated that our method can take input data from times when the system is seemingly dephased in order to predict recurrences by evolving the reconstructed covariance matrix in time, strongly underlining the predictive power of the obtained recovery scheme. Finally, we have studied the occupation numbers of the eigenmodes over time and obtained strong constraints on the rate of their growth. We have reconstructed the contribution of phase and density fluctuations over time and found that their oscillations are damped. We expect that a quantitative experimental assessment of possible reasons of the deviations from the non-interacting effective model will become possible by following the lines of this work.
Our work paves the way towards new intriguing experiments by giving access to quadrature operators which can be used as the basic ingredients for many quantum information processing protocols [42][43][44] . The method presented offers a novel window into quantum simulators, allowing to assess initial states, notions of entanglement and various other quantities previous read-out schemes did not allow for. It is our hope that our new quantum read-out method will enable exciting insights into the physics of ultra-cold superfluids, but also due to its generality that it will become a versatile tool used in state-of-the art quantum technologies allowing to fully use the power of the existing quantum simulation platforms 45 .
Note added. After the completion of the manuscript, we became aware of similar developments in the discrete setting of optical lattices with applications to topological band insulators 22,23,46 -it would be interesting to also include there our ideas of using semidefinite constraints ensuring that the reconstructed covariance matrix is physical and the recovery stable.

Data availability
Data are available upon reasonable request.

Code availability
Code is available upon reasonable request.