# Exact dynamics and decoherence of two cold bosons in a 1D harmonic trap

###### Abstract

We study dynamics of two interacting ultra cold Bose atoms in a harmonic oscillator potential in one spatial dimension. Making use of the exact solution of the eigenvalue problem of a particle in the delta-like potential we study time evolution of initially separable state of two particles. The corresponding time dependent single particle density matrix is obtained and diagonalized and single particle orbitals are found. This allows to study decoherence as well as creation of entanglement during the dynamics. The evolution of the orbital corresponding to the largest eigenvalue is then compared to the evolution according to the Gross-Pitaevskii equation. We show that if initially the center of mass and relative degrees of freedom are entangled then the Gross-Pitaevskii equation fails to reproduce the exact dynamics and entanglement is produced dynamically. We stress that predictions of our study can be verified experimentally in an optical lattice in the low-tunneling limit.

###### pacs:

03.75Kk, 03.75.Gg, 67.85.-d## I Introduction

Theoretical description of a Bose-Einstein condensate of trapped weakly interacting atomic system is traditionally based on a mean field approximation GPS (). By assuming that many-body wave function can be written in a form of N-fold product state, i.e. that all atoms occupy the same single particle orbital, the stationary Gross-Pitaevskii (GP) equation for the order parameter is found. Assuming further that the N-fold product approximation holds also in dynamical situations one arrives at the time dependent Gross-Pitaevskii equation. Under most of experimental conditions there are no strong correlations in the system and the GP equation turned-out to be extremely fruitful in predicting and describing a variety of features of those systems. Soon it occurred that also high energy solutions of the GP equation can be useful in studying Bose systems at finite temperatures. The GP equation has become a work horse of the theory of weakly interacting ultra cold bosons.

On the other hand examples when the mean field description does not reproduce the real dynamics have been studied. For instance direct comparison of the mean field and many body theory of vortex nucleation showed that the GP equation fails to describe this phenomenon Maciek1 (); Maciek2 (); Maciek3 (). Similarly a mean field description of attractive Bose systems encounters difficulties Gajda (); Gajda2 (); Hedelberg (). Due to permanent progress in experimental techniques the physics of ultra cold atomic gases started to penetrate areas traditionally associated with condensed matter physics where correlations play a crucial role. Evidently, in such situations simple mean field description based on the GP equation becomes insufficient. The Mott insulator-superfluid transition Bloch () or the Tonks-Girardeau gas Tonks (); Tonks2 () are some examples.

It is commonly believed that creation of the Mott insulator with a small and controlled number of atoms per lattice site allows for applications of such systems in quantum information. All quantum information processing is based on the quantum correlations which cannot be described by any classical theory based on local realism. States showing these non-local correlations are known as entangled states. Two entangled spin– particles are the essence of the Einstein, Podolsky, and Rosen paradox EPR (). The entangled states vialote Bell inequalities Bell (). It has been realized that also continuous variable entangled systems can be used in quantum computations qcomp1 (); qcomp2 (); qcomp3 (); qcomp4 (); qcomp5 (); qcomp6 (); qcomp7 (). Several authors Entropy1 (); Entropy2 (); Entropy3 () studied recently creation of entanglement during dynamics of two interacting particles in a harmonic trap. In particular it has been shown that this dynamically entangled state violates a Bell-type inequality for a certain choice of observables Entropy1 ().

In this paper we study dynamically created entanglement, and its measure – the von Neuman entropy, for a realistic system of two identical bosonic atoms in a harmonic trap. We consider low energy collisions of the atoms. At such energies the range of van der Waals interactions is smaller than the s-wave scattering length. Therefore, the interaction potential can be approximated by a contact pseudo potential. This approximation occurred to be in excellent agreement with experimental results Esslinger () where binding energy of molecular system has been measured. The molecules were created from atoms in an optical lattice in the limit of small tunneling. This experimental arrangement is perfectly suited for a study of exact dynamics of two trapped atoms. We consider a realistic case of two atoms per lattice site deep in the Mott insulator phase. By applying a Bragg pulses Phillips () one creates a state in which each atom is in a superposition of two counter propagating wave packets. Initially, such a state is a two-fold product state of two identical wave functions, i.e. is separable. Each wave function has two components moving initially with opposite momenta. The center of mass dynamics is generated by a different Hamiltonian than dynamics of the relative coordinate therefore this two particle continuous variable system becomes entangled. We use the von Neuman entropy as a measure of the entanglement and study its behavior in time for various interaction strengths. We also analyze a coherence of the system which is directly related to the largest eigenvalue of the single particle density matrix and compare the exact dynamics to the mean field description based on the GP equation.

