Collision Dynamics of One-Dimensional Bose–Einstein Condensates

Abstract

We study the collision dynamics of two Bose–Einstein condensates, with their dynamical wave functions modeled by a set of coupled, time-dependent Gross–Pitaevskii equations. In an effective one-dimensional system, we identify regimes characterized by the relationship between inter- and intra-atomic interactions and the initial configuration of the system, akin to the equilibrium phase diagram of two interacting Bose condensates. We consider a dynamical setup in which two wave packets are initially at rest, with a small separation about the center of an anisotropic harmonic trap. Upon release, we observe a rapid approach to dynamical equilibrium in the limits of very large and very small inter-particle repulsion, characterized by periodic transmission or reflection of the condensates as distinguishable units, whereas the intermediate, critical regime is characterized by extended transient dynamics, density fracturing, and dynamical mixing.

1. Introduction

The study of Bose–Einstein condensates (BECs) has led to remarkable insights into the behavior of quantum systems at ultra-low temperatures, uncovering phenomena with applications ranging from fundamental physics to quantum technologies. Recent studies have started to investigate the dynamics of interacting BECs. For example, Guan et al. have explored the role of interactions in double-well systems created using dilute rubidium BECs, emphasizing the distinct dynamics arising from Raman and lattice coupling [1]. Their study uses a Ramsey-type interference pulse sequence, crucial for atom interferometers. On the other hand, Burchianti et al. have focused on the effect of interactions in the phase evolution of expanding atomic BECs, challenging the conventional assumption about the modification of phase evolution through an effective force [2]. Both papers contribute to our understanding of the intricate interplay between interactions and the dynamics of BECs, offering valuable insights for applications in areas such as quantum simulations and matter–wave interference. Our work expands on these previous studies by focusing on the creation of dynamic mixing regimes in colliding BECs.

An interesting feature of binary BECs is that they exhibit a range of ground-state configurations depending on the atomic interactions [3,4,5]. Binary BECs exhibit miscibility or immiscibility depending on the relative strength of these interactions. As observed in equilibrium experiments [6,7,8,9,10,11,12], when the inter-species repulsive interaction is large enough, the two components repel each other so that they separate into two distinct clouds with small spatial overlap, corresponding to the immiscible phase. When the effect of inter-species interactions is weak compared to the intra-species interactions, the two components are miscible and overlap with each other at the center of the trapping potential.

It is well known for a binary BEC in equilibrium, where the kinetic energy is negligible compared with the interaction energy, that the transition between the miscible and immiscible phases can be determined by the inter-species and intra-species interactions. The criterion for the equilibrium immiscible ground state is $g_{12}^2>g_{11}{g}_{22}$ [13,14], where $g_{11}$ and $g_{22}$ are the effective intra-species interactions of components 1 and 2, whereas $g_{12}$ is the inter-species interaction strength. One can then control the miscibility of the gases by using the Feshbach resonance technique to adjust the interaction strengths [11,12].

Research on binary BECs has been primarily concerned with static ground-state properties, most likely due to the difficulty in achieving binary mixtures experimentally for sufficiently long durations. There have been a few studies investigating dynamical behavior over relatively short time scales (≈10 ms) and small kinetic energies [7,15,16]. Our studies suggest that interesting behavior may exist in the collisions of binary BECs over long-time evolutions (>1 s) and for relatively large kinetic energies in a system composed of two fully condensed clouds of ultra-cold Rubidium atoms confined to an anisotropic harmonic trap. The primary purpose of this paper is to consider the boundary between miscible and immiscible phases in a dynamical setting. Specifically, we observe that there are certain non-equilibrium phases which can be mapped onto the ground-state phases (miscible and immiscible), as well as an additional phase related to dynamical mixing.

We begin by introducing a theoretical model for condensate evolution in an anisotropic harmonic trap and outline the static ground-state phases for comparison with the dynamical scenario that we study in this paper. Then, we describe our computational modeling of condensate evolution using a split-step Fourier method (SSM) over a finite spatio-temporal grid. Our results are then presented in a way that illustrates the dynamical phases we observe by varying the trap frequency and inter-atomic interaction strength. Our methods are described in more detail in the Appendix A.

2. Theoretical Model

We consider two condensates composed of 1000 ⁸⁷Rb atoms in the hyperfine spin states $|1\rangle =|1,-1\rangle$ and $|2\rangle =|2,+1\rangle$. Let $a_{11}$ be the s-wave scattering length of a BEC composed of atoms in the $|1\rangle$ state, and $a_{22}$ be the s-wave scattering length for $|2\rangle$. There is also an inter-atomic interaction between the condensates parameterized by $a_{12}$. Experimentally derived values of the scattering lengths are given in

