format: text

text:
# Numerical study of the coupled time-dependent Gross-Pitaevskii equation: Application to Bose-Einstein condensation

## I Introduction

The experimental detection  of Bose-Einstein condensation (BEC) at ultralow temperature in dilute bosonic atoms (alkali and hydrogen atoms) employing magnetic traps have spurred intense theoretical activities on various aspects of the condensate . Many properties of the condensate are usually described by the nonlinear mean-field Gross-Pitaevskii (GP) equation . The GP equation in both time-dependent and independent forms is formally similar to the Schrödinger equation including a nonlinear term .

More recently, there has been experimental realization of BEC involving two types of atoms . In the actual experiment 87Rb atoms formed in the $`F=1`$, $`m=1`$ and $`F=2`$, $`m=1`$ states by the use of a laser served as two different species of atoms, where $`F`$ and $`m`$ are the total angular momentum and its projection . In another experiment a coupled BEC was formed with the 87Rb atoms in the $`F=1`$, $`m=1`$ and $`F=2`$, $`m=2`$ states. . Experimentally, it is possible to use the same magnetic trap to confine the atoms in two quantum states, which makes this study easier technically compared to the formation of a BEC with two different types of atoms requiring two different trapping mechanisms. It has also been found in these studies  that the 87Rb atoms have a repulsive interaction in all three states considered above. These experiments initiated theoretical activities in multicomponent BEC described by the coupled GP equation .

A numerical study of the time-dependent coupled GP equation is interesting as this can provide solution to many time-evolution problems involving more than one types of atoms forming a BEC. The solution of the coupled nonlinear GP equation is nontrivial  and here we undertake the challenging task of the numerical study of these time-evolution problems.

In a multi-component BEC the main feature is the coupling between different types of condensates, which can lead to new effects associated with this novel and surely much richer situation inexistent in a single-component BEC. We list some interesting possibilities below which we shall investigate numerically in this study. It is possible to have a distinct trapping frequency for each of the components, one of which can be suddenly altered experimentally. Of course, this change would affect the component of the BEC directly trapped by this field. However, it is more interesting to study how this sudden change affects the other component of BEC not directly trapped by this field. Also it is possible to suddenly vary  the atomic scattering length of one of the species and study its effect on the other component. The above variation of one of the trap frequencies or scattering lengths can be carried through in a periodic fashion and its effect on the other component can be studied. The studies mentioned above are peculiar to a coupled BEC and are of interest as it is now possible to vary both the trap frequencies  as well as the scattering lengths  both abruptly or in a periodic fashion. These effects do not have any analogue in the uncoupled BEC and we motivate the present investigation with special emphasis on these effects in this study. We investigate these time-evolution problems using a set of two coupled GP equations in the purely repulsive case. In addition we study the stationary solution to the coupled GP equation describing a multi-component condensate, where we consider up to four components.

We solve the coupled BEC problem using the time-dependent coupled GP equation in cases of attractive and repulsive atomic interactions by discretization with the Crank-Nicholson-type rule complimented by the known boundary conditions at origin and infinity . This procedure leads to good convergence for both the stationary and time-evolution problems.

First we consider stationary coupled condensates under the action of trap potentials. Stable and converged numerical results are obtained for up to four coupled equations in the repulsive case and two in the attractive case. The time-dependent GP equation is directly solved to obtain the full time-dependent solution in the case of stationary problems, from which a trivial time dependent phase factor is separated and the stationary solution obtained as in the uncoupled case .

We also study the numerical stability of the calculational scheme which is more difficult to obtain when the nonlinearity is large. For this purpose we only consider repulsive interatomic interactions in the time-evolution problems where condensates with large nonlinearity can be formed. In the case of attractive interaction a large nonlinearity leads to the collapse of the condensate .

In Sec. II we present the coupled time-dependent GP equation which we use. In Sec. III we describe the numerical method in some detail. In Sec. IV we report the numerical results for the stationary case and in Sec. V we report a study of three types of time-evolution problems. Finally, in Sec. V we give a summary of our investigation.

## II Nonlinear Coupled Gross-Pitaevskii Equation

The GP equation  for a coupled trapped Bose-Einstein condensate at zero temperature is written as

$`[{\displaystyle \frac{\mathrm{}^2}{2m}}^2`$ $`+`$ $`{\displaystyle \frac{1}{2}}c_im\omega ^2r^2+{\displaystyle \underset{l=1}{\overset{M}{}}}g_{il}N_l|\mathrm{\Psi }_l(𝐫,\tau )|^2`$ (1)
$``$ $`i\mathrm{}{\displaystyle \frac{}{\tau }}]\mathrm{\Psi }_i(𝐫,\tau )=0,`$ (2)
where $`\mathrm{\Psi }_i(𝐫,\tau )`$ at position $`𝐫`$ and time $`\tau `$ is the wave function for the component $`i`$ of the condensate, $`m`$ is the mass of a single bosonic atom, $`N_l`$ the number of condensed atoms of type $`l`$, $`M`$ the number of types of atoms, $`c_im\omega ^2r^2/2`$ the attractive harmonic-oscillator trap potential, $`\omega `$ the oscillator frequency. The parameter $`c_i`$ has been introduced to independently modify the frequency of the harmonic oscillator trap for each type of atoms. Here $`g_{il}=4\pi \mathrm{}^2a_{il}/m`$ is the coupling constant for elastic interaction between atoms of types $`i`$ and $`l`$, where $`a_{il}`$ is the corresponding scattering length. The masses of different types of atoms are taken to be equal, as this is necessary while considering the coupled BEC formed of different spin states of 87Rb $``$ one of the most important realization to date. In this work we shall not allow the transition of one type of atoms of the BEC to the other and take the number of atoms of each component to be constant as in the experiment of Ref. .

Here we shall be interested in the spherically symmetric solution $`\mathrm{\Psi }_i(𝐫,\tau )\psi _i(r,\tau )`$ to Eq. (1), which can be written as

$`[{\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{1}{r}}{\displaystyle \frac{^2}{r^2}}r`$ $`+`$ $`{\displaystyle \frac{1}{2}}c_im\omega ^2r^2+{\displaystyle \underset{l=1}{\overset{M}{}}}g_{il}N_l|\psi _l(r,\tau )|^2`$ (3)
$``$ $`i\mathrm{}{\displaystyle \frac{}{\tau }}]\psi _i(r,\tau )=0,`$ (4)
The above limitation to the spherically symmetric solution (in the zero angular momentum state) reduces the coupled GP equation to a one-dimensional coupled partial differential equation.

As in Refs.  it is convenient to use dimensionless variables defined by $`x=\sqrt{2}r/a_{\text{ho}}`$ , and $`t=\tau \omega ,`$ where $`a_{\text{ho}}\sqrt{\mathrm{}/(m\omega )}`$, and $`\varphi _i(x,t)=x\psi _i(r,\tau )(\sqrt{2}\pi a_{\text{ho}}^3)^{1/2}`$ . In terms of these variables Eq. (3) becomes

$`[{\displaystyle \frac{^2}{x^2}}`$ $`+`$ $`{\displaystyle \frac{1}{4}}c_ix^2+{\displaystyle \underset{l=1}{\overset{M}{}}}n_{il}{\displaystyle \frac{|\varphi _l(x,t)|^2}{x^2}}`$ (5)
$``$ $`i{\displaystyle \frac{}{t}}]\varphi _i(x,t)=0,`$ (6)
where $`n_{il}2\sqrt{2}N_la_{il}/a_{\text{ho}}`$ is the reduced number of particles and this number could be negative (positive) when the corresponding scattering length is negative (positive) representing an attractive (repulsive) interatomic interaction. The normalization of the wave function is $`_0^{\mathrm{}}|\varphi _i(x,t)|^2𝑑x=1`$ and its root-mean-square (rms) radius $`x_{\text{rms}}^{(i)}`$ is given by

$$x_{\text{rms}}^{(i)}=\left[_0^{\mathrm{}}x^2|\varphi _i(x,t)|^2𝑑x\right]^{1/2}.$$
(7)
## III Numerical Method

To solve Eq. (5) numerically one needs the proper boundary conditions at $`x0`$ and $`\mathrm{}`$. For a confined condensate, for a sufficiently large $`x`$, $`\varphi _i(x,t)`$ must vanish asymptotically. Hence the nonlinear term proportional to $`|\varphi _i(x,t)|^3`$ can eventually be neglected in the GP equation for large $`x`$. Consequently the the asymptotic form of the physically acceptable solution is given by $`lim_x\mathrm{}|\varphi _i(x,t)|\mathrm{exp}(x^2/4).`$ Next we consider Eq. (5) as $`x0`$. The nonlinear term approaches a constant in this limit because of the regularity of the wave function at $`x=0`$. Then one has the condition $`|\varphi _i(0,t)|=0.`$

