Entangled Quantum Dynamics of Many-Body Systems using Bohmian Trajectories

TL;DR

This paper introduces a new ab-initio method using Bohmian trajectories to simulate many-body bosonic systems, capturing quantum entanglement effects efficiently.

Contribution

It presents a novel approach that evolves one-particle wavefunctions to model entanglement in many-body quantum systems.

Findings

Successfully modeled breathing dynamics of interacting bosons

Accurately predicted ground state properties

Demonstrated computational efficiency over traditional methods

Abstract

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. In this approach, quantum entanglement effects arise due to the interactions between different configurations of Bohmian particles evolving simultaneously. The method is used to study the breathing dynamics and ground state properties in a system of interacting bosons.

Equations

$$\frac{d{X}_i}{dt}={\frac{\hslash }{m_i}\text{Im}\left\{\frac{{\partial }_{x_i}{\psi }_i(x_i,t)}{{\psi }_i(x_i,t)}\right\}|}_{x_i=X_i(t)},$$

$${\psi }_i^n(x_i,t)\equiv {\frac{{\partial }^n\Psi (x_1,x_2,t)}{\partial x_j^n}|}_{x_j=X_j(t),j\ne i}$$

$$i\hslash \frac{\partial {\psi }_i^n(x_i,t)}{\partial t}=-\frac{{\hslash }^2}{2m_i}\frac{{\partial }^2{\psi }_i^n(x_i,t)}{\partial x_i^2}+{\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right){\psi }_i^{n-k}(x_i,t)\frac{{\partial }^{k}V(x_i,x_j)}{\partial x_j^k}|}_{x_j=X_j(t),j\ne i}-\frac{{\hslash }^2}{2m_j}{\psi }_i^{n+2}(x_i,t)+i\hslash \frac{d{X}_j(t)}{dt}{\psi }_i^{n+1}(x_i,t)$$

$$\Psi (x,y,t)=\sum _{i,j}{c}_{ij}(t){\varphi }_i(x){\varphi }_j(y),$$

$${\psi }^0(x)\equiv {\Psi (x,y)|}_{y=Y}=\sum _i{a}_i{\varphi }_i(x)$$

$${\psi }^1(x)\equiv {\frac{\partial \Psi (x,y)}{\partial y}|}_{y=Y}=\sum _i{b}_i{\varphi }_i(x)$$

$${\psi }^2(x)\equiv {\frac{{\partial }^2\Psi (x,y)}{\partial y^2}|}_{y=Y}=\sum _i{c}_i{\varphi }_i(x)$$

$$\vec{\mathbf{a}}(Y)=\mathbf{C}\vec{\mathbf{\varphi }}(Y),\vec{\mathbf{b}}(Y)={\mathbf{C}\vec{\mathbf{\varphi }}}^{\prime }(Y)\text{and}\vec{\mathbf{c}}(Y)={\mathbf{C}\vec{\mathbf{\varphi }}}^{\prime \prime }(Y)$$

$${\psi }^1(x)=\sum _k{\alpha }_k{\psi }_k^0(x;Y_k)$$

$${\psi }^2(x)=\sum _k{\beta }_k{\psi }_k^0(x;Y_k)$$

$$⟨\hat{A}⟩\approx \sum _w{\Delta }_w\int {\psi }_w^{\ast }(x){\hat{A}}_w{\psi }_w(x)dx,$$

$${\psi }^0(x;Y,Z,\cdots )\equiv {\Psi (x,y,z,\cdots )|}_{y=Y,z=Z,\cdots }$$

$${\psi }^1(x;Y^{\prime },Z,\cdots )\equiv {\frac{\partial \Psi (x,y,z,\cdots )}{\partial y}|}_{y=Y,z=Z,\cdots }$$

$${\psi }^1(z;X,Y^{\prime },\cdots )\equiv {\frac{\partial \Psi (x,y,z,\cdots )}{\partial y}|}_{x=X,y=Y,\cdots },$$

$$i\hslash {\partial }_t{\psi }^0(x;Y,Z,\cdots ;t)=$$

$$\left(-\frac{{\hslash }^2}{2m}{\nabla }_x^2+V(x;Y,Z,\cdots )\right){\psi }^0(x;Y,Z,\cdots ;t)$$

$$+i\hslash \frac{dY}{dt}{\psi }^1(x;Y^{\prime },Z,\cdots ;t)-\frac{{\hslash }^2}{2m}{\psi }^2(x;Y^{\prime \prime },Z,\cdots ;t)$$