Table 1. Table of parameters used in numerical simulations of colliding condensates composed of rubidium atoms.

Spatio-Temporal Grid		Scattering Lengths		Trapping Frequencies		Cloud Composition
Box size	320 μm	$a_{11}$	$100.40a_0$	${\omega }_x$	$2\pi \times 100$	Atomic mass	87 u
Grid points	${10}^4$	$a_{22}$	$95.44a_0$	${\omega }_y$	$2\pi \times 100$	Atom number	${10}^3$
Evolution time	$40\times {\displaystyle \frac{2\pi }{{\omega }_z}}$	$a_{12}$	$98.006a_0$	${\omega }_z$	$2\pi \times (5-50)$
Trotter steps	${10}^6$

from [17]. The condensates reside in an anisotropic harmonic trap,

$$V\left(\mathbf{r}\right)={\displaystyle \frac{1}{2}}m\left({\omega }_x^2{x}^2+{\omega }_y^2{y}^2+{\omega }_z^2{z}^2\right),$$

(1)

where m is the mass of rubidium. The trap is designed with the $x-y$ confinement being strong enough to create a quasi one-dimensional condensate in the z direction, but not strong enough to violate the Born approximation for the pseudo-potential of the atomic interactions. That is, if the radial trap frequency is given by ${\omega }_r={\omega }_x={\omega }_y,$ then the trap energy ${ϵ}_r=\hslash {\omega }_r$ is large enough that the condensate remains in the ground state of the trap in the $x-y$ plane but is free to explore new states in the z direction. Additionally, we assume that if ${\ell }_r=\sqrt{\hslash /m{\omega }_r}$ is the characteristic harmonic oscillator length, and $a_s$ is the s-wave scattering length of the condensate, then the Born approximation requires $a_s/{\ell }_r\ll 1,$ which is monitored and strictly enforced in our simulations.

The use of the 1D Gross–Pitaevskii equation (GPE) to describe the effective dynamics of Bose–Einstein condensates in a strongly anisotropic trapping potential has been justified by several foundational studies. In such systems, where the confinement is much stronger in two dimensions than in the third, the condensate wave function becomes effectively frozen in the tightly confined directions, reducing the dynamics to the remaining weakly confined dimension. This dimensional reduction is supported by the work of Pérez-García et al. [18], who demonstrated that the wave function can be factorized and the dynamics captured in one dimension. Olshanii further established how the effective 1D interaction strength emerges from the 3D scattering length under strong transverse confinement [19], a result that underpins the validity of the 1D GPE. Petrov et al. [20] and Menotti and Stringari [21] provided additional analyses of the quasi-1D regime, emphasizing the conditions under which the 1D approximation holds, particularly when the transverse energy gap is much larger than the interaction energy. These theoretical foundations are comprehensively discussed in Pitaevskii and Stringari’s textbook [13], which outlines the scenarios where the 1D Gross–Pitaevskii equation accurately describes the condensate’s behavior, making it a reliable tool for studying BECs in elongated traps.

The effective one-dimensional dynamics of two rubidium BECs occupying different hyperfine states can thus be described within a mean-field model by the following set of time-dependent, coupled Gross–Pitaevskii equations [13]:

$$\begin{array}{cc}\hfill i\hslash {\displaystyle \frac{\partial {\psi }_1}{\partial t}}& =\left(-{\displaystyle \frac{{\hslash }^2{\nabla }^2}{2m}}+V+N{g}_{11}{\left|{\psi }_1\right|}^2+N{g}_{12}{\left|{\psi }_2\right|}^2\right){\psi }_1,\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill i\hslash {\displaystyle \frac{\partial {\psi }_2}{\partial t}}& =\left(-{\displaystyle \frac{{\hslash }^2{\nabla }^2}{2m}}+V+N{g}_{22}{\left|{\psi }_2\right|}^2+N{g}_{12}{\left|{\psi }_1\right|}^2\right){\psi }_2,\hfill \end{array}$$

(2b)

where N is the atom number, and ${\psi }_1$ and ${\psi }_2$ are the 1D wave functions describing the condensates of atoms in $|1\rangle$ and $|2\rangle$

$$g_{ij}={\displaystyle \frac{2{\hslash }^2{a}_{ij}}{m{\ell }_r^2}},$$

