Unitary dynamics of strongly-interacting Bose gases with time-dependent variational Monte Carlo in continuous space

We introduce time-dependent variational Monte Carlo for continuous-space Bose gases. Our approach is based on the systematic expansion of the many-body wave-function in terms of multi-body correlations and is essentially exact up to adaptive truncation. The method is benchmarked by comparison to exact Bethe-ansatz or existing numerical results for the integrable Lieb-Liniger model. We first show that the many-body wave-function achieves high precision for ground-state properties, including energy and first-order as well as second-order correlation functions. Then, we study the out-of-equilibrium, unitary dynamics induced by a quantum quench in the interaction strength. Our time-dependent variational Monte Carlo results are benchmarked by comparison to exact Bethe ansatz results available for a small number of particles, and also compared to quench action results available for non-interacting initial states. Moreover, our approach allows us to study large particle numbers and general quench protocols, previously inaccessible beyond the mean-field level. Our results suggest that it is possible to find correlated initial states for which the long-term dynamics of local density fluctuations is close to the predictions of a simple Boltzmann ensemble.


I. INTRODUCTION
The study of equilibration and thermalization properties of complex many-body systems is of fundamental interest for many areas of physics and natural sciences [1]. For systems governed by classical physics, an exact solution of Newton's equations of motion is often numerically feasible, using for instance moleculardynamics simulations. For quantum systems, the mathematical structure of the time-dependent Schrödinger equation is instead fundamentally more involved. Quantum Monte Carlo algorithms, the de facto tool for simulating quantum many-body systems at thermal equilibrium [2][3][4], cannot be directly used to study timedependent unitary dynamics. Out-of-equilibrium properties are then often treated on the basis of approximations drastically simplifying the microscopic physics. Irreversibility is either enforced with an explicit breaking of unitarity, e.g. within the quantum Boltzmann approach, or the dynamics is reduced to mean-field description using time-dependent Hartree-Fock and Gross-Pitaevskii approaches. Although these approaches may qualitatively describe thermalization [5,6], their range of validity cannot be assessed because genuine quantum correlations and entanglement are ignored.
For specific systems, exact dynamical results can be derived. This is the case for integrable 1D models, for which Bethe ansatz (BA) solutions exist [7]. However, also in this case many open questions still persist. For example, the exact evaluation of correlation functions for out-ofequilibrium dynamics is at present an unsolved problem. As a result, despite important theoretical and experimental progress [8][9][10][11][12][13][14][15], a complete picture of thermalization (or its absence), e.g. based on general quench protocols [16], is still missing.
Numerical methods for strongly-interacting systems face important challenges as well. Numerical Renormalization Group (NRG) and Density-Matrix Renormalization Group (DMRG) approaches provide an essentially exact description of arbitrary 1D lattice systems in-and out-of-equilibrium [17][18][19][20], but they have less predictive power when applied to continuos-space systems. On one hand, multi-scale extensions of the DMRG optimization scheme to the limit of continuous-space lattices [21] are so-far limited to relatively small system sizes [22]. On the other hand, efficient ground-state optimization schemes for continuous quantum field matrix product states (c-MPS) [23] have been introduced only very recently [24], and applications to quantum dynamics are still to be realized. A further formidable challenge is the efficient ex-tension of these approaches to higher dimensions, which is a fundamentally hard problem.
Another class of methods for strongly-interacting systems is based on variational Monte Carlo (VMC), combining highly-entangled variational states with robust stochastic optimization schemes [25]. Such approaches have been successfully applied to the description of continuous quantum systems, in any dimension and not only in 1D [26,27]. More recently, out-of-equilibrium dynamics has become accessible with the extension of these methods to real-time unitary dynamics, within timedependent variational Monte Carlo (t-VMC) [28,29]. So far, the t-VMC approach has been developed for lattice systems with bosonic [28,29], spin [30][31][32], and fermionic [33] statistics, yielding a description of dynamical properties with an accuracy often comparable with MPS-based approaches.
In this Paper, we extend t-VMC to access dynamical properties of interacting quantum systems in continuous space. Our approach is based on a systematic expansion of the wave function in terms of few-body Jastrow correlation functions. Using the 1D Lieb-Liniger model as a test case, we first show that the inclusion of highorder correlations allows us to systematically approach the exact BA ground state energy. Our results improve by orders of magnitude on previously published VMC and c-MPS results, and are in line with latest state-ofthe-art developments in the field. We further compute single-body and pair correlation functions, hardly accessible by current BA methods. We then calculate the time evolution of the contact pair correlation function following a quench in the interaction strength. For the non-interacting initial state, we benchmark our results to exact BA calculations available for a small number of bosons and further compare to the quench action approach for large systems approaching the thermodynamic limit. We finally apply our method to the study of general quenches from arbitrary initial states, for which no exact results in the thermodynamic limit are currently available.

