Solving the Quantum Many-Body Problem with Bohmian Trajectories

Bohmian mechanics is an interpretation of quantum mechanics that describes the motion of quantum particles with an ensemble of deterministic trajectories. Several attempts have been made to utilize Bohmian trajectories as a computational tool to simulate quantum systems consisting of many particles, a very demanding computational task. In this paper, we present a novel ab-initio approach to solve the many-body problem for bosonic systems by evolving a system of one-particle wavefunctions representing pilot waves that guide the Bohmian trajectories of the quantum particles. The method is used to study the breathing dynamics and ground state properties in a system of several interacting bosons. The precision of our method is benchmarked against the numerically exact multiconfigurational time-dependent Hartree method for bosons.

Bohmian mechanics is an interpretation of quantum mechanics that describes the motion of quantum particles with an ensemble of deterministic trajectories.Several attempts have been made to utilize Bohmian trajectories as a computational tool to simulate quantum systems consisting of many particles, a very demanding computational task.In this paper, we present a novel ab-initio approach to solve the many-body problem for bosonic systems by evolving a system of one-particle wavefunctions representing pilot waves that guide the Bohmian trajectories of the quantum particles.The method is used to study the breathing dynamics and ground state properties in a system of several interacting bosons.The precision of our method is benchmarked against the numerically exact multiconfigurational time-dependent Hartree method for bosons.
Numerical simulation of the quantum dynamics of many-body systems is plagued by the dimension of the Hilbert space which increases exponentially with the number of particles.Much of the progress in theoretical condensed matter, atomic and molecular physics in the past few decades has been achieved by finding new ways to circumvent this problem.Some of the most powerful approaches are density functional theory [1], quantum Monte Carlo [2], density matrix renormalization group [3] and the multi-configuration time-dependent Hartree method [4].Recently, nonconventional approaches based on artificial intelligence [5], Bohmian mechanics [6][7][8][9] and wavelet transforms [10] have been proposed as well.Several methods aim to reduce the complexity of the numerical simulation of many-body systems by resorting to low dimensional objects such as density functions [1,11] and natural orbitals [12] or by using mixed classical-quantum dynamics such as the Ehrenfest approach [13] and the surface hopping method [14,15].In this paper, we use another class of low-dimensional objects, namely singleparticle pilot waves evolved concurrently with Bohmian trajectories to extract all the physical information about the system.Within the de Broglie-Bohm interpretation of quantum mechanics [16][17][18], the quantum mechanical wavefunction is a pilot wave that guides the motion of the particle in the physical space.While this interpretation does not alleviate the need for dealing with manydimensional functions, the prospect of replacing the full many-particle wavefunction by single-particle wavefunctions that guide the Bohmian particles in the physical space was recently explored [19].However, this idea was only explored when the entanglement between the particles could be neglected [19,20], thus ruling out its application to strongly correlated systems.Other approaches to treat many-body wavefunctions with trajectories involve other approximations such as the mean-field ap-proximation or assume a certain form for the wavefunction [21,22].
Here, we introduce a novel approach that simulates multi-particle quantum dynamics in a non-perturbative manner using pilot waves without neglecting the entanglement or relying on particular assumptions about the underlying quantum state.We apply our approach to study the breathing dynamics of many-boson systems in a trap with long-and short-range interactions.We use the same method to compute the ground state energies for an exactly solvable system.We compare our results with the multiconfigurational time-dependent Hartree method for bosons (MCTDHB), which has become a standard method to analyze systems of many interacting bosons beyond the mean-field method [23,24].