A convenient way to solve Eq. (5) numerically is to discretize it in both space and time and reduce it to a set of algebraic equations which could then be solved by using the known asymptotic boundary conditions. The procedure is similar to that in the uncoupled case . We discretize Eq. (5) by using a space step $`h`$ and time step $`\mathrm{\Delta }`$ with a finite difference scheme using the unknown $`(\varphi _i)_p^k`$ which will be approximation of the exact solution $`\varphi _i(x_p,t_k)`$ where $`x_p=ph`$ and $`t_k=k\mathrm{\Delta }`$. The time derivative in Eq. (5) involves the wave function at times $`t_k`$ and $`t_{k+1}`$. As in the uncoupled case we express the wave functions and their derivatives in Eq. (5) in terms of the average over times $`t_k`$ and $`t_{k+1}`$  and the resultant scheme leads to accurate results and good convergence. In practice we use the following Crank-Nicholson-type scheme  to discretize the partial differential equation (5)

$`i{\displaystyle \frac{(\varphi _i)_p^{k+1}(\varphi _i)_p^k}{\mathrm{\Delta }}}`$ $`=`$ $`{\displaystyle \frac{1}{2h^2}}[(\varphi _i)^{k+1}_{p+1}2(\varphi _i)^{k+1}_p+(\varphi _i)^{k+1}_{p1}`$ (8)
$`+`$ $`(\varphi _i)_{p+1}^k2(\varphi _i)_p^k+(\varphi _i)_{p1}^k]`$ (9)
$`+`$ $`{\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{c_ix_p^2}{4}}+{\displaystyle \underset{l=1}{\overset{M}{}}}n_{il}{\displaystyle \frac{|(\varphi _l)_p^k|^2}{x_p^2}}\right]`$ (10)
$`\times `$ $`[(\varphi _i)_p^{k+1}+(\varphi _i)_p^k].`$ (11)
Considering that the wave function components $`\varphi _i`$ are known at time $`t_k`$, Eq. (11) is an equation in the unknowns $``$ $`(\varphi _i)_{p+1}^{k+1},(\varphi _i)_p^{k+1}`$ and $`(\varphi _i)_{p1}^{k+1}`$. In a lattice of $`N`$ points Eq. (11) represents a tridiagonal set for $`p=2,3,\mathrm{},(N1)`$ for a specific component $`\varphi _i`$. This set has a unique solution if the wave functions at the two end points $`(\varphi _i)_1^{k+1}`$ and $`(\varphi _i)_N^{k+1}`$ are known . In the present problem these values at the end points are provided by the known asymptotic conditions.

To find the ground state of the condensate we start with the analytically known properly normalized wave functions of the uncoupled harmonic oscillator problems described by Eq. (5) with $`n_{il}=0`$. We then repeatedly propagate these solutions in time using the Crank-Nicholson-type algorithm (11). Starting with $`n_{il}=0`$, at each time step we increase or decrease the nonlinear parameter $`n_{il}`$ by an amount $`\mathrm{\Delta }_1`$. This procedure is continued until the desired final value of $`n_{il}`$ is reached. The resulting solution is the ground state of the condensate corresponding to the specific nonlinearity.

The time-dependent approach is the most suitable for solving time-evolution problems. In the present study we only consider evolutions problem starting from a stable condensate at $`t=0`$. In these cases the stationary problem is solved first and the wave function so obtained serves as the starting wave function for the time-evolution problem.

## IV Results for the Stationary Problem

First we consider the stationary ground-state solution of Eq. (5) in cases of both attractive and repulsive interactions using two and four coupled equations. The numerical integration was performed up to $`x_{\text{max}}=15`$ with $`h=0.0001`$ using time step $`\mathrm{\Delta }=0.05`$ and the parameter $`\mathrm{\Delta }_1=0.01`$. After some experimentation we find that good convergence is obtained with parameters $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }_1`$ near these values. The convergence is fast for small nonlinearity. The final convergence of the scheme breaks down if nonlinearity is too large. In practice these difficulties start for $`n_{ii}>20`$ for the ground state for a repulsive interaction in a computational analysis in double precision. For an attractive interaction the coupled GP equation does not sustain a large nonlinearity $`|n_{il}|`$ and leads to collapse. Except for values of nonlinearity near collapse, the GP equation in the attractive case leads to good convergence.

### A Repulsive Atomic Interaction

In most of the experimental realization of BEC in trapped atoms, the interatomic interaction is repulsive and we consider this case first. We consider the simple case of two coupled GP equations in the case of repulsive interaction with (a) $`n_{11}=n_{22}=10`$, $`n_{12}=n_{21}=5`$, $`c_1=1`$ and $`c_2=0.25`$, (b) $`n_{11}=n_{22}=5`$, $`n_{12}=n_{21}=100`$, $`c_1=1`$ and $`c_2=0.25`$, and (c) $`n_{11}=n_{22}=10`$, $`n_{12}=n_{21}=5`$, and $`c_1=c_2=1`$. In this case all interactions are repulsive corresponding to the positive sign of all $`n_{il}2\sqrt{2}N_ia_{il}/a_{\text{ho}}`$. Although these parameters are in dimensionless units, it is easy to associate them to an actual physical problem of experimental interest. For example, for the mixture of $`|F=1,m=1`$ and $`|F=2,m=1`$ 87Rb states, the ratio of scattering lengths $`a_{|1,1}/a_{|2,1}=1.062`$ . If we label the $`|1,1`$ state by 1 and the $`|2,1`$ state by 2, and consider $`a_{11}/a_{\text{ho}}a_{22}/a_{\text{ho}}0.002`$, then $`n_{11}=n_{22}=10`$ corresponds to $`N_1N_21770`$. This estimate gives an idea of the actual experimental condition that the present set of parameters simulate. The three models considered above can simulate actual experimental situations composed of two states of 87Rb. The different values of $`n_{12}`$ and $`n_{21}`$ considered above corresponds to different possible unknown repulsive interactions among the two species of condensates. It is realized from the coupled GP equation that in case of model (c) $`\varphi _1=\varphi _2`$. We show results for the two components of the wave function for two sets of values of $`\mathrm{\Delta }_1`$: 0.01 (full line) and 0.015 (dashed line) in Fig. 1. The difference between the two sets of results increases as the nonlinearity of the GP equation given by $`n_{il}`$ increases, e.g., for the case (b) above compared to (a). The difference reduces to zero as the nonlinearity decreases.

Next we consider the more complicated case of four coupled GP equations with repulsive atomic interaction. This is a purely theoretical case with no experimental analogue, as all experiments to date are limited to two coupled condensates only. In this case the numerical method works in the same fashion and good convergence is obtained with moderate values of nonlinearity. Again we consider the case of repulsive interaction between all possible pairs. In this case in the four-component model we take $`n_{11}=4`$, $`n_{22}=5`$, $`n_{33}=6`$, $`n_{44}=8`$, and $`n_{il}=2,il`$; $`c_1=4`$, $`c_2=1`$, $`c_3=4`$, and $`c_4=0.25`$. The solution for the wave-function components obtained with $`\mathrm{\Delta }=0.05`$ and $`\mathrm{\Delta }_1=0.01`$ (full line), and $`\mathrm{\Delta }_1=0.015`$ (dashed line) are shown in Fig. 2. The maximum difference between the two calculations is about 6$`\%`$ for the largest component ($`\varphi _1`$) near $`x=0`$, although the average difference is much smaller as can be seen in Fig. 2.

### B Attractive Atomic Interaction

The case of attractive atomic interaction demands special attention as one can have the phenomenon of collapse in this case. We consider the case of two coupled GP equations with attractive interactions between like atoms $`ii`$, $`i=1,2,`$ and with repulsive interactions between unlike atoms $`ij`$, $`ij.`$ In this case some of the atomic interactions are attractive in nature. Consequently, with large (attractive) nonlinearity the system may undergo collapse and for stable stationary solution of the GP equation the nonlinearity should be maintained small. We consider the following three cases: (a) $`n_{11}=1`$, $`n_{22}=1.5`$, $`n_{12}=n_{21}=4`$, $`c_1=4`$, and $`c_2=1`$; (b) $`n_{11}=1`$, $`n_{22}=1.5`$, $`n_{12}=n_{21}=80`$, $`c_1=4`$, and $`c_2=1`$; (c) $`n_{11}=1`$, $`n_{22}=1.5`$, $`n_{12}=n_{21}=100`$, $`c_1=4`$, and $`c_2=1`$. The only observed case of BEC with attractive interaction is the case of 7Li with $`a/a_{\text{ho}}0.0005`$ . If we label this state by 1 then $`n_{11}=1`$ corresponds to the actual particle number $`N_1=700`$. In an uncoupled condensate of 7Li the BEC collapses for more than 1400 atoms. Although there has been no experimental realization of coupled BEC in the case of attractive interaction, the parameters cited above may simulate the BEC of ground-state atoms of 7Li coupled to one of its excited states, where the atomic interaction is also attractive. If we assume that the excited-state atoms have the same value of $`a/a_{\text{ho}}`$ as in the ground state then $`n_{22}=1.5`$ corresponds to the number of atoms $`N_2=1500`$ in the excited state where the excited state is labeled by the index 2.

The wave-function components in this case are shown in Fig. 3. As for stable stationary solution the nonlinearity in this case has to be smaller than in the purely repulsive case, numerically it is easier to obtain precise solution except for values of the nonlinearity close to (and beyond) collapse. The parameters above in cases (a), (b), and (c) are chosen to illustrate the collapse of the system arising from the divergence of the first component of the wave function ($`\varphi _1`$). The nonlinear couplings $`n_{12}`$ and $`n_{21}`$ are increased as we move from case (a) to (b) and then to (c). Other parameters of the model are kept fixed. Because of the attractive interactions in this case, the system moves towards collapse as we move from case (a) to (c) through (b). The central density corresponding to $`\varphi _1`$ increases (eventually tends to infinity) and the rms radius of the system decreases (eventually tends to zero) with the increase of nonlinearity. This is clear from a comparison of Fig. 3 with Figs. 1 and 2. The very large central value of the wave function $`\varphi _1`$ ($`35`$) and its small radial extension indicates a large central density and a small rms radius.

Finally, we consider the case of two interacting systems with all interactions repulsive. In this case the system is more vulnerable to collapse if the nonlinearity is large. We consider the following three sets of parameters in this case for which we show the solution in Fig. 4: (a) $`n_{11}=1`$, $`n_{22}=1`$, $`n_{12}=n_{21}=0.4`$, $`c_1=4`$, and $`c_2=0.25`$; (b) $`n_{11}=1`$, $`n_{22}=1`$, $`n_{12}=n_{21}=0.5`$, $`c_1=4`$, and $`c_2=0.25`$; (c) $`n_{11}=1`$, $`n_{22}=1`$, $`n_{12}=n_{21}=0.55`$, $`c_1=4`$, and $`c_2=0.25`$. These parameters simulate the possible coupled BEC composed of the attractive ground and excited states of 7Li where the interaction between a ground- and an excited-state atom is also taken to be attractive. It is possible to calculate the number of the two types of atoms as in the discussion related to Fig. 3. Here the nonlinearity increases as we move from case (a) to (b) and then to (c), and consequently, the wave-function components become more and more localized with a large central density and small rms radii signaling the onset of collapse of the system. This is clear from Fig. 4. In case (c) the nonlinearity is the highest and one is closer to collapse. However, there is a difference between the two collapses shown in Figs. 3 and 4. In Fig. 3 the route to collapse is manifested through a singular behavior of component $`\varphi _1`$ of the wave function; in Fig. 4 both components $`\varphi _1`$ and $`\varphi _2`$ exhibit the singular behavior. The spacial extension of the wave function components close to collapse in Figs. 3 and 4 is much smaller than the wave function components in the purely repulsive cases shown in Figs. 1 and 2.

### C Estimate of Numerical Error

It is appropriate to comment quantitatively on the numerical accuracy of the present method. If we iterate the final solution in time without changing the nonlinearity, the numerical result keeps on oscillating with a small amplitude around the converged value. This oscillation gives a good estimate of the numerical error of the method. This error manifests in a different way in Fig. 1, where we have varied $`\mathrm{\Delta }_1`$. The numerical solution of the time-dependent method is independent of the space step $`h`$ provided that a typical value around $`h=0.0001`$ is employed as in the present study. No visible difference in the solution is found if $`h`$ is increased by a factor of two. However, the above oscillations with time iteration are sensitive to the parameters $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }_1`$. The values of these parameters ($`\mathrm{\Delta }=0.05`$ and $`\mathrm{\Delta }_1=0.01`$) are chosen to minimize the oscillation of the results with time iterations. The oscillation increases if larger or smaller values of one or both of these parameters are employed and can really be large if an improper value of step $`\mathrm{\Delta }`$ or $`\mathrm{\Delta }_1`$ is chosen. This oscillation is quite similar to that in the uncoupled case studied in detail in Ref. .

