Regular and chaotic behavior of collective atomic motion in two-component Bose-Einstein condensates

The first experimental observation of interference fringes between two freely expanding Bose-Einstein condensates (BECs) has launched the fascinating possibility to probe new quantum phenomena on macroscopic scales related to the superfluid nature of the Bose condensates Andrews1997 . A chief effect is the Josepshon analog, in which the phase differences between two trapped BECs in a double-well potential can generate Josephson-like current-phase oscillations via collectively atomic tunneling through the barrier Smerzi1997 ; Raghavan1999 . Nevertheless, the nonlinearity of tunneling effects produces novel anharmonic Josephson oscillations and macroscopically quantum self-trapping (MQST) effects. The potentially existing such interesting phenomena indeed triggers a wide variety of theoretic investigations in the settings of two-component BECs systems in the double-well potential JuliaDiaz2009 ; Sun2009 ; Satija2009 ; Messeguer2011 , triple-well Richaud2018 and the single BEC system in the triple-well Liu2007 , four-well Liberato2006 and multiple-well potentials Cataliotti2001 ; Nigro2018 , to cite a few. The effective Josephson dynamics can also be observed in some systems without an external potential aba . More importantly, several experiments have successfully observed coherence oscillations through quantum tunneling of atoms between the condensates trapped in a double-well potential Albiez2005 ; LeBlanc2011 ; Pigneur2018 , whereas similar phenomena were also seen from attractive to repulsive atomic interactions by experimentally changing atomic scattering lengths through Feshbach resonances Spagnolli2017 . In addition to the regular coherent oscillations mentioned above, the studies on the possibility of the appearance of chaotic motions, by applying an externally time-dependent trap potential, are also a fascinating task with findings that deserve further experimental justification Abdullaev2000 ; Lee2001 ; Jiang2014 ; Tomkovic2017 . Chaotic behavior is also studied in the dynamics of of three coupled BECs fra .

In the present work, we consider the two-component BECs system with atoms in two different hyperfine states in a single-well harmonic trap potential. Such a two-component BEC system provides us an ideal platform for the study of intriguing phenomena, for example, mimicking quantum gravitational effects in Ref. Syu2019 . Using Feshbach resonances to experimentally tune the scattering lengths of atoms allows us to construct the diagram in Fig. 1, showing the typical ground-state wave function of the binary condensates Thalhammer2008 ; Tojo2010 ; Catani2008 ; Wacker2015 ; Moses2015 . The effective parameter characterizing the miscibility or immiscible regime of binary condensates is mainly determined by the value of $\Delta=g_{11}g_{22}-g_{12}^{2}$, where $g_{ij}$ denotes the atomic interaction between $i$ and $j$ atoms Riboli2002 ; Papp2008 ; Sasaki2009 . The parameter regime for $\Delta>0$ corresponds to the miscible distribution of two-component condensates, whereas for $\Delta<0$ the different species of atoms repel each other and the distribution becomes immiscible Riboli2002 ; Papp2008 ; Sasaki2009 . In particular, we focus on the situations of the scattering lengths in the immiscible regime, where the condensate of one of the components is located at the center of the potential trap, and can effectively be regarded as an effective potential barrier between bilateral condensates of the second component of the system. Experimentally, spatially asymmetric bilateral condensates can be prepared by means of a magnetic field gradient McCarron2011 ; Eto2015 ; Eto2016 , allowing the possibility to observe collectively coherence oscillations of atoms between them. This is an extension of the works of Smerzi1997 ; Raghavan1999 , where the system of a single BEC in an external double-well potential is studied. The collective oscillations of atoms through quantum tunneling between the condensates centered at each of the potential wells have been observed where the full dynamics of the system can be further simplified, and turned into an analog pendulum’s dynamics in the so-called two modes approximation Burchianti2017 ; Smerzi1997 ; Raghavan1999 ; Ananikian2006 . In our setting, the dynamics of the system can also be approximated by analogous coupled pendulums dynamics using the variational approach as well as two modes approximation. One can then reproduce regular oscillations around the stable states studied in Refs. Smerzi1997 ; Raghavan1999 within some parameter regime by considering the effects from time-averaged trajectories of the central condensate. Additionally, the collective motion of atoms, if departing from the unstable states, could show chaotical behavior, oscillating between bilateral condensates due to the time-dependent driving force given by the dynamics of the condensate centered in the middle Abdullaev2000 ; Lee2001 ; Jiang2014 ; Tomkovic2017 . The scattering lengths in the parameter regime, which give the chaotic dynamics of the system, will be useful as a reference for experimentalists to observe the effects.