## Ii Two bosons in a harmonic trap

We are going to study dynamics of the simplest nontrivial system – two atoms confined in a one dimensional harmonic potential. In fact generalization of our results to two or three spatial dimensions is straightforward. We limit our analysis to the 1D case as this situation captures all features of the dynamics. For simplicity we are using harmonic-oscillator units. It means that all energies are measured in , all lengths in , and all momenta in . Hamiltonian of the system of two interacting bosons in the harmonic trap has the form:

(1) |

where and are positions of atoms interacting via a short range potential modeled by the delta function. This form of the short range interaction is justified in the limit of vanishing relative velocity of colliding atoms, where atomic de Broglie wavelength is much larger then a range of two body potential. In 2D and 3D the corresponding Hamiltonian is not a self adjoint operator. To correct for this fact a regularization is required. In contrast to many dimensions, the regularization of the delta function is not necessary in one dimensional case Se86 (). In 1D the parameter is given by , where is a scattering length Al88 (). It is worth to notice that finite range interactions between particles modeled by the Gaussian function in the context of the dynamics of two bosons was studied in Sascha ().

To demonstrate entanglement formation we study the evolution of two bosons which initially are in a product quantum state

(2) |

Function is a one-particle wave function called the order parameter in the mean field context.

The exact dynamics of the two interacting bosons in the harmonic trap can be found because all eigenstates of the full two-body Hamiltonian (1) are known. They are found in Bu98 (). The two particle problem has to be first brought to a single particle one by introducing the center of mass and the relative coordinates:

(3a) | ||||

(3b) |

In these coordinates Hamiltonian (1) separates into two independent parts – the center of mass part, and the relative part:

(4a) | ||||

(4b) |

As one can see, the dynamics of the center of mass is described by the standard one dimensional harmonic oscillator Hamiltonian (4a). Its eigenstates are well known and have a standard form

(5a) | |||

where are Hermite polynomials. The energy of -th eigenstate in our units is obviously given by | |||

(5b) |

The eigenstates of the Hamiltonian (4b) describing relative dynamics of two particles are given in Bu98 () and for one dimensional problem have a form

(6a) | |||||

(6b) |

where are confluent hypergeometric functions, and are normalization coefficients. Since the wave function of identical bosons must be symmetric under exchange of the two particles, therefore the physical wave function is composed from functions with even only. The energies of these even states are given by a sequence of zeros of the function:

(7) |

The quantum number is equal to . The initial wave function can be easily decomposed to the superposition of the eigenstates of the Hamiltonian:

(8) |

Obviously the evolution of the initial two boson state is given by:

(9) |

The last step is to return to the original coordinates by using relations (3):

(10) |

Standard method of detection of ultra cold trapped atomic systems are destructive. The optical lattice potential is turned off and the system is allowed to expand balistically. Only after expansion a size of the system exceeds a resolution of a CCD camera. The picture of the CCD camera gives therefore direct insight into the initial momentum distribution of atoms. The wave function Eq.(10) written in the momentum space of the two atoms is:

(11) |

In repeated single particle detections preceded by the ballistic expansion of the system one-particle momentum distribution is monitored:

(12) |

where is the reduced one particle density matrix in the momentum representation:

(13) |

By making its spectral decomposition we can determine the number of orbitals and their relative occupations needed for accurate description of the two bosons dynamics. Time dependence of the eigenvalues of the density matrix is discussed below. Let us mention that the largest eigenvalue is a direct measure of the coherence of the system.

We shall compare this exact dynamics with the approximate one governed by the Gross-Pitaevskii equation. The main idea leading to the mean field approximation relies on the assumption that generation of entanglement between bosons during the evolution is negligible and therefore the quantum state of the system remains separable. In other words all correlations between bosons are neglected and the same wave functions of every particle is assumed during the evolution:

(14) |

This assumption leads directly to the Gross-Pitaevskii equation which determines the dynamics of the one-particle wave function :

(15) |

The probability density in momentum space reads:

(16) |

where is the Fourier transform of the time dependent solution of the GP equation, . We compare the exact one-particle momentum distribution with that predicted by the Gross-Pitaevskii approximation (16). Moreover, in the situation when many eigenvalues of the density matrix (13) are of the same order we can also compare the Gross-Pitaevskii momentum distribution (16) with the momentum distribution of the dominant orbital obtained from diagonalization of the exact one-particle density matrix in the momentum space. Obviously, in a general case the GP solution overestimates the coherence of the system.

The GP equation is solved numerically on a spatial grid of points separated by . The time step is equal to . We use the split-operator method which is very stable for the chosen temporal and spatial steps.

## Iii Results

