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, it is the read-out that limits such quantum simulations. In this work, we introduce a new paradigm of experimental read-out exploiting coherent non-interacting dynamics in order to extract otherwise inaccessible observables. Specifically, we present a novel tomographic recovery method allowing to indirectly measure second moments of relative density fluctuations in one-dimensional superfluids which until now eluded direct measurements. We achieve this by relating second moments of relative phase fluctuations which are measured at different evolution times through known dynamical equations arising from unitary non-interacting multi-mode dynamics. Applying methods from signal processing we reconstruct the full matrix of second moments, including the relative density 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 allow to explore static and dynamical properties of otherwise intractable quantum manybody systems. In many instances, however, it is the read-out that limits such quantum simulations. In this work, we introduce a new paradigm of experimental read-out exploiting coherent non-interacting dynamics in order to extract otherwise inaccessible observables. Specifically, we present a novel tomographic recovery method allowing to indirectly measure second moments of relative density fluctuations in one-dimensional superfluids which until now eluded direct measurements. We achieve this by relating second moments of relative phase fluctuations which are measured at different evolution times through known dynamical equations arising from unitary non-interacting multi-mode dynamics. Applying methods from signal processing we reconstruct the full matrix of second moments, including the relative density 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 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 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 non-equilibrium 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 one-dimensional 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 one-dimensional 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 one-dimensional 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.
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 lowenergy 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 )] = iδ(z − z ). 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 SM). 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 SM 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 {ω k , k = 1, 2, . . . } is discrete [28]. We denote eigenmode operators of phase and density fluctuations byφ k and δρ k respectively and use their corresponding wave Note that here we have [δρ k ,φ k ] = iδ k,k where δ k,k 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 SM 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,[30][31][32][33] Φ(z, z , t) = (φ(z, t) −φ(z 0 , t))(φ(z , t) −φ(z 0 , t)) .
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 ,r k ] = iΩ k,k with Ω = 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 Q(V ) = V + 1 2 iΩ 0 is positive-semidefinite, i.e., has non-negative eigenvalues. This condition reflects the Heisenberg uncertainty principle of canonically conjugated observables [34]. It will be convenient to use the notation . . .  (1) is decoupled on the level of momentum space operators {φ k , δρ k } involving multiple modes rotating at different frequencies (3). Measurements of real-space continuum fields {φ(z, t)} yield the referenced two-point correlation functions Φ(z, z , 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 non-local eigenmodes of the Hamiltonian. This approach is general and can be applied to any other system in which quenches to non-interacting multi-mode Hamiltonians are available.
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 realspace resolution in the experiment.
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 least squares recovery problem. For this we collect all measured values of Φ(z, z , t i ) at different points z and z 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 , t i ) sorted as in b via Eq. (8). The timeevolution 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 · 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 optimal solution to this convex quadratic problem with a semi-definite 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 semi-definite programming.
We use the package cvx [35], see the SM 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 ) 0. 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 en- 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 , t) grows typically for increasing spatial separations |z − z | but Φ(z, z , t)/σ[Φ(z, z , 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 box-like 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 , t) defined in Eq. (5) can be measured at points z, z spaced by the pixel size of the camera δ ≈ 2 µm [33]. 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 SM). 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 , t) (left) to the reconstructed ones ΦRec (center). 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 ∆Φ = |ΦData − ΦRec|/σ[Φ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 SM). (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 [36]. The thermal predictions are corrected for the suppression due to the finite imaging resolution, see SM and Ref. [33]. 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 φ 2 k and δρ 2 k . The strong suppression of the higher mode signals by the imaging renders a meaningful comparison impossible.
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., cos(φ) 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 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 δρ j δρ k .
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 SM), 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 Following reconstructions based on input data from a given time window we calculated the phase correlation functions C(z, t) capturing the full spatio-temporal dynamics (left). The cut atzc (right) shows the quantitative comparison to the measured data. The shaded area around the curves, as well as the error bars of the data, indicate 80% confidence intervals obtained from a bootstrap analysis [36]. The upper row corresponds to the propagation of the initial state reconstruction presented in Fig. 2 (red interval). The middle and lower row show the propagation of reconstructions based on seemingly dephased data (green and blue interval, respectively). 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 SM). 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 quantum information perspective in the context of the resource theory of coherence [37,38].
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 φ2 k + δρ 2 k − 1 2 and can study energy expectation values Ĥ = 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 φ2 k and δρ 2 k compared to the thermal theory with the Hamiltonian (11). Due to the finite imaging resolution we are only able to resolve the lowest-lying 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 generalised 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 preand redict the system's dynamics. In Ref. . We observe a gradual decay of the oscillation amplitude reflecting the apparent equilibration observed in Fig. 3. The error bars indicate the 80% confidence intervals obtained from a bootstrap analysis [36]. The lines connecting the data points are a guide to the eye. was visualized and quantified through the correlator The phase-locked initial state corresponds to C ≈ 1 independent of the longitudinal separationz = |z − z |. 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]. Fig. 3b 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 is 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 SM). 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 [39][40][41][42]. 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, φ2 k (t) and δρ 2 k (t) , and calculate the phonon occupation numbers n k (t) = 1 2 φ2 k (t) + δρ 2 k (t) − 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 H. 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 SM). 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' φ2 k (0) / δρ 2 k (0) 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 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 and outlook. 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 one-dimensional 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 noninteracting 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 [43][44][45]. 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 [46].
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 [23,32,47] -it would be interesting to also include there our ideas of using semi-definite constraints ensuring that the reconstructed covariance matrix is physical and the recovery stable.