From Figs. 1 and 2 we find that the error in the wave function component $`|\varphi _i(x,t)/x|`$ as a function of $`x`$ is the largest at $`x=0`$. In Table I we show the percentage error of $`lim_{x0}|\varphi _i(x,0)/x|`$, defined by $`100[|\varphi _i(0,t)||\varphi _i(0,0)|]/|\varphi _i(0,0)|`$ at those iterations where this error is maximum. For illustration we consider the models (a) discussed in Figs. 1 and 4. At positive $``$ there is overshooting and at negative $``$ there is undershooting. Between a positive $``$ and a negative $``$ there is a zero of $``$ denoting zero error. We find that the wave functions oscillate with time around the stationary solution. The maximum reported error is about 6$`\%`$. Considering that we are dealing with coupled nonlinear equations these errors are well within the acceptable limit. The errors shown in Table I would also be the typical errors in time-evolution problems with same nonlinearity which we study in Sec. V. From Table I we see that the period of oscillation of the result varies from one case to another. However, the error increases as the nonlinearity increases or as the system approaches collapse in the case of attractive interaction. For example, in models (c) of Figs. 3 and 4, which are close to collapse corresponding to almost maximum permissible nonlinearity, the error increases quickly with time iteration and a large numerical error could be generated easily.

## V Results for the Time-Evolution Problem

Now we consider three types of time-evolution problems, some of which could possibly be studied experimentally. As the repulsive case leads to more stable configuration of the condensate, in this study we consider only this case in the process of time evolution. The attractive case of coupled BEC is also very interesting from a physical point of view because of the occurrence of collapse. We have performed a study of the dynamics of collapse in coupled BEC using the present numerical method, which will be reported elsewhere. The two types of parameters that can be varied in the time-evolution study are the frequencies of the harmonic oscillator traps and the different scattering lengths. Recently, it has been possible to vary the scattering length experimentally by varying an external field . It is also possible to vary the trap frequency by varying the currents in the magnets responsible for confinement .

In the first type of problems we consider a sudden change of the harmonic oscillator frequencies or scattering lengths at $`t=0`$ and study its effect on a preformed condensate. In the second type we study the effect of a periodic temporal variation of these frequencies on a preformed condensate. Finally, we study the effect of a periodic temporal variation of the scattering lengths on the preformed condensate. In all cases we take the preformed condensate as the one described by the model (a) of Fig. 1. We have commented before that the parameters of this model can simulate the coupled BEC composed of the ground and an excited state of 87Rb, which gives a motivation for this choice. When we implement these time-dependent perturbations, the system starts to oscillate (grow and shrink) with time. The corresponding evolution can be studied best through the rms radii  which execute periodic oscillation with time.

### A Sudden Change of Trap Frequency or Scattering Length

By varying the external fields it is possible to vary the harmonic oscillator trap frequency of the confining trap as well as the atomic interactions (scattering lengths) . First we consider a sudden change of both the trap frequencies on the preformed coupled stationary BEC state corresponding to model (a) of Fig. 1. We set the reduced time $`Tt/0.05=0`$ when we start the time evolution. As the time step $`\mathrm{\Delta }`$ is 0.05, $`T`$ is just the number of iterations. In this model we have two different trap frequencies for the two components given by $`c_1=1`$ and $`c_2=0.25`$. In the first case at $`T=0`$ we suddenly interchange the constants $`c_1`$ and $`c_2`$, i.e., we set $`c_1=0.25`$ and $`c_2=1`$ and study the time evolution. The evolution of the rms radii corresponding to $`\varphi _1`$ and $`\varphi _2`$ are shown in Fig. 5 (a). Both rms radii execute oscillations. However, that corresponding to $`\varphi _1`$ has a much larger amplitude. The periods of oscillation of the two radii are different.