To make the detailed comparison we concentrate on a one particular class of the initial states. We assume that at the beginning each particle is in the state described by the Schrödinger cat like wave function

(17) |

Parameter measures the separation between two wave packets moving in the opposite direction in the relative coordinates space. Such a choice is motivated by the preparation procedure described above, i.e. we assume that Bragg pulses bringing the atoms into the superposition of wave packets moving in opposite directions are applied. When then the initial state is very close to the ground state of the system so we expect that the exact dynamics is almost indistinguishable from the dynamics in the mean filed approximation. For large the initial state is still separable but it is highly delocalized. Relative and center of mass degrees of freedom are entangled in the initial state. They evolve in a different way, therefore we expect that the exact dynamics could be dramatically different than the dynamics predicted by a simple mean field approach.

### iii.1 Dependence on delocalization of one-particle state

First let us discuss situation for generic interaction strength () when , i.e. when the extension of the initial state is equal to the trap length unit. We observe that the single particle density matrix obtained from the exact dynamics develops more then one nonzero eigenvalue, i.e. many one particle orbitals are involved. Fig. 1a shows time dependence of the eigenvalues of the one-particle density matrix (13). Because one of the eigenvalues is incessantly much larger than the others the system coherence is large and the Gross-Pitaevskii description is quite accurate in this case. Time dependence of the momentum distributions deduced from the Gross-Pitaevskii equation and the exact solution are shown in Fig. 1b and they are in agreement with our predictions (whole time dependence of momentum distributions is available on-line movies ()).

Situation changes dramatically when we increase the delocalization parameter. When is large enough then a few orbitals can play the crucial role in the dynamics and the mean field approximation is no longer valid. Fig. 2a shows the time dependence of the eigenvalues of the density matrix for . As we see, the main orbital (its eigenvalue is represented by a thick solid line) initially dominates. But after a few periods of the trap oscillations the other orbital becomes much more important than the first one. The dynamics is obviously much more complicated than it is predicted by the mean field approach. It is clear when we compare the momentum density distribution predicted by the exact and the mean field solutions (Fig. 2b and movie available on-line movies ()).

We see that evidently Gross-Pitaevskii equation properly describes the dynamics of the first orbital rather then the whole system, Fig. 2b. It is the reason why the Gross-Pitaevskii equation gives good predictions when only one eigenvalue of the one particle matrix dominates during the entire evolution.

It is also interesting to study similarities and differences between prediction of the mean field approach and exact solution in the situation when the initial state of each particle is antisymmetric in position space, i.e. is described by the wave function of the form

(18) |

Nevertheless the wave function of the system is still symmetric under particle exchange. Since corresponding GP Hamiltonian is invariant under reflection , therefore the symmetry of the initial state will be preserved. As we observe, it is not true for the exact two body dynamics. The evolution preserves only the symmetry of each orbital separately, but not the symmetry of the whole system. It is clearly demonstrated in Fig. 3 where we compare the momentum distribution predicted by the mean field approach with the single particle density obtained from the exact dynamics. Time dependence of the eigenvalues of one-particle density matrix is identical as the one for the corresponding symmetric case (Fig. 2a). Whole movie is also available on-line movies ().

### iii.2 Dependence on interaction strength

Now we want to show that correctness of the mean field approximation significantly depends on the interaction strength parameter . It is quite obvious that in the situation when the interaction is switched off, the two particles initially in the state which is not entangled (product state) will stay in such a state during the whole evolution even for a highly delocalized state. In this case the mean field approximation naturally leads to the same solution as the exact solution. It is the interparticle interaction which can produce entangled two body states during the evolution.

Time dependence of the eigenvalues for a moderate interaction strength () is presented in Fig. 1 and 2. In those situations only two eigenvalues (i.e. two orbitals) are important for many trap periods. For stronger interactions this picture changes significantly. Time evolution of eigenvalues for strong interaction () and are shown in Fig. 4. After a few trap periods many different orbitals become important. Moreover the orbital which dominates at the beginning becomes unimportant after a very short time. Therefore we do not expect that the Gross-Pitaevskii approximation may give correct predictions in this case.

On the other hand when the interaction is very weak we can expect that the production of entanglement will be very slow even for highly delocalized states and therefore the mean field approximation may be correct for a long evolution time. Time dependence of the eigenvalues of the one-particle density matrix when the interaction is weak but the initial state is highly delocalized is presented in fig. 5.

### iii.3 Revivals of product states