To be concrete, consider a mixture of ${}^{85}\text{Rb}$ $|1\rangle=|F=2,m_{F}=-2\rangle$ and ${}^{87}\text{Rb}$ $|2\rangle=|F=1,m_{F}=-1\rangle$ by ignoring the small mass difference Papp2008 . The scattering lengths we choose are estimated to be $a_{22}=50a_{0}$ and $a_{12}=85a_{0}$ with $a_{0}$ the Bohr radius. The existence of the Feshbach resonance at $160\text{G}$ permits us to tune the scattering length $a_{11}$ of ${}^{85}\text{Rb}$ atoms from $50a_{0}$ to $150a_{0}$ to be specified in each case later Papp2008 . Let us consider two-component elongated BECs in quasi-one-dimensional (1D) geometries with the trap potential Becker2008 $\omega_{x}=2\pi\times 12\text{Hz},\,\omega_{y,z}=2\pi\times 100\text{Hz}$. We take the numbers of atoms $N_{1}=500$, and $N_{2}=1000$ as an example. The experimental realization of the Josephson dynamics by two weakly linked BECS in a double-well potential in a single BEC system is confirmed in Ref. Albiez2005 . Both Josephson oscillations and nonlinear trapping are observed by tuning the external potential. Such phenomena are also explored in two-component BECs system in a double-well potential where apart from their own atomic scattering processes the Rabi coupling between the atoms from each components of the BECs is introduced. Tuning the scattering lengths using the associated Feshbach resonances and/or the Rabi coupling constant allow us to observe the transition from Josephson oscillations to nonlinear self-trappings zib . Our proposal also in the two-component BECs system but in a single-well potential suggests a novel setup, although involving more rich dynamical aspects of condensates, the only tunable parameters are atomic scattering lengths again through Feshbach resonances to achieve the miscible condensates that needs more experimental endeavor as compared with zib .

Our presentation is organized as follows. In Sec. II, we first introduce the model of the binary BECs system, and the corresponding time-dependent Gross-Pitaeviskii (GP) equations for each of condensate wave functions. In this two-component BECs system, we study the evolution of the wave functions within a strong cigar-shape trap potential so that it can effectively treated as a quasi one-dimensional system. We further assume that the condensates start being displaced from their ground state in the immiscible regime where one of the condensates is distributed in the trap center while the other condensate surrounds the central condensate on its two sides. Then we propose the ansatz of the Gaussian wave functions, which can reduce the dynamics of the system into few degrees of freedom, such as the center of the central condensate as well as the population imbalance and the relative phase difference between bilateral condensates in the two-mode approximation. We later obtain their equations of motion with the form of coupled oscillators equations. In Sec. III, we analyze the stationary-state solutions and their stability property. In Sec. IV, we examine the time evolution of the system around the stable stationary state, showing the regular oscillatory motion. In Sec. V, the study is focused on when the system starts from near unstable stationary-state solutions, exhibiting the chaotic behavior by numerical studies. The Melnikov homoclinic method is adopted for analytically studying the existence of chaos. Concluding remarks are in Sec. VI. In Appendix, we provide more detailed derivations or approximations to arrive at coupled oscillators equations from the time-dependent GP equations using the appropriate ansatz of Gaussian wave functions.

We consider the binary BECs for same atoms in two different hyperfine states confined in a strong cigar-shape potential with the size of the trap $L_{x}$ along the axial direction,taken in the $x$ direction, which is much larger than $L_{r}$ along the radial direction. This system can be effectively treated as the pseudo one-dimensional system with the Lagrangian described as Garcia1996 ; Salasnich2002