$$+i\hslash \frac{dZ}{dt}{\psi }^1(x;Y,Z^{\prime },\cdots ;t)-\frac{{\hslash }^2}{2m}{\psi }^2(x;Y,Z^{\prime \prime },\cdots ;t)+$$

$$+\cdots \text{etc.}$$

$$\Psi (x,y,z,...,t)=\sum _{i,j,k}{c}_{ijk...}(t){\varphi }_i(x){\varphi }_j(y){\varphi }_k(z)z...$$

$$\frac{\partial {\psi }_1^n(x_1,t)}{\partial t}=\frac{\partial }{\partial t}{\frac{{\partial }^n\Psi (x_1,x_2)}{\partial x_2^n}|}_{x_2=X_2(t)}+\frac{d{X}_2(t)}{dt}{\frac{{\partial }^n\Psi (x_1,x_2)}{\partial x_2^n}|}_{x_2=X_2(t)}.$$

$$i\hslash \frac{\partial \Psi (x_1,x_2)}{\partial t}=\left(-\frac{{\hslash }^2{\nabla }_1^2}{2m_1}-\frac{{\hslash }^2{\nabla }_2^2}{2m_2}+V(x_1,x_2)\right)\Psi (x_1,x_2),$$

$$i\hslash \frac{\partial {\psi }_1^n(x_1,t)}{\partial t}=-\frac{{\hslash }^2}{2m_1}\frac{{\partial }^2{\psi }_1^n(x_1,t)}{\partial x_1^2}+\frac{{\partial }^n}{\partial x_2^n}{\left(V(x_1,x_2)\Psi (x_1,x_2)\right)|}_{x_2=X_2(t)}-\frac{{\hslash }^2}{2m_2}{\psi }_1^{n+2}(x_1,t)+i\hslash \frac{d{X}_2(t)}{dt}{\psi }_1^{n+1}(x_1,t).$$

$$i\hslash \sum _i{\dot{a}}_i(t){\varphi }_i(x)=-\frac{{\hslash }^2}{2m}\sum _i{a}_i(t)\frac{{\partial }^2{\varphi }_i(x)}{\partial x^2}+V(x,Y)\sum _i{a}_i(t){\varphi }_i(x)-\frac{{\hslash }^2}{2m}\sum _i{c}_i(t){\varphi }_i(x)+i\hslash \frac{dY}{dt}\sum _i{b}_i(t){\varphi }_i(x).$$

$$i\frac{\partial \psi (x,t)}{\partial t}=H\psi (x,t)+W(x,t),$$

$$\psi (x,t)=e^{-iHt}\left[{\int }_0^t{e}^{iH{t}^{\prime }}W(x,t^{\prime })d{t}^{\prime }+\psi (x,t_0)\right]$$

Full text

Introduction

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 kohn1965 , quantum Monte Carlo foulkes2001 , density matrix renormalization group white1992 and the multi-configuration time-dependent Hartree method beck2000 . Recently, nonconventional approaches based on machine learning Carleo2017 , Bohmian mechanics wyatt2006 ; lopreore1999 ; benseny2014 ; hall2014 and wavelet transforms poirier2003 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 kohn1965 ; parr1980 and natural orbitals bauer2014 or by using mixed classical-quantum dynamics such as the Ehrenfest approach meyera1979 and the surface hopping method tully1971 ; tully1990 . In this paper, we use another class of low-dimensional objects, namely single-particle 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 durr2012 ; durr2009 ; holland1995 , 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 many-dimensional 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 norsen2015 . However, this idea was only applied when the entanglement between the particles could be neglected norsen2015 ; oriols2007 , thus ruling out its application to strongly correlated systems. Other approaches to treat many-body wavefunctions with trajectories involve approximations such as the mean-field approximation christov2017 ; christov2006 or the semi-classical approaches norsen2016 ; struyve2015 or assume a wavefunction of a certain form christov2007 .

In the de Broglie-Bohm interpretation, quantum effects are captured by a so-called quantum potential which together with the classical potential governs the motion of the particle, see, e.g., lopreore1999 . It has been shown recently that this term can be computed by modeling quantum phenomena by many interacting classical worlds hall2014 ; herrmann2017 . In this paper, we introduce a novel approach to model quantum phenomena using interacting configurations of quantum particles guided by pilot waves. This approach simulates the multi-particle quantum dynamics in a non-perturbative manner without neglecting the entanglement or relying on particular assumptions about the underlying quantum state. We apply our approach to study the breathing dynamics of few-boson systems in a trap with long- and short-range interactions and compute the ground state energy for an exactly solvable system. The proposed approach in its current stage does not supersede established numerical methods nor overcomes the scaling problem of simulating-many body systems, but offers a new way forward that may be further developed into a full-fledged method.