Looking at fig. 2 and fig. 5 one can observe that as initially only one eigenvalue dominates in the Schmidt decomposition of the single particle density matrix, the other eigenvalues become more important at later times. However, the time dependence of the dominant eigenvalue exhibits some oscillations and a partial revival of the ‘initial’ eigenvalue can be observed. It is interesting to find a physical explanation of this behavior. To this end in fig. 6 we show the spectrum of the two-particle state for the two studied parameter sets and (as in fig. 2) and and (as in fig. 5). In this figure we plot the probability of the given eigenenergy, , resulting from the decomposition (8).

First let us notice that the eigenenergies do appear ‘in pairs’. The effect of ‘pairing’ of eigenenergies can be easily explained. For given and the eigenenergy has two components. is the energy of the center of mass while corresponds to the relative coordinate. As was mentioned before, energies of the relative excitations are very close to the energies of harmonic trap. This is because the potential in the relative coordinate space is the harmonic potential of the trap modified at by the presence of delta function. This delta function shifts the harmonic energy by very small amount. Therefore state labeled by is almost degenerate with the state labeled by and therefore spectrum is paired.

If by we denote the difference between two energies of dominant pair (maximum in the spectral decomposition) it is clear that after time these two the most important amplitudes of the two-atom wavefunction match in phase and the partial revival of the product state can be observed. This is signified by a reappearance of the initial eigenvalue of the single particle matrix. The revival time calculated this way for parameters of Fig. 2 is equal to while for Fig. 4 is which agrees perfectly with predictions of the exact dynamics.

Obviously revival time depends on the initial state as well as eigenmodes of the Hamiltionian (1).

In addition to this large time oscillations of the eigenvalues of the density matrix some small fast oscillations can be also observed. This fast modulations appear every half of the trap period when two wavepackets meet at the trap center and results from interaction between them.

### iii.4 Entanglement of particles

Mutual interactions between particles obviously leads to the quantum correlations between particles. To study them we use the correlation measure introduced in Kazik ():

(19) |

where are the eigenvalues of the one-particle density matrix . This measure has very simple interpretation. It gives an effective number of single particle orbitals occupied in the given many body state. In particular when one-particle density matrix has equal eigenvalues then .

Other commonly used Entropy1 (); Entropy2 (); Entropy3 (); Entropy0 () measure of entanglement in the system is von Neumann entropy defined as

(20) |

This entropy is even more interesting than the number of dominant eigenvalues since it is directly connected with the entropy defined in thermodynamical context. Time dependence of this two measures of entanglement in the system for two different regimes of interaction strength are presented in Fig. 7. Obviously in the beginning, when the system is in separable state, entanglement and von Neumann entropy are equal to and respectively. We observe that correlation and entropy increase in time and seem to saturate for large time. Even though they have different physical interpretation they behave very similarly which might seem quite surprising. They reach ‘stationary regime’ faster for stronger interactions.

Both quantities exhibit fast oscillations modulated by a slowly varying functions. These fast oscillations are related to partial revival of dominant eigenvalue discussed in previous subsection. Every minimum observed in correlation function corresponds to the moment when there is a dominant eigenvalue in the Schmidt decomposition of the one-particle density matrix. Let us remind that this revivals are related to phase matching of two dominant eigenmodes of the two particle state.

Long time modulations of correlation functions are related to the quantum nature of the system and discreetness of the energy spectrum. In such cases evolution is always quasi-periodic and due to the interference of amplitudes long time scale oscillations do appear. In our case the number of modes with no zero amplitudes is relatively small and therefore oscillations of correlation functions appear on a time scale of few hundred trap periods.

## Iv Summary

In this paper we study the exact dynamics of two particles trapped in a harmonic trap and interacting by a contact potential. We assumed that initially each particle is transferred by the Bragg pulses to the state being the superposition of two wave packets moving in opposite directions. We show that the two particle state, although initially being a product state does not preserve the product form during the evolution. The reason is that the initial state entangles the center of mass and relative coordinates of the two particle system. These two degrees of freedom evolve according to different Hamiltonians. As a result the single particle reduced matrix develops many eigenvalues during the evolution what signifies decreasing coherence of the system. This situation cannot be correctly described by the GP equation. Our predictions can be verified in the experiment with deep optical lattices when two atoms occupy each site. We show one-particle momenta distributions for different initial states and compare them to those obtained from the mean-field dynamics. The differences between the two signify the two atom entanglement. The momentum distribution is directly measured by exposure of the system to a resonant light after ballistic expansion and therefore creation of entanglement in the two particle system can be easily traced in time and compared to exact solutions. We monitored the von Neumann entropy which is common measure of entanglement. We show that entanglement is dynamically created during evolution, however it is not very surprising for interacting system. A comment about a system of two fermionic particles would be also in place. As two identical fermions do not interact in the s-wave channel, as long as other partial waves can be neglected, their dynamic is driven by the noninteracting Hamiltonian. On the other hand, if the spatial part of a wave function of two fermions is symmetric and a spin part is responsible for the antisymmetrization of the total wave function, then our exact solution evidently applies to such a situation.