Entangled quantum dynamics using Bohmian trajectories
Let us illustrate the usage of pilot waves in a 1D system consisting of 2 particles.The coordinates are denoted by x 1 and x 2 , the potential by V (x 1 , x 2 ) and the wavefunction describing the full system is Ψ(x 1 , x 2 , t).In order to evolve the Bohmian trajectories X 1 (t) and X 2 (t) for the two particles (we denote the Bohmian trajectories throughout this paper by uppercase letters), we need to evaluate the pilot waves ψ i (x i , t).The pilot waves are the full wavefunction projected on the coordinates of all the particles except one, i.e., ψ i (x i , t) ≡ Ψ(x 1 , x 2 , t)| xj =Xj (t), j =i ; hence, they are also called conditional wavefunctions (CWs) [25].In the absence of gauge fields, the Bohmian velocities are computed in terms of the pilot waves as where m i is the mass of particle i.
It is guaranteed that the density of the Bohmian particles evolved by Eq. ( 1) follows the evolution of the density function as computed by Schrödinger's equation [8].In order to evolve the CWs without having to solve the time-dependent Schrödinger equation (TDSE) for Ψ(x 1 , x 2 , t), we introduce a generalized set of conditional The trajectories of the light particle computed from the exact pilot waves (solid blue) and using pilot waves evolved using Eq.(3) (dashed).The case for N = 7 corresponds to interacting pilot waves evolved with the help of a hierarchy of conditional wavefunctions {ψ n 1 , ψ n 2 } up to n = 7 (see Eq.( 3)), while N = 0 corresponds to noninterating pilot waves (the small entanglement approximation).We notice that the former case is more accurate than the latter.The two particles have mass ratio 1:100.ω is the trap frequency of the light particle.
wavefunctions ψ n i (x i , t) defined as where the pilot waves correspond to ψ 0 i .It can be easily proved (see Methods) that the equation of motion for We see from this equation that the pilot waves corresponding to different particles interact indirectly through the last three terms of Eq. (3).In [19], a similar equation of motion is derived in terms of nonlocal potentials.The price to be paid in order to evolve the interacting pilot waves using Eq. ( 3) in an exact manner instead of evolving the multi-dimensional full wavefunction is that we need to evolve the whole hierarchy of {ψ n i }.One can truncate the hierarchy at a finite order N .The higher the value of N , the longer it takes till the truncation error propagates from ψ N i down to ψ 0 i and affects the evolution of X i (t) and the precision of the results.Since the errors afflict the pilot waves through the last two terms in Eq. (3), this method of evolving the pilot waves is most suitable when we are interested in the dynamics of a very light particle interacting with a much heavier one over a very short time scale.In this case, we can omit the last two terms for the heavy particle while retaining them for the light particle.
To illustrate this method, we consider the entan-gled dynamics of two particles of masses m 1 = 1 and m 2 = 100 subject to the harmonic potential with k = 0.1 (atomic units are used in the rest of this paper).Let us take an initial entangled state (see Methods) and evolve the Bohmian trajectories for the initial conditions X 1 = 1, X 2 = 2.We first truncate the hierarchy at N = 0, thus making Eq.( 3) unitary.This case corresponds to the small entanglement approximation, i.e., noninterating pilot waves.Fig. 1 shows that the Bohmian trajectory evolved by the corresponding pilot wave deviates from the trajectory computed from the exact pilot wave already at half a cycle of the oscillatory motion.Increasing the depth of the hierarchy to N = 7 only extends the range of accurate dynamics for another cycle.
In this example, we see clearly that although the particles are non-interacting, we need to account for the interaction between the pilot waves of the two particles correctly through the higher order CWs even when the ratio of the particles' masses is 1:100.Otherwise, the errors originating from truncating the hierarchy of {ψ n i } propagate very fast to ψ 0 i .Increasing N beyond N ≈ 10 will not help in prolonging the range of accuracy because of the numerical errors in the calculation of the higherorder derivatives of the wavefunction.Is there a more accurate way to propagate the pilot waves to avoid the errors originating from the truncated orders?It turns out, as we will show shortly, that the answer is yes!To this end, we assume an ansatz for the full wavefunction that allows the calculation of the first and second CWs ψ 1 i and ψ 2 i which we subsequently use to evolve the interacting pilot waves ψ 0 1 and ψ 0 2 .In the rest of the paper, we consider only identical particles.
A novel approach to propagate the pilot waves The most general form for the wavefunction of a 2particle system is where {φ i } is a complete basis for the one-body Hilbert space.The conditional wavefunctions of the first particle conditioned on the second particle located at y = Y are expressed as (the time variable and the particle index are dropped to simplify the notation) where ∂y 2 | y=Y .These relations can be written in vector form as: The problem of finding ψ 1 and ψ 2 boils down to finding the coefficients b i and c i constituting the vectors b and c.This is accomplished by making use of an ensemble of Bohmian pairs of coordinates {X, Y } which are selected initially from the one-particle density function ρ(x) at t = 0.Each of these pairs is called a configuration.If we can represent both φ and φ for a certain value of Y as a linear superposition of all { φ} corresponding to all members of the ensemble, i.e.,if φ = k α k φ k and φ = k β k φ k then it follows from the linearity in Eq. ( 5) that b = k α k a k and c = k β k a k .Finding the values of α k and β k is equivalent to solving a system of linear equations.In this way, we can predict ψ 1 , ψ 2 without ever constructing the coefficient matrix C.
It should be noted that ψ 1 and ψ 2 can be obtained without expressing ψ 0 in terms of a basis at all, since {α k } and {β k } depend only on the amplitudes of any complete basis at the location of the Bohmian particles.After {α k } and {β k } are obtained, we can express ψ 1 and ψ 2 as With the CWs at our disposal, we use the equation of motion (3) for n = 0 to evolve the ensemble of CWs for all Bohmian particles as described in the Methods Section.We call this scheme Interacting Pilot Waves (IPW).
Some observables can be computed by averaging over the ensemble of Bohmian configurations {X, Y } such as x 2 .Since we have access to the CWs, we can devise a more accurate method that approximates the exact expression of the expectation value of an operator Â, Â = Ψ * (x, y) ÂΨ(x, y)dxdy by performing the integral over one variable as a Riemann sum over its Bohmian coordinates, i.e., where ψ w (x) is the conditional wavefunction of the first particle conditioned on the coordinate of the second particle belonging to the w th configuration of the ensemble, ∆ w is the distance between adjacent values of Y at the w th configuration and Âw is the operator Â conditioned on X 2w .Similarly, the reduced one-body density can be approximated as ρ(x) ≈ w ∆ w ψ * w (x)ψ w (x).Let us apply this method to study the breathing dynamics of two bosons in a harmonic trap, V (x) = 1 2 k t x 2 .The bosons are initially condensed in the ground state of the trap, and start a breathing motion when a harmonic interaction V (x, y) = 1 2 k i (x−y) 2 is suddenly switched on.A finite fixed set of orbitals are taken to be the lowest set of eigenfunctions of the one-body problem with the effective potential felt by one particle due to the other one, namely V eff (x) = 1 2 k t x 2 + ρ(y)V (x, y)dy.We illustrate in Fig. 2 the behavior of ρ(0) for k t = k i = 1 computed by the IPW method with 6 orbitals compared with the exact dynamics (i.e., from the exact pilot waves obtained by solving the TDSE exactly).
Generalization to many-particle systems Generalizing our algorithm to a many-particle problem consisting of N B bosons is straightforward.Let us denote the coordinates of the particles by x, y, z, • • • etc., while, as before, we denote the Bohmian coordinates by upper case letters.A single configuration of Bohmian walkers is denoted by (X, Y, Z, • • • ).
Let us denote the conditional wavfunctions by