A. Expansion of the Many-Body Wave-Function
Consider a non-relativistic quantum system of N identical bosons in d dimensions, and governed by the firstquantization Hamiltonian where v 1 ( x) and v 2 ( x, y) are, respectively, a one-body external potential and a pair-wise inter-particle interaction [34]. Without loss of generality, a time-dependent N -body state can be written as Φ(X, t) = exp [U (X; t)], where X = x 1 , x 2 . . . x N is the ensemble of particle positions and U is a complex-valued function of the Nparticle coordinates, R N ×d → C. Since the Hamiltonian (1) contains only two-body interactions, it is expected that an expansion of U in terms of few-body Jastrow functions containing at most m-body terms, rapidly converges towards the exact solution. Truncating this expansion up to a certain order M ≤ N , leads to the Bijl-Dingle-Jastrow-Feenberg expansion [35][36][37][38] where u m (r; t) are functions of m particle coordinates, r = x i1 , x i2 . . . x im , and of the time t. A global constraint on the function U (X, t) is given by particle statistics.
In the bosonic case, we demand that U ( . . x σ(N ) ), for all particle permutations σ.
In general, the functions u m (r; t) can have an arbitrarily complex dependence on the m particle coordinates, which can prove problematic for practical applications. Nonetheless, a simplified functional dependence can often be imposed, resulting from the two-body character of the interactions in the original Hamiltonian. For m ≥ 3, u m (r; t) can be conveniently factorized in terms of general two-particle vector and tensor functions following Ref. [26]. Details of this approach and the present implementation are presented in Appendix A and F. An appealing property of the many-body expansion (2) is that it is able to describe intrinsically non-local correlations in space. For instance the two-body u 2 ( x i , x j ; t), as well as any m-body function u m (r; t), can be long range in the particle separation | x i − x j |. This non-local spatial structure allows for a correct description of gapless phases, where a two-body expansion may already capture all the universal features, in the sense of the renormalization group approach [39]. This is in contrast with the MPS decomposition of the wave-function, which is intrinsically local in space.