###### Acknowledgements.

We acknowledge the support of the Polish Ministry of Science and Education grants for 2008-2010 (T.S.), 2009-2011 (M.B., K.R.), and EU project NAME-QUAM (M.G.).## References

- (1) F. Dalfovo, S. Giorgini, L. Pitaevskii, S. Stringari, Rev. Mod. Phys. 71, 463512 (1999)
- (2) D. Dagnino, N. Barberán, M. Lewenstein, and J. Dalibard, Nature Physics 5, 431 (2009)
- (3) D. Dagnino, N. Barberán, and M. Lewenstein, Phys. Rev. A 80, 053611 (2009)
- (4) B. Juliá-Díaz,1 D. Dagnino, M. Lewenstein, J. Martorell, and A. Polls, Phys. Rev. A 81, 023615 (2010)
- (5) M. Gajda, Phys. Rev. A 73, 023603 (2006)
- (6) M. Gajda, M. Zaluska-Kotur, J. Mostowski, Optics Express 8, 106 (2001)
- (7) E. Alon, Alexej I. Streltsov, and Lorenz S. Cederbaum, Phys. Rev. A 77, 033613 (2008)
- (8) M. Greiner, O. Mandel, T. Esslinger, T.W. Hänsch, and I. Bloch, Nature 415, 39 (2002)
- (9) B. Paredes, A. Widera, V. Murg, O. Mandel, S. Fölling, I. Cirac, G. V. Shlyapnikov, T. W. Hänsch, and I. Bloch, Nature 429, 277 (2004)
- (10) T. Kinoshita, T. Wenger, and D. S. Weiss, Science 305, 1125 (2004)
- (11) A. Einstein, B. Podolsky, and N. Rosen, Phys. Rev. 47, 777780 (1935)
- (12) J. S. Bell, Rev. Mod. Phys. 38, 447 (1966)
- (13) L. Vaidman, Phys. Rev. A 49, 1473 (1994)
- (14) S. L. Braunstein, H. J. Kimble, Phys. Rev. Lett. 80, 869 (1998)
- (15) P. van Loock, S. L. Braunstein, H. J. Kimble, Phys. Rev. A 62, 022309 (2000)
- (16) P. van Loock, S.L. Braunstein, Phys. Rev. Lett. 84, 3482 (2000)
- (17) P. van Loock, S.L. Braunstein, Phys. Rev. Lett. 87, 247901 (2001)
- (18) S. L. Braunstein, Nature (London) 394, 47 (1998)
- (19) A. Furusawa, J.L. Sørensen, S.L. Braunstein, C.A. Fuchs, H. J. Kimble and E. S. Polzik, Science 282, 706 (1998)
- (20) H. Mack, M. Freyberger, Phys. Rev. A 66, 042113 (2002)
- (21) B. Sun, D. L. Zhou, L. You, Phys. Rev. A 73, 012336 (2006)
- (22) J. Goold, L. Heaney, T. Busch, V. Vedral, Phys. Rev. A 80, 022338 (2009)
- (23) T. Stöferle, H. Moritz, K. Günter, M. Khöl, T. Esslinger, Phys. Rev. Lett. 96, 030401 (2006)
- (24) M. Kozuma, L. Deng, E. W. Hagley, J. Wen, R. Lutwak, K. Helmerson, S. L. Rolston, and W. D. Phillips, Phys. Rev. Lett. 82, 871 (1999)
- (25) P. Šeba, Czech. J. Phys. 36, 559 (1986)
- (26) S. Albeverio, F. Gesztesy, R. Høegh-Krohn, and H. Holden, Solvable Models in Quantum Mechanics (Springer-Verlag, New York, 1988).
- (27) Ch. Matthies, S. Zöllner, H.-D. Meyer, P. Schmelcher, Phys. Rev. A 76, 023602 (2007)
- (28) T. Busch, B. Englert, K. Rza̧żewski, M. Wilkens, Foundations of Physics 28, 549 (1998)
- (29) See supplementary material at [http://link.aps.org/supplemental/10.1103/PhysRevA.82.053631] for movies presenting time dependence of momentum distributions predicted by exact solution and in the mean-field approximation.
- (30) R. Grobe, K. Rza̧żewski, J. H. Eberly, J. Phys. B 27, L503 (1994)
- (31) R. Paškauskas, L. You, Phys. Rev. A 64, 042310 (2001)