where the external potential in the axial direction is given by $V_{\rm ext}(x)=m\,\omega_{x}^{2}\,x^{2}/2$. The mass of atoms is $m$ and the coupling constants between atoms are given in terms of scattering lengths $a_{ij}$ as $g_{ij}=2\hbar^{2}a_{ij}/mL_{r}^{2}$. The wave function $\psi_{j}$ is subject to the normalization

with the numbers of atoms $N_{j}$ in each hyperfine states $j$. The time-dependent GP equations for each component of the BEC condensates are given by Pethick2008

The full dynamics of the condensates in a harmonic trap potential can be explored by solving the GP equations numerically with given initial wave functions. However, in order to extract their feature analytically, we will propose an ansatz of the wave functions with several relevant parameters. All results from solving the equations of motion for the introduced parameters given by the variational approach will be justified by comparing with those given by numerically solving the GP equations. Moreover, we will rescale the spatial or temporal variables into the dimensionless ones by setting $t\rightarrow t/\hbar\omega_{x}$ and $x\rightarrow x/\sqrt{\hbar/m\omega_{x}}$. Henceforth, the dimensionless $t$ and $x$ will be adopted, and their full dimensional expressions will be recovered, if necessary.

Experimentally, the values of parameters $g_{ij}$, which are related to scattering lengths $a_{ij}$, can be tuned by Feshbach resonances  Thalhammer2008 ; Tojo2010 ; Papp2008 . In fact, we might explore the evolution of condensate wave functions, in which $\psi_{1}$ and $\psi_{2}$ are initially prepared in the immiscible regime, determined by the parameter $\Delta=g_{11}g_{22}-g_{12}^{2}$ when $\Delta<0$ Riboli2002 ; Papp2008 ; Sasaki2009 as illustrated in Fig. 1. Let us now assume the general Gaussian wavefunction $\psi_{2}$ centered at the potential trap center as Garcia1996 ; Salasnich2002 ; Lin2000

where four time-dependent variational parameters are the position of the center $\xi(t)$, the width $w(t)$, and $\alpha(t)$ (slope) and $\beta(t)$ $({\rm(curvature)}^{1/2}$) of the wave function. We also assume that the ground-state configuration of $\psi_{1}(x,t)$ is of the two-mode form Smerzi1997 ; Raghavan1999 ; Burchianti2017 ; Ananikian2006

The involved Gaussian-like spatial functions $\psi_{L}(x)$ and $\psi_{R}(x)$ together with the initial wave function $\psi_{1}(x,t=0)$ with width $\sigma_{0}$ will be determined numerically by the stationary-state solutions of the GP equations (3) and (4) with the forms shown in Fig. 2, which is symmetric with respect to $x=0$. $\psi_{L}(x)$ and $\psi_{R}(x)$ presumably serve as a good basis to express the spatial part of the wave functions of bilateral condensates with normalization conditions $\int dx\,|\psi_{L,R}|^{2}=1$, and also negligible overlap $\int dx\,\psi_{L}\psi_{R}\approx 0$ for the weak link between bilateral condensates in the immiscible regime. Note that the analytical approach based on the wave functions above gives transparent and consistent results as compared with those from numerically solving the full GP equations in the case of Josephson oscillations but for the motion of MQST discrepancy arise as will be seen later where one should look for the wave functions of bilateral condensates being not restricted to be symmetric with respect to $x=0$. Thus, the time-dependent part of the wave function can be parametrized as

with the phase $\phi_{L,R}(t)$ and the amplitude $\sqrt{N_{L,R}(t)}$ given by the number of atoms in each sides of the condensates. The total number of atoms ($j=1$) is $N_{1}=N_{L}+N_{R}$. To realize quantum coherence between left and right condensates, and also atomic number oscillations, it finds more convenient to introduce the population imbalance

and relative phase

When all atoms locate in the left(right) side of the central condensate, $z=1(-1)$.