Next we consider a sudden change in one of the trap frequencies on the preformed condensate at $`T=0`$ corresponding to model (a) of Fig. 1 with $`c_1=1,c_2=0.25`$. At $`T=0`$, we set $`c_1=c_2=1`$, which corresponds to parameters of model (c) of Fig. 1. The evolution of the rms radii is shown in Fig. 5 (b). Although in the stationary configuration of model (c) of Fig. 1, $`\varphi _1=\varphi _2`$, this condition is never attained in this evolution problem. The system keeps on oscillating indefinitely with time. The oscillation shown in Figs. 5 (a) and (b) has nothing to do with the nonlinear or coupled nature of the problem. Similar oscillation also appears in an uncoupled linear oscillator when the trap frequency is suddenly changed. In the present coupled nonlinear problem both rms radii execute oscillations with time. However, when the amplitude of oscillation of one of the components increases, that of the other decreases. This behavior denotes the transfer of kinetic energy from one component to the other.

Now we study the effect of a sudden change of the scattering length(s) on the preformed condensate . We consider the problem when the parameters of model (a) of Fig. 1 are suddenly changed to those of model (b) of Fig. 1 at $`T=0`$. This is achieved by changing the nonlinearities suddenly at $`T=0`$ from $`n_{11}=n_{22}=10,n_{12}=n_{21}=5`$ to $`n_{11}=n_{22}=5,n_{12}=n_{21}=100`$ with a variation of the external field which controls the scattering length(s). In this case the oscillations of the system are shown in Fig. 5 (c) where we plot the time evolution of the rms radii of the two components. Both components of the condensate execute oscillations but with different frequencies and amplitudes. One of the components execute giant oscillations with large amplitude, whereas the amplitude of the other is much smaller.

Finally, we consider the case when one of the trapping potentials is switched off at $`T=0`$ on the preformed condensate of Fig. 1 model (a) by setting $`c_2=0`$. The oscillation in this case is shown in Fig. 5 (d) where we plot the two rms radii. In the absence of the trapping potential the second component of the condensate can no longer remain localized in space. However, it does not expand monotonically before evaporating. It starts to execute giant oscillation and eventually escapes to infinity. Similar oscillation was found in the case of an uncoupled BEC when the trapping potential was removed . The first component essentially remains unchanged during the process under the action of the unchanged trap potential. The minor oscillation of the rms radii of the first component is due to the coupling to the expanding second component.

### B Periodic Oscillation of Trap Frequency

Instead of making a sudden change in the parameters of the model, next we introduce periodic oscillation in some of the parameters of the model for $`T0`$ and study the consequence on the system. We introduce a periodic variation in the parameters $`c_i`$ which are related to the harmonic oscillator trap frequencies. Experimentally, this variation is possible via a variation of the external fields which are controlled by currents.

We again consider at $`T=0`$ the preformed condensate of the model (a) of Fig. 1. First we consider the variation $`c_1=10.5\mathrm{sin}(\pi T/20)`$ which corresponds to varying the frequency of the first trap. The resultant variation of the two rms radii are shown in Fig. 6. The first radius (full line) oscillates more rapidly with larger amplitude and frequency than the second radius (dashed line). This is reasonable as we are directly varying the first frequency in this case. The rms radius of the second wave function feels the effect through its coupling to the first component. We also varied both the papameters $`c_1`$ and $`c_2`$ in a periodic fashion which corresponds to varying both the frequencies. In this case both the rms radii execute oscillation. However, no interesting effect is observed and we do not show the details of that oscillation here.

### C Periodic Oscillation of Scattering Length

Now we study the effect of a periodic variation of the scattering length(s) of the system on a preformed condensate. In our formulation this corresponds to a periodic variation of the parameters $`n_{il}`$. This variation of the atomic interactions or the scattering lengths is now feasible experimentally . We consider the periodic variation in one of the scattering lengths ($`a_{11}`$) by setting $`n_{11}=10.5\mathrm{sin}(\pi T/20)`$ for $`T0`$ on the preformed condensate of model (a) of Fig. 1. The resultant oscillation of the rms radii are shown in Fig. 7. This variation corresponds to a variation of the atomic interaction among atomic states of the first type. Consequently, the rms radii of the first component of the BEC executes pronounced oscillation with moderate amplitude. There is no direct variation in the parameters of the second component. The second component of the condensate feels the effect of variation of $`n_{11}`$ through the coupling to the first component. Because of this secondary effect the second component also executes oscillation as can be seen from its rms radii in Fig. 7, albeit with a much smaller amplitude compared to the first component.

We also considered a periodic variation of the scattering length between one atom of each type ($`a_{12}=a_{21}`$) by setting $`n_{12}=n_{21}=0.50.25\mathrm{sin}(\pi T/20)`$ on the same preformed condensate for $`T0`$ and studied the resultant oscillation of the rms radii. However, no interesting behavior was observed and we do not show details of that oscillation.

## VI Summary

In this paper we present a numerical study of the coupled time-dependent Gross-Pitaevskii equation for BEC in three space dimensions under the action of harmonic oscillator trap potentials with attractive and repulsive interparticle interactions between different types of atoms .

The time-dependent coupled GP equation is solved by discretizing it using a Crank-Nicholson-type scheme . This method leads to good convergence for small nonlinearity. Numerical difficulty appears for large nonlinearity ($`n_{il}>20`$). For medium nonlinearity, the accuracy of the method can be increased by reducing the space step $`h`$.

The ground-state stationary wave functions are found to be sharply peaked near the origin for attractive interatomic interaction for larger nonlinearity (Fig. 4). For a repulsive interatomic interaction the wave function extends over a larger region of space (Figs. 1 and 2). In the case of an attractive potential, the rms radii decrease with an increase of nonlinearity. There could be a collapse for attractive interaction when the nonlinear parameters $`n_{il}`$ are increased as in the uncoupled case . In the purely repulsive case we solved two and four coupled GP equations. In problems involving attraction we solved only the two coupled GP equations.

In addition to the stationary problem we studied three types of evolution problems. A stable coupled condensate is considered at $`T=0`$ on which a time-dependent perturbation is introduced. Two types of perturbations were considered on a two-component condensate with purely repulsive interactions. In the first type a sudden change in the parameters related to the frequencies of the trap and the scattering lengths was introduced. In the second type a periodic variation of the different scattering lengths and the frequencies of the harmonic oscillator trap was introduced. In all cases the condensates execute periodic oscillation which is studied via the time evolution of the rms radii as in the uncoupled case . We conclude that the present time-dependent approach is very suitable for studying both the stationary and time-evolution problems of a coupled BEC.

The work is supported in part by the CNPq and FAPESP of Brazil.

Figure Caption:

1. Wave function components $`\varphi _1`$ (label 1) and $`\varphi _2`$ (label 2) for two coupled GP equations with (a) $`n_{11}=n_{22}=10,`$ $`n_{12}=n_{21}=5,`$ $`c_1=1`$, $`c_2=0.25`$; (b) $`n_{11}=n_{22}=5,`$ $`n_{12}=n_{21}=100,`$ $`c_1=1`$, $`c_2=0.25`$; and (c) $`n_{11}=n_{22}=10,`$ $`n_{12}=n_{21}=5,`$ $`c_1=1`$, $`c_2=1`$ calculated with $`\mathrm{\Delta }=0.05`$ and $`\mathrm{\Delta }_1=0.01`$ (full line); $`\mathrm{\Delta }=0.05`$ and $`\mathrm{\Delta }_1=0.015`$ (dashed line). In case (c) only the results for $`\mathrm{\Delta }_1=0.01`$ are shown.

2. Wave function components $`\varphi _1`$ (label 1), $`\varphi _2`$ (label 2), $`\varphi _3`$ (label 3), and $`\varphi _4`$ (label 4) for four coupled GP equations with $`n_{11}=4`$, $`n_{22}=5`$, $`n_{33}=6`$, $`n_{44}=8`$, and $`n_{il}=2,il`$; $`c_1=4`$, $`c_2=1`$, $`c_3=4`$, and $`c_4=0.25`$ calculated with $`\mathrm{\Delta }=0.05`$ and $`\mathrm{\Delta }_1=0.01`$ (full line); $`\mathrm{\Delta }=0.05`$ and $`\mathrm{\Delta }_1=0.015`$ (dashed line).