(3)

where i and j refer to intra-atomic $i=j$ and inter-atomic $i\ne j$ interaction scattering lengths.

When released at a fixed separation relative to the trap minimum, we expect the condensates to oscillate and collide multiple times over the course of their evolution. We seek to understand how the inter-atomic interaction energy $g_{12}$ and the trap frequency ${\omega }_z$ affect the ensuing dynamics. Before moving to the dynamical simulations, let us first discuss the static ground-state phases, which will serve as a point of reference for the dynamical phases discussed later on. Beginning from the energy functional [13]

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill E& =\hfill \\ \hfill & \int dz\{\sum _{j=1,2}\left[{\displaystyle \frac{N{\hslash }^2}{2m}}{\left|{\partial }_z{\psi }_j\right|}^2+{\displaystyle \frac{1}{2}}Nm{\omega }_z^2{z}^2{\left|{\psi }_j\right|}^2\right]+\hfill \\ \hfill & {\displaystyle \frac{N^2{g}_{11}}{2}}{\left|{\psi }_1\right|}^4+{\displaystyle \frac{N^2{g}_{22}}{2}}{\left|{\psi }_2\right|}^4+N^2{g}_{12}{\left|{\psi }_1\right|}^2{\left|{\psi }_2\right|}^2\},\hfill \end{array}\end{array}$$

(4)

we apply a variational method for finding the approximate ground-state properties. In this work, we consider the two condensates to be initially unmixed and take a form closely resembling the ground state of the trap, which is approximately Gaussian in nature. Assuming sufficiently low particle numbers, we therefore take the following Gaussian ansatz to perform our analysis:

$${\psi }_j={\pi }^{-1/4}{w}_j^{-1/2}{e}^{-{\left(z-z_j\right)}^2/\left(2w_j^2\right)},$$

(5)

where $j=1,2$, $w_j$ is a variational parameter related to the width of the Gaussian distribution, and $z_j$ is another variational parameter representing the central location of the Gaussian distribution. Inserting this ansatz into Equation (4), we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill E=& {\displaystyle \frac{1}{4}}Nm{\omega }_z^2(w_1^2+2z_1^2)+{\displaystyle \frac{1}{4}}Nm{\omega }_z^2(w_2^2+2z_2^2)\hfill \\ \hfill & +{\displaystyle \frac{1}{4}}\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{N^2{g}_{11}}{w_1}}+{\displaystyle \frac{1}{4}}\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{N^2{g}_{22}}{w_2}}\hfill \\ \hfill & +{\displaystyle \frac{{\hslash }^2}{4m}}\left({\displaystyle \frac{N}{w_1^2}}+{\displaystyle \frac{N}{w_2^2}}\right)+{\displaystyle \frac{N^2{g}_{12}}{\sqrt{\pi }}}{\displaystyle \frac{e^{-{\delta }^2/(w_1^2+w_2^2)}}{\sqrt{w_1^2+w_2^2}}},\hfill \end{array}\end{array}$$

(6)

where we have defined $\delta =z_1-z_2.$

For the variational step, we now minimize Equation (6) with respect to $z_i$ and $w_i$ to obtain the following conditions for the ground-state configuration:

$$\begin{array}{c}\begin{array}{cc}\hfill (1-& {\displaystyle \frac{{\hslash }^2}{m^2{\omega }_z^2}}{\displaystyle \frac{1}{w_i^4}}-{\displaystyle \frac{N{g}_{ii}}{\sqrt{2\pi }m{\omega }_z^2}}{\displaystyle \frac{1}{w_i^3}})=\hfill \\ & {\displaystyle \frac{2N{g}_{12}}{\sqrt{\pi }m{\omega }_z^2}}{\displaystyle \frac{w_1^2+w_2^2-2{\delta }^2}{{\left(w_1^2+w_2^2\right)}^{5/2}}}{e}^{-{\displaystyle \frac{{\delta }^2}{w_1^2+w_2^2}}},\hfill \end{array}\end{array}$$

(7a)

$$\begin{array}{c}\delta \left[1-{\displaystyle \frac{4N{g}_{12}}{\sqrt{\pi }m{\omega }_z^2}}{\displaystyle \frac{e^{-{\delta }^2/\left(w_1^2+w_2^2\right)}}{{\left(w_1^2+w_2^2\right)}^{3/2}}}\right]=0,\end{array}$$

(7b)

where $i=1,2.$