References 7
A. Calculating the eigenmodes 9 1. Discretization of fields 9 2. Decoupling of the effective model using symplectic transformations D. Extended data 16 This supplemental material is structured as follows. In Sec. A begin by explaining how to obtain eigenmodes in which the dynamics is decoupled and is a simple rotation of the eigenmodes. In Sec. B we give more details concerning the reconstruction procedure and its implementation. In Sec. C we show the functioning of the recovery procedure on simulated Gaussian data. Finally, in Sec. D we show figures analogous to the main text but based on additional data for systems of different sizes.

Appendix A: Calculating the eigenmodes
In this section we describe in detail how we obtain the eigenmodes of the quench Hamiltonian which can be viewed as a CFT in curved space-time background whenever the GP profile is not homogeneous. In this case, a discretization of fields allows to approximate the low-lying eigenmodes of the continuum Hamiltonian by the eigenmodes of a Hamiltonian involving a finite number of degrees of freedom which are the average fields in a given discretization cell. We then show how to diagonalize the coarse-grained Hamiltonian numerically taking into account that the quench Hamiltonian has a zero-mode. Finally, we describe how to use the numerically obtained wavefunctions for finitely many modes as an approximation to the corresponding eigenmodes in the continuum limit.
The Hamiltonian describing the quench dynamics is functionally parametrized by the GP profile n GP . Due to transverse broadening of the wave functions [8] the density-density interaction is functionally dependent on the GP profile and reads g(z) = ω ⊥ a s (2 + 3a s n GP (z))/(1 + 2a s n GP (z)) 3/2 (A1) where ω ⊥ is the radial trapping frequency and a s is the scattering length [26]. Hence, by knowing the GP profile, we know the Hamiltonian and so we can find the eigenmodes. Here we show how to do this even if the GP profile is not homogeneous n GP = const