3. Wave function components $`\varphi _1`$ (label 1) and $`\varphi _2`$ (label 2) for two coupled GP equations with (a) $`n_{11}=1`$, $`n_{22}=1.5`$, $`n_{12}=n_{21}=4`$, $`c_1=4`$, and $`c_2=1`$ (dashed-dotted line); (b) $`n_{11}=1`$, $`n_{22}=1.5`$, $`n_{12}=n_{21}=80`$, $`c_1=4`$, and $`c_2=1`$ (dashed line); (c) $`n_{11}=1`$, $`n_{22}=1.5`$, $`n_{12}=n_{21}=100`$, $`c_1=4`$ , and $`c_2=1`$ (full line) calculated with $`\mathrm{\Delta }=0.05`$ and $`\mathrm{\Delta }_1=0.01`$.

4. Wave function components $`\varphi _1`$ (label 1) and $`\varphi _2`$ (label 2) for two coupled GP equations with (a) $`n_{11}=1`$, $`n_{22}=1`$, $`n_{12}=n_{21}=0.4`$, $`c_1=4`$, and $`c_2=0.25`$ (dashed-dotted line); (b) $`n_{11}=1`$, $`n_{22}=1`$, $`n_{12}=n_{21}=0.5`$, $`c_1=4`$, and $`c_2=0.25`$ (dashed line); (c) $`n_{11}=1`$, $`n_{22}=1`$, $`n_{12}=n_{21}=0.55`$, $`c_1=4`$, and $`c_2=0.25`$ (full line) calculated with $`\mathrm{\Delta }=0.05`$ and $`\mathrm{\Delta }_1=0.01`$.

5. The rms radii of the two components $`\varphi _1`$ (full line) and $`\varphi _2`$ (dashed line) of the wave function at different reduced times $`Tt/0.05`$ for the oscillating condensate when on the preformed condensate of model (a) of Fig. 1 we suddenly inflict the following changes: (a) $`c_1=0.25`$ and $`c_2=1`$; (b) $`c_1=c_2=1`$; (c) $`n_{11}=n_{22}=5`$; $`n_{21}=n_{12}=100`$; (d) $`c_2=0`$ calculated with $`\mathrm{\Delta }=0.05`$ and $`\mathrm{\Delta }_1=0.01`$. All other parameters are maintained unchanged.

6. The rms radii of the two components $`\varphi _1`$ (full line) and $`\varphi _2`$ (dashed line) of the wave function at different reduced times $`Tt/0.05`$ for the oscillating condensate when on the preformed condensate of model (a) of Fig. 1 we suddenly inflict the following change: $`c_1=10.5\mathrm{sin}(\pi T/20)`$ calculated with $`\mathrm{\Delta }=0.05`$ and $`\mathrm{\Delta }_1=0.01`$. All other parameters are maintained unchanged.

7. The rms radii of the two components $`\varphi _1`$ (full line) and $`\varphi _2`$ (dashed line) of the wave function at different reduced times $`Tt/0.05`$ for the oscillating condensate when on the preformed condensate of model (a) of Fig. 1 we suddenly inflict the following change: $`n_{11}=10.5\mathrm{sin}(\pi T/20)`$ calculated with $`\mathrm{\Delta }=0.05`$ and $`\mathrm{\Delta }_1=0.01`$. All other parameters are maintained unchanged.

Table I: Percentage error $`100[|\varphi _i(0,T)||\varphi _i(0,0)|]/|\varphi _i(0,0)|`$ of $`|\varphi _i(0,T)|`$ ($`i=1,2`$) at successive reduced times $`T(t/0.05)`$, where this error is maximum, calculated with $`\mathrm{\Delta }=0.05`$ and $`\mathrm{\Delta }_1=0.01`$. The cases considered correspond to model (a) of Figs. 1 and 4

| Fig. 1 | $`T`$ | 0 | 19 | 42 | 67 | 101 | 125 | 165 | 190 | 225 | 249 | 271 | 311 | 337 | 371 | 395 | 414 | 458 | 483 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $`\varphi _1`$(a) | $``$ | 0 | $`1.7`$ | 2.5 | $`1.0`$ | 3.5 | $`1.4`$ | 3.2 | $`1.9`$ | $`2.9`$ | $`2.7`$ | 2.1 | $`3.2`$ | 4.0 | $`2.1`$ | $`4.2`$ | $`1.8`$ | 5.4 | $`3.7`$ |
| Fig. 1 | $`T`$ | 0 | 36 | 85 | 163 | 219 | 275 | 353 | 395 | 477 |  |  |  |  |  |  |  |  |  |
| $`\varphi _2`$(a) | $``$ | 0 | 5.4 | $`4.9`$ | 6.2 | $`4.2`$ | $`5.5`$ | $`4.2`$ | $`6.2`$ | $`4.6`$ |  |  |  |  |  |  |  |  |  |
| Fig. 4 | $`T`$ | 0 | 13 | 110 | 192 | 250 | 318 | 376 | 467 |  |  |  |  |  |  |  |  |  |  |
| $`\varphi _1`$(a) | $``$ | 0 | 1.2 | $`0.6`$ | 1.0 | $`0.5`$ | 1.0 | $`0.3`$ | 0.9 |  |  |  |  |  |  |  |  |  |  |
| Fig. 4 | $`T`$ | 0 | 27 | 113 | 183 | 245 | 323 | 386 | 462 |  |  |  |  |  |  |  |  |  |  |
| $`\varphi _2`$(a) | $``$ | 0 | 3.1 | $`2.7`$ | $`3.2`$ | $`2.5`$ | $`3.2`$ | $`1.7`$ | 3.1 |  |  |  |  |  |  |  |  |  |  |


text:
# Bohmian prediction about a two double-slit experiment and its disagreement with standard quantum mechanics

## 1 Introduction

According to the standard quantum mechanics (SQM), the complete description of a system of particles is provided by its wavefunction. The empirical predictions of SQM follow from a mathematical formalism which makes no use of the assumption that matter consists of particles pursuing definite tracks in space-time. It follows that the results of the experiments designed to test the predictions of the theory, do not permit us to infer any statement regarding the particle–not even its independent existence. But, in the Bohmian quantum mechanics (BQM), the additional element that is introduced apart from the wavefunction is the particle position, conceived in the classical sense as pursuing a definite continuous track in space-time \[1-3\]. The detailed predictions made by this causal interpretation explains how the results of quantum experiments come about but it is claimed that they are not tested by them. In fact when Bohm  presented his theory in 1952, experiments could be done with an almost continuous beam of particles, but not with individual particles. Thus, Bohm constructed his theory in such a fashion that it would be impossible to distinguish observable predictions of his theory from SQM. This can be seen from Bell’s comment about empirical equivalence of the two theories when he said:“It \[the de Broglie-Bohm version of non-relativistic quantum mechanics\] is experimentally equivalent to the usual version insofar as the latter is unambiguous. However, could it be that a certain class of phenomena might correspond to a well-posed problem in one theory but to none in the other? Or might the additional particles and definite trajectories of Bohm’s theory lead to a prediction of an observable where SQM would just have no definite prediction to make? To draw discrepancy from experiments involving the particle track, we have to argue in such a way that the observable predictions of the modified theory are in some way functions of the trajectory assumption. The question raised here is whether the de Broglie-Bohm particle law of motion can be made relevant to experiment. At first, it seems that definition of time spent by a particle within a classically forbidden barrier provides a good evidence for the preference of BQM. But, there are difficult technical questions, both theoretically and experimentally, that are still unsolved about this tunnelling times . A recent work indicates that it is not practically feasible to use tunnelling effect to distinguish between the two theories .

On the other hand, Englert et al.  and Scully  have claimed that in some cases Bohm’s approach gives results that disagree with those obtained from SQM and, in consequence, with experiment. Again, at first Dewdney et al.  and then Hiley et al.  showed that the specific objections raised by Englert et al. and Scully cannot be sustained. Furthermore, Hiley believes that no experiment can decide between the standard interpretation and Bohm’s interpretation. However, Vigier , in his recent work, has given a brief list of new experiments which suggest that the U(1) invariant massless photon assumed properties of light within the standard interpretation, are too restrictive and that the O(3) invariant massive photon causal de Broglie-Bohm interpretation of quantum mechanics, is now supported by experiments. Furthermore, in some of the recent investigation, some feasible experiments have been suggested to distinguish between SQM and BQM . In one work, Ghose indicated that although BQM is equivalent to SQM when averages of dynamical variables are taken over a Gibbs ensemble of Bohmian trajectories, the equivalence breaks down for ensembles built over clearly separated short intervals of time in specially entangled two-bosonic particle systems . Another one  is an extension of Ghose’s work to show disagreement between SQM and BQM in a two-particle system with an unentangled wavefunction, particularly at the statistical level333To clarify our discussion it is worth noting that in this paper we have used the following definitions:
* By statistically level we mean our final interference pattern.
\** The individual level refers to our experiment with a pair of particles which are emitted in clearly separated short intervals of time.. Further discussion of this subject can be found in \[13-15\]. In that experiment, to obtain a different interference pattern from SQM, we must deviate the source from its geometrically symmetric location.