B. Time-Dependent Variational Monte Carlo
The time evolution of the variational state (2) is entirely determined by the time-dependence of the Jastrow functions u m (r; t). In order to establish optimal equations of motion for the variational parameters, we start noticing that the functional derivative of U (X; t) with respect to the variational, complex-valued, Jastrow functions u m (r; t), yields the m-body density operators The expectation values of the operators ρ m over the state |Φ(t) give the instantaneous m-body correlations. For instance, ρ 2 (r 1 , , is proportional to the two-point density-density correlation function. We can then express the time derivative of the truncated variational state U (M ) (X; t) using the functional derivatives, Eq. (3), as a sum of the few-body density operators up to the truncation M , i.e.
The exact wave-function satisfies the Schrödinger equa- is the so-called local energy. The optimal time evolution of the truncated Bijl-Dingle-Jastrow-Feenberg expansion (2) can be derived imposing the Dirac-Frenkel time-dependent variational principle [40,41]. In geometrical terms, this amounts to minimizing the Hilbert-space norm of the residuals loc (X; t) , thus yielding a variational many-body state as close as possible to the exact one [42]. The minimization can be performed explicitly and yields a closed set of integro-differential equation for the Jastrow functions u m (r; t) : . (6) In practice, these equations are numerically solved for the time derivatives ∂ t u p (r ; t) at each time step t. The expectation values taken over the time-dependent state, · t , which enter Eq. (6), are found via a stochastic sampling of the probability distribution Π(X, t) = |Φ(X, t)| 2 . This is efficiently achieved by means of the Metropolis-Hastings algorithm, as per conventional Monte Carlo schemes (See Appendix B for details). It then yields the full time evolution of the truncated Bijl-Dingle-Jastrow-Feenberg state (2) after time integration.
The t-VMC approach as formulated here provides, in principle, an exact description of the real-time dynamics of the N -body system. The essential approximation lies however in the truncation of the Bijl-Dingle-Jastrow-Feenberg expansion to the M most relevant terms. In practical applications, the M = 2 or M = 3 truncation is often sufficient. Systematic improvement beyond M = 3 is possible [26], but may require a substantial computational effort.

III. LIEB-LINIGER MODEL
As a first application of the continuous-space t-VMC approach, we consider the Lieb-Liniger model [44]. On one hand, some exact results and numerical data are available, allowing us to benchmark the Jastrow expansion and the t-VMC approach. On the other hand, several aspects of the out-of-equilibrium dynamics of this model are unknown, which we compute here for the first time using t-VMC.
The Lieb-Liniger model describes N interacting bosons in one dimension with contact interactions. It corresponds to the Hamiltonian (1) with v 1 (x) = 0 and v 2 (x, y) = gδ(x − y), where g is the coupling constant. Here, we consider periodic boundary conditions over a ring of length L, and the particle density n = N/L. The density dependence is as usual expressed in terms of the dimensionless parameter γ = mg/ 2 n. The Lieb-Liniger model is the prototypal model of continuous onedimensional strongly-correlated gas exactly solvable by Bethe ansatz [44]. This model is experimentally realized in ultracold atomic gases strongly confined in onedimensional optical traps, and several studies on out-ofequilibrium physics have been already realized [14,[45][46][47][48][49].
As a result of translation invariance, we have u 1 (x; t) = const, and the first non-trivial term in the many-body expansion (2) is the two-body, translation invariant, function u 2 (|x − y| ; t). To compute the functional derivatives of the many-body wave function, we proceed with a projection of the continuous Jastrow fields u m (r; t) onto a finite basis set. Here, we have found convenient to represent both the field Hamiltonian (1) and the Jastrow functions, onto a uniform mesh of spacing a, leading to n var = L 2a variational parameters for the 2-body Jastrow term. In the following our results are extrapolated to the continuous limit, corresponding to a → 0. The finitebasis projection as well as the numerical time-integration of Eq. (6) are detailed in the Appendix C and benchmarked against exact diagonalization results in the three particle case in Appendix E.

A. Ground-state properties
To assess the quality of truncated Bijl-Dingle-Jastrow-Feenberg expansions for the Lieb-Liniger model, we start our analysis considering ground-state properties. An exact solution can be found from the Bethe ansatz (BA) and gives access to exact ground state energies and local properties [44]. Other non-local properties are substantially more difficult to extract from the BA solution, and unbiased results for ground state correlation functions have not been reported so far. In order to determine the best possible variational description of the ground-state within our many-body expansion, different strategies are possible. A first possibility is to consider the imaginarytime evolution |Ψ(τ ) = e −τ H |Ψ 0 which systematically AG is the parametrized 2-body Jastrow state of Ref. [43], c-MPS1 results are from Ref. [23], and c-MPS2 are those very recently reported in [24]. Distances r are in units of the inverse density 1/n. Our variational results have been obtained for N = 100 particles. Finite-size corrections on local quantities are negligible, and very mildly affect the reported large-distances correlations. Overall statistical errors are of the order of symbol sizes, for (a), and line widths, for (b) and (c).
converges to the exact ground-state in the limit τ ∆ 1 , where ∆ 1 = E 1 −E 0 is the gap with the first excited state on a finite system and provided that the trial state Ψ 0 is non-orthogonal to the exact ground state. Imaginarytime evolution in the variational subspace can be implemented considering the formal substitution t → −iτ in the t-VMC equations (6). The resulting equations are equivalent to the stochastic reconfiguration approach [50]. However, direct minimization of the variational energy can be significantly more efficient, in particular for systems becoming gapless in the thermodynamic limit, where ∆ 1 ∼ poly(1/N ). Given the gapless nature of the Lieb-Liniger model, we have found computationally more efficient to adopt a Newton method to minimize the energy variance [51].
For the ground state, the many-body expansion truncated at M = 2 is exact not only in the non-interacting limit γ = 0 but also in the Tonks-Girardeau limit γ → ∞. In this fermionized limit the wave-function can be written as the modulus of a Vandermonde determinant of plane waves, corresponding to the two-body Jastrow function u 2 (r) = log (sin rπ/L) [52]. To assess the overall quality of pair wave functions for ground state properties, we start comparing the variational ground-state energies E obtained for M = 2 with the exact BA result. In Fig. 1(a) we show the relative error ∆E/E as a function of in the interaction strength γ. We find that the relative error is lower than 10 −4 for all values of γ and the accuracy of our two-body Jastrow function is superior to previously published variational results based either on c-MPS [23] or VMC [43] [see Fig. 1(a) for a quantitative comparison]. Notice that the improvement with respect to previous VMC results is due to the larger variational freedom of our u 2 (r) function, which is not restricted to any specific functional form as done in Ref. [43].
Even though the accuracy reached by the two-body Jastrow function may already be sufficient for most practical purposes for all values of γ, we have also considered higher order terms with M = 3, shown in Fig. 2(a). The introduction of the third-order term yields a sizable improvement such that the maximum error is about three orders of magnitude smaller than original c-MPS results [23], and feature a similar accuracy of recently reported c-MPS results [24]. Overall, our approach reaches a precision on a continuous-space system, which is comparable to state-of-the-art MPS/DMRG results for gapless systems on a lattice [53].
Finally, to further assess the quality of our ground state ansatz beyond the total energy, we have also studied non-local properties of the ground state wave function, which are not accessible by existing exact BA methods. In Figs. 1(b) and (c) we show, respectively, our results for the off-diagonal part of the one-body density matrix, g 1 (r) ∝ Ψ † (r)Ψ(0) , and for the pair correlation function, g 2 (r) ∝ Ψ † (r)Ψ(r)Ψ † (0)Ψ(0) , where Ψ(r) is the bosonic field operator. We find an overall excellent agreement with the results that have been obtained with c-MPS in Ref. [24], except for some small deviations at large values of r which we attribute to residual finite-sizeeffects in our approach. We found that the addition of the 3-body terms does not change significantly correlation functions. Already at the 2-body level, the present results are statistically indistinguishable from exact results obtained using our implementation [54] of the worm algorithm [55] and for the same system (not shown).

B. Quench dynamics
Having assessed the quality of the ansatz for local and non-local ground-state properties, we now turn to the study of the out-of-equilibrium properties of the Lieb-Liniger model. We focus on the description of the unitary dynamics induced by a global quantum quench of the interaction strength, from an initial value γ i to a final value γ f . Exact BA results are available only in the case of a non-interacting initial state (γ i = 0). Even in this case, the dynamical BA equations can be exactly solved only for a modest number of particles with further truncation in the number of energy eigen-modes [56,58], N 10, since the complexity of the BA solution increases exponentially with the number of particles. Simplifications in the thermodynamic limit are exploited by the quench action [59], and have been recently applied to quantum quenches starting from a non-interacting initial state [58]. In the following we first compare our t-VMC results to these existing results, and then present new results for quenches following a non-vanishing initial interaction strength.
To assess the quality of the time-dependent wave function we compare our results for the evolution of local density-density correlations, g 2 (0, t), with the truncated BA results obtained in Refs. [56,57] for a small number of particles, N 6. Appendix E provides also further validation of our method accessing g 2 (r, t) at nonvanishing distances. The comparison shown in Fig. 2(a) shows an overall good agreement. The t-VMC and BA results are indistinguishable for weak interactions (γ f = 1). For larger interactions, we notice systematic but small differences between BA and t-VMC with M = 2 or M = 3. These differences amount to a small increase in the amplitude of the oscillations. This effect tends to increase with the interaction strength, being hardly visible for γ = 1 and more pronounced for γ = 10. However, these oscillations result at large times from the discrete mode structure due to the very small number of particles. They vanish in the physical thermodynamic limit. In turn, the comparison at small particle numbers indicates an accuracy better than a few percent for timeaveraged quantities in the asymptotic large time limit up to γ = 10, with results at M = 3 systematically improving on the M = 2 case. Concerning small and intermediate time scales, we do not observe systematic deviations between the t-VMC results and the BA solution. In particular, the relaxation times are remarkably well captured by the t-VMC approach. On the basis of this comparison and of the comparison for three particles with exact diagonalization results presented in Appendix E, we conclude that t-VMC allows accurate quantitative studies of both the relaxation and equilibration dynamics. This careful benchmarking now allows us to confidently apply the t-VMC approach regimes that are inaccessible to exact BA namely large but finite N and long times, as well as the case of non-vanishing initial interactions.
Let us consider relaxation of density correlations for a large number of particles, close to the thermodynamic limit (here we use N = 100). As shown in Fig. 2(b), we notice that the amplitudes of the large-time oscillations, attributed to the discrete mode spectrum, are now drastically suppressed compared to the quenches with N = 6. After an initial relaxation phase, the quantity g 2 (0, t) approaches a stationary value. Comparing our curves with those obtained with the quench action method, we find a qualitatively good agreement, albeit a general tendency to underestimate the quench action predictions is observed.
We now turn to quenches from interacting initial states (γ i = 0) to different interacting final states for which no results have been obtained by means of exact BA nor simplified quench action method so far. In Fig. 3 we show the asymptotic equilibrium values obtained with our t-VMC approach for quantum quenches from γ i = 1 (Left panel), and γ i = 4 (Right panel) to several values of γ f . Since, by the variational theorem, the ground state of H i gives an upper bound for the ground state energy of H f , the system is pushed into a linear combination of excited states of the final hamiltonian. For systems able to thermalize to the Boltzmann ensemble (BE), relaxation to a stationary state described by the density matrix ρ T = e −H f /T , at an effective temperature T , would occur. Comparing the stationary value,ḡ 2 (0), of our t-VMC calculations at long times to the thermal values of the pair correlation functions, g T 2 (0), a necessary condition for simple Boltzmann thermalization is given The effective temperature, T , is determined by imposing the energy expectation value of the final Hamiltonian, H f , in the ground state, Φ 0 (γ i ) of the initial Hamiltonian, Here, the thermal expectation value, H f T at the equilibrium temperature T is computed from the Yang-Yang BA equations [60]. The quantity H f T then depends on a single parameter T , that is fitted to match the value of As shown in Fig.2-(b), Boltzmann thermalization certainly does not occur in the case for the Lieb-Liniger model when quenching from a non-interacting state, γ i = 0, where we findḡ 2 = g T 2 . This can be understood in terms of the existence of dynamically conserved charges (beyond energy and density conservation) which can yield an equilibrium value substantially different from the BE prediction. In particular, it is widely believe that the Generalized Gibbs Ensemble (GGE) is the correct thermal distribution approached after the quench [61,62]. Several constructive approaches for the GGE have been put forward in past years [8,9,58], and the quench action predictions reported in Fig.2-(b) converge to the GGE predictions for the thermal values. In Fig.2-(b) we also show the thermal GGE values g GGE 2 (rightmost dashed lines), and notice that our results are much closer to the GGE predictions than the simple BE. Deviations from the asymptotic GGE results are observed at large γ f , a regime in which the accuracy of our approach is still sufficient to resolve the difference between the BE and the long-term equilibration value.
For correlated initial ground states, γ i = 0, GGE predictions are fundamentally harder to obtain than for the non-interacting initial states, and the BE is the only reference thermal distribution we can compare with at this stage. From our results we observe that the difference betweenḡ 2 and the simple BE prediction, g T 2 , is quantitatively reduced, see Fig. 3. In particular, for γ i = 4, the stationary valuesḡ 2 are quantitatively close to the ones predicted by the Boltzmann thermal distribution at the effective temperature T . Even though this quantitative agreement is likely to be coincidental, the regimes of parameters quenches studied here provide guidance for future experimental studies. In particular, it will be of great interest to understand whether a cross-over from a strongly non-Boltzmann to a close-to-Boltzmann thermal behavior might occur as a function of the initial interaction strength also for other local observables.

IV. CONCLUSIONS
In this paper we have introduced a novel approach to the dynamics of strongly-correlated quantum systems in continuous space. Our method is based on correlated many-body wave-function systematically expanded in terms of reduced m-body Jastrow functions. The unitary dynamics in the subspace of these correlated states was realized using time-dependent variational Monte Carlo. We have demonstrated the possibility or performing calculations up to the three-body level, m ≤ 3, for the Lieb-Liniger model, for both static and dynamical properties. The improvement from m = 2 to m = 3 provides an internal criterium to judge the validity of our results whenever exact results are unavailable. Benchmarking t-VMC with exact or numerical approaches whenever available, we have found a very good agreement with existing results. For static properties, our approach is at the level of state-of-the-art MPS techniques in lattice systems and of latest c-MPS results for interacting gases. For dynamical properties, we have investigated for the first time general interaction quenches which are at the moment unaccessible to Bethe-ansatz approaches. Since the general structure of our t-VMC method does not depend on the dimensionality of the system, it can be directly applied to bosonic systems in higher dimensions with a polynomial increase in computational cost. The methods presented here therefore pave the way to accurate out-of-equilibrium dynamics of two-and three-dimensional quantum gases and fluids beyond mean field approximations.

ACKNOWLEDGMENTS
We acknowledge discussions with J. De Nardis, M. Dolfi, M. Fagotti, T. Osborne, and M. Troyer. We thank M. Ganahl for providing us the c-MPS results in Fig. 1, J. Zill for the Bethe ansatz results in Fig.  2 (a), and J. De Nardis for the quench actions re-sults in Fig. 2 (b). loc (X, t), however, may contain effective interaction terms involving a number of bodies larger than M , which leads to a systematic error in the truncation. However, the structure of these additional terms stemming from the local energy can be systematically used to deduce the functional structure of the higher order terms. For example, the one-body truncated local energy reads, for one-dimensional particles, and contains a 2-body term which cannot be accounted for exactly by u 1 . Introduction of a symmetric two-body Jastrow factor u 2 (x i , x j ; t), then leads to In the latter expression, one recognizes an effective twobody term which can be accounted for by u 2 and an additional three-body term in the form of product of twobody functions. The functional form of the three-body Jastrow can be therefore deduced from this additional term and formed accordingly: with two-body functionsū 3 (x i , x j ; t) containing new variational parameters to be determined. Upon pursuing this approach, the expansion can be systematically pushed to higher orders and the functional structure of the higher order functions inferred. The same constructive approach we have discussed here is also valid for the Schrödinger equation in imaginary-time ∂ τ U (X; τ ) = −E loc (x; τ ), and has been successfully used to infer the functional structure for ground-state properties [63].

Appendix B: Monte Carlo Sampling
In order to solve the t-VMC equations of motion, Eq. (6), expectation values of some given operator O need to be computed over the many-body wave-function Φ(X, t). This is achieved by means of Monte Carlo sampling of the probability distribution Π(X) = |Φ(X)| 2 (In the following we omit explicit reference to the time t, assuming that all expectation values are taken over the wave-function at a given fixed time). An efficient way of sampling the given probability distribution is to devise a Markov chain of configurations X(1), X(2), . . . X(N c − 1), X(N c ) which are distributed according to Π(X). Quantum expectation values of a given operator O can then be obtained as statistical expectation values over the Markov chain as where O loc (X) = X|O|Φ X| Φ , and the equivalence is achieved in the limit N c → ∞.
The Markov chain is realized by the Metropolis-Hastings algorithm. Given the current state of the Markov chain, X(i), a configuration X is generated according to a given transition probability T (X(i) → X ). The proposed configuration is then accepted (i.e. X(i + 1) = X ) with probability With this discretization, a two-body Jastrow factor reads where u 2 (a, b; t) is a time-dependent matrix of size N s × N s which, in 1D and in the presence of translational symmetry, depends only on dist(a − b), i.e. it has N s /2 variational parameters.
Appendix E: Benchmark Study for N = 3 on a Lattice Here, we use exact diagonalization of a Hamiltonian within a given finite basis for a quantitative test of our method. Exact diagonalization is limited to very small systems on a finite basis, and we haven chosen a system containing N = 3 particles on L/a = 40 lattice sites as a simple, but highly non-trivial reference. In contrast to our comparison with BA methods, all observables can be accessed by exact diagonalization and we have used the off-diagonal one body density matrix g 1 (r, t) and the pair correlation function g 2 (r, t) at different distances r = |x 1 − x 2 | of two particles after time t where the system is quenched from the non-interacting initial state, γ i = 0, to a final interaction γ f > 0, to provide a benchmark on a more general observable.
We first benchmark the influence of the time-step lattice size discretization ∆t error on g 2 . From Fig.4(a) we see that the t-VMC dynamics is stable over a long time and the time step error can be brought to convergence. Further, we see that for final interaction γ f = 4 the truncation at the two-body level, U (2) , introduces only a small systematic error, mainly a dephasing effect, which is almost negligible at the scale of the figure. Due to the stochastic noise of the Monte Carlo integration, t-VMC introduces additional high frequency oscillations which are, however, well separable from the deterministic propagation. The amplitude of these high frequency oscillations also quantifies the purely statistical error of our data.
Whereas exact diagonalization is limited to rather small basis sets, we can access much a larger basis within t-VMC. In Fig.4 (b) we show results within the U (2) approximation with L/a = 80 and L/a = 160 with time discretization ∆t ∼ (L/a) 2 . We see that the basis set truncation in general introduces a dephasing at large enough time.
The systematic error of U (2) increases towards quenches to stronger interaction strength and becomes more visible for γ f = 8 shown in Fig.5(a). However, even in this case, the most important effect remains to be a simple dephasing, a small shift of averaged quantities is probable, but difficult to quantify precisely. Introducing a general three-body Jastrow fields, U (3) , described in detail in Appendix F, the systematic error for N = 3 can be fully eliminated.
In figure 5 (b), we also benchmark the possibility of calculating the off-diagonal one-body density matrix after a quench. Here the systematic error of U (2) is more pronounced at smaller γ f in the long range and time regime.
From our study of the three particle problem, we conclude that truncation of the many-body wave function at the level of U 2 may provide an excellent approximation for g 2 (r, t) for quenches involving not too strong interaction strengths, γ 5. The systematic error due to the U (2) truncation is mainly a dephasing at large times involving small relative errors of time averaged quantities. Similar systematic dephasing errors will occur for too large time discretization or basis set truncation.
Since our method provides a parametrization of the full wave function for a given time, many different observables can be evaluated via usual Monte Carlo methods. However, the quality of different observables may vary and depend more sensitive on the inclusion of higher order correlations U (n) with n > 2 as in the case of the single body density matrix. Although these higher order terms are computationally expensive, the scaling is not exponential, and we have explicitly shown that calculations with n = 3 are feasible. We notice that the computational complexity may be further reduced by functional forms adapted to the problem [26]. For a general time dependent wave function, we have to go beyond the usual ground state structure of the threebody Jastrow given in Eq. (A3). Here, we provide details of our three-body term in a general form beyond the present application in one dimension.
Introducing M basis functions, b a (r), a = 1, . . . M , we can introduce many-body vectors [26], B a iα = j r α ij b a (r ij ), where α = 1, . . . D indicates the summation over directions and i = 1, . . . N . The variational parameters of a general three-body structure can then be written in terms of a matrix w ab , such that   Figure 5. (a) Time-dependent expectation value of the two-body correlations, g2(r, t), after a quantum quench from a noninteracting state, γi = 0, to γ f = 8 at three different distance |x1 − x2| = 0, L/10, L/4. Here, the system is on a lattice with L/a = 40 lattice sites, the full line is obtained by exact diagonalization of the Hamiltonian, the other curves are from t-VMC truncated at the level of U (2) or U (3) . In contrast to γ f = 4 shown in Fig.4, systematic differences of U (2) compared to the exact results are more visible here, the exact dynamics is recovered by inclusion of three-body terms, U (3) , into the tVMC wave function. (b) shows the off-diagonal single particle density matrix, g1(r, t), at three different distances, r = L/10, L/4 and L/2.