Discretization of fields
We want to find approximations to the wave functions and eigenmodes discussed above by discretizing the interval [−R 1D , R 1D ] into N pixels, each of size 2R 1D /N . Fixing N , for l = 1, . . . , N + 1 the coordinates of the discretization lattice read z l = −R 1D + 2R 1D l−1 N and we define discretization pixels which are the closed intervals p l = [z l , z l+1 ] for l = 1, . . . , N . We then introduce the discretized operators as the integration of the field operators viâ with ∆z := |p l | = 2R 1D /N . Following Refs. [43,44], these discretized operators yields a vector of canonical coordinateŝ satisfying the bosonic canonical commutation relations [Q j ,Q k ] = iΩ j,k /∆z where Ω = 0 1 1 N −1 1 N 0 . as can be verified easily. Observe that the right-hand side will yield a Dirac delta in the continuum limit N → ∞. The discretization of the effective model will be a quadratic operator in the discretized modesφ (N ) l and δˆ (N ) l which can be efficiently diagonalized using single particle transformations only as we want to explain below in the next section.
2. Decoupling of the effective model using symplectic transformations Using the general notation of quadraturesQ, we consider quadratic Hamiltonians of the form where H = H ∈ R 2N ×2N are the couplings We will assume that H is positive semi-definite, i.e., H 0 and that there is no coupling between the phases and densities in the effective model and all Hamiltonians considered in this work will have this property. In this case the couplings H will be block diagonal and we will use the decomposition To discretize the integral we define the geometric mean η l = (n GP (z l )n GP (z l+1 )) 1/2 for l = 1 . . . , N which giveŝ From this we read off Depending on the detail of the simulation g(z) ≈ g(0) can be taken constant or H ρ may include the pressure term as discussed above with a similar discretization scheme. With this notation, we obtain Starting from a set of canonical coordinatesQ thenr = SQ for S ∈ R 2N ×2N will again denote a vector of canonically commuting operators if S is symplectic, i.e., it fulfills S Ω S T = Ω (A13) which can be seen by explicitly checking thatr again fulfills [r j ,r k ] = iΩ j,k /∆z. In view of diagonalizing the Hamiltonians of interest, it is important to note that matrices of the form S = Q ⊕ Q for any orthogonal Q ∈ O(N ) as well as S = A ⊕ A −1 for any invertible A ∈ GL(N, R) that is symmetric, i.e., A T = A are both symplectic matrices and that the inverse as well as the product of symplectic matrices are again symplectic. We can then diagonalize Hamiltonians of the form as given in Eq. (A5) under the assumption that H ρ is invertible. This property allows us to define a symplectic matrix The matrix of the phase couplings in the new coordinates readsH φ = (H T ρ ) 1/2 H φ (H ρ ) 1/2 and is again real and symmetric such that it can be diagonalized by an orthogonal transformation Q ∈ O(N ) withH φ = QΣQ T . Here, Σ is diagonal and we assume that all zero eigenvalues are sorted to the first N 0 ≥ 0 positions, i.e., Σ = 0 N0 ⊕Σ withΣ 0 diagonal and we define the eigenfrequencies ω viaΣ 1/2 = diag(ω N0+1 , . . . , ω N ). With the diagonal matrix Σ φ = 1 1 N0 ⊕Σ and the transformation That is, in the canonical coordinates (φ we have that the Hamiltonian in Eq. (A5) takes the formĤ such thatφ (N ) j ≈φ j and δρ (N ) j ≈ δρ j asĤ ≈Ĥ N . We will therefore not distinguish betweenφ and δρ j outside of this section.

Discrete approximations
With this, we can read off the discrete approximation to the wave functions f φ k and f ρ k relatingφ (N ) andφ (N ) or correspondingly δρ (N ) δˆ (N ) as the rows of S = S 1 S 2 which is of block structure, i.e.,S = S φ ⊕ S ρ . Specifically we find and Note that when relatingφ (N ) andφ (N ) or δρ (N ) and δˆ (N ) we included a factor √ ∆z. The inclusion of this factor allows to change the commutation relations from [δˆ k ] = iδ j,k → iδ j,k as one would expect from the discrete canonical eigenmodes of the system. Let us furthermore observe that the relation to the real-space correlators are given by where we exploited that S has a block-diagonal structure and the inverse scaling in the discretization step ∆z should be noted.  In this section, we describe the data analysis and formulate the reconstruction procedure with additional details. As described in Refs. [7,8,31] through matter-wave interferometry phase profiles ϕ of the superfluid can be measured. After the two superfluids were coupled with tunneling strength J ≈ 3.5 Hz for sufficiently long times, the separation potential is increased rapidly in about 1 ms. The end of the ramp defines the initial time t 0 = 0 ms for the quench evolution that follows. The gas can then be held for a specific time t which defines the time during which the system evolves under the quench Hamiltonian in Eq. (1). Here we focus on the referenced second moments obtained from the profiles, for which the corresponding physical observable is where z 0 denotes a fixed reference point in the system. For each hold time t about n Sample ≈ 200 phase profiles were obtained, this number is limited by the stability of the setup, but can be increased if fewer hold times are considered in total. Fig. 6 shows an example of two of many phase profiles used in the analysis. In the data analysis, data from N p = 19 central pixels is used which corresponds to about 60% of the cloud as each pixel has size = 1.95 µm, so the total size observed is about 37 µm. The positions of the pixels are z a = (a − 10) (B2) and the reference point is chosen as z 0 = 0 µm, i.e., in the middle of the cloud. We consider referenced second moments here, as by this we are able to consistently remove any offset phase between two different measured profiles. In practice we subtract at each pixel a the central phase profile value ϕ (i) (z 0 ). The experimental estimate for Eq. (B1) is then In terms of the eigenmodes this reads Using the expressions for the time evolution, we get It must be noted that the measured value at a pixel z a does not exactly reflect the value of the fieldφ(z a ) but rather a convolution of the field with a Gaussian function, i.e., it probes the value averaged over a patch of specific characteristic length. More precisely, the experiment allows us only access to measurements of and we define the correspondingly convoluted wave functioñ For the considered experimental setup we find the estimation σ ≈ 3 µm. In order to include the convolution in the reconstruction it then suffices to use the convoluted wave functions and set

Recovery procedure
In the implementation, we vectorize the covariance matrix V such that each block is a vector, i.e., v φφ = vec(V φφ ) etc. and define v = v φφ ⊕ v φρ ⊕ v ρρ ∈ R 3N 2 and use the notation v = vec(V ). Formula (B5) shows that for each input (z a , z b , t) we can find a vector w ∈ R 3N 2 such that Φ(z a , z b , t) = w T v. For a fixed time t i k we then collect the from the measured data extracted second moments Φ est (z a , z b , t i ) in a vector b k ∈ R N 2 p and construct the corresponding vectors w and collect them as rows in a matrix A k ∈ R N 2 p ×3N 2 . Doing this for all n times t i1 , . . . , t in which are used for the reconstruction as input, we then stack all b k and A k into a large vector b matrix A correspondingly, i.e., We furthermore define at each time step a diagonal matrix W ∈ R N 2 p ×N 2 p which contains the inverse statistical errors of the experimental measurement of the second moments W (k) (a,b),(a,b) = 1/σ est (Φ(z a , z b , t k )) and collect all W (k) in on large blockdiagonal matrix W in order to define a more uniform target function for the optimization. With this definition we aim at minimizing the vector Hilbert Schmidt-norm subject to the semi-definite constraint where in the main text we have introduced the notation A(V ) = A vec(V ). The numerical reconstruction has been implemented with use of the cvx package. The standard theory of semi-definite programming shows that there is always a unique solution v Opt to this optimization problem. Unfolding the vectorization yields the reconstructed covariance matrix V Opt . As a final remark, it is interesting to note that positivity constraints (imposing that the density operator is positive semi-definite) similar to the Heisenberg constraint (reflecting the Heisenberg uncertainty principle as a semi-definite constraint) characterizing bosonic covariance matrices can significantly increase stability of least squares reconstructions [48]. In fact, wide classes of recoveries with a positivity constraint [48] can be interpreted as compressed sensing schemes [49]. In this context, is important to stress that the semi-definite constraint V + 1 2 iΩ 0 readily implies that V > 0, so that V is strictly positive, so that the constraint of Ref. [48] is readily enforced. Hence, it is interesting to see that much of the intuition on the positive cone for density operators carries over to the Heisenberg cone for covariance matrices. Further explorations of seeing our scheme as a compressed sensing scheme will be left to future work.
Appendix C: Simulation of the reconstruction procedure Various aspects of the reconstruction can be modeled by considering a thermal state of the effective Hamiltonian of the strongly coupled condensates as discussed in the main text. Thermal correlations in the discretized model can be obtained either by considering the exact formulas from Refs. [50,51] or by classical phase approximation [52]. In the following we consider the latter and study the effects of finite sample size and finite measurement resolution. We denote the real-space phase fluctuation functions by and in classical field approximation we can calculate these via Γ φφ ≈ (H φ (J = 0)) −1 /k B T [52]. Together with the correlations for density fluctuations, we can propagate these under the quench Hamiltonian. Fig. 7 shows the correlations Γ φφ and additionally the effect of the referencing which removes the running phase, of the convolution which comes from the measurement resolution and finally discretization due to a finite amount of pixels.

Influence of finite sample size
At each time the phase correlations can be used as a positive matrix parametrizing a (classical) Gaussian distribution from which single-shot random profiles can be sampled. After convolution these can be seen to correspond to the direct observable of the interferometry which is useful to assess the systematic imperfections of our procedure. Here we study two possible aspects. By resampling with n Sample = 200 and n Sample = 2000 we study the sensitivity of our reconstructions to statistical fluctuations in the experiment. Secondly, we study the real-space correlations and the convolution of these to see how finite measurement resolution impacts the information that can be obtained from our procedure.
We have simulated theĤ N for N = 400 and J/ = 2π × 3.5 Hz obtaining the correlation functions of the thermal state at T = 40 nK. Using the phase-phase correlation functions at different times we have resampled the profiles. After referencing the profiles the resampled phase-phase correlation functions were obtained. We then perform a recovery at the times indicated in Fig. 8 which shows that finite sample size of about 200 experimental runs for each time constrains the possibility of recovering the phase locking in its full extent. From left to right, the first plot shows the rotation of the real space covariance matrix Γ φφ j,k = φ(zj)φ(z k ) defined at discretization pixels rotated into the eigenmode space S φ Init of the initial Hamiltonian. On all plots we show the lowest-lying eigenmodes and do not show the zero-energy mode considering only the eigenmodes relevant in the experiment. Second plot presents the Gaussian convolution of Γ φφ denoted byΓ φφ which takes into account the finite measurement resolution in the eigenmode space of the initial Hamiltonian. Observe, that in both cases we obtain diagonal matrices, but the convolution introduces a cut-off for resolving the occupation of the higher modes. The last two plots show the same comparison of Γ φφ andΓ φφ but now rotated with the eigenmodes S φ Quench of the quench Hamiltonian. Observe, that the eigenmode occupations are rearranged due to a different mode transformation (nGP = const) and minor coherences are introduced by the convolution.

Influence of finite measurement resolution
The measured phase fluctuations at a given pixel are in fact a convolution which spreads into the neighboring pixels too. Its primary effect is introducing a frequency cutoff, as the higher modes oscillate quickly and are averaged out. The secondary effect is introducing non-universal additional artifacts into the correlations that come from the convolution and are a feature of the measurement setup. We consider the thermal real-space covariance matrix. It is diagonal in the eigenmodes of the initial strongly coupled Hamiltonian, see Fig. 10 left. The convolution introduces however new correlations that are not present in the state and are an artifact of the coarse-graining, Fig. 10 second from the left. If we consider the post-quench Hamiltonian, the modes change slightly due to the non-homegenous GP profile and the covariance matrix is slightly off-diagonal in these modes Fig. 10 second from the right. After the convolution again a cut-off is introduced, but also additional stray artifacts, Fig. 10 right.
Thus, the convolution will introduce in an uncontrolled way additional artifacts at different times and hence the measured real-space second moments of phases will have a discrepancy incorporated by the finite measurement resolution. This explains why in the data analysis, and the resampling simulation above, a perfect reconstruction is impossible. We conclude that the sample size is a smaller limitation than the finite experimental resolution.

Appendix D: Extended data
Here we give further results on additional experimental scans that were performed in the study of revivals in Ref. [8]. In total we consider 5 systems (one of them already presented in the main text) with varying system size and particle number -the corresponding values are listed in Tab I. Table of the experimental scans performed characterized by the system size specified via the external trap and the average particle number per well. The fluctuations of the particle number is about δN Avg Well ≈ 50 independent of the system size. These two parameters together with the harmonic longitudinal ω l = 2π × 7 Hz and radial ω ⊥ = 2π × 1400 Hz trapping frequencies allows to calculate the Gross-Pitaevskii profile nGP and hence parametrize the effective model. The first scan is the one presented in the main text.
In Fig. 11 we show the reconstructed covariance matrices of the initial state as described in the main text for all 5 experimental scans listed in Tab. I. Consistently with the result presented in the main text, we find that the reconstructed covariance matrices are close to being diagonal with significant squeezing suppressing phase fluctuations and enhancing density fluctuations. Furthermore, in Fig. 12 we show the covariance matrix obtained for the first scan considering wave functions convoluted with a Gaussian distribution (B7). Here we find similar V j,k to the covariance matrix shown in the main text for low lying modes but it is noticeable that the higher modes are populated with no clear decay tendency. This can be explained by overfitting noise as the convolution of wave functions for large k vanishes f φ k (z) ≈ 0. Indeed, examining (B10) we find that the optimizer is not sensitive to changes of V j,k with j, k above the effective cut-off, or in other words the least squares recovery becomes numerically ill-conditioned. We did not observe any significant improvement in quantitatively predicting the revivals using convoluted modes.
Secondly, we investigate in Figs. 13 and 14 the correlator C defined in the main text based on data obtained in scans 2 to 5. We show the values extracted from the experimental measurement as well as the results obtained from three reconstructions with different input intervals. The results are consistent with the ones presented and discussed in the main text. The experimental data shows a slow dephasing and weakening of the initial phase locking. The reconstructions are able to recover and predict the signal well if the reconstruction interval includes a recurrence. Reconstructions from dephased data (reconstruction regions a) and c) ) are able to qualitatively describe the system but fail to predict quantitatively for instance the strength of the recurrence. FIG. 11. We show as in Fig. 2 in the main text the blocks of the covariance matrix of the initial state reconstructed with N = 10 modes for all experimental scans. The system sizes and particle number per well are as given in Table I. FIG. 12. We show as in Fig. 11 the reconstructed covariance matrix of the initial state of the first scan, but using convoluted eigenfunctions. As modes with increasing energy display more and more oscillations, the convoluted modes with higher energy become smaller in amplitude once the convolution starts to average over a full oscillation. This renders the least squares less stable and higher modes can have large occupation numbers without changing the real-space correlations because the mode functions are suppressed by the convolution.  Fig. 3 from the main text, we show the correlator C calculated from the experimental data and based on three reconstructions with varying input windows (indicated with dashed boxes) for scan number 2 (left) of the largest system size and number 3 (right) of the smallest system size. For the scan number 2 we have moved input window a) to an earlier time which results in very accurate reconstruction of the revival which would not be the case after moving the input window to a later time by one unit ∆t. Note, however, that the extrapolation works well which indicates that our method given enough input can yield very good results even with relatively small values of the dynamical phase. For the scan number 3 we can reconstruct reliably in the regions between the revivals. Note that in both cases input window c) does not yield strong reconstructed revivals but they are timed well and also in the experimental data the second revivals are not pronounced.  Table I and bootstrap error-bars by resampling the phase profiles nBootstrap = 500 times. Note that the fluctuations of the reconstructed occupation of the first modes increase with increasing system sizes which shows that it is important that all the modes acquire enough dynamical phase. Note that often the jumps in the occupation numbers coincide with the input intervals being placed in regions between the revivals where the reconstruction is difficult because of enhanced phase fluctuations due to the in-rotated density fluctuations. We have checked that taking a larger number of input times I does smoothen the occupation numbers but then the size of the input window is large enough so that interaction effects may start playing a role and the value of the occupation numbers need not be accurate.