In this investigation, we are offering a new thought experiment which can decide between SQM and BQM. Here, the deviation of the source from its geometrical symmetric location is not necessary and we have used a system consisting two correlated particles with an entangled wavefunction.

In the following section, we have introduced a two double-slit experimental set-up. In section 3, Bohm’s interpretation is used to find some observable results about our suggested experiment. Predictions of the standard interpretation and their comparison with Bohmian predictions is examined in section 4. In section 5, we have used selective detection and have compared SQM and BQM with our thought experiment at the ensemble level of particles, and we state our conclusion in section 6.

## 2 Two double-slit experiment presentation

To distinguish between SQM and BQM we consider the following experimental set-up. A pair of identical non-relativistic particles with total momentum zero labelled by 1 and 2, originate from a point source S that is placed exactly in the middle of a two double-slit screens as shown in Fig. 1. We assume that the intensity of the beam is so low that during any individual experiment we have only a single pair of particles passing through the slits and the detectors have the opportunity to relate together for performing selective detection process. In addition, we assume that the detection screens $`S_1`$ and $`S_2`$ register only those pairs of particles that reach the two screens simultaneously. Thus, we are sure that the registration of single particles is eliminated from final interference pattern. The detection process at the screens $`S_1`$ and $`S_2`$ may be nontrivial but they play no causal role in the basic phenomenon of the interference of particles waves . In the two-dimensional system of coordinates $`(x,y)`$ whose origin $`S`$ is shown, the center of slits lie at the points $`(\pm d,\pm Y)`$. The wave incident on the slits will be taken as a plane of the form

$$\psi _{in}(x_1,y_1;x_2,y_2;t)=ae^{i[k_x(x_1x_2)+k_y(y_1y_2)]}e^{iEt/\mathrm{}}$$
(1)
where $`a`$ is a constant and $`E=E_1+E_2=\mathrm{}^2(k_x^2+k_y^2)/m`$ is the total energy of the system of the two particles. The plane wave assumption comes from large distance between source $`S`$ and double-slit screens. To avoid the mathematical complexity of Fresnel diffraction at a sharp-edge slit, we suppose the slits have soft edges that generate waves having identical Gaussian profiles in the $`y`$-direction while the plane wave in the $`x`$-direction is unaffected . The instant at which the packets are formed will be taken as our zero of time. Therefore, the four waves emerging from the slits $`A`$, $`B`$, $`A^{^{}}`$ and $`B^{^{}}`$ are initially

$$\psi _{A,B}(x,y)=a(2\pi \sigma _0^2)^{1/4}e^{(\pm yY)^2/4\sigma _0^2}e^{i[k_x(xd)+k_y(\pm yY)]}$$
(2)
$$\psi _{A^{^{}},B^{^{}}}(x,y)=a(2\pi \sigma _0^2)^{1/4}e^{(\pm y+Y)^2/4\sigma _0^2}e^{i[k_x(x+d)+k_y(\pm y+Y)]}$$
(3)
where $`\sigma _0`$ is the half-width of each slit. At time $`t`$ the general total wavefunction at a space point $`(x,y)`$ of our considered system for bosonic and fermionic particles is given by

$`\psi (x_1,y_1;x_2,y_2;t)=N[\psi _A(x_1,y_1,t)\psi _B^{^{}}(x_2,y_2,t)\pm \psi _A(x_2,y_2,t)\psi _B^{^{}}(x_1,y_1,t)`$
$`+\psi _B(x_1,y_1,t)\psi _A^{^{}}(x_2,y_2,t)\pm \psi _B(x_2,y_2,t)\psi _A^{^{}}(x_1,y_1,t)]`$ (4)
with

$$\psi _{A,B}(x,y,t)=a(2\pi \sigma _t^2)^{1/4}e^{(\pm yYu_yt)^2/4\sigma _0\sigma _t}e^{i[k_x(xd)+k_y(\pm yYu_yt/2)E_xt/\mathrm{}]}$$
(5)
$$\psi _{A^{^{}},B^{^{}}}(x,y,t)=a(2\pi \sigma _t^2)^{1/4}e^{(\pm yYu_yt)^2/4\sigma _0\sigma _t}e^{i[k_x(x+d)+k_y(\pm yYu_yt/2)E_xt/\mathrm{}]}$$
(6)
where $`N`$ is a reparametrization constant that its value is unimportant in this paper and

$`\sigma _t=\sigma _0(1+{\displaystyle \frac{i\mathrm{}t}{2m\sigma _0^2}})`$ (7)
$`u_y={\displaystyle \frac{\mathrm{}k_y}{m}}`$
$`E_x={\displaystyle \frac{1}{2}}mu_x^2`$ (8)
where $`u_x`$ and $`u_y`$, according to BQM, are initial group velocities corresponding to each particle in the $`x`$\- and $`y`$-directions, respectively. In addition, the upper and lower sings in the total wavefunction refer to symmetric and anti-symmetric wavefunction under exchange of particle 1 to particle 2, corresponding to bosonic and fermionic property, while in the equations (5) and (6) they refer to upper and lower slits, respectively. In the next section, we have used BQM to drive some results of this experiment.

## 3 Bohmian predictions about the suggested experiment

In BQM, the complete description of a system is given by specifying the position of the particles in addition to their wavefunction which has the role of guiding the particles according to following guidance condition for $`n`$ particles, with masses $`m_1,m_2,\mathrm{},m_n`$

$$\stackrel{}{\dot{x}_i}(\stackrel{}{x},t)=\frac{1}{m_i}\stackrel{}{_i}S(\stackrel{}{x},t)=\frac{\mathrm{}}{m_i}Im\left(\frac{\stackrel{}{_i}\psi (\stackrel{}{x},t)}{\psi (\stackrel{}{x},t)}\right)$$
(9)
where $`\stackrel{}{x}=(\stackrel{}{x_1},\stackrel{}{x_2},\mathrm{},\stackrel{}{x_n})`$ and

$$\psi (\stackrel{}{x_1},\stackrel{}{x_2},\mathrm{},\stackrel{}{x_n};t)=R(\stackrel{}{x_1},\stackrel{}{x_2},\mathrm{},\stackrel{}{x_n};t)e^{iS(\stackrel{}{x_1},\stackrel{}{x_2},\mathrm{},\stackrel{}{x_n};t)/\mathrm{}}$$
(10)
is a solution of Schr$`\ddot{o}`$dinger’s wave equation. Thus, instead of SQM with indistinguishable particles, in BQM the path of particles or their individual histories distinguishes them and each one of them can be studied separately . In addition, Belousek  in his recent work, concluded that the problem of Bohmian mechanical particles being statistically (in)distinguishable is a matter of theory choice underdetermined by logic and experiment, and that such particles are in any case physically distinguishable. For our considered experiment, the speed of the particles 1 and 2 in the $`y`$-direction is given , respectively, by

$$\dot{y}_1(x_1,y_1;x_2,y_2;t)=\frac{\mathrm{}}{m}Im\frac{_{y_1}\psi (x_1,y_1;x_2,y_2;t)}{\psi (x_1,y_1;x_2,y_2;t)}$$
(11)
$$\dot{y}_2(x_1,y_1;x_2,y_2;t)=\frac{\mathrm{}}{m}Im\frac{_{y_2}\psi (x_1,y_1;x_2,y_2;t)}{\psi (x_1,y_1;x_2,y_2;t)}.$$
(12)
With the replacement of $`\psi (x_1,y_1;x_2,y_2;t)`$ from (4), (5) and (6), we have