Results

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 $\Psi(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 $\psi_{i}(x_{i},t)$. The pilot waves are the full wavefunction projected on the coordinates of all the particles except one, i.e., $\psi_{i}(x_{i},t)\equiv\Psi(x_{1},x_{2},t)|_{x_{j}=X_{j}(t),\ j\neq i}$; hence, they are also called conditional wavefunctions (CWs) durr1992 . In the absence of gauge fields, the Bohmian velocities are computed in terms of the pilot waves as

[TABLE]

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 benseny2014 . In order to evolve the CWs without having to solve the time-dependent Schrödinger equation (TDSE) for $\Psi(x_{1},x_{2},t)$, we introduce a generalized set of conditional wavefunctions $\psi_{i}^{n}(x_{i},t)$ defined as

[TABLE]

where the pilot waves correspond to $\psi_{i}^{0}$. The equation of motion for $\psi_{i}^{n}(x_{i},t)$ is given by norsen2010

[TABLE]

We see from this equation that the pilot waves corresponding to different particles interact indirectly through the last three terms of Eq. (3). In norsen2015 , a similar equation of motion is derived in terms of nonlocal potentials.

In order to evolve the pilot waves in an exact manner using Eq. (3) instead of evolving the multi-dimensional full wavefunction, we need to evolve the whole hierarchy of $\{\psi_{i}^{n}\}$. We illustrate in the Methods section that truncating this hierarchy at a finite order $N$ is not an efficient method to obtain the correct dynamics of an entangled system as the truncation errors propagate very quickly to $\psi_{i}^{0}$. Is there a way to avoid the errors originating from the truncated orders? It turns out 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 $\psi_{i}^{1}$ and $\psi_{i}^{2}$ which we subsequently use to evolve the interacting pilot waves $\psi_{i}^{0}$.

Interacting Pilot Waves.

The most general form for the wavefunction of a 2-particle system is

[TABLE]

where $\{\phi_{i}\}$ is a complete basis for the one-body Hilbert space, also referred to as ’orbitals’ later. Let us assume that a finite number of basis states $M$ is sufficient to capture all the important features of the wavefunction. 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 omitted to simplify the notation)

[TABLE]

where $a_{i}=\sum_{j}c_{ij}\phi_{j}(Y)$, $b_{i}=\sum_{j}c_{ij}\frac{\partial\phi_{j}(y)}{\partial y}|_{y=Y}$ and $c_{i}=\sum_{j}c_{ij}\frac{\partial^{2}\phi_{j}(y)}{\partial y^{2}}|_{y=Y}$. These relations can be written in vector form as:

[TABLE]

where $\bm{\vec{a}}(Y)=\{a_{1}(Y),\ a_{2}(Y),...\}$, $\bm{\vec{\phi}}(Y)=\{\phi_{1}(Y),\phi_{2}(Y),...\}$,…etc.

