Bosonic Quantum Dynamics Following Colliding Potential Wells
BBosonic Quantum Dynamics Following Colliding Potential Wells
Fabian K¨ohler ∗ Center for Optical Quantum Technologies,Department of Physics, University of Hamburg,Luruper Chaussee 149, 22761 Hamburg Germany
Peter Schmelcher † Center for Optical Quantum Technologies,Department of Physics, University of Hamburg,Luruper Chaussee 149, 22761 Hamburg Germany andThe Hamburg Centre for Ultrafast Imaging, University of Hamburg,Luruper Chaussee 149, 22761 Hamburg, Germany (Dated: February 3, 2021) a r X i v : . [ qu a n t - ph ] F e b bstract We employ the multi-configuration time-dependent Hartree method for bosons (MCTDHB) inorder to investigate the correlated non-equilibrium quantum dynamics of two bosons confined intwo colliding and uniformly accelerated Gaussian wells. As the wells approach each other aneffective, transient double-well structure is formed. This induces a transient and oscillatory over-barrier transport. We monitor both the amplitude of the intra-well dipole mode in the courseof the dynamics as well as the final distribution of the particles between the two wells. For fastcollisions we observe an emission process which we attribute to two distinct mechanisms. Energytransfer processes lead to an untrapped fraction of bosons and a resonant enhancement of thedeconfinement for certain kinematic configurations can be observed. Despite the comparativelyweak interaction strengths employed in this work, we identify strong inter-particle correlations byanalyzing the corresponding Von Neumann entropy.
I. INTRODUCTION
Ever since the first realizations of Bose-Einstein condensates (BEC)[1–3], ultracold quan-tum gases were in the focus of experimental and theoretical research in quantum physics.Their nearly perfect isolation from the environment as well as their excellent tunabilityrender them ideal platforms to simulate a wide variety of quantum many-body systems[4–6]in order to unravel their fundamental physical properties. Experimental advancements inrecent years enabled the study of ensembles of ultracold atoms with a controlled num-ber of particles[7, 8] confined in almost arbitrarily shaped external potentials[9] like op-tical lattices[10, 11], harmonic traps[12] and ring traps[13]. By varying the confinementit is possible to realize effectively three-dimensional[14, 15], two-dimensional[16, 17] andone-dimensional[18, 19] systems. Magnetic Feshbach[20, 21] and confinement-inducedresonances[22–25] provide fine-grained control of the inter-particle interaction. Recentstudies have employed this versatile toolbox of ultracold atoms to establish links to solid-state systems[26, 27], the electronic structure of molecules[28], light-matter interaction[29],topological matter[30, 31], and even black-hole analogs[32]. ∗ [email protected] † [email protected]
2n recent years, optical tweezers have become important instruments to confine andmove microscopic objects by exerting small forces via highly focused laser beams. Thistool was originally developed to manipulate micrometer-sized particles[33, 34] but has beenlater refined to manipulate objects on many different length scales ranging from individualatoms[35, 36] to bacteria and viruses[37]. These advancements sparked strong interest touse optical tweezers for the precise manipulation of ensembles of ultracold neutral atoms[38]including Rydberg atoms[39–41]. A very interesting direction of research is to use multipleoptical tweezers to accelerate atomic clouds[42] which allows to set up optical colliders[43–45]. In these experiments, fundamental properties of quantum scattering processes wereobserved such as partial wave interference or the loss of particles on resonant collisions. Inthis light, colliding ultracold atoms could be used to mimic electrons during atom-atomcollisions. Since the dynamics of ultracold atoms take place on much larger time scales, theusually very fast electronic processes could be slowed down[29, 46, 47], potentially providingin depth insights into the fundamental processes of atom-atom or atom-ion collisions suchas projectile ionization[48, 49] or charge transfer[50, 51].Another interesting application of ultracold atoms is quantum information processing[52].In this context, time-dependent colliding trap potentials have been proposed for the real-ization of two-qubit quantum gates as well as the efficient creation of highly entangledstates[53, 54] which are two essential features required for a quantum computer.In the present investigation two bosonic particles are confined in two colliding Gaussianpotential wells. We solve this time-dependent problem using the ab-initio multi-configurationtime-dependent Hartree method for bosons (MCTDHB) which provides an exact descriptioncapturing all correlations[55, 56]. This allows us to compute the time evolution of the two-body wave function across a wide range of kinematic parameters in contrast to the othertheoretical investigations of colliding potentials in the literature[53, 54] which relied on em-ploying effective models and were limited to adiabatic movements of the traps. We show thatduring the time evolution of this system an effective, time-dependent double-well structureforms that drives an oscillatory over-barrier bosonic transport between the wells. This pro-cess terminates when the wells have been separated sufficiently after penetrating each other.During the collision process the displacement of the bosons from the well trajectories inducesan intra-well dipole mode and determines the final distribution of the particles between thewells. For fast collisions this setup exhibits deconfinement of the particles which we can at-3ribute to two different mechanisms. Firstly, for very fast accelerations an increase in kineticenergy leads to a positive total energy of the system towards the end of the time evolutionthereby causing an untrapping of particles. Secondly, we observe a resonant enhancement ofthe emission for certain kinematic parameters similar to the ionization processes that takeplace in atom-atom collisions.Our work is structured as follows. In Section II we introduce the physical setup anddescribe the computational approach used to solve the time-dependent problem. We proceedby presenting the results for the dynamics of two interacting bosonic particles in Section IIIand discuss suitable observables to unravel the properties of the system. We summarizeour findings in Section IV and provide an outlook on possible future studies. Finally, wecomment on the convergence of our variational multi-configuration time-dependent Hartreemethod for bosons (MCTDHB) approach in Appendix A.
II. PHYSICAL SETUP AND COMPUTATIONAL APPROACH
In the present work, we investigate the non-equilibrium quantum dynamics of a closedsystem of N = 2 interacting bosons. We employ MCTDHB[55, 57, 58] to solve the time-dependent many-body Schr¨odinger equation and gain access to the correlated quantumdynamics of the particles. This approach employs a time-dependent, variationally optimalbasis { φ i ( x, t ) } Mi =1 of M single particle functions (SPF). The many-body wave function | Ψ( t ) (cid:105) is then expanded as a superposition | Ψ( t ) (cid:105) = (cid:88) (cid:126)n | N C (cid:126)n ( t ) | (cid:126)n ; t (cid:105) (1)of all (cid:0) N + M − N (cid:1) time-dependent N -particle number states | (cid:126)n ; t (cid:105) that can be built from the M SPFs using time-dependent coefficients C (cid:126)n ( t ). Finally, the Lagrangian formulation of thetime-dependent variational principle[59, 60] yields equation of motions (EoM) for the SPFsand the coefficients[55, 57] are then solved numerically. MCTDHB provides access to thetime evolution of complete full many-body wave function which allows us to compute allrelevant characteristics of the underlying system.We consider N = 2 bosons of mass m interacting repulsively with a contact interaction4f strength of g [61, 62]. The Hamiltonian of the system reads H ( { x i } , t ) = N (cid:88) i =1 h ( x i , t ) + g N (cid:88) i,j =1 i In the scope of the present work we limit ourselves to N = 2 particles when investigatingthe setup described in Section II in order to unravel the main signatures of the dynamicsof the system. This provides an ideal starting point for future works addressing the case oflarger particle numbers. We choose wells of equal width, i.e. α = 1 and depth V = V (cid:48) =20 E G which are deep enough to support 10 trapped states of the one-body Hamiltonian (3).Initially, the wells are located at µ (0) = − . l G and µ (cid:48) (0) = 3 . l G which corresponds toan initial separation of d (0) = 7 l G . For the interaction strength we choose a value of g =0 . E G l G which is comparable to an interaction strength of g HO ≈ . 199 in harmonic oscillatorunits. We find that for this value of g , M = 6 SPFs are sufficient for the convergence of ourMCTDHB simulations (see Section A). We solve the time-dependent problem for varyingvalues of the acceleration a chosen such that the corresponding inverse final speeds v − are equally spaced in the interval (cid:2) . v − , . v − (cid:3) . The reason for this choice will becomeapparent during the analysis since many quantities scale with the inverse speed. A. Time Evolution of the One-Body Density In order to analyze the dynamics of the system and to guide our further analysis approach,we inspect the one-body density[64, 65] ρ (1) ( x, t ) = N (cid:90) | Ψ( x, x , . . . , x N , t ) | d x . . . d x N . (8)This quantity provides insight into the temporal evolution of the spatial distribution of theparticles since ρ (1) ( x, t ) corresponds to the probability density of finding a particle at theposition x at the time t .Figures 2 (a)–(f) show the time evolution of ρ (1) ( x, t ) for various values of the accelerationwhich correspond to different inverse final speeds v − . If the acceleration is not too fast (seeFigure 2 (a)–(e)), we can identify three distinct stages of the dynamics indicated by (I)–(III).7 10 20 30 − x / l G (a) v − v G = 2 . (b) v − v G ≈ . − x / l G (c) v − v G ≈ . (d) v − v G ≈ . t/t G − x / l G (e) v − v G ≈ . . . . t/t G (f) v − v G = 0 . FIG. 2. Time evolution of the one-body density ρ (1) ( x, t ) (see Equation 8) for different inverse finalspeeds v − ∝ a − / . The dashed white lines indicate the trajectories of the well centers while thedotted white lines indicate the positions of the FWHM of the Gaussian wells. The particles are initially localized in the well centered at µ (0) = − . l G and follow itsparabolic trajectory µ ( t ) during stage (I) of the dynamics while wells approach each other.No effect of the presence of the second well centered around µ (cid:48) ( t ) is visible during thisphase of the dynamics. During stage (II) the wells are in close proximity and they evenpenetrate each other. Hence, an effective double-well structure forms (see Figure 1) thatchanges its shape over time and we observe a collective oscillatory particle transport overthe central barrier. Towards the end of the propagation, during stage (III), we find severaleffects depending on the acceleration and hence v − . In general the particles are delocalizedover both wells with varying ratios. For certain values of v − however, the bosons arealmost completely localized in one of the wells. Additionally, we observe a sloshing motionof the particles within each well. We characterize this motion as a dipole mode[61, 62]since the center of mass (CoM) position of the particles oscillates around the center of thewells in which they are confined. This collective excitation is accompanied by a breathingmode which manifests in a periodic widening and contraction of the atomic cloud in each8ell. However, the breathing is much less pronounced compared to the dipole oscillationsuch that we refer to the sloshing motion as a dipole mode in the following. Generally, weobserve that the one-body density is well contained within one full width at half maximum(FWHM) around the well centers as indicated by the white lines in Figure 2. However, forfast collisions (see Figure 2 (e)) we notice a faint density halo in the region between thewells, which indicates an untrapped fraction of particles. When moving towards even fasteraccelerations we also observe effects of the inertia of the bosons (see Figure 2 (f)) which seemto move more slowly than the left well and leave the FWHM region before finally catchingup with the well towards the end of the dynamics. B. Center of Mass Position In order to analyze the transport of particles, we introduce the CoM position (cid:104) X (cid:105) ( t ) = 1 N N (cid:88) i =1 (cid:104) x i (cid:105) ( t ) (9)which measures the average position of the particles. In Figure 3 (a) and (b) we showtwo examples for the time evolution of this quantity. We can clearly make out the threeaforementioned phases (I)–(III) of the dynamics. During stage (I) of the time evolution, (cid:104) X (cid:105) ( t ) matches the trajectory of the left well µ ( t ) as the particles simply follow the motionof the potential. In part (II) we observe an oscillation of (cid:104) X (cid:105) ( t ) around 0 which indicatesthe oscillatory particle transport in the effective double-well structure from the left to theright well and vice versa. During stage (III) we notice that the evolution of (cid:104) X (cid:105) ( t ) stronglydepends on the kinematic parameters. For some values of v − v G , (cid:104) X (cid:105) ( t ) closely follows oneof the trajectories µ ( t ) and µ (cid:48) ( t ) and the dipole mode vanishes (see Figure 3 (b)). In othercases (see Figure 3 (a)) lies in the region between µ ( t ) and µ (cid:48) ( t ) and the dipole mode iswell-pronounced. The amplitude of the dipole mode varies depending on a and is maximalwhen (cid:104) X (cid:105) ( t ) oscillates close to zero.As a next step, we quantify the number of transport processes during phase (II) of thedynamics by determining the number of zero crossings N (II)ZC of the signal (cid:104) X (cid:105) ( t ) for eachvalue of v − during this stage (see Figure 3 (d)). N (II)ZC increases monotonically with v − sincethe effective double-well structure persists for a longer time period and more oscillations cantake place. Since the number of zero crossings has to be a non-negative integer, N (II)ZC is a9tep function of v − . We find the step width to be approximately equal for all steps with anaverage width of 0 . v − . t/t G − − (cid:104) X (cid:105) ( t ) / l G (a) (I) (II) (III) t/t G (b) (I) (II) (III) v − v G µ (cid:48) ( t f )0 µ ( t f ) (cid:104) X (cid:105) ( t f ) / l G (c) v − v G N ( II ) Z C (d) FIG. 3. (a–b) Time-evolution of the CoM position (solid blue line) as a function of time for v f ≈ . 355 (a) and v f ≈ . 247 (b). The orange dashed line indicates the trajectory µ ( t ) while thegreen dotted line visualizes µ (cid:48) ( t ). (c) Expectation value of the CoM position of the particles in thefinal state as a function of v − . The dashed orange line corresponds to a cosine fit of the signal.(d) Number of zero crossings N (II)ZC of (cid:104) X (cid:105) ( t ) in the region (II) as a function of v − . As mentioned before, the final location of the particles strongly depends on the acceler-ation a . Figure 3 (c) shows the final CoM position of the particles (cid:104) X (cid:105) ( t f ) as a functionof v − which resembles a cosine-like structure. Using a least squares fit we can extract theperiod ∆ v − = 0 . v − and the amplitude 3 . l G of the signal. From the amplitude ofthe oscillation, we can deduce that indeed for certain values of v − the density is almostcompletely located in one of the wells. A value of (cid:104) X (cid:105) ( t f ) = ± . l G would indicate that theaverage position of the particles coincides with the final position of one of the well centers.For most values of v − however, the final center of mass position lies somewhere betweenthese extreme cases and indicates that the particles are delocalized across both wells.A further analysis of the center of mass motion shows that the final distribution of theparticles as well as the amplitude of the dipole mode depend on the displacement of theCoM position from the trajectories of the wells at the transition from stage (II) to (III) ofthe dynamics. If the CoM position (cid:104) X (cid:105) ( t ) is close to one of the well centers at this transitionpoint, the particles get pinned in that particular well. A small deflection of (cid:104) X (cid:105) ( t f ) from thewell center leads then to small amplitudes of the corresponding dipole mode in this well. For10ost values of v − however, the separation of the wells splits the one-body density into twoparts and the particles are delocalized across both wells. As emphasized the displacementof the particles within the wells induces an intra-well dipole mode, the amplitude of whichis maximal if (cid:104) X (cid:105) ( t ) is close to 0 at the transition from stage (II) to (III) which correspondsto the maximal deflection of the particles from the well center.In order to distinguish between the intra-well dynamics different wells, we introduce thetruncated CoM observables (cid:104) X ± (cid:105) ( t ) = 1 N N (cid:88) i =1 (cid:104) x i Θ( ± x i ) (cid:105) ( t ) . (10)which measure the average position of particles on either the positive or the negative sidewith respect to x = 0. During part (III) of the dynamics, the dipole motion of the particlesin the initially left (right) well manifests itself in an oscillatory modulation of (cid:104) X + (cid:105) ( t )( (cid:104) X − (cid:105) ( t )). By analyzing the turning points of these modulations, we determine a phase of π / between the two oscillations. Furthermore, we notice that the oscillation period of bothobservables lies in the range 0 . . . . . t G and is approximately constant across all values of a which is to be expected since the frequency of the dipole mode only depends on the shapeof the potential well. C. Nature of the Particle Transport In order to classify the transport process that takes place, we analyze the two-body wavefunction | Ψ( t ) (cid:105) with respect to the time-dependent one-body Hamiltonian h ( x, t ) (Equa-tion (3)). We consider the instantaneous eigenbasis of h ( x, t ) spanned by the time-dependenteigenstates {| Φ i ( t ) (cid:105)} with the corresponding eigenenergies ε i ( t ), i.e. h ( x, t )Φ i ( x, t ) = ε i ( t )Φ i ( x, t ),while assuming an energetic ordering ε i ( t ) ≤ ε i +1 ( t ) for all times. Figure 4 shows theeigenenergies of the ten energetically lowest eigenstates as a function of the well separation d ( t ) = d (0) − at . At the initial ( d (0)) and final ( d ( t f )) separation, the external potentialis able to support ten trapped eigenstates, i.e. states with negative eigenenergies, which arepairwise degenerate. It should be noted that for positive energies the system exhibits adiscrete spectrum of untrapped states instead of a continuous spectrum of extended contin-uum states since we employ a finite grid for the numerical treatment of the problem whichimposes periodic boundary conditions (see Appendix A). However, this does not impact our11nalysis of the trapped fraction or the occupation of the trapped states. If the wells reachclose proximity, an effective double-well structure forms (see Figure 1), where V ( x = 0)determines the height of the barrier and the energetic degeneracies are lifted. In the vicinityof d ( t ) = 0 the central barrier vanishes and the external potential is a single Gaussian wellcentered around x = 0 with a depth V ( x = 0) = − V . Here, the eigenenergies ε ( t ), ε ( t )and ε ( t ) cross zero and reach positive value such that the associated eigenstates becomeuntrapped. d (0) 0 d ( t f ) d ( t ) − V − V ε i ( d ( t )) FIG. 4. Spectrum of the one-body Hamiltonian h ( x, t ) ((3)) as a function of the well separation d ( t ). We show the 10 energetically lowest eigenenergies (solid colored lines) and the values of thecentral potential V ( x = 0) (black dashed line). We proceed with our analysis by defining the operator P j ( t ) = 1 N N (cid:88) i =1 | Φ ij ( t ) (cid:105) (cid:104) Φ ij ( t ) | (11)where | Φ ij ( t ) (cid:105) (cid:104) Φ ij ( t ) | projects the i -th particle onto the j -th one-body eigenstate | Φ j ( t ) (cid:105) .Computing the expectation value of this projector with respect to the many-body wavefunction yields the probability p j ( t ) = (cid:104) Ψ( t ) | P j ( t ) | Ψ( t ) (cid:105) of finding a particle in the j -thone-body eigenstate.In order to unravel the nature of the particle transport, it is instructive to subdivide theset of one-body eigenstates into two categories. Firstly, we introduce the set B A ( t ) thatcontains all states that lie below the central barrier, i.e. all states | Φ i ( t ) (cid:105) with eigenenergies ε i ( t ) < V ( x = 0 , t ). Secondly, B B ( t ) captures all remaining trapped states, i.e. all states | Φ i ( t ) (cid:105) with eigenenergies V ( x = 0 , t ) ≤ ε i ( t ) < 0. It should be noted that both theeigenenergies as well as the central potential and consequently also the sets B σ ( t ) changeover time. 12s a next step we construct the operators O σ ( t ) = (cid:88) j such that | Φ j ( t ) (cid:105)∈B σ ( t ) P j ( t ) , σ ∈ { A , B } (12)that project the many-body wave function onto the states in the respective basis sets. Theexpectation values (cid:104) O σ ( t ) (cid:105) can be understood as the probabilities of a particle to occupyany of the states included in the corresponding basis set B σ ( t ). Additionally, we define theoperator O C ( t ) = 1 − O A ( t ) − O B ( t ) that projects the wave function onto the orthogonal spaceof all untrapped eigenstates. Consequently, the expectation value (cid:104) O C ( t ) (cid:105) correctly capturesthe occupation of untrapped continuum which is discretized due to our finite numerical grid.Figure 5 shows examples for the time evolution of these quantities. In the initial state,only under-barrier states are occupied and hence (cid:104) O A ( t ) (cid:105) ≈ (cid:104) O A ( t ) (cid:105) drops to zero while the occupation (cid:104) O B ( t ) (cid:105) of the trapped over-barrier states rises to approximately one. Consequently we classify theparticle transport that occurs during this stage of the time evolution as an over-barrierprocess. A deeper analysis shows that the start of transport coincides with the crossing of V ( x = 0 , t ) of the eigenenergies ε ( t ) and ε ( t ) (see Figure 4). The corresponding states | Φ ( t ) (cid:105) and | Φ ( t ) (cid:105) are predominantly occupied (see Figure 6). Consequently, the particletransport occurs when these states lie above the central barrier. Towards the end of thepropagation, the over-barrier states become under-barrier states again such that (cid:104) O A ( t ) (cid:105) → (cid:104) O B (cid:105) ( t ) → t → t f .For fast collisions (see Figure 5 (c) and (d)) untrapped states come into play as can beseen in an increase of (cid:104) O C ( t ) (cid:105) towards the end of the dynamics. We analyze this phenomenonfurther in Section III D where we investigate the emission of particles. D. Deconfinement of Particles As a next step in our analysis, we investigate the origin of the faint density halo betweenthe wells that we observe for fast collisions (see Figure 2 (e)), indicating a deonfinement ofparticles. The increase of (cid:104) O C ( t ) (cid:105) > h ( x, t ) (see Equation 3) come into play.In order to understand how the occupation of the individual eigenstates evolves over time,13 t/t G . . . (cid:104) O σ ( t ) (cid:105) (a) v G / v f = 2 . t/t G (b) v G / v f ≈ . 691 0 2 t/t G (c) v G / v f ≈ . 221 0 1 t/t G (d) v G / v f = 0 . FIG. 5. Time evolution of the projections (cid:104) O A ( t ) (cid:105) (solid blue line), (cid:104) O B ( t ) (cid:105) (solid orange line) and (cid:104) O C ( t ) (cid:105) (solid green line) for different final speeds v − . In (c) and (d) we also show the evolutionof (cid:104) O C ( t ) (cid:105) if the initially right well is absent during the propagation ( V (cid:48) = 0) in order to highlightthe influence of the second well on the deconfinement of the particles (see Section III D). we analyze the probabilities p j ( t ) = (cid:104) P j ( t ) (cid:105) of finding a particle in a specific one-bodyeigenstate. Figures 6 (a)–(d) show the time evolution of these quantities for specific valuesof v − . For slow collisions (see Figure 6 (a)) we observe that the eigenstates | Φ ( t ) (cid:105) and V (cid:48) = E G j (a) v G / v f = 2 . (b) v G / v f ≈ . (c) v G / v f ≈ . (d) v G / v f = 0 . − − − − − − l og ( p j ( t )) t/t G V (cid:48) = j (e) t/t G (f) t/t G (g) t/t G (h) − − − − − − l og ( p j ( t )) FIG. 6. Time evolution of the occupations log ( p j ( t )) of the 40 energetically lowest, instantaneouseigenfunctions of the one-body Hamiltonian (3). The top row shows the occupation under thepresence of the well centered around µ (cid:48) ( t ) while the bottom row shows the case V (cid:48) = 0. All statesbelow the red dashed line are trapped states while the states below the orange line are under-barrierstates. | Φ ( t ) (cid:105) are predominantly occupied while the other excited trapped states play a minor roleand no occupation of the untrapped states takes place. When increasing the acceleration andhence the collision speed, we observe a higher occupation of the excited trapped states anda minor population of several untrapped ones (see Figure 6 (b)). For the fastest collisions14nder consideration (see Figures 6 (c) and (d)) all 40 depicted eigenstates play a significantrole and we even observe an equal population of all eigenstates towards the end of thesimulation.We remark that the occupation of untrapped states occurs at different stages of thedynamics when comparing Figures 6 (b)–(d). In Figure 6 (b) the population of untrappedstates increases abruptly towards the end of the considered dynamics while still remainingsmall overall (cid:104) O B ( t ) (cid:105) (cid:28) (cid:104) O C ( t f ) (cid:105) ≈ . (cid:29) (cid:104) O A ( t f ) (cid:105) + (cid:104) O B ( t f ) (cid:105) . Here, we alsoobserve an additional steady increase in the population of untrapped states that alreadystarts in part (I) of the time evolution. Even though this is a small effect, it still suggests theexistence of two distinct mechanisms of the particle deconfinement. For very fast collisions(see Figure 6 (d)) the steady increase of the untrapped population becomes dominant. Thisenhancement for faster collisions suggests that it is a kinematic effect of the particle whichget spilled out of the potential wells due to the fast acceleration.In order to distinguish between the two effects leading to deconfinement and to unraveltheir origins, it is instructive to compare the results in Figures 6 (a)–(d) with simulationswhere the second, initially empty well is not present, i.e. for V (cid:48) = 0 (see Figure 6 (e)–(h)).The first striking difference is the absence of a sudden jump in the occupation of untrappedstates towards the end of the time evolution (see Figure 6 (b) and (c) vs. Figure 6 (f) and(g)). This contribution to the deconfinement can only be explained due to the presenceof the second well. However, the steady increase in the occupation of untrapped one-bodystates is still present (see Figure 6 (c) and (d) vs. Figure6 (g) and (h)). In Figure 5 theseobservations become even clearer when comparing the evolution of (cid:104) O C ( t ) (cid:105) with and withoutthe presence of the initially empty well (see Figure 5). For very fast collisions (see Figure 5(d)) the curves match for the biggest part of the dynamics and only deviate slightly towardsthe end of the time evolution. Consequently, the presence of the second well plays only aminor role concerning the emission of particles. For other parameters however (see Figure 5(c)), the differences are striking and the occupation of untrapped states is greatly enhanceddue to the presence of the second well.As mentioned before, the emission process during early times of the dynamics is of kine-matic origin. We employ the energy of the system as well as its composition to study this15henomenon further. Figure 7 (a) shows the total energy E ( t ) as a function of t for variousinverse final speeds v − . Since we prepare the system in the ground state all energy curvesstart at the ground state energy E ( t = 0) = E ≈ − . E G . When focussing on a very slowmotion of the wells (see curve for v − v G = 2 . t ≈ . t f where it starts to drop as the particles are now impacted by the second potential well. Asthe wells separate, the energy increases back to its initial value. The behavior of the totalenergy changes gradually as we turn towards faster accelerations. First, the dip of the energybecomes less deep and a modulation of the energy becomes visible towards the end of thesimulated dynamics. For v − v G ≈ . t ≈ . t f for v − v G ≈ . . . . t/t f − E ( t ) / E G (a) v − v G = 2 . v − v G ≈ . v − v G ≈ . v − v G = 0 . E ( t f ) / E G (b) V (cid:48) = 20 E G V (cid:48) = 0 050 E k i n ( t f ) / E G (c) v − v G − − − E p o t ( t f ) / E G ∆ v − (d) v − v G . . E i n t ( t f ) / E G ∆ v − (e) FIG. 7. (a) Time evolution of the total energy of the two bosons during the collision dynamics forvarious inverse final speeds v − . The other panels show (b) the total, (c) kinetic, (d) potential and(e) interaction energy of the final state as a function of v − . The orange dotted lines in (b) and(c) correspond to computations performed in the absence of the second, initially right, well, i.e. V (cid:48) = 0, thereby highlighting the impact of this well on the total and kinetic energies. we analyze the energy composition of the final state to get an overview of all simulations.Figure 7 (b)–(d) show the total, kinetic and potential energies of the final state as a function16f the final inverse speed v − . We notice a drastic increase of the kinetic (see Figure 7 (c))and hence the total energy (see Figure 7 (b)) towards large final speeds, i.e. small 1 /v f . For v − v G < . 266 with V (cid:48) = V as well as for v − v G < . 170 with V (cid:48) = 0 the total energyexceeds zero indicating that untrapping takes place solely kinetic energy reasons. The po-tential energy (see Figure 7 (d)) exhibits equidistant peaks whose height increases towardssmall values of v − as the particles become less deeply trapped. As indicated in the figure,the difference between neighboring peaks is equal to half of the period ∆ v − = 0 . v G thatwe introduced in our discussion of the final CoM position of the particles. The same char-acteristics and effects can be seen for the interaction energy (see Figure 7 (e)). The maximaof the interaction energy coincide with the extrema of (cid:104) X (cid:105) ( t f ) since the interaction energyis higher when both particles reside in the same well. The potential energy, on the otherhand side, becomes maximal if where (cid:104) X (cid:105) ( t f ) is zero.So far, our discussion of the particle untrapping relied on the projection onto one-bodyeigenstates. We conclude our analysis of this phenomenon using a two-body or in generalmany-body analysis that relies on projecting the many-body wave function onto number-states built from the instantaneous eigenbasis of the one-body Hamiltonian. Let N ( t ) bethe time-dependent set of all N = 2-particle number states that can be constructed fromall trapped eigenstates of the instantaneous one-body Hamiltonian. We then define themagnitude M B ( t ) = (cid:80) | (cid:126)n (cid:105)∈N ( t ) |(cid:104) (cid:126)n | Ψ( t ) (cid:105)| which captures the total overlap of the many-bodywave function with the number state basis N ( t ). The maximal possible value of M B ( t ) = 1indicates that the many-body wave function lies completely in the Hilbert space spannedby the basis N ( t ) while a value of zero would indicate that | Ψ( t ) (cid:105) is orthogonal to thisspace. Consequently the quantity M U ( t ) = 1 − M B ( t ) can then be used to quantify theuntrapped fraction, i.e. the projection of the many-body function onto the orthogonal spaceof untrapped eigenstates.Figures 8 (a)–(d) show the time evolution of M U ( t ) for different values of v − . For slowto moderately fast collisions (see Figure 8 (a) and (b)), no deconfinement of particles isvisible in the absence of the second well, i.e. for V (cid:48) = 0. As discussed previously, only the‘kinematic emission’ of particles takes place when only a single well is present. This processis enhanced by the collisional speed and we only observe untrapping for the fastest collisionsunder consideration (see Figure 8 (c) and (d)). When comparing these results with thesimulations with V (cid:48) = V , the importance of the presence of both wells becomes evident.17or certain values of v − a drastic increase in the untrapped fraction is noticeable that stemsfrom the final stage of the dynamics (see Figure 8 (a) and (c)). At very high speeds however,the kinematic untrapping is the dominant contribution to the emission of particles such thatthe two curves for M U ( t ) (single and two well dynamics) match each other.The logarithmic representation of the one-body density in Figures 8 (e)–(h) increases thevisibility of the density halo outside of the wells in contrast to the earlier discussion (seeFigure 2). For very fast collisions (see Figure 8 (h)), we notice a density halo on the left sideof the initially occupied well due as a fraction of the density gets spilled out of the potentialwells due to the inertia of the particles. Furthermore, we observe that in the case of theresonant emission of particles at certain values of v − , the density halo is located in thespace between the two well trajectories (see Figures 8 (e) and (g)). At other values, wherealmost no deconfinement takes place, this halo is vanishingly small (see Figure 8 (f)).Figure 8 (i) shows the value of M U ( t ) for the final state. In the absence of the secondwell, i.e. for V (cid:48) = 0, the curve of M U ( t f ) is flat and close to a value of zero for v − v G (cid:39) . V (cid:48) = V ), M U ( t f )exhibits peaks in the parameter regime v − v G (cid:39) . 39 that are not present for V (cid:48) = 0.Figure 8 (j) shows the difference ∆ M U ( t f ) between the simulations with V (cid:48) = V and V (cid:48) =0. This removes all contributions to the untrapping process that exclusively stem fromthe acceleration and not from the influence of the second well. We are able to identifythree distinct peaks at 0 . v − , 0 . v − and 0 . v − where the emission of particles isresonantly enhanced. ∆ M U ( t f ) as a function of v − is reminiscent of an ionization spectrum. E. Inter-Particle Correlations and Entanglement We now analyze the emergence of correlations and entanglement during the collisiondynamics by employing the Von Neumann entropy[66] which reads S (1) ( t ) = − Tr (cid:2) ˆ ρ (1) ( t ) ln (cid:0) ˆ ρ (1) ( t ) (cid:1)(cid:3) = − M (cid:88) i =1 λ i ( t ) ln( λ i ( t )) . (13)Here ˆ ρ (1) ( t ) refers to the one-body density matrix[64] with eigenvalues λ i ( t ). It should benoted that the natural populations λ i ( t ) possess the property 0 ≤ λ i ( t ) ≤ . . . M U ( t ) (a) v G / v f ≈ . (b) v G / v f ≈ . (c) v G / v f ≈ . (d) v G / v f = 0 . 10 5 t/t G − x / l G (e) t/t G (f) . . t/t G (g) t/t G (h) − − . . . . v − v G . . . M U ( t f ) (i) V (cid:48) = V V (cid:48) = 0 0 . . . . v − v G ∆ M U ( t f ) (j) FIG. 8. (a)–(d) Time evolution of the untrapped fraction M U ( t ) for varying v − (blue lines). Theorange lines indicate the evolution of M U ( t ) in the absence of the second, initially empty well (i.e. V (cid:48) = 0) highlighting its importance for the untrapping process for certain values of v − . (e)–(h) Time evolution of the one-body density log ( ρ (1) ( x, t )) (see Equation (8)) for V = V (cid:48) in alogarithmic representation which increases the visibility of the density halo outside the potentialwells in comparison to Figure 2. (i) Untrapped magnitude M U ( t f ) of the final state as a functionof v − . The dotted vertical lines indicate the values of v − that have been used for panels (a)–(d)and (e)–(h). (j) Untrapped magnitude ∆ M U ( t f ) due to the presence of the second well (see maintext). relation (cid:80) Mi =1 λ i ( t ) = 1.A value of S (1) ( t ) = 0 indicates a mean-field state and implies the absence of any correla-tions between the two particles. In the same light, a finite value of S (1) ( t ) (cid:54) = 0 correspondsto inter-particle correlations and hence a deviation from the mean-field product state. For amaximally entanled state within our simulations using six SPFs, the Von Neumann entropyreaches the maximal value of S (1)max = ln( M ) = ln(6) ≈ . 79 (14)which is here solely determined by the dimensionality of the one-body Hilbert space M = 6.Figure 9 shows the entropy of the final state as a function of the final inverse speednormalized to the maximal possible value. We observe a structure of equidistant peaks ofvarying height indicating large values of S (1) ( t f ). The spacing is approximately equal to the19eriod ∆ v − = 0 . v − obtained during the CoM analysis suggesting a relation to the finallocation of the particles. This hypothesis can be easily confirmed by analyzing the one-body density and the CoM observable which show that the maxima of the Von Neumannentropy correspond to situations where the particles are distributed uniformly over bothwells in the final state. Furthermore, we notice that the entropy reaches its largest valueof S (1) ( t f ) ≈ . S max for v − v G ≈ . 21 indicating a highly entangled state for which thetwo largest natural populations are almost equal ( λ ( t f ) ≈ . 517 and λ ( t f ) ≈ . v − where the particles are localizedin one of the wells, i.e. extrema of the CoM position. Here, the first natural population isdominant λ ( t f ) ≈ 1. We notice that the height of the local maxima decreases towards faster . . . . . . v − v G . . . . S ( ) ( t f ) / S m a x ∆ v − FIG. 9. Von Neumann entropy of the final state S (1) ( t f ) normalized by the maximal possible value S (1)max as a function of the inverse final speed v − . collisions and the entropy drops to zero indicating a mean-field product state. The reasonfor this behavior is that for v − → λ ≈ v − v G (cid:39) S (1) ( t f ) vanishes butthe entropy does not drop to zero. This indicates that still measurable correlations betweenthe two particles exist. IV. CONCLUSIONS AND OUTLOOK We have investigated the collisional non-equilibrium quantum dynamics of ultracoldbosons confined in two colliding potential wells. We were able to subdivide the dynam-ics into three distinct stages by identifying the underlying physical processes. Initially, theparticles follow the trajectories of the wells closely. When the well separation falls belowa certain threshold, a periodic collective particle transport takes place in an effective time-dependent double-well structure. By analyzing the population of SPF states we were ableto classify this transport as an over-barrier process. Using the CoM position of particles20e have been able to quantify the number of oscillatory transitions that occur during thedynamics. During the separation of the wells in the third part of the time evolution, wenotice a mode motion of the particles within each well. The amplitude of this motion de-pends on the location of the particles with respect to the well centers at the end of thecollision process. We determine a phase of π / between the dipole modes of both wells whilethe frequency of this motion is independent of the acceleration. Furthermore, we observethat for certain final speeds the particles are strongly localized in one of the wells whilethey are generally delocalized. This phenomenon resembles the charge transfer that takesplace during atom-atom collisions. Another important feature of our time-dependent setupis the untrapping of particles which we characterize in detail using a SPF, number state andenergetic analysis. We have been able to quantify the untrapped fraction unraveling twodifferent contributions to it. During fast collisions, the kinetic energy grows continuouslywhich leads to a positive total energy and consequently to a particle untrapping. However,we also observe a resonant untrapping effect for certain kinematic parameters leading toa rapid emission of particles as the wells separate. We have been able to determine thedependence of this second mechanism on the kinematic parameters which is reminiscent ofan ionization spectrum.Our findings serve as a promising starting point for further studies in different directions.By increasing the inter-particle interaction strength one could enhance the amount of cor-relations that arises during the dynamics and it is interesting to explore the correspondingimpact on the resonant particle untrapping. A variation of the potential wells for exampleby decreasing the depth or introducing an asymmetry between the two Gaussians couldmodify the particle transport. In the light of atom-atom collisions, a particularly intriguingprospect is to employ different initial states. Employing an initial state that incorporatesparticles in both wells, could lead to an enhancement of the emission due to opposite mo-menta of the bosons. Furthermore it would be interesting to investigate the impact of thetrajectories of the wells. Lastly, the multi-configuration time-dependent Hartree method forfermions[67, 68] allows to study the non-equilibrium dynamics of fermions in a similar setup.It would be instructive to analyze the role of the particle statistics and how the phenomenadescribed in this work might be modified.Another exciting route would be the investigation of mixtures of different compo-nents which is of particular interest for ultracold atom research. Such ensembles can21e composed of different elements[69, 70], isotopes[71] or hyperfine states[72] and exhibita plethora of exciting and unique properties like relative phase evolution[73], compositefermionization[74], non-linear[75] and collective excitations[76] as well as miscible-immisciblephase transitions[77, 78]. Depending on the particle statistics this allows for the realizationof Bose-Bose[79, 80], Fermi-Fermi[81, 82] and Bose-Fermi mixtures[83–86]. The multi-layermulti-configuration time-dependent Hartree method for mixtures[56] is a powerful numericalapproach to treat the correlated non-equilibrium dynamics of such systems which allows toextend the setup presented in the present work to such mixtures. The role of the inter-speciesinteraction as well as a possible mass-imbalance between the constituents are particularlyof interest. ACKNOWLEDGMENTS The authors acknowledge fruitful discussions with K. Keiler. This work has been fundedby the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) SFB 925project 170620586. Appendix A: Technical Aspects and Convergence In the present work, we employ the fast Fourier transform (FFT)[87–89] to obtain aspatially discretized representation of the operators and the SPFs. This scheme allows theefficient numerical treatment of large grids consisting of n (cid:39) 100 grid points compared toanother approaches relying on discrete variable representations (DVR)[89]. We use n = 675grid points that are equally spaced in the interval ( − l G , l G ]. It should be noted that theFFT scheme implies periodic boundary conditions for the physical system. We repeat thesame set of simulations presented in the main text using a sine DVR[89] which incorporateshard-wall boundary conditions. Thereby we are able to confirm that spacing between thepotential wells and the edges of the grid is large enough such that no influence of theboundary conditions is visible in the observables discussed in the present work.The underlying time-dependent variational principle used to derive the MCTDHB equa-tions of motion guarantees that the SPF basis is rotated such that the many-body wavefunction optimally captures the state of the physical system. However, care has to be taken22n order to ensure that the number M of SPFs is sufficiently large and thereby the numericalconvergence of the method is guaranteed[55, 89]. We compare the results presented in themain text with simulations that include an additional, seventh SPF and observe that theobservables discussed in the main text do not change significantly. The ground state energyexhibits a relative change of the order of 10 − and the energy of the final state of 10 − in theworst case. We observe that the untrapped fraction of the final state ∆ M U ( t f ) determinedchanges at most by an absolute value of 4 · − when including the additional orbital. Theabsolute change in the relative entropy S (1) / S (1)max of the final state is limited by 0 . 03. Thecenter of mass position of the particles at the end of the time evolution changes at most byone percent.Additionally, the spectral representation of the one-body density matrix is importantto judge the convergence of the approach. The eigenvalues of ρ (1) ( t ), the so-called naturalpopulations, should exhibit a rapidly decreasing hierarchy. This indicates that any naturalorbitals (eigenstates of the one-body density matrix) that are neglected due to the truncationof the single particle Hilbert space play a negligible role. We find that this is the case forall parameters considered in the present work and that the least occupied orbital taken intoaccount shows a population of λ < − for all simulations. Therefore, we consider M = 6SPFs sufficient to describe the time evolution of the physical system accurately. [1] K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, andW. Ketterle, Phys. Rev. Lett. , 3969 (1995).[2] C. C. Bradley, C. A. Sackett, J. J. Tollett, and R. G. Hulet, Phys. Rev. Lett. , 1687 (1995).[3] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science , 198 (1995).[4] I. Bloch, Nature , 1016 (2008).[5] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. , 863 (2011).[6] I. Bloch, J. Dalibard, and S. Nascimb`ene, Nat. Phys. , 267 (2012).[7] F. Serwane, G. Z¨urn, T. Lompe, T. B. Ottenstein, A. N. Wenz, and S. Jochim, Science ,336 (2011).[8] A. M. Kaufman, B. J. Lester, C. M. Reynolds, M. L. Wall, M. Foss-Feig, K. R. A. Hazzard, . M. Rey, and C. A. Regal, Science , 306 (2014).[9] K. Henderson, C. Ryu, C. MacCormick, and M. G. Boshier, New J. Phys. , 043030 (2009).[10] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. , 3108(1998).[11] I. Bloch, Nat. Phys. , 23 (2005).[12] S. Chu, J. E. Bjorkholm, A. Ashkin, J. P. Gordon, and L. W. Hollberg, Opt. Lett. , 73(1986).[13] O. Morizot, Y. Colombe, V. Lorent, H. Perrin, and B. M. Garraway, Phys. Rev. A , 023617(2006).[14] M. Greiner, O. Mandel, T. Esslinger, T. W. H¨ansch, and I. Bloch, Nature , 39 (2002).[15] L.-M. Duan, E. Demler, and M. D. Lukin, Phys. Rev. Lett. , 090402 (2003).[16] O. Zobay and B. M. Garraway, Phys. Rev. Lett. , 1195 (2001).[17] Y. Colombe, E. Knyazchyan, O. Morizot, B. Mercier, V. Lorent, and H. Perrin, Europhys.Lett. , 593 (2004).[18] C. Orzel, A. K. Tuchman, M. L. Fenselau, M. Yasuda, and M. A. Kasevich, Science ,2386 (2001).[19] B. Paredes, A. Widera, V. Murg, O. Mandel, S. F¨olling, I. Cirac, G. V. Shlyapnikov, T. W.H¨ansch, and I. Bloch, Nature , 277 (2004).[20] T. K¨ohler, K. G´oral, and P. S. Julienne, Rev. Mod. Phys. , 1311 (2006).[21] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. , 1225 (2010).[22] M. Olshanii, Phys. Rev. Lett. , 938 (1998).[23] J. I. Kim, V. S. Melezhik, and P. Schmelcher, Phys. Rev. Lett. , 193203 (2006).[24] P. Giannakeas, F. K. Diakonos, and P. Schmelcher, Phys. Rev. A , 042703 (2012).[25] P. Giannakeas, V. S. Melezhik, and P. Schmelcher, Phys. Rev. Lett. , 183201 (2013).[26] B. P. Anderson and M. A. Kasevich, Science , 1686 (1998).[27] G.-B. Jo, Y.-R. Lee, J.-H. Choi, C. A. Christensen, T. H. Kim, J. H. Thywissen, D. E.Pritchard, and W. Ketterle, Science , 1521 (2009).[28] D.-S. L¨uhmann, C. Weitenberg, and K. Sengstock, Phys. Rev. X , 031016 (2015).[29] S. Sala, J. F¨orster, and A. Saenz, Phys. Rev. A , 011403 (2017).[30] G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif, and T. Esslinger,Nature , 237 (2014). 31] N. Goldman, J. C. Budich, and P. Zoller, Nat. Phys. , 639 (2016).[32] J. Steinhauer, Nat. Phys. , 959 (2016).[33] A. Ashkin, Phys. Rev. Lett. , 156 (1970).[34] A. Ashkin, J. M. Dziedzic, J. E. Bjorkholm, and S. Chu, Opt. Lett. , 288 (1986).[35] A. Ga¨etan, Y. Miroshnychenko, T. Wilk, A. Chotia, M. Viteau, D. Comparat, P. Pillet,A. Browaeys, and P. Grangier, Nat. Phys. , 115 (2009).[36] S. Saskin, J. T. Wilson, B. Grinkemeyer, and J. D. Thompson, Phys. Rev. Lett. , 143002(2019).[37] A. Ashkin and J. M. Dziedzic, Science , 1517 (1987).[38] K. O. Roberts, T. McKellar, J. Fekete, A. Rakonjac, A. B. Deb, and N. Kjærgaard, Opt.Lett. , 2012 (2014).[39] K.-N. Schymik, V. Lienhard, D. Barredo, P. Scholl, H. Williams, A. Browaeys, and T. Lahaye,Phys. Rev. A , 063107 (2020).[40] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lien-hard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. L¨auchli, and A. Browaeys, “Programmablequantum simulation of 2D antiferromagnets with hundreds of Rydberg atoms,” (2020),arXiv:2012.12268.[41] M. Morgado and S. Whitlock, “Quantum simulation and computing with Rydberg-interactingqubits,” (2020), arXiv:2011.03031.[42] A. Rakonjac, A. B. Deb, S. Hoinka, D. Hudson, B. J. Sawyer, and N. Kjærgaard, Opt. Lett. , 1085 (2012).[43] N. Kjærgaard, A. S. Mellish, and A. C. Wilson, New J. Phys. , 146 (2004).[44] R. Thomas, M. Chilcott, C. Chisholm, A. B. Deb, M. Horvath, B. J. Sawyer, and NielsKjærgaard, J. Phys.: Conf. Ser. , 012007 (2017).[45] R. Thomas, M. Chilcott, E. Tiesinga, A. B. Deb, and N. Kjærgaard, Nat. Commun. , 4895(2018).[46] S. V. Rajagopal, K. M. Fujiwara, R. Senaratne, K. Singh, Z. A. Geiger, and D. M. Weld,Ann. Phys. , 1700008 (2017).[47] R. Senaratne, S. V. Rajagopal, T. Shimasaki, P. E. Dotti, K. M. Fujiwara, K. Singh, Z. A.Geiger, and D. M. Weld, Nat. Commun. , 2065 (2018).[48] Y. D. Wang, J. H. McGuire, T. J. M. Zouros, D. H. Lee, J. M. Sanders, and P. Richard, ucl. Instrum. Meth. B , 124 (1993).[49] E. C. Montenegro, W. E. Meyerhof, and J. H. McGuire, Adv. Atom. Mol. Opt. Phy. , 249(1994).[50] W. L. Fite, R. F. Stebbings, D. G. Hummer, and R. T. Brackmann, Phys. Rev. , 663(1960).[51] R. E. Olson and A. Salop, Phys. Rev. A , 531 (1977).[52] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. , 885 (2008).[53] D. Jaksch, H.-J. Briegel, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. , 1975(1999).[54] T. Calarco, E. A. Hinds, D. Jaksch, J. Schmiedmayer, J. I. Cirac, and P. Zoller, Phys. Rev.A , 022304 (2000).[55] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Phys. Rev. A , 033613 (2008).[56] L. Cao, V. Bolsinger, S. I. Mistakidis, G. M. Koutentakis, S. Kr¨onke, J. M. Schurer, andP. Schmelcher, J. Chem. Phys. , 044106 (2017).[57] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, J. Chem. Phys. , 154103 (2007).[58] A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Phys. Rev. Lett. , 030402 (2007).[59] P. Kramer and M. Saraceno, Geometry of the Time-Dependent Variational Principle in Quan-tum Mechanics , Lecture Notes in Physics, Vol. 140 (Springer-Verlag Berlin Heidelberg, 1981).[60] J. Broeckhove, L. Lathouwers, E. Kesteloot, and P. Van Leuven, Chem. Phys. Lett. , 547(1988).[61] L. Pitaevskii and S. Stringari, Bose-Einstein Condensation , 1st ed., International Series ofMonographs on Physics No. 116 (Oxford University Press, Oxford, New York, 2003).[62] C. J. Pethick and H. Smith, Bose–Einstein Condensation in Dilute Gases , 2nd ed. (CambridgeUniversity Press, 2008).[63] H.-D. Meyer and G. A. Worth, Theor. Chem. Acc. , 251 (2003).[64] K. Sakmann, A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Phys. Rev. A , 023615(2008).[65] P. A. M. Dirac, Math. Proc. Cambridge , 376 (1930).[66] J. von Neumann, Nachr. Akad. Wiss. Gorringen Math.-Phys. Kl. , 273 (1927).[67] J. Zanghellini, M. Kitzler, C. Fabian, T. Brabec, and A. Scrinzi, Laser Phys. , 1064 (2003).[68] J. Caillat, J. Zanghellini, M. Kitzler, O. Koch, W. Kreuzer, and A. Scrinzi, Phys. Rev. A , , 053202 (2002).[70] D. A. Brue and J. M. Hutson, Phys. Rev. Lett. , 043201 (2012).[71] N. Poli, R. E. Drullinger, G. Ferrari, J. L´eonard, F. Sorrentino, and G. M. Tino, Phys. Rev.A , 061403 (2005).[72] C. J. Myatt, E. A. Burt, R. W. Ghrist, E. A. Cornell, and C. E. Wieman, Phys. Rev. Lett. , 586 (1997).[73] R. P. Anderson, C. Ticknor, A. I. Sidorov, and B. V. Hall, Phys. Rev. A , 023603 (2009).[74] M. A. Garcia-March, B. Juli´a-D´ıaz, G. E. Astrakharchik, T. Busch, J. Boronat, and A. Polls,Phys. Rev. A , 063604 (2013).[75] A. M. Kamchatnov, Y. V. Kartashov, P.-´E. Larr´e, and N. Pavloff, Phys. Rev. A , 033618(2014).[76] K. M. Mertes, J. W. Merrill, R. Carretero-Gonz´alez, D. J. Frantzeskakis, P. G. Kevrekidis,and D. S. Hall, Phys. Rev. Lett. , 190402 (2007).[77] C. Ticknor, Phys. Rev. A , 013623 (2013).[78] E. Nicklas, W. Muessel, H. Strobel, P. G. Kevrekidis, and M. K. Oberthaler, Phys. Rev. A , 053614 (2015).[79] G. Modugno, M. Modugno, F. Riboli, G. Roati, and M. Inguscio, Phys. Rev. Lett. , 190404(2002).[80] J. Catani, L. De Sarlo, G. Barontini, F. Minardi, and M. Inguscio, Phys. Rev. A , 011603(2008).[81] E. Wille, F. M. Spiegelhalder, G. Kerner, D. Naik, A. Trenkwalder, G. Hendl, F. Schreck,R. Grimm, T. G. Tiecke, J. T. M. Walraven, S. J. J. M. F. Kokkelmans, E. Tiesinga, andP. S. Julienne, Phys. Rev. Lett. , 053201 (2008).[82] C. Kohstall, M. Zaccanti, M. Jag, A. Trenkwalder, P. Massignan, G. M. Bruun, F. Schreck,and R. Grimm, Nature , 615 (2012).[83] F. Schreck, L. Khaykovich, K. L. Corwin, G. Ferrari, T. Bourdel, J. Cubizolles, and C. Sa-lomon, Phys. Rev. Lett. , 080403 (2001).[84] Z. Hadzibabic, C. A. Stan, K. Dieckmann, S. Gupta, M. W. Zwierlein, A. G¨orlitz, andW. Ketterle, Phys. Rev. Lett. , 160401 (2002). 85] C. Ospelkaus, S. Ospelkaus, K. Sengstock, and K. Bongs, Phys. Rev. Lett. , 020401 (2006).[86] J. Heinze, S. G¨otze, J. S. Krauser, B. Hundt, N. Fl¨aschner, D.-S. L¨uhmann, C. Becker, andK. Sengstock, Phys. Rev. Lett. , 135303 (2011).[87] D. Kosloff and R. Kosloff, J. Comput. Phys. , 35 (1983).[88] R. Kosloff, J. Phys. Chem. , 2087 (1988).[89] M. Beck, A. J¨ackle, G. Worth, and H.-D. Meyer, Phys. Rep. , 1 (2000)., 1 (2000).