$`\dot{y}_1`$ $`=`$ $`N{\displaystyle \frac{\mathrm{}}{m}}Im\{{\displaystyle \frac{1}{\psi }}[[2(y_1Yu_yt)/4\sigma _0\sigma _t+ik_y]\psi _{A_1}\psi _{B_2^{^{}}}`$ (13)
$`\pm `$ $`[2(y_1+Y+u_yt)/4\sigma _0\sigma _tik_y]\psi _{A_2}\psi _{B_1^{^{}}}`$ (14)
$`+`$ $`[2(y_1+Y+u_yt)/4\sigma _0\sigma _tik_y]\psi _{B_1}\psi _{A_2^{^{}}}`$ (15)
$`\pm `$ $`[2(y_1Yu_yt)/4\sigma _0\sigma _t+ik_y]\psi _{B_2}\psi _{A_1^{^{}}}]\}`$ (16)
$`\dot{y}_2`$ $`=`$ $`N{\displaystyle \frac{\mathrm{}}{m}}Im\{{\displaystyle \frac{1}{\psi }}[[2(y_2+Y+u_yt)/4\sigma _0\sigma _tik_y]\psi _{A_1}\psi _{B_2^{^{}}}`$ (17)
$`\pm `$ $`[2(y_2Yu_yt)/4\sigma _0\sigma _t+ik_y]\psi _{A_2}\psi _{B_1^{^{}}}`$ (18)
$`+`$ $`[2(y_2Yu_yt)/4\sigma _0\sigma _t+ik_y]\psi _{B_1}\psi _{A_2^{^{}}}`$ (19)
$`\pm `$ $`[2(y_2+Y+u_yt)/4\sigma _0\sigma _tik_y]\psi _{B_2}\psi _{A_1^{^{}}}]\}`$ (20)
where, for example, the short notation $`\psi _A(x_1,y_1,t)=\psi _{A_1}`$ is used. Furthermore, from (5) and (6) it is clear that

$`\psi _A(x_1,y_1,t)=\psi _B(x_1,y_1,t)`$ (21)
$`\psi _A(x_2,y_2,t)=\psi _B(x_2,y_2,t)`$ (22)
$`\psi _B^{^{}}(x_1,y_1,t)=\psi _A^{^{}}(x_1,y_1,t)`$ (23)
$`\psi _B^{^{}}(x_2,y_2,t)=\psi _A^{^{}}(x_2,y_2,t)`$ (24)
which indicates the reflection symmetry of $`\psi (x_1,y_1;x_2,y_2;t)`$ with respect to the $`x`$–axis. Utilizing this symmetry in (13) and (14), we can see that

$`\dot{y}_1(x_1,y_1;x_2,y_2;t)=\dot{y}_1(x_1,y_1;x_2,y_2;t)`$ (25)
$`\dot{y}_2(x_1,y_1;x_2,y_2;t)=\dot{y}_2(x_1,y_1;x_2,y_2;t)`$ (26)
which are valid for both bosonic and fermionic particles. Relations (16) show that if $`y_1(t)=y_2(t)=0`$, then the speed of each particles in the $`y`$-direction is zero along the symmetry axis $`x`$. This means that none of the particles can cross the $`x`$-axis nor are they tangent to it, provided both of them are simultaneously on this axis. Similar conclusions can be found in some other works \[for example, 8, 9, 11-13\]. It can be seen that there is the same symmetry of the velocity about the $`x`$-axis as for an ordinary double-slit experiment .

If we consider $`y=(y_1+y_2)/2`$ to be the vertical coordinate of the center of mass of the two particles, then we can write

$`\dot{y}`$ $`=`$ $`(\dot{y}_1+\dot{y}_2)/2`$ (27)
$`=`$ $`N{\displaystyle \frac{\mathrm{}}{2m}}Im\{{\displaystyle \frac{1}{\psi }}({\displaystyle \frac{y_1+y_2}{2\sigma _0\sigma _t}})(\psi _{A_1}\psi _{B_1^{^{}}}\pm \psi _{A_2}\psi _{B_1^{^{}}}+\psi _{B_1}\psi _{A_2^{^{}}}\pm \psi _{B_2}\psi _{A_1^{^{}}})\}`$ (28)
$`=`$ $`{\displaystyle \frac{(\mathrm{}/2m\sigma _0^2)^2}{1+(\mathrm{}/2m\sigma _0^2)^2t^2}}yt.`$ (29)
Solving the equation of motion (17), we obtain the path of the $`y`$-coordinate of the center of mass

$$y=y_0\sqrt{1+(\mathrm{}/2m\sigma _0^2)^2t^2}.$$
(30)
If at $`t=0`$ the center of mass of the two particles is exactly on the $`x`$-axis, then $`y_0=0`$, and the center of mass of the particles will always remain on the $`x`$-axis. Thus, according to BQM, the two particles will be detected at points symmetric with respect to the $`x`$-axis, as shown in Fig. 1.

It seems that calculation of quantum potential can give us another perspective of this experiment. As we know, to see the connection between the wave and particle, the Schr$`\ddot{o}`$dinger equation can be rewritten in the form of a generalized Hamilton-jacobi equation that has the form of the classical equation, apart from the extra term

$$Q(\stackrel{}{x},t)=\frac{\mathrm{}^2}{2m}\frac{^2R(\stackrel{}{x},t)}{R(\stackrel{}{x},t)}$$
(31)
where function Q has been called the quantum potential . But, it is clear that the calculation and analysis of Q, by using our total wavefunction (4), is not very simple. On the other hand, we can use the form of Newton’s second law, in which the particle is subject to a quantum force $`(\stackrel{}{}Q)`$ in addition to the classical force $`(\stackrel{}{}V)`$ , namely

$$\stackrel{}{F}=\stackrel{}{}(Q+V).$$
(32)
Now, if we utilize the equation of motion of the center of mass $`y`$-coordinate (18) and the equation (20), we shall obtain the quantum potential for the center of mass motion $`(Q_{cm})`$. Thus, we can write

$$\frac{Q}{x}=m\ddot{x}=0$$
(33)
$$\frac{Q}{y}=m\ddot{y}=\frac{my_0(\mathrm{}/2m\sigma _0^2)^2}{(1+(\mathrm{}t/2m\sigma _0^2)^2)^{3/2}}=\frac{my_0^4}{y^3}(\frac{\mathrm{}}{2m\sigma _0^2})^2$$
(34)
where the result of equation (21) is clearly due to motion of plane wave in the $`x`$-direction. In addition, we assume that $`V=0`$ in our experiment. Thus, our effective quantum potential is only a function of $`y`$-variable and it has the form

$$Q=\frac{my_0^4}{2y^2}(\frac{\mathrm{}}{2m\sigma _0^2})^2=\frac{1}{2}my_0^2\frac{\mathrm{}/2m\sigma _0^2)^2}{1+(\mathrm{}t/2m\sigma _0^2)^2}.$$
(35)
If $`y_0=0`$, the quantum potential for the center of mass of the two particles is zero at all times and it remains on the $`x`$-axis. However, if $`y_00`$, then the center of mass cannot touch or cross the $`x`$-axis. These conclusions are consistent with our earlier result (eq. (18)).

## 4 SQM forecast and its comparison with BQM

So far, we have been studying the results obtained from BQM at the individual level. But it is well known from SQM that the probability of simultaneous detection of two particles at $`y_M`$ and $`y_N`$, at the screen $`S_1`$ and $`S_2`$ is equal to

$$P_{12}(y_M,y_N)=_{y_M}^{y_M+\mathrm{}}𝑑y_1_{y_N}^{y_N+\mathrm{}}𝑑y_2|\psi (x_1,y_1;x_2,y_2;t)|^2.$$
(36)
The parameter $`\mathrm{\Delta }`$, which is taken to be small, is a measure of the size of the detectors. It is clear that the probabilistic prediction of SQM is in disagreement with the symmetrical prediction of BQM, because SQM predicts that probability of asymmetrical detection at the individual level of pair of particles can be different from zero, in opposition to BQM’s symmetrical predictions. In addition, based on SQM’s prediction, the probability of finding two particles at one side of the $`x`$-axis can be nonzero while we showed that BQM’s prediction forbids such events in our experiment. In other words, its probability must be exactly zero. Thus, if we provide necessary arrangements to perform this experiment, we must abandon one of the two theories or even both as a complete description of the universe.

Now the question arises as to whether this difference persists if we deal with an ensemble of pair of particles? To answer this question, we consider an ensemble of pair of particles that have arrived at the detection screens $`S_1`$ and $`S_2`$ at different times $`t_i`$. It is well known that, in order to ensure the compatibility between SQM and BQM for ensemble of particles, Bohm added a further postulate to his three basic and consistent postulates . Based on this further postulate, the probability that a particle in the ensemble lies between $`\stackrel{}{x}`$ and $`\stackrel{}{x}+d\stackrel{}{x}`$ at time $`t`$ is given by

$$P(\stackrel{}{x},t)=R^2(\stackrel{}{x},t).$$
(37)
Thus, using BQM, the probability of simultaneous detection for all pairs of particles of the ensemble arriving at the two screens at different instant of time $`t_i`$, with $`y_0=0`$, is