FIG. 2:
The breathing dynamics of two interacting bosons in a harmonic trap.The two bosons are initially in the ground state of the trap before the harmonic interaction between them is suddenly switched on at t = 0.The value of the reduced one-body density function ρ(x) at the origin is monitored during the evolution.The strength of the harmonic interaction is the same as the trap strength.We plot ρ(0) computed by the Interacting Pilot Waves method (IPW) using 6 orbitals and by the exact two-body wavefunction (Exact).The IPW results were obtained by averaging over 5000 Bohmian configurations of the two particles.The two panels depict different time slices to illustrate the accuracy in the long-time regime.
and so on.The equation of motion for ψ 0 (x; Y, Z, • • • ) is a simple generalization of Eq. ( 3) for the case of many particles: Similar equations can be written for all the CWs corresponding to all particles in every configuration.A generic ansatz for the many-body wavefunction similar to Eq. ( 4) reads Ψ(x, y, z, ..., t) = i,j,k c ijk... (t)φ i (x)φ j (y)φ k (z).... (8) In order to compute ψ 1 (x; Y , Z, This can be done with the existing numerical techniques by rearranging all the M N B −1 terms of the tensors, where M is the number of orbitals, in vector forms and solving a linear system of equations. In order to compute the expectation value of an operator, we can treat the collection of CWs belonging to all the particles as if they describe normal single particle wavefunctions (after renormalizing them).Âx can then be computed as Âx ≈ 1 Nw w ψ * w (x) Âx ψw (x)dx, where ψw (x) is the normalized CW of particle x belonging to the w th configuration.If Â is a two-body operator such as the interaction potential between two particles, we first compute a mean-field operator, and then treat it as a one-body operator.For example, the mean-field interaction potential felt by one particle is computed as Ṽ (x) ≈ 1 Nw w ψ * w (y)V (x − y) ψw (y)dy.The expectation value of Ṽ (x) is then computed as a one-body operator.
From the normalized collection of all the CWs, we can also compute the reduced density matrix of one particle ρ(x , x) ≈ 1 Nw w ψ * w (x ) ψw (x).This matrix can be used to compute a set of natural orbitals in terms of the finite basis set used in the postulated ansatz as the eigenstates of ρ(x , x).
Many particles in a harmonic trap We apply this generalization to the dynamics of 3 bosons and 5 bosons in a harmonic trap V (x) = 1 2 x 2 for two cases of interparticle interactions: long-range attractive harmonic interaction and short-range repulsive interaction.As in the two-particle case, all the bosons (a,b) 5 bosons and 3 bosons are initially in the ground state of a harmonic trap before a harmonic interaction of strength ki = 0.1 and ki = 1 respectively is switched on at t = 0.The breathing dynamics is computed by the Interacting Pilot Wave (IPW) method with 4 and 3 orbitals in the two cases respectively and compared with the numerically exact dynamics computed by MCTDHB and with the Hermitian limit (HL) solution of Eq. ( 7) (see text).In (c) we depict the time evolution for an ensemble of 1000 Bohmian trajectories representing the first particle, and in (d) the 5 trajectories corresponding to a single configuration for the case (a).
initially reside in the ground state of the harmonic trap before the interaction is suddenly switched on at t = 0.
We study the breathing dynamics by computing x 2 as a function of time.
For a harmonic interaction of the form V (x, y) = 1 2 k i (x − y) 2 we consider two cases of 5 bosons with weak interactions (k i = 0.1) and 3 bosons with strong interaction (k i = 1) and we use 3 and 4 orbitals in the two cases, respectively.In both cases we compare the results with the numerically exact simulation using MCT-DHB [23,[26][27][28][29] and with the Hermitian limit (HL) of Eq. ( 7) where all the non-hermitian terms in Eq. ( 7) are dropped out.The Hermitian limit is equivalent to the time-dependent quantum Monte-Carlo (TDQMC) of Ref. [22] which does not take entanglement into consideration.
In Fig. 3, we show the results of computing x 2 by av-eraging over the Bohmian coordinates of all the particles using a single ensemble containing 1000 configurations.We notice that the IPW method is more accurate in the weak interaction regime than in the strong interaction regime.The Bohmian trajectories of the first particle in all configurations are show in Fig. 3(c) for k i = 0.1 while the Bohmian trajectories for all the 5 particles in a single configuration are shown in Fig. 3(d).The few constant trajectories appearing in Fig. 3(c) correspond to the cases where we manually set the Bohmian velocities to be zero when the denominator is Eq. ( 1) is below a certain threshold.
In Fig. 4, we plot x 2 after switching on a gaussian with k i = 0.1, σ = 0.25 and compare the results with MCTDHB simulation and the HL of Eq. (7).In this calculation, we with ki = 0.1, σ = 0.25.Three orbitals are used in the IPW calculation.
2.5 Ground state energy for a system of 5 interacting bosons in a harmonic trap with harmonic interaction strength ki = 0.1 computed by the IPW method (solid line) using 3 orbitals and compared with the exact ground state energy (dashed).The system is initialized in the ground state of the trap before the interaction is switched on adiabatically as ki = 0.1 × (1 − e −0.02t 2 ) .The energy is computed with respect to the instantaneous value of ki.compute x 2 from the expectation value of x 2 using the conditional wavefunctions rather than from the Bohmian trajectories.
In order to compute the ground state energy for an interacting system of particles using the IPW scheme, we initialize the CWs and the Bohmian trajectories in the ground state of the noninteracting Hamiltonian.Afterwards, we switch on the interaction adiabatically.According to the adiabatic theory [30], the system remains in the ground state of the instantaneous Hamiltonian.
In Fig. 5, we plot the evolution of the energy of the instantaneous Hamiltonian of a 5-particle system as we switch on the harmonic interaction V (x, y) = 1 2 k i (x − y) 2 adiabatically and compare it with the exact ground state energy E 0 = N B −1 2 √ 1 + k i N B +0.5 [31] for k i = 0.1.The ground state energy computed by MCTDHB [32] is more accurate than the IPW calculation for the same number of orbitals by several significant digits.Perhaps a better method to compute ground state energy is to propagate Eq. ( 7) in complex time, while evolving the Bohmian trajectories in real time [22].The optimal relation between the real and complex time evolution constitutes an interesting topic of research.

Conclusion and discussion
We have presented a promising approach to analyze the dynamic and static properties of systems consisting of several bosons by evolving a system of nonunitary equations that goes beyond the small entanglement approximation and the mean-field approximation.
The accuracy of this new approach is outperformed by the state-of-the-art MCTDHB algorithm.In the MCT-DHB method, increasing the number of orbitals (M ) is confronted with the exponentially large number of configurations of permanents that needs to be taken into account.We have a similar scaling problem in our approach; the number of configurations of Bohmian particles has to be larger than M N B −1 in order to avoid having an undetermined linear system of equations when solving for α w .So, the complexity of our approach still increases exponentially with the number of particles but with a much smaller base for the exponent than the barebones TDSE approach which scales as L N B , where L is the number of grid points per dimension.
It is still an open question whether the non-Hermitian terms in the equations of motion can be replaced by an effective entanglement potential that makes the equations unitary and at the same time captures the entanglement in the system.For systems consisting of many particles, i.e., N B > 20, the entanglement of the ground state is so small [33] that even the Hermitian limit [22] can be efficient for simulating the dynamics involving a small number of excited states.
In principle, generalizing the IPW approach to fermionic systems is straightforward, as long as we choose the initial state with the proper symmetry requirements.However, for fermions two problems arise.First, due to the Pauli exclusion principle, we need a large number of orbitals to describe a fermionic state and therefore, the number of fermions that can be analyzed is small compared to bosons.Second, the conditional wavefunctions for fermions will have nodes, that complicate computing the velocity of the Bohmian walkers around those nodes.Since the node problem is a well known problem for simulating quantum dynamics with Bohmian trajectories [6], the methods developed in this regard in the literature [34,35]  Author contributions TAE conceived the research idea, derived Eq. ( 3), designed the IPW algorithm for two and many particles and implemented it on the computer.LBM and KM provided helpful insights and feedback during the development of the method.All authors participated in the data analysis and the troubleshooting.The manuscript was written by TAE with extensive feedback from LBM and KM.
Competing financial interests The authors declare no competing financial interests.

Derivation of Equation (3)
Let us derive Eq. ( 3) with respect to ψ n , we find by the chain rule that By exchanging the time and spatial derivatives in the first term on the R.H.S. and using the TDSE: we obtain after substituting back in Eq. ( 9) By applying the chain rule to the second term on the R.H.S. we obtain which after substituting in Eq. ( 10) recovers Eq. ( 3).

Initial entangled state
We take the initial entangled state in Fig. 1 for a two-particle system as the ground state of the Hamiltonian with the potential function V (x 1 , x 2 ) = 1 2 k 1 x 2 1 + 1 2 k 2 x 2 2 + 1 2 k 3 (x 1 − x 2 ) 2 with k 1 = k 2 = 0.1, k 3 = 1.0 and the masses of the particle m 1 = 1 and m 2 = 2.
Interacting Pilot Waves for two particles For the two-particle IPW simulations in Fig. 2 we express all the CWs in terms of the basis set {φ i }, and evolve the expansion coefficients {a i }.Equation 3 is then expressed as By taking the inner product with each of the orbitals {φ i (x)} we obtain the time derivative of the expansion coefficients { ȧi }.This system of equations is then solved using a fourth-order Runge-Kutta method.
Propagating the conditional wavefunctions Equations ( 3) and (7) represent a system of coupled nonlinear and nonunitary differential equations that needs to be solved together.Each of these equations can be cast in the form i ∂ψ(x, t) ∂t = Hψ(x, t) + W (x, t), where the first term on the RHS represents the unitary part of the equation and W (x, t) represents the nonunitary part which is a function of all other CWs.If, e.g., H is a constant Hamiltonian, a general solution for this equation takes the form: ψ(x, t) = e −iHt t 0 e iHt W (x, t )dt + ψ(x, t 0 ) In order to propagate ψ(x, t) for a single time step from t = 0 to t = δt using this solution, both the operator e −iHδt and e −iHt are performed using a split-operator method [36].The integral is performed using the trapezoidal rule δt 0 e iHt W (x, t )dt ∼ 1 2 e iHδt W (x, δt) + W (x, 0) .

Miscellaneous numerical techniques
• The evaluation of the conditional wavefunctions at the position of the Bohmian particles was performed by FFT-based interpolation.
• Solving the linear system of equations in Eq. ( 5) was performed by the LAPACK routine gelsd() which uses singular value decomposition and a divide and conquer method to compute the minimum-norm solution to a linear least squares problem [37].

FIG. 1 :
FIG.1: Entangled quantum dynamics of two particles in a harmonic trap.(a) A cartoon representing the pilot waves guiding Bohmian particles moving in 2D.Although the two particles do not interact, their entanglement implies a coupling between the two pilot waves guiding the two Bohmian trajectories (the dashed line).(b) The trajectories of the light particle computed from the exact pilot waves (solid blue) and using pilot waves evolved using Eq.(3) (dashed).The case for N = 7 corresponds to interacting pilot waves evolved with the help of a hierarchy of conditional wavefunctions {ψ n 1 , ψ n 2 } up to n = 7 (see Eq.(3)), while N = 0 corresponds to noninterating pilot waves (the small entanglement approximation).We notice that the former case is more accurate than the latter.The two particles have mass ratio 1:100.ω is the trap frequency of the light particle.

FIG. 3 :
FIG.3: Breathing dynamics of 5 bosons and 3 bosons using IPW compared with numerically exact dynamics.(a,b) 5 bosons and 3 bosons are initially in the ground state of a harmonic trap before a harmonic interaction of strength ki = 0.1 and ki = 1 respectively is switched on at t = 0.The breathing dynamics is computed by the Interacting Pilot Wave (IPW) method with 4 and 3 orbitals in the two cases respectively and compared with the numerically exact dynamics computed by MCTDHB and with the Hermitian limit (HL) solution of Eq.(7) (see text).In (c) we depict the time evolution for an ensemble of 1000 Bohmian trajectories representing the first particle, and in (d) the 5 trajectories corresponding to a single configuration for the case (a).

FIG. 5 :
FIG.5: Ground state energy computation using IPW.Ground state energy for a system of 5 interacting bosons in a harmonic trap with harmonic interaction strength ki = 0.1 computed by the IPW method (solid line) using 3 orbitals and compared with the exact ground state energy (dashed).The system is initialized in the ground state of the trap before the interaction is switched on adiabatically as ki = 0.1 × (1 − e −0.02t 2) .The energy is computed with respect to the instantaneous value of ki.

∂ 2
φ i (x) ∂x 2 + V (x, Y ) i a i (t)φ i (x) − 2 2m i c i (t)φ i (x) + i dY dt i b i (t)φ i (x).(11) may benefit the solution of this problem.Acknowledgments T. A. Elsayed thanks D. A. Deckert, H. Miyagi, A. I. Streltsov, L. F. Buchmann, C. Leveque and Andreas Caranti for discussions and I. Christov for his comments on an earlier version of the manuscript.The authors acknowledge financial support by Villum Foundation.