We can make a few simplifying assumptions to arrive at a more intuitive understanding by setting $g_{11}=g_{22}=g.$ In this limit, $w_1=w_2=w$, and if we further restrict our attention to small $\delta$,

$$w^3=\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{N\left(g+g_{12}\right)}{2m{\omega }_z^2}},$$

(8)

and

$$\delta =\sqrt{2}{\left({\displaystyle \frac{2N}{\pi }}{\displaystyle \frac{g+g_{12}}{2}}\right)}^{1/3}\sqrt{log\left[{\displaystyle \frac{2g_{12}}{g+g_{12}}}\right]}.$$

(9)

We see evidence of a transition at $g_{12}=g.$ For $g_{12}>g,$ the solution is given approximately by Equation (9), and for $g_{12}<g,$ $\delta =0.$ This transition is continuous in the sense that $\delta$ approaches the same value from either side of the transition. More sophisticated numerical techniques are required to arrive at a full solution to Equation (7), but we have gleaned a few key insights from the analysis presented thus far. From Equation (8), it is clear that the widths of the density clouds should increase with increasing atom number and atomic interactions, and should decrease with increasing trap frequency. It is interesting to note that Equation (9) does not appear to depend on the linear trap frequency. This may be the result of our simplifying assumptions, or it could be the result of the symmetry inherent to our harmonic trapping potential. It is well known that asymmetry, whether it be in atomic masses, atom number, or applied potential, can alter the traditional criteria for the immiscibility transition.

3. Computational Model

Now, we turn our attention to the dynamical evolution of the two condensates in the trap. Using a split-step Fourier method (SSM) (see Appendix A) [22], we solve Equation (2) over a finite spatio-temporal grid with the parameters outlined in Table 1. The system is initialized, with each condensate taking a Gaussian form with a width equal to the oscillator length ${\ell }_z=\sqrt{\hslash /m{\omega }_z},$ and the central locations displaced symmetrically about the trap center by a distance of $8{\ell }_z.$ This separation is chosen to ensure minimal initial overlap without introducing excessive kinetic energies throughout the evolution. We choose to scale the separation with the trap frequency in this way to have approximately equivalent kinetic energies upon impact at the trap center. A plot of the initial configuration for a trap frequency of $2\pi \times 10$ Hz is shown in Figure 1. After being released, the condensates are drawn toward the center of the trap and deform slightly in response to the intra-atomic interactions. In the limit of negligible kinetic energy, the condensates would take on a shape given by the Thomas–Fermi approximation [14]. As they approach the center of the trap, the inter-atomic interactions increase with increasing overlap. We expect the condensates to partially transmit and/or reflect in a way that depends on the magnitude of the inter-particle interactions.

The SSM can be used to solve for the evolution of the condensates in the appropriate limit for which the Gross–Pitaevski equation can be applied. From the initial configuration, we use the SSM to propagate the wave function according to the Equation (2). As the condensates evolve, we track their central location and the second and third moments relating to the width and skewedness of their density distributions, respectively. We vary the interaction constant over the expected critical regime and the trapping frequency from 5 to 50 Hz. Furthermore, we have numerically verified that the specifics of the initial state have minimal impact on the stationary state after a few collisions have occurred.

4. Results

We now use the numerically computed wave functions ${\psi }_1(x,t)$ and ${\psi }_2(x,t)$ to calculate the mean position and higher order moments of $\psi$. Let $\delta$ denote the difference in position between the right and left condensate; then,

$$\delta =\int dz\left(\right|{\psi }_1{|}^2-\left|{\psi }_2{|}^2\right)z,$$

(10)

and we can further define

$${\mu }_m\left({\psi }_i\right)=\int dz{\left|{\psi }_i\right|}^2{\left(z-\langle z\rangle \right)}^m$$

(11)

as the higher order moments. Figure 2 illustrates some of the raw data for $\delta$ that we collected in three simulations using a trap with a frequency of 30 Hz. In the case of relatively small or large interactions, we observe a rapid approach to a dynamical equilibrium characterized by a well-defined frequency of motion.

In the regime with small interactions, the condensates mostly transmit upon collision, with minor slowing in response to an energy exchange between kinetic and interaction terms. This is evidenced in Figure 2a,b. The primary mode of oscillation can be extracted from the Fourier transform of $\delta \left(t\right)$, and we find for non-zero interactions that the oscillation frequency is generally smaller than the trap frequency in this regime.