The time evolution of coherent oscillations between bilateral condensates surrounding around the central condensate is mainly determined not only by the self-interaction of atoms Smerzi1997 ; Raghavan1999 , but also the interspecies coupling constant $g_{12}$ with the wave functions $\psi_{1}(x,t)$ and $\psi_{2}(x,t)$, which give an additional effect to influence the dynamics of coherent oscillations due to the wave-function overlap between them. Also, notice that although the initial wave function of the central condensate is centered at $x=0$, which is the minimum of the external trapped potential, the later evolution of the wave function will find its stationary state with the center located at $x=\xi_{0}$, moving around that position with nonzero kinetic energy. Thus, in the presence of $\psi_{2}(x,t)$ shown in Fig. 2 together with the external trapped potential, the resulting effective potential experienced by the wave function $\psi_{1}(x,t)$ in the GP equation in (3) can be seen in the shape of the double-well potential. However, the dynamics of the central condensate leads to the effective potential time dependent, contributing its effects to quantum coherence oscillations between bilateral condensates. This will be the main purpose of studies in Secs. IV and V. This is an extension of the work in Refs. Smerzi1997 ; Raghavan1999 with additional degrees of freedom of the central condensate. When the initial conditions of the central or bilateral condensates are chosen near their stable-state configurations, the dynamics of coherence oscillations modifies slightly the behavior found in Refs. Smerzi1997 ; Raghavan1999 . Nevertheless, as their initial conditions are chosen near the unstable-state configurations, the chaos occurs. Although the analogous coupled pendulums dynamics, which is reduced from the full dynamics of the system and later gives essential information on the time evolution of the condensate wave functions, may produce relatively different trajectories from the solutions of the time-dependent GP equations, they are still very useful to provide us analytic analysis on the chaos dynamics by adopting the Melnikov Homoclinic method.

Substituting the ansatz of the wave function (5) and (6) into the Lagrangian (1), and carrying out the integration over space, the effective action then becomes the functional of time-dependent variables $\alpha$, $\beta$, $\xi$, $w$, $z$, and $\phi$. Their time evolutions follow the Lagrange equations of motion derived from this effective Lagrangian, discussed in detail in the Appendix. When the wave function $\psi_{2}(x,t=0)$ is slight away from the stationary state, the dynamics of the condensate just undergoes small amplitude oscillations around that state. We will further assume that the displacement $\xi$ and the wave function width, deviated from $\sigma_{0}$ and defined by $w=\sigma_{0}+\sigma$, are small as compared with $\sigma_{0}$, namely $\xi\ll\sigma_{0}$ and $\sigma\ll\sigma_{0}$.

Retaining terms up to linear order in $\xi$ and $\sigma$ of interest, we have found the Lagrange equations for the parameters by the variational approach as

The quantities that appear in the equations above are the integrals involving the condensate wave functions, namely

The expressions of $k$ (16) and $\kappa$ (17) are determined by the wave function overlap between the left and right condensates, and thus $k\sim\kappa\ll 1$ by considering the weak link between bilateral condensates that depends not only on atomic scattering lengths but also on the number of the condensates with $N_{2}>N_{1}$. Nevertheless, it will be seen that since the frequency of collectively atomic oscillations $\omega_{J\xi}\propto\omega_{x}\sqrt{k}$ with $\omega_{x}=2\pi\times 12\text{Hz}$ we adopt, the corresponding oscillation timescale $t\propto 1/(\omega_{x}\sqrt{k})$ is constrained to be not larger than the typical lifetime of this type of the condensates of order $\sim 10\text{s}$ Becker2008 , giving $k>10^{-4}$ where $N_{2}$ cannot be too large. The variation of the Gaussian width induces a time-dependent modification to the parameter $k+\kappa\sigma/\sigma_{0}$ in (11) and (12). In the single BEC experiment Smerzi1997 ; Raghavan1999 , the similar modification to $k$ can be achieved by the time-dependent laser intensity or the magnetic field with the possibility to induce the Shapiro effect Grond2011 , analog Kaptiza pendulum Boukobza2010 and even chaotic motions Abdullaev2000 ; Lee2001 ; Jiang2014 ; Tomkovic2017 . Here, since $\sigma\ll\sigma_{0}$, the time-dependent modification is relatively small compared with $k$ to be ignored in this work. Also, for simplicity, we consider that the spatial part of the wave function $\psi_{1}(x,t)$ is symmetric with respect to $x=0$ shown in Fig. 2, thus giving $\Delta E=0$ in (12). The coupling of the population imbalance $z$ to the relative phase $\phi$ in (12) is given by the strength $\Lambda$ in (19), a tunable value by changing the scattering lengths of atoms in bilateral condensates.

