Unitary Dynamics of Strongly Interacting Bose Gases with the Time-Dependent Variational Monte Carlo Method in Continuous Space

Unitary Dynamics of Strongly Interacting Bose Gases with the Time-Dependent Variational... PHYSICAL REVIEW X 7, 031026 (2017) Unitary Dynamics of Strongly Interacting Bose Gases with the Time-Dependent Variational Monte Carlo Method in Continuous Space Giuseppe Carleo Institute for Theoretical Physics, ETH Zurich, Wolfgang-Pauli-Strasse 27, 8093 Zurich, Switzerland Lorenzo Cevolani Laboratoire Charles Fabry, Institut d’Optique, CNRS, Univ Paris Sud 11, 2 avenue Augustin Fresnel, F-91127 Palaiseau cedex, France Laurent Sanchez-Palencia Laboratoire Charles Fabry, Institut d’Optique, CNRS, Univ Paris Sud 11, 2 avenue Augustin Fresnel, F-91127 Palaiseau cedex, France and Centre de Physique Théorique, Ecole Polytechnique, CNRS, Univ Paris-Saclay, F-91128 Palaiseau, France Markus Holzmann LPMMC, UMR 5493 of CNRS, Université Grenoble Alpes, F-38042 Grenoble, France and Institut Laue-Langevin, BP 156, F-38042 Grenoble Cedex 9, France (Received 19 December 2016; revised manuscript received 6 April 2017; published 8 August 2017) We introduce the time-dependent variational Monte Carlo method for continuous-space Bose gases. Our approach is based on the systematic expansion of the many-body wave function in terms of multibody correlations and is essentially exact up to adaptive truncation. The method is benchmarked by comparison to an 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 are also compared to quench action results available for noninteracting 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. DOI: 10.1103/PhysRevX.7.031026 Subject Areas: Atomic and Molecular Physics, Computational Physics, Quantum Physics I. INTRODUCTION more involved. Quantum Monte Carlo algorithms, the de facto tool for simulating quantum many-body systems at The study of equilibration and thermalization properties thermal equilibrium [2–4], cannot be directly used to study of complex many-body systems is of fundamental interest time-dependent unitary dynamics. Out-of-equilibrium for many areas of physics and natural sciences [1].For properties are then often treated on the basis of approx- systems governed by classical physics, an exact solution of imations drastically simplifying the microscopic physics. Newton’s equations of motion is often numerically feasible, Irreversibility is either enforced with an explicit breaking of using, for instance, molecular-dynamics simulations. For unitarity, e.g., within the quantum Boltzmann approach, or quantum systems, the mathematical structure of the time- the dynamics is reduced to mean-field description using dependent Schrödinger equation is instead fundamentally 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 Published by the American Physical Society under the terms of entanglement are ignored. the Creative Commons Attribution 4.0 International license. For specific systems, exact dynamical results can be Further distribution of this work must maintain attribution to derived. This is the case for integrable 1D models, for the author(s) and the published article’s title, journal citation, and DOI. which Bethe ansatz (BA) solutions exist [7]. However, also 2160-3308=17=7(3)=031026(12) 031026-1 Published by the American Physical Society GIUSEPPE CARLEO et al. PHYS. REV. X 7, 031026 (2017) in this case many open questions still persist. For example, further compare to the quench action approach for large the exact evaluation of correlation functions for out-of- systems approaching the thermodynamic limit. Finally, we equilibrium dynamics is at present an unsolved problem. apply our method to the study of general quenches from As a result, despite important theoretical and experimental arbitrary initial states, for which no exact results in the progress [8–15], a complete picture of thermalization (or its thermodynamic limit are currently available. absence), e.g., based on general quench protocols [16],is still missing. II. METHOD Numerical methods for strongly interacting systems face A. Expansion of the many-body wave function important challenges as well. Numerical renormalization group and density-matrix renormalization group (DMRG) Consider a nonrelativistic quantum system of N identical approaches provide an essentially exact description of bosons in d dimensions, and governed by the first- arbitrary 1D lattice systems in and out of equilibrium quantization Hamiltonian [17–20], but they have less predictive power when applied N N X X X to continuous-space systems. On one hand, multiscale 1 1 H ¼ − ∇ þ v ðx⃗ Þþ v ðx⃗ ; x⃗ Þ; ð1Þ 1 i 2 i j extensions of the DMRG optimization scheme to the limit 2 2 i¼1 i¼1 i≠j of continuous-space lattices [21] are, to date, limited to relatively small system sizes [22]. On the other hand, where v ðx⃗ Þ and v ðx;⃗ y⃗ Þ are, respectively, a one-body 1 2 efficient ground-state optimization schemes for continuous external potential and a pairwise interparticle interaction quantum field matrix product states (cMPS) [23] have been [34]. Without loss of generality, a time-dependent N-body introduced only very recently [24], and applications to state can be written as ΦðX;tÞ¼ exp ½UðX; tÞ, where quantum dynamics are still to be realized. A further X ¼ x⃗ ; x⃗ ; …; x⃗ is the ensemble of particle positions 1 2 N formidable challenge is the efficient extension of these and U is a complex-valued function of the N-particle approaches to higher dimensions, which is a fundamentally N×d coordinates, R → C. Since the Hamiltonian Eq. (1) hard problem. contains only two-body interactions, it is expected that an Another class of methods for strongly interacting systems expansion of U in terms of few-body Jastrow functions is based on variational Monte Carlo (VMC) methods, containing at most m-body terms rapidly converges combining highly entangled variational states with robust towards the exact solution. Truncating this expansion up stochastic optimization schemes [25]. Such approaches have to a certain order, M ≤ N, leads to the Bijl-Dingle-Jastrow- been successfully applied to the description of continuous Feenberg expansion [35–38], quantum systems, in any dimension and not only in 1D X X [26,27]. More recently, out-of-equilibrium dynamics has ðMÞ U ðX;tÞ¼ u ðx⃗ ;tÞþ u ðx⃗ ;x⃗ ;tÞ 1 i 2 i j become accessible with the extension of these methods to 2! i¼1 i≠j real-time unitary dynamics, within the time-dependent varia- tional Monte Carlo (tVMC) method [28,29].So far,the þ  u ðx⃗ ;x⃗ ;…;x⃗ ;tÞ; ð2Þ M i i i 1 2 M M! tVMC approach has been developed for lattice systems with i ≠i ≠…i 1 2 M bosonic [28,29], spin [30–32], and fermionic [33] statistics, yielding a description of dynamical properties with an where u ðr; tÞ are functions of m particle coordinates, accuracy often comparable with MPS-based approaches. r ¼ x⃗ ; x⃗ ; …; x⃗ , and of the time t. A global constraint on i i i 1 2 m In this paper, we extend the tVMC approach to access the function UðX;tÞ is given by particle statistics. In dynamical properties of interacting quantum systems in the bosonic case, we demand that Uðx⃗ ; x⃗ ; …; x⃗ Þ¼ 1 2 N continuous space. Our approach is based on a systematic Uðx⃗ ; x⃗ ; …; x⃗ Þ, for all particle permutations σ. σð1Þ σð2Þ σðNÞ expansion of the wave function in terms of few- In general, the functions u ðr; tÞ can have an arbitrarily body Jastrow correlation functions. Using the 1D Lieb- complex dependence on the m particle coordinates, which Liniger model as a test case, we first show that the inclusion can prove problematic for practical applications. of high-order correlations allows us to systematically Nonetheless, a simplified functional dependence can often approach the exact BA ground-state energy. Our results be imposed, resulting from the two-body character of the improve by orders of magnitude on previously published interactions in the original Hamiltonian. For m ≥ 3, VMC and cMPS results, and are in line with the latest state- u ðr; tÞ can be conveniently factorized in terms of general of-the-art developments in the field. We further compute two-particle vector and tensor functions following single-body and pair correlation functions, hardly acces- Ref. [26]. Details of this approach and the present imple- sible by current BA methods. We then calculate the time mentation are presented in Appendixes A and F. evolution of the contact pair correlation function following An appealing property of the many-body expansion a quench in the interaction strength. For the noninteracting Eq. (2) is that it is able to describe intrinsically nonlocal initial state, we benchmark our results to exact BA correlations in space. For instance, the two-body ⃗ ⃗ calculations available for a small number of bosons and u ðx ; x ; tÞ, as well as any m-body function u ðr; tÞ, 2 i j m 031026-2 UNITARY DYNAMICS OF STRONGLY-INTERACTING BOSE … PHYS. REV. X 7, 031026 (2017) can be long range in the particle separation jx⃗ − x⃗ j. This In practice, these equations are numerically solved for the i j time derivatives ∂ u ðr ; tÞ at each time step t. The expect- nonlocal spatial structure allows for a correct description of t p ation values taken over the time-dependent state h·i , which gapless phases, where a two-body expansion may already enter Eq. (6), are found via a stochastic sampling of the capture all the universal features, in the sense of the renormalization group approach [39]. This is in contrast probability distribution ΠðX;tÞ¼jΦðX;tÞj . This is effi- with the MPS decomposition of the wave function, which is ciently achieved by means of the Metropolis-Hastings intrinsically local in space. 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 B. Time-dependent variational Monte Carlo method state Eq. (2) after time integration. The time evolution of the variational state Eq. (2) is The tVMC approach as formulated here provides, in entirely determined by the time dependence of the Jastrow principle, an exact description of the real-time dynamics of functions u ðr; tÞ. In order to establish optimal equations the N-body system. The essential approximation lies, of motion for the variational parameters, we note that the however, in the truncation of the Bijl-Dingle-Jastrow- functional derivative of UðX; tÞ with respect to the varia- Feenberg expansion to the M most relevant terms. In tional, complex-valued, Jastrow functions u ðr; tÞ, practical applications, the M ¼ 2 or M ¼ 3 truncation is often sufficient. Systematic improvement beyond M ¼ 3 is δUðX; tÞ ≡ ρ ðrÞ; ð3Þ possible [26], but may require substantial computational δu ðr; tÞ effort. yields the m-body density operators III. LIEB-LINIGER MODEL X Y As a first application of the continuous-space tVMC ðrÞ¼ δðx − r Þ: ð4Þ m i j m! approach, we consider the Lieb-Liniger model [43]. On one i ≠i ≠i j 1 2 m hand, some exact results and numerical data are available, allowing us to benchmark the Jastrow expansion and the The expectation values of the operators ρ over the tVMC approach. On the other hand, several aspects of the state jΦðtÞi give the instantaneous m-body correlations. out-of-equilibrium dynamics of this model are unknown, For instance, hρ ðr ;r Þi ¼ hδðx − r Þδðx − r Þi , 2 1 2 t i<j i 1 j 2 t which we compute here for the first time using the tVMC where h  i ¼hΦðtÞj…jΦðtÞi=hΦðtÞjΦðtÞi, is propor- approach. tional to the two-point density-density correlation function. The Lieb-Liniger model describes N interacting bosons We can then express the time derivative of the truncated ðMÞ in one dimension with contact interactions. It corresponds variational state U ðX; tÞ using the functional deriva- to the Hamiltonian Eq. (1) with v ðxÞ¼ 0 and v ðx; yÞ¼ 1 2 tives, Eq. (3), as a sum of the few-body density operators up gδðx − yÞ, where g is the coupling constant. Here, we to the truncation M; i.e., consider periodic boundary conditions over a ring of length L and the particle density n ¼ N=L. The density depend- ðMÞ ence is as usual expressed in terms of the dimensionless ∂ U ðX;tÞ¼ drρ ðrÞ∂ u ðr; tÞ: ð5Þ t m t m m¼1 parameter γ ¼ mg=ℏ n. The Lieb-Liniger model is the prototypal model of continuous one-dimensional strongly The exact wave function satisfies the Schrödinger correlated gas exactly solvable by the Bethe ansatz [43]. equation i∂ UðX; tÞ¼ E ðX;tÞ, where E ðX; tÞ¼ t loc loc This model is experimentally realized in ultracold atomic ðhXjHjΦðtÞi=hXjΦðtÞiÞ is the so-called local energy. gases strongly confined in one-dimensional optical traps, The optimal time evolution of the truncated Bijl-Dingle- and several studies on out-of-equilibrium physics have Jastrow-Feenberg expansion Eq. (2) can be derived already been realized [14,44–48]. imposing the Dirac-Frenkel time-dependent variational As a result of translation invariance, we have principle [40,41]. In geometrical terms, this amounts u ðx; tÞ¼ const, and the first nontrivial term in the to minimizing the Hilbert-space norm of the residuals many-body expansion Eq. (2) is the two-body, translation ðMÞ ðMÞ ðMÞ invariant, function u ðjx − yj; tÞ. To compute the functional R ðX; tÞ ≡ ∥i∂ U ðX; tÞ − E ðX; tÞ∥, thus yielding t 2 loc derivatives of the many-body wave function, we proceed a variational many-body state as close as possible to the with a projection of the continuous Jastrow fields u ðr; tÞ exact one [42]. The minimization can be performed onto a finite-basis set. Here, we find it convenient to explicitly and yields a closed set of integro-differential represent both the field Hamiltonian Eq. (1) and the Jastrow equation for the Jastrow functions u ðr; tÞ: functions on a uniform mesh of spacing a, leading to n ¼ var ½L=ð2aÞ variational parameters for the two-body Jastrow δhρ ðrÞi δhHi m t t 0 0 dr ∂ u ðr ; tÞ¼ −i : ð6Þ t p term. In the following, our results are extrapolated to the δu ðr ; tÞ δu ðr; tÞ p m p¼1 continuous limit, corresponding to a → 0. The finite-basis 031026-3 GIUSEPPE CARLEO et al. PHYS. REV. X 7, 031026 (2017) projection as well as the numerical time integration of Eq. (6) the modulus of a Vandermonde determinant of plane are detailed in Appendix C and benchmarked against exact waves, corresponding to the two-body Jastrow function diagonalization results for N ¼ 3 in Appendix E. u ðrÞ¼ log ðsin rπ=LÞ [51]. To assess the overall quality of pair wave functions for ground-state properties, we start comparing the variational ground-state energies E obtained A. Ground-state properties for M ¼ 2 with the exact BA result. In Fig. 1(a), we show To assess the quality of truncated Bijl-Dingle-Jastrow- the relative error ΔE=E as a function of the interaction Feenberg expansions for the Lieb-Liniger model, we start −4 strength γ. We find that the relative error is lower than 10 our analysis considering ground-state properties. An exact for all values of γ and the accuracy of our two-body Jastrow solution can be found from the Bethe ansatz and gives function is superior to previously published variational access to exact ground-state energies and local properties results based on either cMPS [23] or VMC [52] methods [43]. Other nonlocal properties are substantially more [see Fig. 1(a) for a quantitative comparison]. Notice that the difficult to extract from the BA solution, and unbiased improvement with respect to previous VMC results is due results for ground-state correlation functions have not been to the larger variational freedom of our u ðrÞ function, reported so far. In order to determine the best possible which is not restricted to any specific functional form as variational description of the ground state within our many- done in Ref. [52]. body expansion, different strategies are possible. A first Even though the accuracy reached by the two-body possibility is to consider the imaginary-time evolution −τH Jastrow function may already be sufficient for most jΨðτÞi ¼ e jΨ i, which systematically converges to practical purposes for all values of γ, we also consider the exact ground state in the limit τ ≫ Δ , where Δ ¼ 1 1 higher-order terms with M ¼ 3; see Fig. (2)(a). The E − E is the gap with the first excited state on a finite 1 0 introduction of the third-order term yields a sizable system and provided that the trial state Ψ is nonorthogonal improvement such that the maximum error is about 3 to the exact ground state. Imaginary-time evolution in the orders of magnitude smaller than original cMPS results variational subspace can be implemented considering the [23], and features a similar accuracy of recently reported formal substitution t → −iτ in the tVMC equations cMPS results [24]. Overall, our approach reaches a pre- [Eq. (6)]. The resulting equations are equivalent to the cision on a continuous-space system, which is comparable stochastic reconfiguration approach [49]. However, direct to state-of-the-art MPS or DMRG results for gapless minimization of the variational energy can be significantly systems on a lattice [53]. more efficient, in particular, for systems becoming gapless Finally, to further assess the quality of our ground-state in the thermodynamic limit, where Δ ∼ polyð1=NÞ.Given ansatz beyond the total energy, we also study nonlocal the gapless nature of the Lieb-Liniger model, we find it properties of the ground-state wave function, which are computationally more efficient to adopt a Newton method not accessible by existing exact BA methods. In Figs. 1(b) to minimize the energy variance [50]. and 1(c), we show, respectively, our results for the For the ground state, the many-body expansion truncated off-diagonal part of the one-body density matrix, at M ¼ 2 is exact not only in the noninteracting limit γ ¼ 0 g ðrÞ ∝ hΨ ðrÞΨð0Þi, and for the pair correlation function, but also in the Tonks-Girardeau limit γ → ∞. In this † † fermionized limit, the wave function can be written as g ðrÞ ∝ hΨ ðrÞΨðrÞΨ ð0ÞΨð0Þi, where ΨðrÞ is the bosonic (a) (b) (c) FIG. 1. Ground-state properties of the Lieb-Liniger model as obtained from different variational approaches: (a) relative accuracy of ð2Þ ð3Þ the ground-state energy, (b) one-body density matrix, (c) density-density pair correlation functions. U and U denote results for the ð2Þ two- and three-body expansion (present work), U is the parametrized two-body Jastrow state of Ref. [52], cMPS results are from AG Ref. [23], and cMPS are those very recently reported in Ref. [24]. Distances r are in units of the inverse density 1=n. Our variational results are obtained for N ¼ 100 particles. Finite-size corrections on local quantities are negligible, and very mildly affect the reported large-distance correlations. Overall statistical errors are of the order of symbol sizes for (a) and line widths for (b) and (c). 031026-4 UNITARY DYNAMICS OF STRONGLY-INTERACTING BOSE … PHYS. REV. X 7, 031026 (2017) field operator. We find an overall excellent agreement with validation of our method accessing g ðr; tÞ at nonvanishing the results that have been obtained with cMPS in Ref. [24], distances. The comparison shown in Fig. 2(a) shows an except for some small deviations at large values of r, which overall good agreement. The tVMC and BA results are we attribute to residual finite-size effects in our approach. We indistinguishable for weak interactions (γ ¼ 1) at the scale find that the addition of the three-body terms does not of the figure; similarly, for this interaction strength, tVMC ð2Þ ð3Þ significantly change correlation functions. Already at the calculations based on U and U agree with each other two-body level, the present results are statistically indistin- up to possible dephasing effects not visible within our noise guishable from exact results obtained using our implemen- level. For larger interactions, we notice systematic but small tation [54] of the worm algorithm [55] and for the same differences between BA and tVMC with M ¼ 2 or M ¼ 3 system (not shown). results. These differences amount to a small increase in the amplitude of the oscillations. This effect tends to increase B. Quench dynamics with the interaction strength, being hardly visible for γ ¼ 1 and more pronounced for γ ¼ 10. However, these oscil- Having assessed the quality of the ansatz for local and lations result at large times from the discrete mode structure nonlocal ground-state properties, we now turn to the study due to the very small number of particles. They vanish in of the out-of-equilibrium properties of the Lieb-Liniger the physical thermodynamic limit. In turn, the comparison model. We focus on the description of the unitary dynamics at small particle numbers indicates an accuracy better induced by a global quantum quench of the interaction than a few percent for time-averaged quantities in the strength, from an initial value γ to a final value γ . Exact i f asymptotic large time limit up to γ ¼ 10, with results at BA results are available only in the case of a noninteracting M ¼ 3 systematically improving on the M ¼ 2 case. initial state (γ ¼ 0). Even in this case, the dynamical BA Concerning small and intermediate time scales, we do equations can be exactly solved only for a modest number not observe systematic deviations between the tVMC of particles with further truncation in the number of energy results and the BA solution. In particular, the relaxation eigenmodes [56,57], N ≲ 10, since the complexity of the times are remarkably well captured by the tVMC approach. BA solution increases exponentially with the number of On the basis of this comparison and of the comparison for particles. Simplifications in the thermodynamic limit are three particles with exact diagonalization results presented exploited by the quench action [58], and have recently been in Appendix E, we conclude that tVMC approach allows applied to quantum quenches starting from a noninteracting accurate quantitative studies of both the relaxation and initial state [57]. In the following, we first compare our equilibration dynamics. This careful benchmarking now tVMC results to these existing results, and then present new allows us to confidently apply the tVMC approach regimes results for quenches following a nonvanishing initial that are inaccessible to exact BA, namely large but finite N interaction strength. and long times, as well as the case of nonvanishing initial To assess the quality of the time-dependent wave interactions. function, we compare our results for the evolution of local Let us consider relaxation of density correlations for a density-density correlations g ð0;tÞ with the truncated BA large number of particles, close to the thermodynamic limit results obtained in Refs. [56,59] for a small number of (here, we use N ¼ 100). As we show in Fig. 2(b), we notice particles, N ≃ 6. Appendix E also provides further (a) (b) FIG. 2. Time-dependent expectation value of local two-body correlations after a quantum quench from a noninteracting state, γ ¼ 0, to γ . (a) tVMC results are compared with BA results obtained for a small number of particles [56,59]. The correlation function is rescaled to have g ð0; 0Þ¼ 1. (b) tVMC results for N ¼ 100 particles and γ ¼ 1, 2, 4, 8 (from top to bottom) compared to the quench 2 f action (QA) predictions from Ref. [57] (dashed lines), to the Boltzmann thermal averages at the effective temperature T (dot-dashed lines), and to the GGE thermal averages prediction (rightmost dashed lines). Statistical error bars on tVMC data are of the order of the line widths. 031026-5 GIUSEPPE CARLEO et al. PHYS. REV. X 7, 031026 (2017) As we show in Fig. 2(b), Boltzmann thermalization certainly does not occur in the case for the Lieb-Liniger model when quenching from a noninteracting state, γ ¼ 0, where we find g¯ ≠ g . 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 believed that the generalized Gibbs ensemble (GGE) is the correct thermal (a) (b) distribution approached after the quench [61,62]. Several constructive approaches for the GGE have been put forward in past years [8,9,57], and the quench action predictions reported in Fig. 2(b) converge to the GGE FIG. 3. Time-dependent expectation value of local two-body predictions for the thermal values. In Fig. 2(b), we also correlations after a quantum quench from the interacting ground GGE show the thermal GGE values g (rightmost dashed state at γ ¼ 1 (a) and γ ¼ 4 (b). Long-term dynamical averages i i lines), and note that our results are much closer to the GGE (red continuous lines) are compared to thermal averages at the predictions than the simple BE. Deviations from the effective temperature set by energy conservation (black dashed asymptotic GGE results are observed at large γ , a regime lines). Uncertainties on the thermal averages are of the order of the line width, and are larger for small values of γ . 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 γ ≠ 0, GGE predic- that the amplitudes of the large-time oscillations, attributed tions are fundamentally harder to obtain than for the to the discrete mode spectrum, are now drastically sup- noninteracting initial states, and the BE is the only pressed compared to the quenches with N ¼ 6. After an reference thermal distribution we can compare with at this initial relaxation phase, the quantity g ð0;tÞ approaches a stage. From our results we observe that the difference stationary value. Comparing our curves with those obtained between g¯ and the simple BE prediction g is quantita- with the quench action method, we find a qualitatively good 2 2 tively reduced; see Fig. 3. In particular, for γ ¼ 4, the agreement, albeit a general tendency to underestimate the stationary values g¯ are quantitatively close to the ones quench action predictions is observed. 2 predicted by the Boltzmann thermal distribution at the We now turn to quenches from interacting initial states effective temperature T . Even though this quantitative (γ ≠ 0) to different interacting final states for which no agreement is likely to be coincidental, the regimes of results have been obtained by means of exact BA nor the parameter quenches we study here provide guidance for simplified quench action method so far. In Fig. 3, we show future experimental studies. In particular, it will be of great the asymptotic equilibrium values obtained with our tVMC interest to understand whether a crossover from a strongly approach for quantum quenches from γ ¼ 1 (left-hand non-Boltzmann to a close-to-Boltzmann thermal behavior panel) and γ ¼ 4 (right-hand panel) to several values of might occur as a function of the initial interaction strength γ . Since, by the variational theorem, the ground state of H f i also for other local observables. gives an upper bound for the ground-state energy of H ,the system is pushed into a linear combination of excited states of the final Hamiltonian. For systems able to thermalize to the IV. CONCLUSIONS Boltzmann ensemble (BE), relaxation to a stationary state ⋆ In this paper, we introduce a novel approach to the −H =T described by the density matrix ρ ⋆ ¼ e ,atan dynamics of strongly correlated quantum systems in effective temperature T , would occur. Comparing the sta- continuous space. Our method is based on a correlated tionary value g¯ ð0Þ of our tVMC calculations at long times to many-body wave function systematically expanded in the thermal values of the pair correlation functions g ð0Þ,a terms of reduced m-body Jastrow functions. The unitary necessary condition for simple Boltzmann thermalization is dynamics in the subspace of these correlated states is T ⋆ given by g¯ ¼ g . The effective temperature T is deter- realized using the time-dependent variational Monte Carlo mined by imposing the energy expectation value of the final method. We demonstrate the possibility of performing Hamiltonian H in the ground state Φ ðγ Þ of the initial f 0 i calculations up to the three-body level, m ≤ 3, for the Hamiltonian, hH i ⋆ ¼hΦ ðγ ÞjH jΦ ðγ Þi. Here, the ther- f 0 i f 0 i T Lieb-Liniger model, for both static and dynamical proper- mal expectationvalue, hH i ⋆ at the equilibrium temperature f ties. The improvement from m ¼ 2 to m ¼ 3 provides an T is computed from the Yang-Yang BA equations [60].The internal criterion to judge the validity of our results quantity hH i ⋆ then depends on a single parameter T , whenever exact results are unavailable. Benchmarking f T which is fitted to match the value of hΦ ðγ ÞjH jΦ ðγ Þi. the tVMC approach with exact or numerical approaches 0 i f 0 i 031026-6 UNITARY DYNAMICS OF STRONGLY-INTERACTING BOSE … PHYS. REV. X 7, 031026 (2017) whenever available, we find a very good agreement with and contains a two-body term that cannot be accounted for existing results. For static properties, our approach is at the exactly by u . Introduction of a symmetric two-body level of state-of-the-art MPS techniques in lattice systems Jastrow factor u ðx ;x ; tÞ then leads to 2 i j and of the latest cMPS results for interacting gases. For dynamical properties, we investigate, for the first time, ð2Þ ð1Þ E ðX; tÞ¼ E ðX; tÞþ loc loc general interaction quenches, which are at the moment unaccessible to Bethe ansatz approaches. Since the general − ½∂ u ðx ; tÞ∂ u ðx ;x ; tÞþ x 1 i x 2 i j i i structure of our t-VMC method does not depend on the i≠j dimensionality of the system, it can be directly applied to bosonic systems in higher dimensions with a polynomial − ∂ u ðx ;x ; tÞþ x 2 i j i≠j increase in computational cost. The methods we present XX here therefore pave the way to accurate out-of-equilibrium − ∂ u ðx ;x ; tÞ∂ u ðx ;x ; tÞ: x 2 i j x 2 i k dynamics of two- and three-dimensional quantum gases i i i≠j k≠i and fluids beyond mean-field approximations. ðA2Þ ACKNOWLEDGMENTS In the latter expression, one recognizes an effective two- We acknowledge discussions with M. Dolfi, J. De body term which can be accounted for by u and an Nardis, M. Fagotti, T. Osborne, and M. Troyer. We thank additional three-body term in the form of a product of two- M. Ganahl for providing the cMPS results in Fig. 1, J. Zill body functions. The functional form of the three-body for the Bethe ansatz results in Fig. 2(a), and J. De Nardis for Jastrow term can therefore be deduced from this additional the quench actions results in Fig. 2(b). This research was term and formed accordingly: supported by the Marie Curie IEF Program (FP7/2007– 2013 Grant Agreement No. 327143), the European u ðx ;x ;x ; tÞ¼ u¯ ðx ;x ; tÞu¯ ðx ;x ; tÞ; ðA3Þ 3 i j k 3 i j 3 j k Research Council Starting Grant “ALoGlaDis” (FP7/ 2007-2013 Grant Agreement No. 256294) and Advanced with two-body functions u¯ ðx ;x ; tÞ containing new varia- 3 i j Grant “SIMCOFE” (FP7/2007-2013 Grant Agreement tional parameters to be determined. Upon pursuing this No. 290464), the European Commission FET-Proactive approach, the expansion can be systematically pushed to QUIC (H2020 Grant No. 641122), the French ANR-16- higher orders and the functional structure of the higher- CE30-0023-03 (THERMOLOC), and the Swiss National order functions inferred. The same constructive approach Science Foundation through NCCR QSIT. It was per- we discuss here is also valid for the Schrödinger equation in formed using HPC resources from GENCI-CCRT/ imaginary time, ∂ UðX; τÞ¼ −E ðx; τÞ, and has been τ loc CINES (Grant No. c2015056853). successfully used to infer the functional structure for ground-state properties [63]. APPENDIX A: FUNCTIONAL STRUCTURE OF MANY-BODY TERMS APPENDIX B: MONTE CARLO SAMPLING ðMÞ ðMÞ The local residuals R ðX; tÞ¼ i∂ U ðX; tÞ − ðMÞ In order to solve the tVMC equations of motion, Eq. (6), E ðX; tÞ are vanishing if the Schrödinger equation is loc expectation values of some given operator O need to be exactly satisfied by the many-body wave function truncated computed over the many-body wave function ΦðX;tÞ. This ðMÞ at some order M. The local energy E ðX;tÞ, however, loc is achieved by means of Monte Carlo sampling of the may contain effective interaction terms involving a number probability distribution ΠðXÞ¼jΦðXÞj (in the following, of bodies larger than M, which leads to a systematic error in we omit explicit reference to the time t, assuming that all the truncation. However, the structure of these additional expectation values are taken over the wave function at a terms stemming from the local energy can be systemati- given fixed time). An efficient way of sampling the given cally used to deduce the functional structure of the higher- probability distribution is to devise a Markov chain of order terms. For example, the one-body truncated local configurations Xð1Þ; Xð2Þ; …; XðN − 1Þ; XðN Þ, which c c energy reads, for one-dimensional particles, are distribute according to ΠðXÞ. Quantum expectation values of a given operator O can then be obtained as ð1Þ 2 2 E ðX; tÞ¼ − f½∂ u ðx ; tÞ þ ∂ u ðx ; tÞg statistical expectation values over the Markov chain as x 1 i x 1 i loc i i N N X X c hΦjOjΦi 1 þ v ðx Þþ v ðx ;x Þ; ðA1Þ ≃ O (XðiÞ); ðB1Þ 1 i 2 i j loc hΦjΦi N i i≠j i¼1 031026-7 GIUSEPPE CARLEO et al. PHYS. REV. X 7, 031026 (2017) where O ðXÞ¼ðhXjOjΦi=hXjΦiÞ, and the equivalence where we introduce the Hermitian correlation matrix: loc is achieved in the limit N → ∞. ∂hρ ðtÞi The Markov chain is realized by the Metropolis-Hastings K t S ðtÞ¼ K;K algorithm. Given the current state of the Markov chain 0 ∂u XðiÞ, a configuration X is generated according to a given 0 0 ¼hρ ðtÞρ ðtÞi − hρ ðtÞi hρ ðtÞi : ðC2Þ K K K K 0 t t t transition probability T(XðiÞ → X ). The proposed con- figuration is then accepted [i.e., Xði þ 1Þ¼ X ] with At a given time, all the expectations values in Eq. (C1) can probability be explicitly computed with the stochastic approach described in Appendix B. We are therefore left with a 0 0 ΠðX Þ T(X → XðiÞ) linear system in the n unknowns u_ ðtÞ, which needs to A(XðiÞ → X ) ¼ min 1; ; var K Π(XðiÞ) T(XðiÞ → X ) be solved at each time t. In the presence of a large number of variational parameters ðB2Þ n , the solution of the linear system can be achieved using var iterative solvers, e.g., conjugate gradient methods, which do otherwise, it is rejected and Xði þ 1Þ¼ XðiÞ. not need to explicitly form the matrix S. Calling n the In the present tVMC calculations, we use simple iter number of iterations needed to obtain a solution for the linear transition probabilities in which a single particle is dis- system, the computational cost to solve Eq. (C1) is OðM × placed, while leaving all the other particle positions unchanged. In particular, a particle index p is chosen with n × n Þ as opposed to the OðM × n Þ operations var iter var uniform probability 1=N and the position of particle p is needed by a standard solver in which the matrix S is formed then displaced according to x ¼ x þ η , where η is a explicitly. In the present work, we resort to the minimal p p Δ Δ residual method, which is a variant of the Lanczos method, random number uniformly distributed in ½−ðΔ=2Þ; ðΔ=2Þ. working in the Krylov subspace spanned by the repeated The amplitude Δ is an adjustable parameter, and it can action of the matrix S onto an initial vector. In typical typically be chosen to be of the order of the average applications we obtain that n ≪ n , and several thou- interparticle distance. With this choice, the transition iter var sands of variational parameters can be efficiently treated. probability is simply This is of fundamental importance when the continuous (infinite-basis) limit must be taken, for which n → ∞. var Tðx → x Þ¼ ; ðB3Þ p p Once the unknowns u_ ðtÞ are determined, we can solve NΔ K numerically the first-order differential equations given in and the acceptance probability is therefore given by the mere Eq. (6) for given initial conditions u ð0Þ. In the present ratio of the probability distributions, ΠðX Þ=Π(XðiÞ). work, we adopt an adaptive fourth-order Runge-Kutta scheme for the integration of the differential equations. APPENDIX C: FINITE-BASIS PROJECTION APPENDIX D: LATTICE REGULARIZATION The numerical solution of the equations of motion FOR THE LIEB-LINIGER MODEL [Eq. (6)] requires the projection of the Jastrow fields We consider a general wave function Ψðx ;x ; …;x Þ u ðr; tÞ onto a finite basis. The continuous variable r is 1 2 N for N one-dimensional particles, governed by the Lieb- reduced to a finite set of P values for each order m, Liniger Hamiltonian. By means of the variational ðm; rÞ → ðr ;r ; …;r Þ. We introduce a super-index 1;m 2;m P;m K spanning all possible values of the discrete variables r . Monte Carlo method, we want to sample jΦðXÞj . This i;m is achieved via a lattice regularization, i.e., The complete set of variational parameters resulting from the projection on the finite basis can then be written as u ðtÞ and the associated functional derivatives read ρ ðtÞ. 2 2 K K dXjΦðXÞj ≃ jΦðl ;l ; …;l Þj ; 1 2 N The integro-differential equations [Eq. (6)] are then l ;l …l 1 2 N brought to the algebraic form, where l ¼f0;a; …;L − ag are discrete particle positions, a S 0u_ 0ðtÞ¼ −ihE ðtÞρ ðtÞi; ðC1Þ the lattice spacing, L the box size and N ¼ 1 þ L=a the K;K K loc K number of lattice sites. As a discretized Hamiltonian, we take 2 X ℏ 4 H Φðl …l Þ¼ − ½Φðl …l − a; …;l Þþ Φðl …l þ a; …;l Þ a 1 N 1 i N 1 i N 2ma 3 1 5 g − ½Φðl …l − 2a; …;l Þþ Φðl …l þ 2a; …;l Þ þ − Φðl …l ; …;l Þ þ Φðl …l Þ δðl ;l Þ: 1 i N 1 i N 1 i N 1 N i j 12 2 a i<j 031026-8 UNITARY DYNAMICS OF STRONGLY-INTERACTING BOSE … PHYS. REV. X 7, 031026 (2017) The first terms constitute just the fourth-order approximation separable from the deterministic propagation. The ampli- of the Laplacianvia central finite differences, whereas the last tude of these high-frequency oscillations also quantifies the term corresponds to the two-body delta interaction part. purely statistical error of our data. With this discretization, a two-body Jastrow factor reads Whereas exact diagonalization is limited to rather small basis sets, we can access a much larger basis within tVMC ð2Þ u ðx ; tÞ¼ u ðl ;l ; tÞ; method. In Fig. 4(b), we show results within the U 2 ij 2 i j approximation with L=a ¼ 80 and L=a ¼ 160 with time where u ða; b; tÞ is a time-dependent matrix of size 2 discretization Δt ∼ ðL=aÞ . We see that the basis set N × N , which, in 1D and in the presence of translational s s truncation in general introduces a dephasing at large symmetry, depends only on distða − bÞ; i.e., it has N =2 s enough time. ð2Þ variational parameters. The systematic error of U increases towards quenches to stronger interaction strength and becomes more visible APPENDIX E: BENCHMARK STUDY for γ ¼ 8; see Fig. 5(a). However, even in this case, the FOR N = 3 ON A LATTICE most important effect remains to be a simple dephasing; a small shift of averaged quantities is probable, but difficult Here, we use exact diagonalization of a Hamiltonian to quantify precisely. Introducing general three-body within a given finite basis for a quantitative test of our ð3Þ Jastrow fields U , described in detail in Appendix F, method. Exact diagonalization is limited to very small the systematic error for N ¼ 3 can be fully eliminated. systems on a finite basis, and we choose a system containing In Fig. 5(b), we also benchmark the possibility of N ¼ 3 particles on L=a ¼ 40 lattice sites as a simple, but calculating the off-diagonal one-body density matrix after highly nontrivial reference. In contrast to our comparison ð2Þ a quench. Here, the systematic error of U is more with BA methods, all observables can be accessed by exact pronounced at smaller γ in the long range and long time diagonalization, and we use the off-diagonal one-body f regime. density matrix g ðr; tÞ and the pair correlation function From our study of the three-particle problem, we con- g ðr; tÞ at different distances r ¼jx − x j of two particles 2 1 2 clude that truncation of the many-body wave function at the after time t, where the system is quenched from the non- interacting initial state, γ ¼ 0, to a final interaction, γ > 0, level of U may provide an excellent approximation for i f 2 g ðr; tÞ for quenches involving not too strong interaction to provide a benchmark on a more general observable. ð2Þ We first benchmark the influence of the time-step lattice strengths, γ ≲ 5. The systematic error due to the U size discretization Δt error on g . From Fig. 4(a), we see truncation is mainly a dephasing at large times involving that the tVMC dynamics is stable over a long time and the small relative errors of time-averaged quantities. Similar time step error can be brought to convergence. Further, we systematic dephasing errors will occur for too large time see that for final interaction γ ¼ 4, the truncation at the discretization or basis set truncation. ð2Þ Since our method provides a parametrization of the full two-body level U introduces only a small systematic wave function for a given time, many different observables error, mainly a dephasing effect, which is almost negligible can be evaluated via usual Monte Carlo methods. However, at the scale of the figure. Because of the stochastic noise of the quality of different observables may vary and depend the Monte Carlo integration, tVMC introduces additional more sensitively on the inclusion of higher-order high-frequency oscillations which are, however, well (a) (b) FIG. 4. Time-dependent expectation value of the two-body correlations after a quantum quench from a noninteracting state, γ ¼ 0,to γ ¼ 4 at three different distances, jx − x j¼ 0, L=10, L=4. Here, the system is on a lattice with L=a ¼ 40 lattice sites. The full line is f 1 2 ð2Þ obtained by exact diagonalization (ED) of the Hamiltonian; the other curves are from tVMC results truncated at the level of U . In (a), we show the convergence with different time-step discretization. In (b), we show the approach to the continuum for tVMC simulations ð2Þ using U for discretizations L=a ¼ 40, 80, and 160. 031026-9 GIUSEPPE CARLEO et al. PHYS. REV. X 7, 031026 (2017) (a) (b) FIG. 5. (a) Time-dependent expectation value of the two-body correlations g ðr; tÞ after a quantum quench from a noninteracting state, γ ¼ 0,to γ ¼ 8 at three different distances, jx − x j¼ 0, L=10, L=4. Here, the system is on a lattice with L=a ¼ 40 lattice sites. The i f 1 2 ð2Þ full line is obtained by exact diagonalization of the Hamiltonian; the other curves are from tVMC results truncated at the level of U or ð3Þ ð2Þ U . In contrast to γ ¼ 4 shown in Fig. 4, systematic differences of U compared to the exact results are more visible here; the exact ð3Þ dynamics is recovered by inclusion of three-body terms U into the tVMC wave function. (b) Off-diagonal single-particle density matrix g ðr; tÞ at three different distances, r ¼ L=10, L=4, and L=2. The exact results (ED–black dashed line) are for all purposes ð3Þ indistinguishable from the U curves. ðnÞ [2] D. Ceperley, Path Integrals in the Theory of Condensed correlations U with n> 2, as in the case of the single- Helium, Rev. Mod. Phys. 67, 279 (1995). body density matrix. Although these higher-order terms are [3] L. Pollet, Recent Developments in Quantum Monte Carlo computationally expensive, the scaling is not exponential, Simulations with Applications for Cold Gases, Rep. Prog. and we explicitly show that calculations with n ¼ 3 are Phys. 75, 094501 (2012). feasible. We note that the computational complexity may [4] T. Plisson, B. Allard, M. Holzmann, G. Salomon, A. Aspect, be further reduced by functional forms adapted to the P. Bouyer, and T. Bourdel, Coherence Properties of a Two- problem [26]. Dimensional Trapped Bose Gas around the Superfluid Transition, Phys. Rev. A 84, 061606 (2011). [5] N. G. Berloff and B. V. Svistunov, Scenario of Strongly ð3Þ APPENDIX F: GENERAL STRUCTURE OF U Nonequilibrated Bose-Einstein Condensation, Phys. Rev. A 66, 013603 (2002). For a general time-dependent wave function, we have to [6] N. Navon, A. L. Gaunt, R. P. Smith, and Z. Hadzibabic, go beyond the usual ground-state structure of the three- Emergence of a Turbulent Cascade in a Quantum Gas, body Jastrow term given in Eq. (A3). Here, we provide Nature (London) 539, 72 (2016). details of our three-body term in a general form beyond the [7] M. Gaudin, La Fonction d’Onde de Bethe (Masson, Paris, present application in one dimension. 1983). Introducing M basis functions, b ðrÞ, a ¼ 1; …;M,we [8] J.-S. Caux and R. M. Konik, Constructing the Generalized a α a can introduce many-body vectors [26], B ¼ r b ðr Þ, ij iα j ij Gibbs Ensemble after a Quantum Quench, Phys. Rev. Lett. where α ¼ 1; …;D indicates the summation over direc- 109, 175301 (2012). tions and i ¼ 1; …;N. The variational parameters of a [9] J.-S. Caux and F. H. L. Essler, Time Evolution of Local general three-body structure can then be written in terms of Observables After Quenching to an Integrable Model, Phys. Rev. Lett. 110, 257203 (2013). a matrix w , such that ab [10] M. Collura, S. Sotiriadis, and P. Calabrese, Equilibration of X X X a Tonks-Girardeau Gas Following a Trap Release, Phys. ab ab a b u ðr ; r ; r Þ¼ w W ;W ¼ B B : 3 1 2 3 ab iα iα Rev. Lett. 110, 245301 (2013). iα i≠j≠k ab [11] S. Trotzky, Y-A. Chen, A. Flesch, I. P. McCulloch, U. ðF1Þ Schollwöck, J. Eisert, and I. Bloch, Probing the Relaxation Towards Equilibrium in an Isolated Strongly Correlated One-Dimensional Bose Gas, Nat. Phys. 8, 325 (2012). In order to reduce the variational parameters (∼M ), we [12] T. Langen, R. Geiger, M. Kuhnert, B. Rauer, and J. may perform a singular value decomposition of the matrix Schmiedmayer, Local Emergence of Thermal Correlations w to reduce the effective degrees of freedom. ab in an Isolated Quantum Many-Body System, Nat. Phys. 9, 640 (2013). [13] D. Greif, G. Jotzu, M. Messer, R. Desbuquois, and T. Esslinger, Formation and Dynamics of Antiferromagnetic Correlations in Tunable Optical Lattices, Phys. Rev. Lett. [1] Nature Physics Insight on Non-equilibrium physics, Nat. Phys. 11, 103 (2015). 115, 260401 (2015). 031026-10 UNITARY DYNAMICS OF STRONGLY-INTERACTING BOSE … PHYS. REV. X 7, 031026 (2017) [14] M. Gring, M. Kuhnert, T. Langen, T. Kitagawa, B. [32] G. Carleo and M. Troyer, Solving the Quantum Many-Body Rauer, M. Schreitl, I. Mazets, D. A. Smith, E. Demler, Problem with Artificial Neural Networks, Science 355, 602 and J. Schmiedmayer, Relaxation and Prethermalization in (2017). an Isolated Quantum System, Science 337, 1318 (2012). [33] K. Ido, T. Ohgoe, and M. Imada, Time-Dependent Many- [15] M. Cominotti, D. Rossini, M. Rizzi, F. Hekking, and A. Variable Variational Monte Carlo Method for Nonequili- Minguzzi, Optimal Persistent Currents for Interacting brium Strongly Correlated Electron Systems, Phys. Rev. B Bosons on a Ring with a Gauge Field, Phys. Rev. Lett. 92, 245106 (2015). 113, 025301 (2014). [34] We have conveniently set the particle mass and the reduced [16] P. Calabrese and J. Cardy, Time Dependence of Correlation Planck constant to unity. Functions Following a Quantum Quench, Phys. Rev. Lett. [35] A. Bijl, The Lowest Wave Function of the Symmetrical Many 96, 136801 (2006). Particles System, Physica (Amsterdam) 7, 869 (1940). [17] S. R. White, Density Matrix Formulation for Quantum [36] R. B. Dingle, The Zero-Point Energy of a System of Renormalization Groups, Phys. Rev. Lett. 69, 2863 (1992). Particles, Philos. Mag. 40, 573 (1949). [18] S. R. White and A. E. Feiguin, Real-Time Evolution Using [37] R. Jastrow, Many-Body Problem with Strong Forces, Phys. the Density Matrix Renormalization Group, Phys. Rev. Lett. Rev. 98, 1479 (1955). 93, 076401 (2004). [38] E. Feenberg, Theory of Quantum Fluids, Pure and Applied [19] A. J. Daley, C. Kollath, U. Schollwock, and G. Vidal, Time- Physics (Academic Press, New York, 1969). Dependent Density-Matrix Renormalization-Group Using [39] C. L. Kane, S. Kivelson, D. H. Lee, and S. C. Zhang, Adaptive Effective Hilbert Spaces, J. Stat. Mech. (2004) General Validity of Jastrow-Laughlin Wave Functions, Phys. Rev. B 43, 3255 (1991). P04005. [40] P. A. M. Dirac, Note on Exchange Phenomena in the Thomas [20] F. B. Anders and A. Schiller, Real-Time Dynamics in Quantum-Impurity Systems: A Time-Dependent Numerical Atom, Math. Proc. Cambridge Philos. Soc. 26, 376 (1930). Renormalization-Group Approach, Phys. Rev. Lett. 95, [41] I. Frenkel, Wave Mechanics: Advanced General Theory,in 196801 (2005). The International Series of Monographs on Nuclear [21] M. Dolfi, B. Bauer, M. Troyer, and Z. Ristivojevic, Energy: Reactor Design Physics Vol. 2 (Clarendon Press, Multigrid Algorithms for Tensor Network States, Phys. Oxford, 1934). [42] The natural norm induced by a quantum Hilbert space is the Rev. Lett. 109, 020604 (2012). Fubini-Study norm, which is gauge invariant and therefore [22] M. Dolfi, A. Kantian, B. Bauer, and M. Troyer, Minimizing Nonadiabaticities in Optical-Lattice Loading, Phys. Rev. A insensitive to the unknown normalizations of the quantum 91, 033407 (2015). states we are dealing with here. [23] F. Verstraete and J. I. Cirac, Continuous Matrix Product States [43] E. H. Lieb and W. Liniger, Exact Analysis of an Interacting for Quantum Fields, Phys. Rev. Lett. 104, 190405 (2010). Bose Gas. I. The General Solution and the Ground State, [24] M. Ganahl, J. Rincon, and G. Vidal, Continuous Matrix Phys. Rev. 130, 1605 (1963). Product States for Quantum Fields: An Energy Minimiza- [44] H. Moritz, T. Stoferle, M. Kohl, and T. Esslinger, Exciting tion Algorithm., Phys. Rev. Lett. 118, 220402 (2017). Collective Oscillations in a Trapped 1D Gas, Phys. Rev. [25] C. J. Umrigar, J. Toulouse, C. Filippi, S. Sorella, and R. G. Lett. 91, 250402 (2003). Hennig, Alleviation of the Fermion-Sign Problem by Opti- [45] T. Kinoshita, T. Wenger, and D. S. Weiss, A Quantum mization of Many-Body Wave Functions, Phys. Rev. Lett. Newton’s Cradle, Nature (London) 440, 900 (2006). 98, 110201 (2007). [46] J. P. Ronzheimer, M. Schreiber, S. Braun, S. S. Hodgman, [26] M. Holzmann, B. Bernu, and D. M. Ceperley, Many-Body S. Langer, I. P. McCulloch, F. Heidrich-Meisner, I. Bloch, Wavefunctions for Normal Liquid He, Phys. Rev. B 74, and U. Schneider, Expansion Dynamics of Interacting 104510 (2006). Bosons in Homogeneous Lattices in One and Two Dimen- [27] M. Taddei, M. Ruggeri, S. Moroni, and M. Holzmann, sions, Phys. Rev. Lett. 110, 205301 (2013). Iterative Backflow Renormalization Procedure for Many- [47] B. Fang, G. Carleo, A. Johnson, and I. Bouchoule, Quench- Body Ground-State Wave Functions of Strongly Interacting Induced Breathing Mode of One-Dimensional Bose Gases, Normal Fermi Liquids, Phys. Rev. B 91, 115106 Phys. Rev. Lett. 113, 035301 (2014). [48] G. Boéris et al., Mott Transition for Strongly Interacting (2015). [28] G. Carleo, F. Becca, M. Schiro, and M. Fabrizio, Locali- One-Dimensional Bosons in a Shallow Periodic Potential, zation and Glassy Dynamics Of Many-Body Quantum Phys. Rev. A 93, 011601 (2016). Systems., Sci. Rep. 2, 243 (2012). [49] S. Sorella, Generalized Lanczos Algorithm for Variational [29] G. Carleo, F. Becca, L. Sanchez-Palencia, S. Sorella, and M. Quantum Monte Carlo, Phys. Rev. B 64, 024512 (2001). Fabrizio, Light-Cone Effect and Supersonic Correlations in [50] C. J. Umrigar and C. Filippi, Energy and Variance Opti- One- and Two-Dimensional Bosonic Superfluids, Phys. mization of Many-Body Wave Functions, Phys. Rev. Lett. Rev. A 89, 031602 (2014). 94, 150201 (2005). [30] L. Cevolani, G. Carleo, and L. Sanchez-Palencia, Protected [51] M. Girardeau, Relationship between Systems of Impen- Quasilocality in Quantum Systems with Long-Range Inter- etrable Bosons and Fermions in One Dimension, J. Math. actions, Phys. Rev. A 92, 041603 (2015). Phys. (N.Y.) 1, 516 (1960). [31] B. Blaß and H. Rieger, Test of Quantum Thermalization in [52] G. E. Astrakharchik and S. Giorgini, Correlation Functions the Two-Dimensional Transverse-Field Ising Model, Sci. and Momentum Distribution of One-Dimensional Bose Rep. 6, 38185 (2016). Systems, Phys. Rev. A 68, 031602 (2003). 031026-11 GIUSEPPE CARLEO et al. PHYS. REV. X 7, 031026 (2017) [53] P. Pippan, S. R. White, and H. G. Evertz, Efficient Matrix- [59] J. C. Zill, T. M. Wright, K. V. Kheruntsyan, T. Gasenzer, and Product State Method for Periodic Boundary Conditions, M. J. Davis, A Coordinate Bethe Ansatz Approach to the Phys. Rev. B 81, 081103 (2010). Calculation of Equilibrium and Nonequilibrium Correla- [54] G. Carleo, G. Boéris, M. Holzmann, and L. Sanchez- tions of the One-Dimensional Bose Gas, New J. Phys. 18, Palencia, Universal Superfluid Transition and Transport 045010 (2016). Properties of Two-Dimensional Dirty Bosons, Phys. Rev. [60] C. N. Yang and C. P. Yang, Thermodynamics of a Lett. 111, 050406 (2013). One-Dimensional System of Bosons with Repulsive Delta- [55] M. Boninsegni, N. Prokof’ev, and B. Svistunov, Worm Function Interaction, J. Math. Phys. (N.Y.) 10, 1115 Algorithm for Continuous-Space Path Integral Monte Carlo (1969). Simulations, Phys. Rev. Lett. 96, 070601 (2006). [61] M. Rigol, Fundamental Asymmetry in Quenches between [56] J. C. Zill, T. M. Wright, K. V. Kheruntsyan, T. Gasenzer, Integrable and Nonintegrable Systems, Phys. Rev. Lett. and M. J. Davis, Relaxation Dynamics of the Lieb-Liniger 116, 100601 (2016). Gas Following an Interaction Quench: A Coordinate [62] M. Kollar and M. Eckstein, Relaxation of a One- Bethe-Ansatz Analysis, Phys. Rev. A 91, 023611 (2015). Dimensional Mott Insulator after an Interaction Quench, [57] J. De Nardis, L. Piroli, and J.-S. Caux, Relaxation Dynamics Phys. Rev. A 78, 013626 (2008). of Local Observables in Integrable Systems., J. Phys. A 48, [63] M. Holzmann, D. M. Ceperley, C. Pierleoni, and K. Esler, 43FT01 (2015). Backflow Correlations for the Electron Gas and Metallic [58] J.-S. Caux and F. H. L. Essler, Time Evolution of Local Hydrogen, Phys. Rev. E 68, 046707 (2003). Observables After Quenching to an Integrable Model, Phys. Rev. Lett. 110, 257203 (2013). 031026-12 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Physical Review X American Physical Society (APS)