As we move closer to the critical regime, we observe diverging transient times. Figure 2c,d illustrate a typical example of a system with an interaction strength near the expected critical value. There is no clear sign of approach to dynamical equilibrium over the specified time evolution. It is possible that at a later time, this system may settle into a mode of oscillation, but it is unclear from the present work if this is the case. We observe diverging transient times in this regime over a broad range of trap frequencies, justifying the claim that the system does indeed exhibit critical behavior.

For even larger interactions, we observe yet another rapid approach to a different kind of oscillation characterized by periodic reflections of the condensates, or “bouncing”. The dominant feature in the corresponding Fourier spectrum (Figure 2f) indicates the existence of an average, and there appears to be a sub-leading term which captures the frequency of the bouncing mode. This frequency is generally larger than the trap frequency, although never quite double, as would be expected in the case of infinite repulsion at the trap center.

We perform similar analyses for a range of interactions and trapping frequencies. By taking the time average of $\delta \left(t\right)$ for each simulation, we can compare our dynamical results to the corresponding static case, i.e., the ground state of two interacting Bose–Einstein condensates. The results are summarized in Figure 3. There are a few key similarities between the two plots. We see in both systems that there is zero average separation for interactions below the critical value, and the magnitude of the separations follow similar trends generally speaking. As the interaction increases, so do the separations in both cases. The separation also decreases in response to stronger trapping frequencies. The first key difference to point out is that there is an alteration in the critical line for the dynamic setting. Specifically, we observe a dependence on the trapping frequency which is not evident in the static case. For stronger traps, it appears that there is an extended dynamical mixing regime, with the transition to the bouncing regime moved to higher critical repulsions, $g_{12}>\sqrt{g_{11}{g}_{22}}$, as indicated in Figure 3b.

Please note that the 1D approximation for the Gross–Pitaevskii equation becomes less accurate as the trapping frequencies in the tightly confined directions increase. Specifically, when the energy associated with the interactions within the condensate becomes comparable to the energy gap of the transverse confinement, the condensate may start to occupy higher transverse modes, an effect that is not captured by the effective 1D approach.

Another method of analysis involving the higher order moments leads us to further insights regarding the configuration of the binary mixture as a result of dynamical evolution in the trap. Figure 4 contains plots of the first three moments (time-averaged) versus interaction strength for a selection of trap frequencies. While Figure 4a contains mostly the same information discussed in Figure 3, perhaps we can more clearly see the dependence of the transition on trapping frequency. Again, for larger trapping frequencies, it would appear that the transition shifts to larger values. The third moment pictured in Figure 4c also shows signs of the transition. This can be understood by noting that the third moment addresses the asymmetry of the distribution, and for large interactions, the condensates are repelling each other with a force great enough to induce deviations from the symmetric Gaussian structure they began with. The second moment, on the other hand, shows no sign of a transition in $g_{12}.$ This can be understood as the second moment capturing the width of the distribution, which is more strongly dependent on the intra-atomic interactions and the trapping frequency. As we increase the trapping frequency, ${\ell }_z$ becomes smaller, and the condensates are initialized with smaller widths. It is also interesting to note the slight decrease in width corresponding to an increase in interaction strength. This stands in opposition to the analysis presented for the static case, and the resolution to this discrepancy has not yet been uncovered.

5. Discussion

Collision dynamics of condensate mixtures has only been given limited attention in experimental research to date [15,23]. From our numerical simulations, we have reason to believe that non-trivial behavior exists during such collisions as a result of the nonlinear interactions, and that these effects may be observable in a suitable laboratory environment. Specifically, we have identified an extended critical regime close to the miscible–immicscible phase boundary, where the time scale of equilibration diverges, akin to what was observed in Ref. [15]. These results call for experiments that enable observation of long evolution times, the ability to suspend disparate atomic clouds in proximity to one another, and the ability to tune their interactions via Feshbach resonances. Collision dynamics of bosonic mixtures should be observable in laboratories with sufficiently wide time windows of up to 100ms, given that the relevant time scales of the transient phenomena under consideration fall within a range of 15–50 ms.

Our numerical simulations of colliding rubidium BECs in an anisotropic harmonic trap indicate dynamical signatures of a possibly modified miscibility phase transition. The transition seems to occur at values of $g_{12}>\sqrt{g_{11}{g}_{22}}$ which increase with greater trap frequency. Three dynamical phases have been characterized by the relative strength of atomic interactions. In future work, we plan to continue working toward expanding the current capabilities of our SSM to tackle systems in two and three dimensions for the study of new phases resulting from rotational degrees of freedom, and more exotic condensate geometries.