It is known that the dynamical equations of pair variables $(z,\phi)$ in the single BEC system show an analogy to that of a nonrigid pendulum model, whose length depends on the angular momentum $z$ Raghavan1999 . Nevertheless, in our binary BECs system, the pair variables $(z,\phi)$ are also coupled to the displacement of the wave function $\psi_{2}$ with the coupling constant $\eta$ in (20) depending on the interspecies coupling constant $g_{12}$ of atoms, and also the wave function overlap between the bilateral and central condensates. With the parameters shown in Fig. 1 and the chosen scattering lengths obeying $\Delta<0$, which lead to the spatial parts of the wave functions schematically shown in Fig. 1, $\eta$ is of order one $\eta\sim 1.0$, and $\varepsilon\sim 0.1$ due to cancellation between the contributions from the left and right wave functions. In the case of the initial deviation of the wave-function width $\sigma(0)=\dot{\sigma}(0)=0$, the time-dependent $\sigma$ is driven by $\xi$, giving $\sigma\sim\varepsilon\xi$, which in turn leads to the term $\varepsilon\sigma$ in the equation of $\xi$ (13) of the order of $\varepsilon^{2}\xi$ with $\varepsilon^{2}\ll\eta$ to be safely ignored. From the viewpoint of the wave function $\psi_{2}$, the presence of $\psi_{R}$ and $\psi_{L}$ will make trap potential narrower and affects slightly the natural vibration frequency in $\xi$ and give a modified frequency $\omega_{\xi}$ in (22). Since we just focus on various types of coherent oscillations between bilateral condensates with the dynamical variables $z,\,\phi$, apart from their mutual coupling, they are also coupled to the displacement of the central condensate $\xi$. Ignoring the $\varepsilon\sigma$ term in (13), the small deviation from the wave-function width $\sigma$ is decoupled. In fact, in this approximation, once the time-dependence of $\xi$ is found, the equation of $\sigma$ in (14) can be solved, and then plugging all solutions to the equations (86) in Appendix can find $\alpha$ and $\beta$. Also notice that taking into account the time-dependent $\sigma$ by involving (14) certainly can improve the agreement with the full numerical results. Thus, after the further replacement $\xi\rightarrow\sqrt{2N_{2}/N_{1}}\xi$, the set of the key equations can be simplified as

These are the main results of this section and we will term them the coupled pendulums (CP) dynamics. The corresponding Hamiltonian then reads as

In our binary BECs case, the above Hamiltonian also shows an analogy between the dynamics of coupled pendulums and the BEC dynamics as in Ref.  Raghavan1999 .

In the two-component BEC system, the relevant equations of motion for describing quantum coherent oscillations between bilateral condensates in the presence of the central condensate are effectively described by the dynamics of two-coupled pendulums. In this case, the stationary states with the vanishing time-derivative terms now obey

As a result, the stationary relative phase is given by

Also, Eq. (30) leads to a relation between the population imbalance $z_{0}$ and the displacement of the central condensate $\xi_{0}$

Inserting the relation (32) into Eq. (29) together with $\phi_{0}=0$ or $\pm\pi$, the population imbalance is found to be