Unitary Dynamics of Strongly Interacting Bose Gases with the Time-Dependent Variational Monte Carlo Method in Continuous Space

Free
12 pages

Loading next page...
 
/lp/aps_physical/unitary-dynamics-of-strongly-interacting-bose-gases-with-the-time-K3xcjtxcGD
Publisher
The American Physical Society
Copyright
Copyright © Published by the American Physical Society
eISSN
2160-3308
D.O.I.
10.1103/PhysRevX.7.031026
Publisher site
See Article on Publisher Site

Abstract

PHYSICAL REVIEW X 7, 031026 (2017) Unitary Dynamics of Strongly Interacting Bose Gases with the Time-Dependent Variational Monte Carlo Method in Continuous Space Giuseppe Carleo Institute for Theoretical Physics, ETH Zurich, Wolfgang-Pauli-Strasse 27, 8093 Zurich, Switzerland Lorenzo Cevolani Laboratoire Charles Fabry, Institut d’Optique, CNRS, Univ Paris Sud 11, 2 avenue Augustin Fresnel, F-91127 Palaiseau cedex, France Laurent Sanchez-Palencia Laboratoire Charles Fabry, Institut d’Optique, CNRS, Univ Paris Sud 11, 2 avenue Augustin Fresnel, F-91127 Palaiseau cedex, France and Centre de Physique Théorique, Ecole Polytechnique, CNRS, Univ Paris-Saclay, F-91128 Palaiseau, France Markus Holzmann LPMMC, UMR 5493 of CNRS, Université Grenoble Alpes, F-38042 Grenoble, France and Institut Laue-Langevin, BP 156, F-38042 Grenoble Cedex 9, France (Received 19 December 2016; revised manuscript received 6 April 2017; published 8 August 2017) We introduce the time-dependent variational Monte Carlo method for continuous-space Bose gases. Our approach is based on the systematic expansion of the many-body wave function in terms of multibody correlations and is essentially exact up to adaptive truncation. The method is benchmarked by comparison to an 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 are also compared to quench action results available for noninteracting 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. DOI: 10.1103/PhysRevX.7.031026 Subject Areas: Atomic and Molecular Physics, Computational Physics, Quantum Physics I. INTRODUCTION more involved. Quantum Monte Carlo algorithms, the de facto tool for simulating quantum many-body systems at The study of equilibration and thermalization properties thermal equilibrium [2–4], cannot be directly used to study of complex many-body systems is of fundamental interest time-dependent unitary dynamics. Out-of-equilibrium for many areas of physics and natural sciences [1].For properties are then often treated on the basis of approx- systems governed by classical physics, an exact solution of imations drastically simplifying the microscopic physics. Newton’s equations of motion is often numerically feasible, Irreversibility is either enforced with an explicit breaking of using, for instance, molecular-dynamics simulations. For unitarity, e.g., within the quantum Boltzmann approach, or quantum systems, the mathematical structure of the time- the dynamics is reduced to mean-field description using dependent Schrödinger equation is instead fundamentally 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 Published by the American Physical Society under the terms of entanglement are ignored. the Creative Commons Attribution 4.0 International license. For specific systems, exact dynamical results can be Further distribution of this work must maintain attribution to derived. This is the case for integrable 1D models, for the author(s) and the published article’s title, journal citation, and DOI. which Bethe ansatz (BA) solutions exist [7]. However, also 2160-3308=17=7(3)=031026(12) 031026-1 Published by the American Physical Society GIUSEPPE CARLEO et al. PHYS. REV. X 7, 031026 (2017) in this case many open questions still persist. For example, further compare to the quench action approach for large the exact evaluation of correlation functions for out-of- systems approaching the thermodynamic limit. Finally, we equilibrium dynamics is at present an unsolved problem. apply our method to the study of general quenches from As a result, despite important theoretical and experimental arbitrary initial states, for which no exact results in the progress [8–15], a complete picture of thermalization (or its thermodynamic limit are currently available. absence), e.g., based on general quench protocols [16],is still missing. II. METHOD Numerical methods for strongly interacting systems face A. Expansion of the many-body wave function important challenges as well. Numerical renormalization group and density-matrix renormalization group (DMRG) Consider a nonrelativistic quantum system of N identical approaches provide an essentially exact description of bosons in d dimensions, and governed by the first- arbitrary 1D lattice systems in and out of equilibrium quantization Hamiltonian [17–20], but they have less predictive power when applied N N X X X to continuous-space systems. On one hand, multiscale 1 1 H ¼ − ∇ þ v ðx⃗ Þþ v ðx⃗ ; x⃗ Þ; ð1Þ 1 i 2 i j extensions of the DMRG optimization scheme to the limit 2 2 i¼1 i¼1 i≠j of continuous-space lattices [21] are, to date, limited to relatively small system sizes [22]. On the other hand, where v ðx⃗ Þ and v ðx;⃗ y⃗ Þ are, respectively, a one-body 1 2 efficient ground-state optimization schemes for continuous external potential and a pairwise interparticle interaction quantum field matrix product states (cMPS) [23] have been [34]. Without loss of generality, a time-dependent N-body introduced only very recently [24], and applications to state can be written as ΦðX;tÞ¼ exp ½UðX; tÞ, where quantum dynamics are still to be realized. A further X ¼ x⃗ ; x⃗ ; …; x⃗ is the ensemble of particle positions 1 2 N formidable challenge is the efficient extension of these and U is a complex-valued function of the N-particle approaches to higher dimensions, which is a fundamentally N×d coordinates, R → C. Since the Hamiltonian Eq. (1) hard problem. contains only two-body interactions, it is expected that an Another class of methods for strongly interacting systems expansion of U in terms of few-body Jastrow functions is based on variational Monte Carlo (VMC) methods, containing at most m-body terms rapidly converges combining highly entangled variational states with robust towards the exact solution. Truncating this expansion up stochastic optimization schemes [25]. Such approaches have to a certain order, M ≤ N, leads to the Bijl-Dingle-Jastrow- been successfully applied to the description of continuous Feenberg expansion [35–38], quantum systems, in any dimension and not only in 1D X X [26,27]. More recently, out-of-equilibrium dynamics has ðMÞ U ðX;tÞ¼ u ðx⃗ ;tÞþ u ðx⃗ ;x⃗ ;tÞ 1 i 2 i j become accessible with the extension of these methods to 2! i¼1 i≠j real-time unitary dynamics, within the time-dependent varia- tional Monte Carlo (tVMC) method [28,29].So far,the þ  u ðx⃗ ;x⃗ ;…;x⃗ ;tÞ; ð2Þ M i i i 1 2 M M! tVMC approach has been developed for lattice systems with i ≠i ≠…i 1 2 M bosonic [28,29], spin [30–32], and fermionic [33] statistics, yielding a description of dynamical properties with an where u ðr; tÞ are functions of m particle coordinates, accuracy often comparable with MPS-based approaches. r ¼ x⃗ ; x⃗ ; …; x⃗ , and of the time t. A global constraint on i i i 1 2 m In this paper, we extend the tVMC approach to access the function UðX;tÞ is given by particle statistics. In dynamical properties of interacting quantum systems in the bosonic case, we demand that Uðx⃗ ; x⃗ ; …; x⃗ Þ¼ 1 2 N continuous space. Our approach is based on a systematic Uðx⃗ ; x⃗ ; …; x⃗ Þ, for all particle permutations σ. σð1Þ σð2Þ σðNÞ expansion of the wave function in terms of few- In general, the functions u ðr; tÞ can have an arbitrarily body Jastrow correlation functions. Using the 1D Lieb- complex dependence on the m particle coordinates, which Liniger model as a test case, we first show that the inclusion can prove problematic for practical applications. of high-order correlations allows us to systematically Nonetheless, a simplified functional dependence can often approach the exact BA ground-state energy. Our results be imposed, resulting from the two-body character of the improve by orders of magnitude on previously published interactions in the original Hamiltonian. For m ≥ 3, VMC and cMPS results, and are in line with the latest state- u ðr; tÞ can be conveniently factorized in terms of general of-the-art developments in the field. We further compute two-particle vector and tensor functions following single-body and pair correlation functions, hardly acces- Ref. [26]. Details of this approach and the present imple- sible by current BA methods. We then calculate the time mentation are presented in Appendixes A and F. evolution of the contact pair correlation function following An appealing property of the many-body expansion a quench in the interaction strength. For the noninteracting Eq. (2) is that it is able to describe intrinsically nonlocal initial state, we benchmark our results to exact BA correlations in space. For instance, the two-body ⃗ ⃗ calculations available for a small number of bosons and u ðx ; x ; tÞ, as well as any m-body function u ðr; tÞ, 2 i j m 031026-2 UNITARY DYNAMICS OF STRONGLY-INTERACTING BOSE … PHYS. REV. X 7, 031026 (2017) can be long range in the particle separation jx⃗ − x⃗ j. This In practice, these equations are numerically solved for the i j time derivatives ∂ u ðr ; tÞ at each time step t. The expect- nonlocal spatial structure allows for a correct description of t p ation values taken over the time-dependent state h·i , which gapless phases, where a two-body expansion may already enter Eq. (6), are found via a stochastic sampling of the capture all the universal features, in the sense of the renormalization group approach [39]. This is in contrast probability distribution ΠðX;tÞ¼jΦðX;tÞj . This is effi- with the MPS decomposition of the wave function, which is ciently achieved by means of the Metropolis-Hastings intrinsically local in space. 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 B. Time-dependent variational Monte Carlo method state Eq. (2) after time integration. The time evolution of the variational state Eq. (2) is The tVMC approach as formulated here provides, in entirely determined by the time dependence of the Jastrow principle, an exact description of the real-time dynamics of functions u ðr; tÞ. In order to establish optimal equations the N-body system. The essential approximation lies, of motion for the variational parameters, we note that the however, in the truncation of the Bijl-Dingle-Jastrow- functional derivative of UðX; tÞ with respect to the varia- Feenberg expansion to the M most relevant terms. In tional, complex-valued, Jastrow functions u ðr; tÞ, practical applications, the M ¼ 2 or M ¼ 3 truncation is often sufficient. Systematic improvement beyond M ¼ 3 is δUðX; tÞ ≡ ρ ðrÞ; ð3Þ possible [26], but may require substantial computational δu ðr; tÞ effort. yields the m-body density operators III. LIEB-LINIGER MODEL X Y As a first application of the continuous-space tVMC ðrÞ¼ δðx − r Þ: ð4Þ m i j m! approach, we consider the Lieb-Liniger model [43]. On one i ≠i ≠i j 1 2 m hand, some exact results and numerical data are available, allowing us to benchmark the Jastrow expansion and the The expectation values of the operators ρ over the tVMC approach. On the other hand, several aspects of the state jΦðtÞi give the instantaneous m-body correlations. out-of-equilibrium dynamics of this model are unknown, For instance, hρ ðr ;r Þi ¼ hδðx − r Þδðx − r Þi , 2 1 2 t i<j i 1 j 2 t which we compute here for the first time using the tVMC where h  i ¼hΦðtÞj…jΦðtÞi=hΦðtÞjΦðtÞi, is propor- approach. tional to the two-point density-density correlation function. The Lieb-Liniger model describes N interacting bosons We can then express the time derivative of the truncated ðMÞ in one dimension with contact interactions. It corresponds variational state U ðX; tÞ using the functional deriva- to the Hamiltonian Eq. (1) with v ðxÞ¼ 0 and v ðx; yÞ¼ 1 2 tives, Eq. (3), as a sum of the few-body density operators up gδðx − yÞ, where g is the coupling constant. Here, we to the truncation M; i.e., consider periodic boundary conditions over a ring of length L and the particle density n ¼ N=L. The density depend- ðMÞ ence is as usual expressed in terms of the dimensionless ∂ U ðX;tÞ¼ drρ ðrÞ∂ u ðr; tÞ: ð5Þ t m t m m¼1 parameter γ ¼ mg=ℏ n. The Lieb-Liniger model is the prototypal model of continuous one-dimensional strongly The exact wave function satisfies the Schrödinger correlated gas exactly solvable by the Bethe ansatz [43]. equation i∂ UðX; tÞ¼ E ðX;tÞ, where E ðX; tÞ¼ t loc loc This model is experimentally realized in ultracold atomic ðhXjHjΦðtÞi=hXjΦðtÞiÞ is the so-called local energy. gases strongly confined in one-dimensional optical traps, The optimal time evolution of the truncated Bijl-Dingle- and several studies on out-of-equilibrium physics have Jastrow-Feenberg expansion Eq. (2) can be derived already been realized [14,44–48]. imposing the Dirac-Frenkel time-dependent variational As a result of translation invariance, we have principle [40,41]. In geometrical terms, this amounts u ðx; tÞ¼ const, and the first nontrivial term in the to minimizing the Hilbert-space norm of the residuals many-body expansion Eq. (2) is the two-body, translation ðMÞ ðMÞ ðMÞ invariant, function u ðjx − yj; tÞ. To compute the functional R ðX; tÞ ≡ ∥i∂ U ðX; tÞ − E ðX; tÞ∥, thus yielding t 2 loc derivatives of the many-body wave function, we proceed a variational many-body state as close as possible to the with a projection of the continuous Jastrow fields u ðr; tÞ exact one [42]. The minimization can be performed onto a finite-basis set. Here, we find it convenient to explicitly and yields a closed set of integro-differential represent both the field Hamiltonian Eq. (1) and the Jastrow equation for the Jastrow functions u ðr; tÞ: functions on a uniform mesh of spacing a, leading to n ¼ var ½L=ð2aÞ variational parameters for the two-body Jastrow δhρ ðrÞi δhHi m t t 0 0 dr ∂ u ðr ; tÞ¼ −i : ð6Þ t p term. In the following, our results are extrapolated to the δu ðr ; tÞ δu ðr; tÞ p m p¼1 continuous limit, corresponding to a → 0. The finite-basis 031026-3 GIUSEPPE CARLEO et al. PHYS. REV. X 7, 031026 (2017) projection as well as the numerical time integration of Eq. (6) the modulus of a Vandermonde determinant of plane are detailed in Appendix C and benchmarked against exact waves, corresponding to the two-body Jastrow function diagonalization results for N ¼ 3 in Appendix E. u ðrÞ¼ log ðsin rπ=LÞ [51]. To assess the overall quality of pair wave functions for ground-state properties, we start comparing the variational ground-state energies E obtained A. Ground-state properties for M ¼ 2 with the exact BA result. In Fig. 1(a), we show To assess the quality of truncated Bijl-Dingle-Jastrow- the relative error ΔE=E as a function of the interaction Feenberg expansions for the Lieb-Liniger model, we start −4 strength γ. We find that the relative error is lower than 10 our analysis considering ground-state properties. An exact for all values of γ and the accuracy of our two-body Jastrow solution can be found from the Bethe ansatz and gives function is superior to previously published variational access to exact ground-state energies and local properties results based on either cMPS [23] or VMC [52] methods [43]. Other nonlocal properties are substantially more [see Fig. 1(a) for a quantitative comparison]. Notice that the difficult to extract from the BA solution, and unbiased improvement with respect to previous VMC results is due results for ground-state correlation functions have not been to the larger variational freedom of our u ðrÞ function, reported so far. In order to determine the best possible which is not restricted to any specific functional form as variational description of the ground state within our many- done in Ref. [52]. body expansion, different strategies are possible. A first Even though the accuracy reached by the two-body possibility is to consider the imaginary-time evolution −τH Jastrow function may already be sufficient for most jΨðτÞi ¼ e jΨ i, which systematically converges to practical purposes for all values of γ, we also consider the exact ground state in the limit τ ≫ Δ , where Δ ¼ 1 1 higher-order terms with M ¼ 3; see Fig. (2)(a). The E − E is the gap with the first excited state on a finite 1 0 introduction of the third-order term yields a sizable system and provided that the trial state Ψ is nonorthogonal improvement such that the maximum error is about 3 to the exact ground state. Imaginary-time evolution in the orders of magnitude smaller than original cMPS results variational subspace can be implemented considering the [23], and features a similar accuracy of recently reported formal substitution t → −iτ in the tVMC equations cMPS results [24]. Overall, our approach reaches a pre- [Eq. (6)]. The resulting equations are equivalent to the cision on a continuous-space system, which is comparable stochastic reconfiguration approach [49]. However, direct to state-of-the-art MPS or DMRG results for gapless minimization of the variational energy can be significantly systems on a lattice [53]. more efficient, in particular, for systems becoming gapless Finally, to further assess the quality of our ground-state in the thermodynamic limit, where Δ ∼ polyð1=NÞ.Given ansatz beyond the total energy, we also study nonlocal the gapless nature of the Lieb-Liniger model, we find it properties of the ground-state wave function, which are computationally more efficient to adopt a Newton method not accessible by existing exact BA methods. In Figs. 1(b) to minimize the energy variance [50]. and 1(c), we show, respectively, our results for the For the ground state, the many-body expansion truncated off-diagonal part of the one-body density matrix, at M ¼ 2 is exact not only in the noninteracting limit γ ¼ 0 g ðrÞ ∝ hΨ ðrÞΨð0Þi, and for the pair correlation function, but also in the Tonks-Girardeau limit γ → ∞. In this † † fermionized limit, the wave function can be written as g ðrÞ ∝ hΨ ðrÞΨðrÞΨ ð0ÞΨð0Þi, where ΨðrÞ is the bosonic (a) (b) (c) FIG. 1. Ground-state properties of the Lieb-Liniger model as obtained from different variational approaches: (a) relative accuracy of ð2Þ ð3Þ the ground-state energy, (b) one-body density matrix, (c) density-density pair correlation functions. U and U denote results for the ð2Þ two- and three-body expansion (present work), U is the parametrized two-body Jastrow state of Ref. [52], cMPS results are from AG Ref. [23], and cMPS are those very recently reported in Ref. [24]. Distances r are in units of the inverse density 1=n. Our variational results are obtained for N ¼ 100 particles. Finite-size corrections on local quantities are negligible, and very mildly affect the reported large-distance correlations. Overall statistical errors are of the order of symbol sizes for (a) and line widths for (b) and (c). 031026-4 UNITARY DYNAMICS OF STRONGLY-INTERACTING BOSE … PHYS. REV. X 7, 031026 (2017) field operator. We find an overall excellent agreement with validation of our method accessing g ðr; tÞ at nonvanishing the results that have been obtained with cMPS in Ref. [24], distances. The comparison shown in Fig. 2(a) shows an except for some small deviations at large values of r, which overall good agreement. The tVMC and BA results are we attribute to residual finite-size effects in our approach. We indistinguishable for weak interactions (γ ¼ 1) at the scale find that the addition of the three-body terms does not of the figure; similarly, for this interaction strength, tVMC ð2Þ ð3Þ significantly change correlation functions. Already at the calculations based on U and U agree with each other two-body level, the present results are statistically indistin- up to possible dephasing effects not visible within our noise guishable from exact results obtained using our implemen- level. For larger interactions, we notice systematic but small tation [54] of the worm algorithm [55] and for the same differences between BA and tVMC with M ¼ 2 or M ¼ 3 system (not shown). results. These differences amount to a small increase in the amplitude of the oscillations. This effect tends to increase B. Quench dynamics with the interaction strength, being hardly visible for γ ¼ 1 and more pronounced for γ ¼ 10. However, these oscil- Having assessed the quality of the ansatz for local and lations result at large times from the discrete mode structure nonlocal ground-state properties, we now turn to the study due to the very small number of particles. They vanish in of the out-of-equilibrium properties of the Lieb-Liniger the physical thermodynamic limit. In turn, the comparison model. We focus on the description of the unitary dynamics at small particle numbers indicates an accuracy better induced by a global quantum quench of the interaction than a few percent for time-averaged quantities in the strength, from an initial value γ to a final value γ . Exact i f asymptotic large time limit up to γ ¼ 10, with results at BA results are available only in the case of a noninteracting M ¼ 3 systematically improving on the M ¼ 2 case. initial state (γ ¼ 0). Even in this case, the dynamical BA Concerning small and intermediate time scales, we do equations can be exactly solved only for a modest number not observe systematic deviations between the tVMC of particles with further truncation in the number of energy results and the BA solution. In particular, the relaxation eigenmodes [56,57], N ≲ 10, since the complexity of the times are remarkably well captured by the tVMC approach. BA solution increases exponentially with the number of On the basis of this comparison and of the comparison for particles. Simplifications in the thermodynamic limit are three particles with exact diagonalization results presented exploited by the quench action [58], and have recently been in Appendix E, we conclude that tVMC approach allows applied to quantum quenches starting from a noninteracting accurate quantitative studies of both the relaxation and initial state [57]. In the following, we first compare our equilibration dynamics. This careful benchmarking now tVMC results to these existing results, and then present new allows us to confidently apply the tVMC approach regimes results for quenches following a nonvanishing initial that are inaccessible to exact BA, namely large but finite N interaction strength. and long times, as well as the case of nonvanishing initial To assess the quality of the time-dependent wave interactions. function, we compare our results for the evolution of local Let us consider relaxation of density correlations for a density-density correlations g ð0;tÞ with the truncated BA large number of particles, close to the thermodynamic limit results obtained in Refs. [56,59] for a small number of (here, we use N ¼ 100). As we show in Fig. 2(b), we notice particles, N ≃ 6. Appendix E also provides further (a) (b) FIG. 2. Time-dependent expectation value of local two-body correlations after a quantum quench from a noninteracting state, γ ¼ 0, to γ . (a) tVMC results are compared with BA results obtained for a small number of particles [56,59]. The correlation function is rescaled to have g ð0; 0Þ¼ 1. (b) tVMC results for N ¼ 100 particles and γ ¼ 1, 2, 4, 8 (from top to bottom) compared to the quench 2 f action (QA) predictions from Ref. [57] (dashed lines), to the Boltzmann thermal averages at the effective temperature T (dot-dashed lines), and to the GGE thermal averages prediction (rightmost dashed lines). Statistical error bars on tVMC data are of the order of the line widths. 031026-5 GIUSEPPE CARLEO et al. PHYS. REV. X 7, 031026 (2017) As we show in Fig. 2(b), Boltzmann thermalization certainly does not occur in the case for the Lieb-Liniger model when quenching from a noninteracting state, γ ¼ 0, where we find g¯ ≠ g . 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 believed that the generalized Gibbs ensemble (GGE) is the correct thermal (a) (b) distribution approached after the quench [61,62]. Several constructive approaches for the GGE have been put forward in past years [8,9,57], and the quench action predictions reported in Fig. 2(b) converge to the GGE FIG. 3. Time-dependent expectation value of local two-body predictions for the thermal values. In Fig. 2(b), we also correlations after a quantum quench from the interacting ground GGE show the thermal GGE values g (rightmost dashed state at γ ¼ 1 (a) and γ ¼ 4 (b). Long-term dynamical averages i i lines), and note that our results are much closer to the GGE (red continuous lines) are compared to thermal averages at the predictions than the simple BE. Deviations from the effective temperature set by energy conservation (black dashed asymptotic GGE results are observed at large γ , a regime lines). Uncertainties on the thermal averages are of the order of the line width, and are larger for small values of γ . 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 γ ≠ 0, GGE predic- that the amplitudes of the large-time oscillations, attributed tions are fundamentally harder to obtain than for the to the discrete mode spectrum, are now drastically sup- noninteracting initial states, and the BE is the only pressed compared to the quenches with N ¼ 6. After an reference thermal distribution we can compare with at this initial relaxation phase, the quantity g ð0;tÞ approaches a stage. From our results we observe that the difference stationary value. Comparing our curves with those obtained between g¯ and the simple BE prediction g is quantita- with the quench action method, we find a qualitatively good 2 2 tively reduced; see Fig. 3. In particular, for γ ¼ 4, the agreement, albeit a general tendency to underestimate the stationary values g¯ are quantitatively close to the ones quench action predictions is observed. 2 predicted by the Boltzmann thermal distribution at the We now turn to quenches from interacting initial states effective temperature T . Even though this quantitative (γ ≠ 0) to different interacting final states for which no agreement is likely to be coincidental, the regimes of results have been obtained by means of exact BA nor the parameter quenches we study here provide guidance for simplified quench action method so far. In Fig. 3, we show future experimental studies. In particular, it will be of great the asymptotic equilibrium values obtained with our tVMC interest to understand whether a crossover from a strongly approach for quantum quenches from γ ¼ 1 (left-hand non-Boltzmann to a close-to-Boltzmann thermal behavior panel) and γ ¼ 4 (right-hand panel) to several values of might occur as a function of the initial interaction strength γ . Since, by the variational theorem, the ground state of H f i also for other local observables. gives an upper bound for the ground-state energy of H ,the system is pushed into a linear combination of excited states of the final Hamiltonian. For systems able to thermalize to the IV. CONCLUSIONS Boltzmann ensemble (BE), relaxation to a stationary state ⋆ In this paper, we introduce a novel approach to the −H =T described by the density matrix ρ ⋆ ¼ e ,atan dynamics of strongly correlated quantum systems in effective temperature T , would occur. Comparing the sta- continuous space. Our method is based on a correlated tionary value g¯ ð0Þ of our tVMC calculations at long times to many-body wave function systematically expanded in the thermal values of the pair correlation functions g ð0Þ,a terms of reduced m-body Jastrow functions. The unitary necessary condition for simple Boltzmann thermalization is dynamics in the subspace of these correlated states is T ⋆ given by g¯ ¼ g . The effective temperature T is deter- realized using the time-dependent variational Monte Carlo mined by imposing the energy expectation value of the final method. We demonstrate the possibility of performing Hamiltonian H in the ground state Φ ðγ Þ of the initial f 0 i calculations up to the three-body level, m ≤ 3, for the Hamiltonian, hH i ⋆ ¼hΦ ðγ ÞjH jΦ ðγ Þi. Here, the ther- f 0 i f 0 i T Lieb-Liniger model, for both static and dynamical proper- mal expectationvalue, hH i ⋆ at the equilibrium temperature f ties. The improvement from m ¼ 2 to m ¼ 3 provides an T is computed from the Yang-Yang BA equations [60].The internal criterion to judge the validity of our results quantity hH i ⋆ then depends on a single parameter T , whenever exact results are unavailable. Benchmarking f T which is fitted to match the value of hΦ ðγ ÞjH jΦ ðγ Þi. the tVMC approach with exact or numerical approaches 0 i f 0 i 031026-6 UNITARY DYNAMICS OF STRONGLY-INTERACTING BOSE … PHYS. REV. X 7, 031026 (2017) whenever available, we find a very good agreement with and contains a two-body term that cannot be accounted for existing results. For static properties, our approach is at the exactly by u . Introduction of a symmetric two-body level of state-of-the-art MPS techniques in lattice systems Jastrow factor u ðx ;x ; tÞ then leads to 2 i j and of the latest cMPS results for interacting gases. For dynamical properties, we investigate, for the first time, ð2Þ ð1Þ E ðX; tÞ¼ E ðX; tÞþ loc loc general interaction quenches, which are at the moment unaccessible to Bethe ansatz approaches. Since the general − ½∂ u ðx ; tÞ∂ u ðx ;x ; tÞþ x 1 i x 2 i j i i structure of our t-VMC method does not depend on the i≠j dimensionality of the system, it can be directly applied to bosonic systems in higher dimensions with a polynomial − ∂ u ðx ;x ; tÞþ x 2 i j i≠j increase in computational cost. The methods we present XX here therefore pave the way to accurate out-of-equilibrium − ∂ u ðx ;x ; tÞ∂ u ðx ;x ; tÞ: x 2 i j x 2 i k dynamics of two- and three-dimensional quantum gases i i i≠j k≠i and fluids beyond mean-field approximations. ðA2Þ ACKNOWLEDGMENTS In the latter expression, one recognizes an effective two- We acknowledge discussions with M. Dolfi, J. De body term which can be accounted for by u and an Nardis, M. Fagotti, T. Osborne, and M. Troyer. We thank additional three-body term in the form of a product of two- M. Ganahl for providing the cMPS results in Fig. 1, J. Zill body functions. The functional form of the three-body for the Bethe ansatz results in Fig. 2(a), and J. De Nardis for Jastrow term can therefore be deduced from this additional the quench actions results in Fig. 2(b). This research was term and formed accordingly: supported by the Marie Curie IEF Program (FP7/2007– 2013 Grant Agreement No. 327143), the European u ðx ;x ;x ; tÞ¼ u¯ ðx ;x ; tÞu¯ ðx ;x ; tÞ; ðA3Þ 3 i j k 3 i j 3 j k Research Council Starting Grant “ALoGlaDis” (FP7/ 2007-2013 Grant Agreement No. 256294) and Advanced with two-body functions u¯ ðx ;x ; tÞ containing new varia- 3 i j Grant “SIMCOFE” (FP7/2007-2013 Grant Agreement tional parameters to be determined. Upon pursuing this No. 290464), the European Commission FET-Proactive approach, the expansion can be systematically pushed to QUIC (H2020 Grant No. 641122), the French ANR-16- higher orders and the functional structure of the higher- CE30-0023-03 (THERMOLOC), and the Swiss National order functions inferred. The same constructive approach Science Foundation through NCCR QSIT. It was per- we discuss here is also valid for the Schrödinger equation in formed using HPC resources from GENCI-CCRT/ imaginary time, ∂ UðX; τÞ¼ −E ðx; τÞ, and has been τ loc CINES (Grant No. c2015056853). successfully used to infer the functional structure for ground-state properties [63]. APPENDIX A: FUNCTIONAL STRUCTURE OF MANY-BODY TERMS APPENDIX B: MONTE CARLO SAMPLING ðMÞ ðMÞ The local residuals R ðX; tÞ¼ i∂ U ðX; tÞ − ðMÞ In order to solve the tVMC equations of motion, Eq. (6), E ðX; tÞ are vanishing if the Schrödinger equation is loc expectation values of some given operator O need to be exactly satisfied by the many-body wave function truncated computed over the many-body wave function ΦðX;tÞ. This ðMÞ at some order M. The local energy E ðX;tÞ, however, loc is achieved by means of Monte Carlo sampling of the may contain effective interaction terms involving a number probability distribution ΠðXÞ¼jΦðXÞj (in the following, of bodies larger than M, which leads to a systematic error in we omit explicit reference to the time t, assuming that all the truncation. However, the structure of these additional expectation values are taken over the wave function at a terms stemming from the local energy can be systemati- given fixed time). An efficient way of sampling the given cally used to deduce the functional structure of the higher- probability distribution is to devise a Markov chain of order terms. For example, the one-body truncated local configurations Xð1Þ; Xð2Þ; …; XðN − 1Þ; XðN Þ, which c c energy reads, for one-dimensional particles, are distribute according to ΠðXÞ. Quantum expectation values of a given operator O can then be obtained as ð1Þ 2 2 E ðX; tÞ¼ − f½∂ u ðx ; tÞ þ ∂ u ðx ; tÞg statistical expectation values over the Markov chain as x 1 i x 1 i loc i i N N X X c hΦjOjΦi 1 þ v ðx Þþ v ðx ;x Þ; ðA1Þ ≃ O (XðiÞ); ðB1Þ 1 i 2 i j loc hΦjΦi N i i≠j i¼1 031026-7 GIUSEPPE CARLEO et al. PHYS. REV. X 7, 031026 (2017) where O ðXÞ¼ðhXjOjΦi=hXjΦiÞ, and the equivalence where we introduce the Hermitian correlation matrix: loc is achieved in the limit N → ∞. ∂hρ ðtÞi The Markov chain is realized by the Metropolis-Hastings K t S ðtÞ¼ K;K algorithm. Given the current state of the Markov chain 0 ∂u XðiÞ, a configuration X is generated according to a given 0 0 ¼hρ ðtÞρ ðtÞi − hρ ðtÞi hρ ðtÞi : ðC2Þ K K K K 0 t t t transition probability T(XðiÞ → X ). The proposed con- figuration is then accepted [i.e., Xði þ 1Þ¼ X ] with At a given time, all the expectations values in Eq. (C1) can probability be explicitly computed with the stochastic approach described in Appendix B. We are therefore left with a 0 0 ΠðX Þ T(X → XðiÞ) linear system in the n unknowns u_ ðtÞ, which needs to A(XðiÞ → X ) ¼ min 1; ; var K Π(XðiÞ) T(XðiÞ → X ) be solved at each time t. In the presence of a large number of variational parameters ðB2Þ n , the solution of the linear system can be achieved using var iterative solvers, e.g., conjugate gradient methods, which do otherwise, it is rejected and Xði þ 1Þ¼ XðiÞ. not need to explicitly form the matrix S. Calling n the In the present tVMC calculations, we use simple iter number of iterations needed to obtain a solution for the linear transition probabilities in which a single particle is dis- system, the computational cost to solve Eq. (C1) is OðM × placed, while leaving all the other particle positions unchanged. In particular, a particle index p is chosen with n × n Þ as opposed to the OðM × n Þ operations var iter var uniform probability 1=N and the position of particle p is needed by a standard solver in which the matrix S is formed then displaced according to x ¼ x þ η , where η is a explicitly. In the present work, we resort to the minimal p p Δ Δ residual method, which is a variant of the Lanczos method, random number uniformly distributed in ½−ðΔ=2Þ; ðΔ=2Þ. working in the Krylov subspace spanned by the repeated The amplitude Δ is an adjustable parameter, and it can action of the matrix S onto an initial vector. In typical typically be chosen to be of the order of the average applications we obtain that n ≪ n , and several thou- interparticle distance. With this choice, the transition iter var sands of variational parameters can be efficiently treated. probability is simply This is of fundamental importance when the continuous (infinite-basis) limit must be taken, for which n → ∞. var Tðx → x Þ¼ ; ðB3Þ p p Once the unknowns u_ ðtÞ are determined, we can solve NΔ K numerically the first-order differential equations given in and the acceptance probability is therefore given by the mere Eq. (6) for given initial conditions u ð0Þ. In the present ratio of the probability distributions, ΠðX Þ=Π(XðiÞ). work, we adopt an adaptive fourth-order Runge-Kutta scheme for the integration of the differential equations. APPENDIX C: FINITE-BASIS PROJECTION APPENDIX D: LATTICE REGULARIZATION The numerical solution of the equations of motion FOR THE LIEB-LINIGER MODEL [Eq. (6)] requires the projection of the Jastrow fields We consider a general wave function Ψðx ;x ; …;x Þ u ðr; tÞ onto a finite basis. The continuous variable r is 1 2 N for N one-dimensional particles, governed by the Lieb- reduced to a finite set of P values for each order m, Liniger Hamiltonian. By means of the variational ðm; rÞ → ðr ;r ; …;r Þ. We introduce a super-index 1;m 2;m P;m K spanning all possible values of the discrete variables r . Monte Carlo method, we want to sample jΦðXÞj . This i;m is achieved via a lattice regularization, i.e., The complete set of variational parameters resulting from the projection on the finite basis can then be written as u ðtÞ and the associated functional derivatives read ρ ðtÞ. 2 2 K K dXjΦðXÞj ≃ jΦðl ;l ; …;l Þj ; 1 2 N The integro-differential equations [Eq. (6)] are then l ;l …l 1 2 N brought to the algebraic form, where l ¼f0;a; …;L − ag are discrete particle positions, a S 0u_ 0ðtÞ¼ −ihE ðtÞρ ðtÞi; ðC1Þ the lattice spacing, L the box size and N ¼ 1 þ L=a the K;K K loc K number of lattice sites. As a discretized Hamiltonian, we take 2 X ℏ 4 H Φðl …l Þ¼ − ½Φðl …l − a; …;l Þþ Φðl …l þ a; …;l Þ a 1 N 1 i N 1 i N 2ma 3 1 5 g − ½Φðl …l − 2a; …;l Þþ Φðl …l þ 2a; …;l Þ þ − Φðl …l ; …;l Þ þ Φðl …l Þ δðl ;l Þ: 1 i N 1 i N 1 i N 1 N i j 12 2 a i<j 031026-8 UNITARY DYNAMICS OF STRONGLY-INTERACTING BOSE … PHYS. REV. X 7, 031026 (2017) The first terms constitute just the fourth-order approximation separable from the deterministic propagation. The ampli- of the Laplacianvia central finite differences, whereas the last tude of these high-frequency oscillations also quantifies the term corresponds to the two-body delta interaction part. purely statistical error of our data. With this discretization, a two-body Jastrow factor reads Whereas exact diagonalization is limited to rather small basis sets, we can access a much larger basis within tVMC ð2Þ u ðx ; tÞ¼ u ðl ;l ; tÞ; method. In Fig. 4(b), we show results within the U 2 ij 2 i j approximation with L=a ¼ 80 and L=a ¼ 160 with time where u ða; b; tÞ is a time-dependent matrix of size 2 discretization Δt ∼ ðL=aÞ . We see that the basis set N × N , which, in 1D and in the presence of translational s s truncation in general introduces a dephasing at large symmetry, depends only on distða − bÞ; i.e., it has N =2 s enough time. ð2Þ variational parameters. The systematic error of U increases towards quenches to stronger interaction strength and becomes more visible APPENDIX E: BENCHMARK STUDY for γ ¼ 8; see Fig. 5(a). However, even in this case, the FOR N = 3 ON A LATTICE most important effect remains to be a simple dephasing; a small shift of averaged quantities is probable, but difficult Here, we use exact diagonalization of a Hamiltonian to quantify precisely. Introducing general three-body within a given finite basis for a quantitative test of our ð3Þ Jastrow fields U , described in detail in Appendix F, method. Exact diagonalization is limited to very small the systematic error for N ¼ 3 can be fully eliminated. systems on a finite basis, and we choose a system containing In Fig. 5(b), we also benchmark the possibility of N ¼ 3 particles on L=a ¼ 40 lattice sites as a simple, but calculating the off-diagonal one-body density matrix after highly nontrivial reference. In contrast to our comparison ð2Þ a quench. Here, the systematic error of U is more with BA methods, all observables can be accessed by exact pronounced at smaller γ in the long range and long time diagonalization, and we use the off-diagonal one-body f regime. density matrix g ðr; tÞ and the pair correlation function From our study of the three-particle problem, we con- g ðr; tÞ at different distances r ¼jx − x j of two particles 2 1 2 clude that truncation of the many-body wave function at the after time t, where the system is quenched from the non- interacting initial state, γ ¼ 0, to a final interaction, γ > 0, level of U may provide an excellent approximation for i f 2 g ðr; tÞ for quenches involving not too strong interaction to provide a benchmark on a more general observable. ð2Þ We first benchmark the influence of the time-step lattice strengths, γ ≲ 5. The systematic error due to the U size discretization Δt error on g . From Fig. 4(a), we see truncation is mainly a dephasing at large times involving that the tVMC dynamics is stable over a long time and the small relative errors of time-averaged quantities. Similar time step error can be brought to convergence. Further, we systematic dephasing errors will occur for too large time see that for final interaction γ ¼ 4, the truncation at the discretization or basis set truncation. ð2Þ Since our method provides a parametrization of the full two-body level U introduces only a small systematic wave function for a given time, many different observables error, mainly a dephasing effect, which is almost negligible can be evaluated via usual Monte Carlo methods. However, at the scale of the figure. Because of the stochastic noise of the quality of different observables may vary and depend the Monte Carlo integration, tVMC introduces additional more sensitively on the inclusion of higher-order high-frequency oscillations which are, however, well (a) (b) FIG. 4. Time-dependent expectation value of the two-body correlations after a quantum quench from a noninteracting state, γ ¼ 0,to γ ¼ 4 at three different distances, jx − x j¼ 0, L=10, L=4. Here, the system is on a lattice with L=a ¼ 40 lattice sites. The full line is f 1 2 ð2Þ obtained by exact diagonalization (ED) of the Hamiltonian; the other curves are from tVMC results truncated at the level of U . In (a), we show the convergence with different time-step discretization. In (b), we show the approach to the continuum for tVMC simulations ð2Þ using U for discretizations L=a ¼ 40, 80, and 160. 031026-9 GIUSEPPE CARLEO et al. PHYS. REV. X 7, 031026 (2017) (a) (b) FIG. 5. (a) Time-dependent expectation value of the two-body correlations g ðr; tÞ after a quantum quench from a noninteracting state, γ ¼ 0,to γ ¼ 8 at three different distances, jx − x j¼ 0, L=10, L=4. Here, the system is on a lattice with L=a ¼ 40 lattice sites. The i f 1 2 ð2Þ full line is obtained by exact diagonalization of the Hamiltonian; the other curves are from tVMC results truncated at the level of U or ð3Þ ð2Þ U . In contrast to γ ¼ 4 shown in Fig. 4, systematic differences of U compared to the exact results are more visible here; the exact ð3Þ dynamics is recovered by inclusion of three-body terms U into the tVMC wave function. (b) Off-diagonal single-particle density matrix g ðr; tÞ at three different distances, r ¼ L=10, L=4, and L=2. The exact results (ED–black dashed line) are for all purposes ð3Þ indistinguishable from the U curves. ðnÞ [2] D. Ceperley, Path Integrals in the Theory of Condensed correlations U with n> 2, as in the case of the single- Helium, Rev. Mod. Phys. 67, 279 (1995). body density matrix. Although these higher-order terms are [3] L. Pollet, Recent Developments in Quantum Monte Carlo computationally expensive, the scaling is not exponential, Simulations with Applications for Cold Gases, Rep. Prog. and we explicitly show that calculations with n ¼ 3 are Phys. 75, 094501 (2012). feasible. We note that the computational complexity may [4] T. Plisson, B. Allard, M. Holzmann, G. Salomon, A. Aspect, be further reduced by functional forms adapted to the P. Bouyer, and T. Bourdel, Coherence Properties of a Two- problem [26]. Dimensional Trapped Bose Gas around the Superfluid Transition, Phys. Rev. A 84, 061606 (2011). [5] N. G. Berloff and B. V. Svistunov, Scenario of Strongly ð3Þ APPENDIX F: GENERAL STRUCTURE OF U Nonequilibrated Bose-Einstein Condensation, Phys. Rev. A 66, 013603 (2002). For a general time-dependent wave function, we have to [6] N. Navon, A. L. Gaunt, R. P. Smith, and Z. Hadzibabic, go beyond the usual ground-state structure of the three- Emergence of a Turbulent Cascade in a Quantum Gas, body Jastrow term given in Eq. (A3). Here, we provide Nature (London) 539, 72 (2016). details of our three-body term in a general form beyond the [7] M. Gaudin, La Fonction d’Onde de Bethe (Masson, Paris, present application in one dimension. 1983). Introducing M basis functions, b ðrÞ, a ¼ 1; …;M,we [8] J.-S. Caux and R. M. Konik, Constructing the Generalized a α a can introduce many-body vectors [26], B ¼ r b ðr Þ, ij iα j ij Gibbs Ensemble after a Quantum Quench, Phys. Rev. Lett. where α ¼ 1; …;D indicates the summation over direc- 109, 175301 (2012). tions and i ¼ 1; …;N. The variational parameters of a [9] J.-S. Caux and F. H. L. Essler, Time Evolution of Local general three-body structure can then be written in terms of Observables After Quenching to an Integrable Model, Phys. Rev. Lett. 110, 257203 (2013). a matrix w , such that ab [10] M. Collura, S. Sotiriadis, and P. Calabrese, Equilibration of X X X a Tonks-Girardeau Gas Following a Trap Release, Phys. ab ab a b u ðr ; r ; r Þ¼ w W ;W ¼ B B : 3 1 2 3 ab iα iα Rev. Lett. 110, 245301 (2013). iα i≠j≠k ab [11] S. Trotzky, Y-A. Chen, A. Flesch, I. P. McCulloch, U. ðF1Þ Schollwöck, J. Eisert, and I. Bloch, Probing the Relaxation Towards Equilibrium in an Isolated Strongly Correlated One-Dimensional Bose Gas, Nat. Phys. 8, 325 (2012). In order to reduce the variational parameters (∼M ), we [12] T. Langen, R. Geiger, M. Kuhnert, B. Rauer, and J. may perform a singular value decomposition of the matrix Schmiedmayer, Local Emergence of Thermal Correlations w to reduce the effective degrees of freedom. ab in an Isolated Quantum Many-Body System, Nat. Phys. 9, 640 (2013). [13] D. Greif, G. Jotzu, M. Messer, R. Desbuquois, and T. Esslinger, Formation and Dynamics of Antiferromagnetic Correlations in Tunable Optical Lattices, Phys. Rev. Lett. [1] Nature Physics Insight on Non-equilibrium physics, Nat. Phys. 11, 103 (2015). 115, 260401 (2015). 031026-10 UNITARY DYNAMICS OF STRONGLY-INTERACTING BOSE … PHYS. REV. X 7, 031026 (2017) [14] M. Gring, M. Kuhnert, T. Langen, T. Kitagawa, B. [32] G. Carleo and M. Troyer, Solving the Quantum Many-Body Rauer, M. Schreitl, I. Mazets, D. A. Smith, E. Demler, Problem with Artificial Neural Networks, Science 355, 602 and J. Schmiedmayer, Relaxation and Prethermalization in (2017). an Isolated Quantum System, Science 337, 1318 (2012). [33] K. Ido, T. Ohgoe, and M. Imada, Time-Dependent Many- [15] M. Cominotti, D. Rossini, M. Rizzi, F. Hekking, and A. Variable Variational Monte Carlo Method for Nonequili- Minguzzi, Optimal Persistent Currents for Interacting brium Strongly Correlated Electron Systems, Phys. Rev. B Bosons on a Ring with a Gauge Field, Phys. Rev. Lett. 92, 245106 (2015). 113, 025301 (2014). [34] We have conveniently set the particle mass and the reduced [16] P. Calabrese and J. Cardy, Time Dependence of Correlation Planck constant to unity. Functions Following a Quantum Quench, Phys. Rev. Lett. [35] A. Bijl, The Lowest Wave Function of the Symmetrical Many 96, 136801 (2006). Particles System, Physica (Amsterdam) 7, 869 (1940). [17] S. R. White, Density Matrix Formulation for Quantum [36] R. B. Dingle, The Zero-Point Energy of a System of Renormalization Groups, Phys. Rev. Lett. 69, 2863 (1992). Particles, Philos. Mag. 40, 573 (1949). [18] S. R. White and A. E. Feiguin, Real-Time Evolution Using [37] R. Jastrow, Many-Body Problem with Strong Forces, Phys. the Density Matrix Renormalization Group, Phys. Rev. Lett. Rev. 98, 1479 (1955). 93, 076401 (2004). [38] E. Feenberg, Theory of Quantum Fluids, Pure and Applied [19] A. J. Daley, C. Kollath, U. Schollwock, and G. Vidal, Time- Physics (Academic Press, New York, 1969). Dependent Density-Matrix Renormalization-Group Using [39] C. L. Kane, S. Kivelson, D. H. Lee, and S. C. Zhang, Adaptive Effective Hilbert Spaces, J. Stat. Mech. (2004) General Validity of Jastrow-Laughlin Wave Functions, Phys. Rev. B 43, 3255 (1991). P04005. [40] P. A. M. Dirac, Note on Exchange Phenomena in the Thomas [20] F. B. Anders and A. Schiller, Real-Time Dynamics in Quantum-Impurity Systems: A Time-Dependent Numerical Atom, Math. Proc. Cambridge Philos. Soc. 26, 376 (1930). Renormalization-Group Approach, Phys. Rev. Lett. 95, [41] I. Frenkel, Wave Mechanics: Advanced General Theory,in 196801 (2005). The International Series of Monographs on Nuclear [21] M. Dolfi, B. Bauer, M. Troyer, and Z. Ristivojevic, Energy: Reactor Design Physics Vol. 2 (Clarendon Press, Multigrid Algorithms for Tensor Network States, Phys. Oxford, 1934). [42] The natural norm induced by a quantum Hilbert space is the Rev. Lett. 109, 020604 (2012). Fubini-Study norm, which is gauge invariant and therefore [22] M. Dolfi, A. Kantian, B. Bauer, and M. Troyer, Minimizing Nonadiabaticities in Optical-Lattice Loading, Phys. Rev. A insensitive to the unknown normalizations of the quantum 91, 033407 (2015). states we are dealing with here. [23] F. Verstraete and J. I. Cirac, Continuous Matrix Product States [43] E. H. Lieb and W. Liniger, Exact Analysis of an Interacting for Quantum Fields, Phys. Rev. Lett. 104, 190405 (2010). Bose Gas. I. The General Solution and the Ground State, [24] M. Ganahl, J. Rincon, and G. Vidal, Continuous Matrix Phys. Rev. 130, 1605 (1963). Product States for Quantum Fields: An Energy Minimiza- [44] H. Moritz, T. Stoferle, M. Kohl, and T. Esslinger, Exciting tion Algorithm., Phys. Rev. Lett. 118, 220402 (2017). Collective Oscillations in a Trapped 1D Gas, Phys. Rev. [25] C. J. Umrigar, J. Toulouse, C. Filippi, S. Sorella, and R. G. Lett. 91, 250402 (2003). Hennig, Alleviation of the Fermion-Sign Problem by Opti- [45] T. Kinoshita, T. Wenger, and D. S. Weiss, A Quantum mization of Many-Body Wave Functions, Phys. Rev. Lett. Newton’s Cradle, Nature (London) 440, 900 (2006). 98, 110201 (2007). [46] J. P. Ronzheimer, M. Schreiber, S. Braun, S. S. Hodgman, [26] M. Holzmann, B. Bernu, and D. M. Ceperley, Many-Body S. Langer, I. P. McCulloch, F. Heidrich-Meisner, I. Bloch, Wavefunctions for Normal Liquid He, Phys. Rev. B 74, and U. Schneider, Expansion Dynamics of Interacting 104510 (2006). Bosons in Homogeneous Lattices in One and Two Dimen- [27] M. Taddei, M. Ruggeri, S. Moroni, and M. Holzmann, sions, Phys. Rev. Lett. 110, 205301 (2013). Iterative Backflow Renormalization Procedure for Many- [47] B. Fang, G. Carleo, A. Johnson, and I. Bouchoule, Quench- Body Ground-State Wave Functions of Strongly Interacting Induced Breathing Mode of One-Dimensional Bose Gases, Normal Fermi Liquids, Phys. Rev. B 91, 115106 Phys. Rev. Lett. 113, 035301 (2014). [48] G. Boéris et al., Mott Transition for Strongly Interacting (2015). [28] G. Carleo, F. Becca, M. Schiro, and M. Fabrizio, Locali- One-Dimensional Bosons in a Shallow Periodic Potential, zation and Glassy Dynamics Of Many-Body Quantum Phys. Rev. A 93, 011601 (2016). Systems., Sci. Rep. 2, 243 (2012). [49] S. Sorella, Generalized Lanczos Algorithm for Variational [29] G. Carleo, F. Becca, L. Sanchez-Palencia, S. Sorella, and M. Quantum Monte Carlo, Phys. Rev. B 64, 024512 (2001). Fabrizio, Light-Cone Effect and Supersonic Correlations in [50] C. J. Umrigar and C. Filippi, Energy and Variance Opti- One- and Two-Dimensional Bosonic Superfluids, Phys. mization of Many-Body Wave Functions, Phys. Rev. Lett. Rev. A 89, 031602 (2014). 94, 150201 (2005). [30] L. Cevolani, G. Carleo, and L. Sanchez-Palencia, Protected [51] M. Girardeau, Relationship between Systems of Impen- Quasilocality in Quantum Systems with Long-Range Inter- etrable Bosons and Fermions in One Dimension, J. Math. actions, Phys. Rev. A 92, 041603 (2015). Phys. (N.Y.) 1, 516 (1960). [31] B. Blaß and H. Rieger, Test of Quantum Thermalization in [52] G. E. Astrakharchik and S. Giorgini, Correlation Functions the Two-Dimensional Transverse-Field Ising Model, Sci. and Momentum Distribution of One-Dimensional Bose Rep. 6, 38185 (2016). Systems, Phys. Rev. A 68, 031602 (2003). 031026-11 GIUSEPPE CARLEO et al. PHYS. REV. X 7, 031026 (2017) [53] P. Pippan, S. R. White, and H. G. Evertz, Efficient Matrix- [59] J. C. Zill, T. M. Wright, K. V. Kheruntsyan, T. Gasenzer, and Product State Method for Periodic Boundary Conditions, M. J. Davis, A Coordinate Bethe Ansatz Approach to the Phys. Rev. B 81, 081103 (2010). Calculation of Equilibrium and Nonequilibrium Correla- [54] G. Carleo, G. Boéris, M. Holzmann, and L. Sanchez- tions of the One-Dimensional Bose Gas, New J. Phys. 18, Palencia, Universal Superfluid Transition and Transport 045010 (2016). Properties of Two-Dimensional Dirty Bosons, Phys. Rev. [60] C. N. Yang and C. P. Yang, Thermodynamics of a Lett. 111, 050406 (2013). One-Dimensional System of Bosons with Repulsive Delta- [55] M. Boninsegni, N. Prokof’ev, and B. Svistunov, Worm Function Interaction, J. Math. Phys. (N.Y.) 10, 1115 Algorithm for Continuous-Space Path Integral Monte Carlo (1969). Simulations, Phys. Rev. Lett. 96, 070601 (2006). [61] M. Rigol, Fundamental Asymmetry in Quenches between [56] J. C. Zill, T. M. Wright, K. V. Kheruntsyan, T. Gasenzer, Integrable and Nonintegrable Systems, Phys. Rev. Lett. and M. J. Davis, Relaxation Dynamics of the Lieb-Liniger 116, 100601 (2016). Gas Following an Interaction Quench: A Coordinate [62] M. Kollar and M. Eckstein, Relaxation of a One- Bethe-Ansatz Analysis, Phys. Rev. A 91, 023611 (2015). Dimensional Mott Insulator after an Interaction Quench, [57] J. De Nardis, L. Piroli, and J.-S. Caux, Relaxation Dynamics Phys. Rev. A 78, 013626 (2008). of Local Observables in Integrable Systems., J. Phys. A 48, [63] M. Holzmann, D. M. Ceperley, C. Pierleoni, and K. Esler, 43FT01 (2015). Backflow Correlations for the Electron Gas and Metallic [58] J.-S. Caux and F. H. L. Essler, Time Evolution of Local Hydrogen, Phys. Rev. E 68, 046707 (2003). Observables After Quenching to an Integrable Model, Phys. Rev. Lett. 110, 257203 (2013). 031026-12

Journal

Physical Review XAmerican Physical Society (APS)

Published: Jul 1, 2017

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

Monthly Plan

  • Read unlimited articles
  • Personalized recommendations
  • No expiration
  • Print 20 pages per month
  • 20% off on PDF purchases
  • Organize your research
  • Get updates on your journals and topic searches

$49/month

Start Free Trial

14-day Free Trial

Best Deal — 39% off

Annual Plan

  • All the features of the Professional Plan, but for 39% off!
  • Billed annually
  • No expiration
  • For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588

$360/year

billed annually
Start Free Trial

14-day Free Trial