The problem of finding $\psi^{1}$ and $\psi^{2}$ boils down to finding the coefficients $b_{i}$ and $c_{i}$ constituting the vectors $\bm{\vec{b}}$ and $\bm{\vec{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 $\rho(x)$ at $t=0$ note . Each of these pairs is called a configuration. If we can represent both $\bm{\vec{\phi}}^{\prime}$ and $\bm{\vec{\phi}}^{\prime\prime}$ for a certain value of $Y$ as a linear superposition of all $\{\bm{\vec{\phi}}(Y_{k})\}$ corresponding to all members of the ensemble, i.e., if $\bm{\vec{\phi}}^{\prime}(Y)=\sum_{k}\alpha_{k}\bm{\vec{\phi}}(Y_{k})$ where $\bm{\vec{\phi}}(Y_{k})$ corresponds to the $k^{\text{th}}$ member of the ensemble and $\bm{\vec{\phi}}^{\prime\prime}=\sum_{k}\beta_{k}\bm{\vec{\phi}}(Y_{k})$ then it follows from the linearity in Eq. (5) that $\bm{\vec{b}}(Y)=\sum_{k}\alpha_{k}\bm{\vec{a}}(Y_{k})$ and $\bm{\vec{c}}(Y)=\sum_{k}\beta_{k}\bm{\vec{a}}(Y_{k})$. Finding the values of $\alpha_{k}$ and $\beta_{k}$ is equivalent to solving a system of linear equations. In this way, we can obtain $\psi^{1}$, $\psi^{2}$ without ever constructing the coefficient matrix $\bm{C}$.

It should be noted that $\psi^{1}$ and $\psi^{2}$ can be determined without expressing $\psi^{0}$ in terms of a basis at all, since $\{\alpha_{k}\}$ and $\{\beta_{k}\}$ depend only on the amplitudes of any complete basis at the location of the Bohmian particles. After $\{\alpha_{k}\}$ and $\{\beta_{k}\}$ are obtained, we can express $\psi^{1}$ and $\psi^{2}$ as

[TABLE]

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 $\langle x^{2}\rangle$. 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 $\hat{A}$, $\langle\hat{A}\rangle=\int\Psi^{*}(x,y)\hat{A}\Psi(x,y)dxdy$ by performing the integral over one variable as a Riemann sum over its Bohmian coordinates, i.e.,

[TABLE]

where $\psi_{w}(x)$ is the conditional wavefunction of the first particle conditioned on the coordinate of the second particle belonging to the $w^{\text{th}}$ configuration of the ensemble, $\Delta_{w}$ is the distance between adjacent values of $Y$ at the $w^{\text{th}}$ configuration and $\hat{A}_{w}$ is the operator $\hat{A}$ conditioned on $Y_{w}$. For two-body operators such as the interaction potential $V(x,y)$, $\hat{A}_{w}(x)$ is given by $V(x;Y_{w})$. Similarly, the reduced one-body density can be approximated as $\rho(x)\approx\sum_{w}\Delta_{w}\psi_{w}^{*}(x)\psi_{w}(x)$.

Let us apply this method to study the breathing dynamics of two bosons in a harmonic trap, $V(x)=\frac{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)=\frac{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_{\text{eff}}(x)=\frac{1}{2}k_{t}x^{2}+\int\rho(y)V(x,y)dy$. We illustrate in Fig. 1 the behavior of $\rho(0)$ for $k_{t}=k_{i}=1$ computed by the IPW method with 6 orbitals compared with the exact dynamics (atomic units are used in the rest of this paper).

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,\cdots$ etc., while, as before, we denote the Bohmian coordinates by upper case letters. A single configuration of Bohmian walkers is denoted by $\left(X,Y,Z,\cdots\right)$.

Let us denote the conditional wavfunctions by

[TABLE]

and so on. The equation of motion for $\psi^{0}(x;Y,Z,\cdots)$ is a simple generalization of Eq. (3) for the case of many particles:

[TABLE]

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

[TABLE]

In order to compute $\psi^{1}(x;Y^{\prime},Z,\cdots;t)$ from $\psi^{0}_{w}(x;Y_{w},Z_{w},\cdots;t)$ belonging to all configurations, we need to express the tensor $[\phi_{j}^{\prime}(Y)\phi_{k}(Z)...]$ as a linear superpositions of all the $[\phi_{j}(Y_{w})\phi_{k}(Z_{w})...]$ tensors belonging to all configurations; i.e., $[\phi_{j}^{\prime}(Y)\phi_{k}(Z)...]=\sum_{w}\alpha_{w}[\phi_{j}(Y_{w})\phi_{k}(Z_{w})...]$. 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. Since the size of the vectors now becomes exponentially bigger as the number of particles becomes larger, the bottleneck of this method would be to take a sufficiently large number of configurations that ensures having a complete linear system. Therefore, an immediate room for improvement here would be to find smart tactics to overcome this problem.

**Computing observables from conditional wavefunctions **

In order to compute the expectation value of an operator, we can treat the collection of normalized CWs belonging to all the particles as if they describe normal single particle wavefunctions. $\langle\hat{A}_{x}\rangle$ can then be computed as $\langle\hat{A_{x}}\rangle\approx\frac{1}{N_{w}}\sum_{w}\int\tilde{\psi}_{w}^{*}(x)\hat{A}_{x}\tilde{\psi}_{w}(x)dx$, where $\tilde{\psi}_{w}(x)$ is the normalized CW of particle $x$ belonging to the $w^{\text{th}}$ configuration and $N_{w}$ is the number of configurations. If $\hat{A}$ 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 $\tilde{V}(x)\approx\frac{1}{N_{w}}\sum_{w}\int\tilde{\psi}_{w}^{*}(y)V(x-y)\tilde{\psi}_{w}(y)dy$. The expectation value of $\tilde{V}(x)$ is then computed as a one-body operator. From the normalized collection of all the CWs, we can also get an approximation for the reduced density matrix of one particle $\rho(x^{\prime},x)\approx\frac{1}{N_{w}}\sum_{w}\tilde{\psi}_{w}(x^{\prime})\tilde{\psi}_{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 $\rho(x^{\prime},x)$.

**Many particles in a harmonic trap **

Let us apply this generalization to the dynamics of 3 bosons and 5 bosons in a harmonic trap $V(x)=\frac{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 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 $\langle x^{2}\rangle$ as a function of time.

For a harmonic interaction of the form $V(x,y)=\frac{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 multiconfigurational time-dependent Hartree method for bosons (MCTDHB) alon2008 ; lode2016 ; fasshauer2016 ; mctdh-x ; mctdh-b and with the Hermitian limit (HL) of Eq. (7) (also referred to as small entanglement approximation) 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. christov2007 which does not take entanglement into consideration. It was also employed recently in albareda2014 ; albareda2015 in order to devise an approximate solution for electron-nuclear dynamics in molecular systems.

In Fig. 2, we show the results of computing $\langle x^{2}\rangle$ by averaging 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 shown in Fig. 2(c) for $k_{i}=0.1$ while the Bohmian trajectories for all the 5 particles in a single configuration are shown in Fig. 2(d). The few constant trajectories appearing in Fig. 2(c) correspond to the cases where we manually set the Bohmian velocities to be zero when the denominator in Eq. (1) is below a certain threshold.

In Fig. 3, we plot $\langle x^{2}\rangle$ after switching on a gaussian interaction $V(x-y)=k_{i}/\sqrt{2\pi\sigma^{2}}\times\text{e}^{\frac{-(x-y)^{2}}{2\sigma^{2}}}$ with $k_{i}=0.1$, $\sigma=0.25$ and compare the results with MCTDHB simulation and the HL of Eq. (7). In this calculation, we compute $\langle x^{2}\rangle$ 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 Born1928 , the system remains in the ground state of the instantaneous Hamiltonian. In Fig. 4, 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)=\frac{1}{2}k_{i}(x-y)^{2}$ adiabatically and compare it with the exact ground state energy $E_{0}=\frac{N_{B}-1}{2}\sqrt{1+k_{i}N_{B}}+0.5$ yan2003 for $k_{i}=0.1$. The ground state energy computed by MCTDHB lode2012 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 christov2007 . The optimal relation between the real and complex time evolution constitutes an interesting topic of research.

**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. Our method builds on the formal expansion (Eq. 3), but as we find that the truncation of this set of equations quickly leads to errors, we introduce and apply a truncation-free method that provides the lowest order pilot waves in a self-consistent manner.

The accuracy of this new approach is confirmed but also outperformed by the state-of-the-art MCTDHB algorithm. In the MCTDHB 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 $\alpha_{w}$. So, the complexity of our approach still increases exponentially with the number of particles. It is worth mentioning also that MCTDHB is much faster than our algorithm. While a typical result in the previous figures takes a few hours to compute, it takes much less time by the well developed MCTDHB.

Improvements on our method may come from: (i) strategies that minimize the number of configurations, and hence the computational power, required to evolve the pilot waves without having an underdetermined system of equations, and (ii) optimal choice of the basis functions (possibly an adaptive set of orbitals) that lead to the most compact representation of the full wavefunction, and hence the smallest number of orbitals to capture the dynamics of the many-particle wavefunction.

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 walmsley2017 that even the Hermitian limit christov2007 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 wyatt2006 , the methods developed in this regard in the literature goldfarb2007 ; Bittner2012 may benefit the solution of this problem.

As a final comment, we note that in order to describe entangled dynamics, we need to consider many interacting configurations of Bohmian particles in the present approach. A similar situation arises in hall2014 ; herrmann2017 where the concept of interacting classical worlds was introduced. This similarity between the two approaches may be worth further attention in discussions of the foundations of quantum mechanics.”

**Methods **

**Derivation of Equation (3) **

Let us derive Eq. (3) with respect to $\psi_{1}^{n}(x_{1},t)$. Since $\psi_{1}^{n}(x_{1},t)\equiv\frac{\partial^{n}\Psi(x_{1},x_{2},t)}{\partial x_{2}^{n}}|_{x_{2}=X_{2}(t)}$, we find by the chain rule that

[TABLE]

By exchanging the time and spatial derivatives in the first term on the R.H.S. and using the TDSE:

[TABLE]

we obtain after substituting back in Eq. (9)

[TABLE]

By applying the chain rule to the second term on the R.H.S. we obtain $\frac{\partial^{n}}{\partial x_{2}^{n}}\left.\left(V(x_{1},x_{2})\Psi(x_{1},x_{2})\right)\right|_{x_{2}=X_{2}(t)}=\sum_{k=0}^{n}\binom{n}{k}\psi_{1}^{n-k}(x_{1},t)\frac{\partial^{k}V(x_{1},x_{2})}{\partial x_{2}^{k}}|_{x_{2}=X_{2}(t)}$ which after substituting in Eq. (10) recovers Eq. (3).

**Evolving a hierarchy of pilot waves **

Let us illustrate the inefficiency of evolving a truncated hierarchy of $\psi_{i}^{n}$ using Eq. (3) in order to compute the dynamics of an entangled system. We consider the entangled dynamics of two particles of masses $m_{1}=1$ and $m_{2}=100$ subject to the harmonic potential $V(x_{1},x_{2})=\frac{1}{2}kx_{1}^{2}+\frac{1}{2}kx_{2}^{2}$ with $k=0.1$. Let us take the initial state to be the entangled ground state of the Hamiltonian with the potential function $V(x_{1},x_{2})=\frac{1}{2}k_{1}x_{1}^{2}+\frac{1}{2}k_{2}x_{2}^{2}+\frac{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$. This is an entangled state. We 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 Hermitian limit, i.e., noninterating pilot waves. Fig. 5 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 $\{\psi_{i}^{n}\}$ propagate very fast to $\psi_{i}^{0}$. Increasing $N$ beyond $N\approx 10$ will not help in prolonging the range of accuracy because of the numerical errors in the calculation of the higher-order derivatives of the wavefunction.

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, i.e., do a semiclassical approximation struyve2015 for one particle only.

**Interacting Pilot Waves for two particles **

For the two-particle IPW simulations in Fig. 1 we express all the CWs in terms of the basis set $\{\phi_{i}\}$, and evolve the expansion coefficients $\{a_{i}\}$. Equation (3) is then expressed as

[TABLE]

By taking the inner product with each of the orbitals $\{\phi_{i}(x)\}$ we obtain the time derivative of the expansion coefficients $\{\dot{a}_{i}\}$. This system of equations is then solved using a fourth-order Runge-Kutta method.

**Propagating the conditional wavefunctions **

Each of equations (3) and (7) represents a system of coupled nonlinear and nonunitary differential equations that can be cast in the form

[TABLE]

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:

[TABLE]

In order to propagate $\psi(x,t)$ for a single time step from $t=0$ to $t=\delta t$ using this solution, both the operator $e^{-iH\delta t}$ and $e^{-iHt^{\prime}}$ are performed using a split-operator method kosloff1988 . The integral is performed using the trapezoidal rule $\int_{0}^{\delta t}e^{iHt^{\prime}}W(x,t^{\prime})dt^{\prime}\sim\frac{1}{2}\left[e^{iH\delta t}W(x,\delta t)+W(x,0)\right]$.

**Miscellaneous numerical techniques **

• •

We use the split operator technique [36] in order to propagate the exact two-particle wavefunction in imaginary time (to generate the ground state) and in real time (as in Fig. 1) or to generate the exact Bohmian trajectories (as in Fig. 5).

• •

The evaluation of the conditional wavefunctions at the position of the Bohmian particles was performed by FFT-based interpolation.

• •

The spatial derivatives of the wavefunction to compute $\{\psi_{1}^{n},\psi_{2}^{n}\}$ are computed using the Fast Fourier Transform with the FFTW package.

• •

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 mkl .

**Acknowledgments ** T. A. Elsayed thanks D. A. Deckert, H. Miyagi, A. I. Streltsov, L. F. Buchmann, C. Leveque and Andreas Caranti for brief discussions, and I. Christov for his comments on an earlier version of the manuscript. The authors acknowledge financial support by Villum Foundation. Part of this paper was written while T. A. Elsayed was visiting Hunter College of the City University of New York.

**Author contributions ** TAE conceived the research idea, 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 interests ** The authors declare no competing interests.