$`P_{12}=\underset{N\mathrm{}}{lim}{\displaystyle \underset{i=1}{\overset{N}{}}}R^2(y_1(t_i),y_1(t_i),t_i){\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑y_1{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑y_2|\psi (y_1,y_2,t)|^2=1`$ (38)
where every term in the sum shows only one pair arriving on the screens $`S_1`$ and $`S_2`$ at the point $`(y_1(t_i),y_1(t_i))`$ and time $`t_i`$, weighted by the corresponding density $`R^2`$. If all times of $`t_i`$ are taken to be $`t`$, then the summation on $`i`$ can be changed to an integral over all paths that cross the screens $`S_1`$ and $`S_2`$ at that time. Now, we can consider the joint probability of two points $`y_M`$ and $`y_N`$ on the two screens at time $`t`$ that are not symmetric about the $`x`$-axis, but we know that they are not detected simultaneously. Then, one can obtain the probability of detecting two particles at two arbitrary points $`y_M`$ and $`y_N`$

$$P_{12}(y_M,y_N)=_{y_M}^{y_M+\mathrm{\Delta }}𝑑y_1_{y_N}^{y_N+\mathrm{\Delta }}𝑑y_2|\psi (y_1,y_2,t)|^2$$
(39)
which is similar to the prediction of SQM (eq. (24)) but obtained in a Bohmian way . Thus, it appears that for a geometrically symmetric arrangement, the possibility of distinguishing the two theories at the statistical level is impossible, as was expected \[1-3, 9, 17\].

## 5 Selective detection and comparison of SQM with BQM at the statistical level

In the previous section, we have shown that SQM and BQM have different predictions for our suggested experiment, at the individual level. Since SQM talks about individual events in probabilistic terms, the existence of different predictions by the two theories at the individual level is not a strange result. On the other hand, we have seen that the two theories, for a geometrically symmetric arrangement, are consistent at the ensemble level. Here, one can ask whether the individual level is the only area to distinguish between the two theories and whether the disagreement between them cannot appear at the ensemble level. In this section, we answer this question in negative, and we shall provide conditions under which SQM can be interpreted as a vague theory at the ensemble level.

### 5.1 The case where $`y_0`$ is exactly zero

We have seen that the assumption of $`y_0=0`$ is not in contradiction with the statistically results of SQM and, in consequence, with experiment. Thus, we can assume that initially each particle in the source is statistically distributed according to the absolute square of the wavefunction, but, this distribution is completely symmetric so that the $`y`$-coordinate of the center of mass is on the $`x`$-axis. If we can prepare such a special source with two correlated particles, we can try to do our experiment in the following fashion: particles emitted from the source $`S`$ into the right hand side of the experimental set-up can pass through slits $`A`$ or $`B`$. Since we have assumed that the total momentum of the pair of particles is zero, if one of the particles goes through the slit $`A`$ for instance, the other particle must go through the slit $`B^{^{}}`$ on the left hand side of the experimental set-up. Based on BQM and using equation (16), we infer that the particle passing through the slit $`A`$ must be detected on the upper half plane of the $`x`$-axis on the screen $`S_1`$. The same thing must occur for the other particle that passes through $`B^{^{}}`$, but in the lower half plane of the $`x`$-axis on the $`S_2`$ screen. Using this prediction, we assume that only those particles arriving at $`S_2`$ are recorded for which there is a simultaneous detection of the other particle at the upper side on $`S_1`$. We called this special detection, in which some of the selected particles are recorded, a selective detection. Thus, based on the prediction power of BQM, we will record two particles symmetric with respect to the $`x`$-axis, for each emitted pairs of particles. If we wait to record an ensemble of particles, we will see an interference pattern of particles on the lower half plane of the $`S_2`$ screen. On the other hand, based on SQM, the probability of finding a particle at any point on the $`S_2`$ screen, even at the upper side, is nonzero and there is no compulsion to detect pairs of particles symmetrically on the two sides of the $`x`$-axis, as it can be seen from equation (24) and is depicted in Fig. 1. Therefore, if we accept that SQM is still efficient and unambiguous for the selective detection, the interference pattern will be seen on the whole screen $`S_2`$, particularly at the upper side of it, at the ensemble level. Consequently, we shall have observable results to distinguish the two theories, SQM and BQM.

### 5.2 The case where $`y_0`$ is statistically distributed

One can argue that $`y_0`$ cannot have a well-defined position and it must be distributed according to Born’s principle. However, we will show that this objection cannot alter our obtained results. Assume that, $`y_0=0`$ but $`\mathrm{}y_00`$. If we provide conditions in which $`\mathrm{}y_0`$ is very small and $`\mathrm{}t/2m\sigma _0^21`$, we can still detect particles symmetrically with respect to the $`x`$-axis, with a good approximation. To obtain symmetrical detection about the $`x`$-axis with reasonable approximation, the center of mass variation from the $`x`$-axis must be smaller than the distance between any two neighboring maxima, that is

$$y\frac{\lambda D}{2Y}\frac{\pi \mathrm{}t}{Ym}$$
(40)
where $`\lambda `$ is the de Broglie wavelength. For conditions $`\mathrm{}t/2m\sigma _0^21`$, $`Y\sigma _0`$ and using equation (18), one can obtain

$$y_0\frac{\pi \mathrm{}t}{Ym}\sigma _0.$$
(41)
Therefore, if we use a source with $`\mathrm{}y_0\sigma _0`$, we shall obtain $`yy_0\sigma _0`$ for each individual observation, and our symmetrical detection can be maintained with a good approximation. It is evident that, if one considers $`\mathrm{}y_0\sigma _0`$, as was done in , the incompatibility between the two theories will be disappeared. But, we believe that, instead of the usual one-particle two-slit experiment with $`\mathrm{}y_0\sigma _0`$, our correlated two-particle system provides a new situation in which we can adjust $`y_0`$ independent of $`\sigma _0`$, so that

$$y_0=\frac{1}{2}(y_1+y_2)_{t=0}\sigma _0.$$
(42)
Although it is obvious that $`(\mathrm{}y_1)_{t=0}=(\mathrm{}y_2)_{t=0}\sigma _0`$, but position correlation between the two entangled particles cause that they always satisfy equation (30). Furthermore, if it is assumed that $`y_0`$ is statistically distributed, another problem can be raised, which is mentioned by Marchildon . We have shown that, if both particles are simultaneously on the $`x`$-axis, both velocities in the $`y`$-direction vanish, and neither particle could cross or be tangent to the $`x`$-axis. However, under $`\mathrm{}y_00`$ condition, pairs of particles cannot be simultaneously on the $`x`$-axis and we have not the aforementioned constraint on the motion of particles (relations (16)). However, using our selective detection, we can still obtain our last result, because the center of mass of the two particles are on the $`x`$-axis with reasonable approximation. It is clear that, under such a condition, one cannot claim that the particle detected at the upper (lower) side must have passed through upper (lower) slit, in spite of $`y_0=0`$ condition. To confirm these results, it is worth noting that, Durr et al.  argue that the selective detection can alter the statistical predictions of the two theories:“note that by selectively forgetting results we can dramatically alter the statistics of those that we have not forgotten. This is a striking illustration of the way in which Bohmian mechanics does not merely agree with the quantum formalism, but, eliminating ambiguities, clarifies, and sharpens it.”. Elsewhere , we have utilized another kind of selective detection by which we could alter statistical prediction of SQM, using BQM for an interference device that contains two unentangled particles.

## 6 Conclusion

In this investigation, we have suggested an experiment to distinguish between SQM and BQM. In fact, we believe that some particular experiments that for one reason or another have not yet been performed can decide between them. Thus, it has been shown that a two double-slit experiment set-up, along with a source of two identical non-relativistic particles with total momentum of zero, emitted at suitable time intervals, has the following characteristics:
1) The suggested experiment will yield different observable predictions for SQM and BQM, at the individual level.
2) The two theories yield the same interference pattern at the ensemble level without using a selective detection, as is expected.
3) Since in BQM the particles are distinguishable and their past history are known, using selective detection, it has been shown that either the two theories will predict different results at the statistical level, or that BQM has more predictive power than SQM. It is shown that selective detection can be considered as a tool for arriving to a new realm in which trajectory interpretation is sharply formulated, while the standard interpretation is ambiguous and silent even at the ensemble level.

Therefore, it seems possible to distinguish between the two theories and to see whether BQM is a worthy successor to SQM.