where $\Lambda_{\text{eff}}=\Lambda-\omega_{\xi}^{-2}\eta^{2}$. The solution $z_{0}\neq 0$ exists as long as $\Lambda_{\text{eff}}/2k\geq 1$ or $\Lambda_{\text{eff}}/2k\leq-1$ in the MQST state for the atom system. Notice that $\Lambda_{\text{eff}}$ can be positive or negative, and it is tunable by changing either $g_{11}$ or $g_{12}$ to vary $k$ (16), $\Lambda$ (19), and $\eta$ (20) . There are four different types of stationary states, $(z_{0},\,\phi_{0})=$ $(0,\,0)-$I, $(0,\,\pm\pi)-$II, $(\pm\sqrt{1-4k^{2}/\Lambda_{\text{eff}}^{2}},\,\,0)-$III and $(\pm\sqrt{1-4k^{2}/\Lambda_{\text{eff}}^{2}},\,\,\pm\pi)-$IV, summarized in Table 1. Furthermore, the corresponding energy for each stationary states can be computed through (27) using the relation of stationary-state solutions (32), where the results are listed in Table 2. In the case $\Lambda_{\text{eff}}/2k>1$, the ground state is in configuration I where the bilateral condensates have the same populations $z_{0}=0$, thus the central condensate is centered at $\xi_{0}=0$, consistent with Fig. 1. As for $\Lambda_{\text{eff}}/2k<-1$, the ground state shows the population imbalance for $z_{0}\neq 0$ and the wavefunction $\psi_{2}$ is centered at $\xi_{0}\neq 0$ due to (32) for a given $z_{0}$ also seen in Fig. 1.

We next study the stability of the system around the equilibrium solutions. To do so, we consider the small perturbations around $\xi_{0}$, $z_{0}$, and $\phi_{0}$ defined by

Substituting (34)–(36) into Eqs. (24)–(26), and using $\xi_{0}=\omega_{\xi}^{-2}\eta z_{0}$, one obtains the linearized equations of motion for $\delta\xi(t)$ and $\delta z(t)$ as

Notice that the linear term in $\delta\phi$ vanishes where $\phi_{0}=0$ or $\pm\pi$ for the stationary-state configurations is evaluated in obtaining the linearized equations of motion. Considering the oscillation motion of $\delta\xi(t)$ and $\delta z(t)$ with $\delta\xi(t),\delta z(t)\propto\exp(i\omega t)$, we obtain the eigenfrequencies

In addition to the nature frequency $\omega_{\xi}$ for the central condensate, there exists an oscillating frequency

with

which manifests the nature frequency for quantum oscillations between two bilateral condensates $\omega_{J}$ with the effects from the dynamics of the central condensate Smerzi1997 . However, the true eigenfrequencies are the mixture of these two fundamental frequencies $\omega_{\xi},\,\omega_{J\xi}$ as they couple. With the parameters in the immiscible regime, typically for small $k$ where the wave-function overlap between bilateral condensates is small, this gives $\omega_{\xi}\gg\omega_{J\xi}$. In the two-component BECs system, the existing oscillation frequency $\omega_{\xi}$ of the central condensate will drive the dynamics of the pair variables $(z,\phi)$ into the oscillatory motions with frequencies $\omega_{-}$ and $\omega_{+}$ where $\omega_{+}\gg\omega_{-}$.

Later we will manipulate the initial conditions of $z(0),\,\phi(0)$ as well as $\xi(0),\,\dot{\xi}(0)$ to effectively switch off the fast varying mode of frequency $\omega_{+}$, giving the results presumably to those in Ref. Smerzi1997 . Then, the presence of the fast varying mode for more general initial conditions will also be studied to see its effects on quantum coherent phenomena $(z,\phi)$. We will learn that, in the case of relatively small $k\sim 10^{-3}$ with very small wavefunction overlap, the time average over the evolution of the fast varying mode is adopted to find the deviation from the results in Ref. Smerzi1997 . Those trajectories are regular motion moving around the stable states of the system listed in Table 2. However, we also probe the parameter regime with relatively large wave-function overlap with relatively large value $k\sim 10^{-2}$ for amplifying the influence of the fast varying mode that can even drive the dynamics of $z$ and $\phi$ to the chaotic behavior, if the system starts from the unstable state also seen in Table 2. In our system, the dynamics of quantum coherence oscillations $(z,\phi)$ is very sensitive to the tunneling energy $k$. For even larger $k$, the proposed wave-function ansatz fails to describe the evolution of $(z,\phi)$ where the dynamics of the system can only be explored by directly solving time-dependent GP equations, which is beyond the scope of the present work.