Manifestation of Superfluidity in Atom-Number-Imbalanced Two-Component Bose–Einstein Condensates

Abstract

Superfluid and dissipative regimes in the dynamics of a two-component quasi-one-dimensional Bose–Einstein condensate (BEC) with unequal atom numbers in the two components have been explored. The system supports localized waves of the symbiotic type owing to the same-species repulsion and cross-species attraction. The minority BEC component moves through the majority component and creates excitations. To quantify the emerging excitations, we introduce a time-dependent function called disturbance. Through numerical simulations of the coupled Gross–Pitaevskii equations with periodic boundary conditions, we have identified a critical velocity of the localized wave, above which a transition from the superfluid to dissipative regime occurs, as evidenced by a sharp increase in the disturbance function. The factors responsible for the discrepancy between the actual critical velocity and the speed of sound, expected from theoretical arguments, have been discussed.

1. Introduction

Superfluidity is the ability of a liquid to flow through tight tubes or narrow slits without dissipation. The experimental discovery [1,2] and theoretical explanation [3] of this macroscopic quantum phenomenon were among the most important achievements of 20th-century physics. While initial studies were restricted to liquid helium, superfluid properties were recently observed and explored in a dilute gas of bosons, called Bose–Einstein condensates (BECs) [4,5]. The advantage of ultracold gases in studies of superfluidity over liquid helium is that they are analytically tractable, weakly interacting systems, and free of surface tension effects.

The phenomenon of superfluidity in quantum liquids and gases has been studied in different contexts. Among them, the breakdown of superfluidity when the flow velocity reaches some critical value is particularly intriguing. The physics behind the suppression of superfluidity and transition to a usual dissipative regime above the critical velocity is linked to the generation of excitations in the system such as phonons, rotons, and vortices when the liquid flows at velocities exceeding some critical value $v_c$, which was predicted by Landau [3,6] $v_c=\min \left(ϵ\left(p\right)/p\right)$, with $ϵ\left(p\right)$ being the energy of the excitation with momentum p. Measurements of the critical velocity in BECs by laser stirring have confirmed that the superfluidity indeed breaks down if the laser beam moves faster than a certain finite speed [7,8,9]. However, the measured values of the critical velocity (∼1.6 mm/s) in the experiment with sodium BEC [7] appeared to be smaller than the speed of sound (∼6.2 mm/s), which were expected from the theoretical arguments. The diameter of the laser spot in the experiment was macroscopic ($\sim 13\mathsf{\mu }\mathrm{m}$), exceeding the intrinsic length scale of the condensate, the so called healing length $\xi =1/\sqrt{8\pi a_{s}n}\sim 0.3\mathsf{\mu }\mathrm{m}$, where $a_s=52a_B$ is the two-body s wave scattering length in units of Bohr radius and $n=1.5\times {10}^{14}$ cm⁻³ is the peak density of the condensate. Under these conditions, vortex shedding by the moving object [10,11], rather than the emission of phonons, was a more realistic phenomenon. Later investigations revealed several factors that could be responsible for the observed discrepancy, such as the size and strength of the probe object [12,13,14], the intrinsic nonlinearity of BEC [15,16,17,18], the system’s dimensionality [19], etc. All of the arguments mentioned indicate the need for further research on the origin of the critical velocity, leading to the breakdown of superfluidity in quantum gases. This is necessary in order to gain a better understanding of this macroscopic quantum phenomenon.

This work aims to explore the manifestation of superfluid properties in the dynamics of quasi-1D binary quantum gases. In our setting, the mixture condensate can support stable localized waves of the symbiotic type [20,21] due to same-species repulsion and cross-species attraction. We consider a significant disparity in atom numbers between the components and a toroidal trap potential for the binary condensate. Under these conditions, the smaller localized component appears immersed in the larger extended component, distributed over the whole integration domain. The minority component is propelled by applying a potential that only affects this component. The majority component is subjected to a uniform (flat) potential, while the minority component experiences a harmonic potential. This causes the minority component to oscillate, similar to a pendulum, with zero velocity at the turning points and maximum velocity at the center. By adjusting the amplitude of these oscillations, we can observe different scenarios where the majority component is either in a superfluid state or not. In cases where the majority component is not in a superfluid state, a sound wave is generated, and the energy of the oscillating “pendulum” is transformed into sound. In fact, suppression of superfluidity and transition to a dissipative regime occurs when the velocity of the probe object (localized smaller component) reaches a critical value, leading to the emergence of density waves superimposed on the uniform background of the larger component. We introduce a measure called disturbance to quantify the emerging excitations in the majority component, whose superfluidity is being tested.

2. Model and Governing Equations

We use the following coupled dimensionless Gross–Pitaevskii equations (GPE) to describe the dynamics of a two-component BEC in a quasi-1D trap potential

$$i{\displaystyle \frac{\partial {\psi }_i}{\partial t}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{{\partial }^2{\psi }_i}{\partial x^2}}+(g_i|{\psi }_i{|}^2+g_{ij}|{\psi }_j{|}^2){\psi }_i=V_i{\psi }_i,i,j=1,2,i\ne j,$$

(1)

where ${\psi }_i(x,t)$ are the mean field wave functions of the condensate components; $g_i,g_{ij}$ are the same-species and cross-species atomic interactions coefficients, respectively; and $V_i$ is the external potential for the i-th component. Specifically, we use a uniform potential for the majority component and a harmonic trap for the minority component. Such a component-selective potential was employed in studies of two-vortex collisions in particle-imbalanced heteronuclear BEC mixtures [22]. The dimensionless quantities in these equations were scaled using the frequency of the radial confinement ${\omega }_{\perp }$, atomic mass m, and radial harmonic oscillator length $l_{\perp }=\sqrt{\hslash /m{\omega }_{\perp }}$ as follows: time $t\to t{\omega }_{\perp }$, space $x\to x/l_{\perp }$, and wave function ${\psi }_i\to \sqrt{2|a_s|}{\psi }_i$, with $a_s$ being the atomic s-wave scattering length. When considering imbalanced settings with different atoms in the components $N_{1,2}=\int {\left|{\psi }_{1,2}\left(x\right)\right|}^{2}dx$, we denote ${\psi }_1(x,t)$ and ${\psi }_2(x,t)$ as the minority and majority components, respectively. The system of coupled GPE (1) is suitable for the description of homonuclear BEC mixtures, such as two internal states of ⁸⁷Rb [23] or ³⁹K [24]. In these mixtures, all of the atoms belong to the same element but can occupy two different internal spin states. We use the periodic boundary conditions for the GPE (1), implying the condensate is held in a toroidal trap potential, similar to that reported in [25,26].

3. Numerical Simulations

For numerical simulations of the coupled GPE (1), we use the split-step Fourier method, which is the most commonly employed numerical scheme for exploring the propagation of light pulses in nonlinear optical media and the evolution of BEC. In this approach, dispersive and nonlinear terms are integrated separately, and then the results are combined to form the complete solution. The linear dispersive term is evaluated in the frequency domain using fast Fourier transform (FFT), while the nonlinear term is evaluated in the spatial domain.

Over the years, several enhanced versions of the split-step FFT method have been developed featuring a higher accuracy and speed. One of the efficient and accurate algorithms closely related to the split-step Fourier method is the fourth-order Runge–Kutta in the interaction picture (RK4IP) method [27,28]. The algorithm transforms the problem into an interaction picture, allowing for the use of conventional techniques to advance the solution a step forward. To explain the idea of the method, we write the governing equation in the form

$${\displaystyle \frac{\partial \psi }{\partial t}}=-i\left[\hat{D}+\hat{N}\right]\psi ,$$

(2)

where $\hat{D}=-{\partial }^2/2\partial x^2$ is the linear dispersion operator, while $\hat{N}=\left[g_i|{\psi }_i{|}^2+g_{ij}{\left|{\psi }_j\right|}^2+V_i\right]$ is the nonlinear operator containing non-derivative terms. In the following, we omit the indices and consider all the steps for a single component, as the extension for the coupled GPE is straightforward.

The condensate wave function in the interaction picture can be presented in the form

$${\psi }_I=e^{i(t-t^{\prime })\hat{D}}\psi ,$$

(3)

where $t^{\prime }$ is the time separating the interaction picture from the normal one. The corresponding evolution equation

$${\displaystyle \frac{\partial {\psi }_I}{\partial t}}=-i{\hat{N}}_I{\psi }_I,{\hat{N}}_I=e^{i(t-t^{\prime })\hat{D}}\hat{N}{e}^{-i(t-t^{\prime })\hat{D}}$$

(4)

can be solved using Runge–Kutta methods, whose efficiency can be greatly enhanced by choosing the separation time as $t^{\prime }=t+\Delta t/2$, with $\Delta t$ being the time step. Then, the algorithm for a one step advancing of the wave function (${\psi }_j\to {\psi }_{j+1}$) is presented as follows [27,28]

$$\begin{array}{ccc}\hfill {\psi }_I& =& e^{-i(\Delta t/2)\hat{D}}{\psi }_j,\hfill \\ \hfill k_1& =& e^{-i(\Delta t/2)\hat{D}}\left[-i\Delta t\hat{N}\left(t_j\right)\right]{\psi }_j,\hfill \\ \hfill k_2& =& -i\Delta t\hat{N}(t_j+\Delta t/2)({\psi }_I+k_1/2),\hfill \\ \hfill k_3& =& -i\Delta t\hat{N}(t_j+\Delta t/2)({\psi }_I+k_2/2),\hfill \\ \hfill k_4& =& -i\Delta t\hat{N}(t_j+\Delta t)e^{-i(\Delta t/2)\hat{D}}({\psi }_I+k_3),\hfill \\ \hfill {\psi }_{j+1}& =& e^{-i(\Delta t/2)\hat{D}}\left\{{\psi }_I+\left[k_1+2(k_2+k_3)\right]/6\right\}+k_4/6.\hfill \end{array}$$

Thus, forwarding the solution one step requires four evaluations of the nonlinear operator $\hat{N}$ and four assessments of the exponential dispersion operator $e^{-i(\Delta t/2)\hat{D}}$, which requires eight FFTs. Overall, the scheme has a global fourth-order accuracy $O(\Delta t^4)$.

We obtain the ground state for the dynamical Equation (1) using the phenomenological damping procedure, first introduced by L. P. Pitaevskii in the context of superfluid helium. Later, the method was applied to numerically find the ground state of a BEC trapped in different potentials [29]. In this approach, the governing equation is presented in a form that includes the terms responsible for the relaxation to equilibrium

$$i{\displaystyle \frac{\partial {\psi }_i}{\partial t}}=(1+i\Lambda )\left[-{\displaystyle \frac{1}{2}}{\displaystyle \frac{{\partial }^2}{\partial x^2}}-(g_i|{\psi }_i{|}^2+g_{ij}|{\psi }_j{|}^2)-V_i-{\mathsf{\mu }}_i\right]{\psi }_i,i,j=1,2,i\ne j,$$

(5)

where $\Lambda <0$ denotes a phenomenological damping term that is assumed to be the same for both components, ${\mathsf{\mu }}_i$ is a chemical potential to be adjusted during the evolution towards the equilibrium state to preserve the number of condensate atoms. A constant value of $\Lambda$ is usually selected based on experimental observations matching the dissipative process. In BEC applications [29], appropriate values are typically $\Lambda \sim [-0.03,-0.5]$.

Numerical simulations start with the creation of a localized state in the particle imbalanced ($N_1\ne N_2$) binary condensate with repulsive intra-component ($g_1<0,g_2<0$) and attractive inter-component ($g_{12}>0$) interactions, as shown in Figure 1a. The density of the majority component $|{\psi }_2{|}^2$ beyond the space occupied by the localized wave $|{\psi }_1{|}^2$ is relatively small and uniform. Here, we are interested in the superfluid and dissipative behaviors exhibited by the majority component.

In the next step, we displace the minority component by $\delta x$ from its initial position at $x=0$ and let it oscillate in the harmonic potential $V_1\left(x\right)=\beta x^2$, as illustrated in Figure 1b. The center-of-mass position of the minority component is evaluated according to its definition

$$\eta \left(t\right)={\displaystyle \frac{1}{N_1}}{\int }_{-L/2}^{L/2}x{\left|{\psi }_1(x,t)\right|}^{2}dx,$$

(6)

where ${\psi }_1(x,t)$ is the solution for GPE (1) at the given time instance. By analyzing the oscillation dynamics of the minority component for different initial displacements, we reveal a critical value $\delta x_{cr}\sim 2.5\pi$ below which the dynamics are fully conservative, while larger displacements demonstrate damped oscillations, as shown in Figure 1b for $\delta x=\pi$ and $\delta x=3\pi$, respecitvely. These numerical experiments emulate the shuttle motion of a probe object or laser beam in a superfluid. It is natural to expect that at sufficiently large oscillation amplitudes, the probe object acquires a velocity greater than the superfluid critical velocity ($v>v_c$) near the minimum of the harmonic trap at $x=0$, producing excitations in the condensate.

Table 1. The maximal velocity v acquired by the minority component in the harmonic trap for different initial shifts $\delta x$, according to numerical solution of the GPE (1).

Initial shift, $\delta x$	$\pi$	$1.5\pi$	$2\pi$	$2.5\pi$	$3\pi$
Max. velocity, v	0.10	0.15	0.19	0.24	0.29

presents the velocity of the localized minority component when it passes through the minimum of the harmonic trap for different initial displacements.

We introduce a time-dependent function called disturbance to quantify the generation and growth of excitations in the background condensate of the majority component [30]

$$D\left(t\right)={\int }_C\left(\right|{\psi }_2{(x,t)|}^2-|{\psi }_2(x,0){{|}^2)}^{2}dx,$$

(7)

where the integration is performed over the spatial domain with the vanishing amplitude of the smaller component $C\in \left[|{\psi }_1(x,t)|\simeq 0\right]$. The latter condition ensures we only consider the excitations on top of the uniform part of the repulsive condensate ${\psi }_2$.

Figure 2 shows the growth rate of excitations for various initial displacements of the minority component in relation to the minimum of the harmonic trap.

As can be seen from Figure 2a, the movement of the minority component does not produce excitations in the majority component, thus the disturbance function Equation (7) does not grow with time, if its initial position is not sufficiently far from the origin ($\delta x<2.5\pi$). In contrast, at a greater initial shift ($\delta x=3\pi$), the disturbance function rapidly grows with time, evidencing the vigorous generation of excitations. These excitations in the majority component appear as density modulations (sound waves) on the initially uniform domain of ${\psi }_2$. It is reasonable to suggest that the generation rate of quasiparticles or density waves starts to grow exponentially once the velocity of the probe object exceeds the critical value. Figure 2b shows the behavior of the integral disturbance $\mathcal{D}$ for different maximal velocities of the minority component moving through the majority component. As can be seen from this figure, the numerically obtained data (red points) nicely fit with the exponential model $\mathcal{D}\left(v\right)=\alpha e^{\gamma v}$.

Different dynamical regimes can be explored by altering either the initial position or the initial velocity of the probe object in the harmonic trap. The simulation results presented in Figure 3 demonstrate the transition from the superfluid regime to the dissipative regime when the initial displacement of the minority component exceeds the critical value. The generation of excitations occurs during time intervals when the localized minority component repeatedly passes the minimum point of the harmonic potential with a supercritical velocity. It is evident that the density modulations on the background condensate $|{\psi }_2{(x,t)|}^2$ do not emerge when the minority component undergoes small-amplitude oscillations and reaches subcritical velocity near the bottom of the parabolic trap, as shown in Figure 3a,c. However, during large-amplitude oscillations, the probe object reaches supercritical velocity, causing the density modulations on $|{\psi }_2{(x,t)|}^2$ to strongly amplify (see Figure 3b,d).

The existence of a definite superfluid critical velocity $v_c$ can be demonstrated by examining the long-time evolution of the oscillating minority component whose initial condition was in the dissipative domain. This implies that it is far shifted from the minimum of the harmonic trap and initially performs large amplitude oscillations, thus attaining a supercritical velocity ($v>v_c$) at the origin. It is reasonable to expect that, after gradually losing its kinetic energy due to generation of excitations in the larger component, the probe object slows down to attain the subcritical velocity ($v<v_c$), and the superfluid regime is reestablished. To illustrate this idea, in Figure 4. we show the long-time evolution of the center-of-mass position and velocity of the minority component for two initial displacements, corresponding to the superfluid and dissipative regimes. As can be seen from this figure, at a smaller displacement ($\delta x=\pi$), the localized wave undergoes free oscillations with a constant amplitude, indicating the presence of the superfluid regime (blue dashed line). On the other hand, a larger displacement ($\delta x=3\pi$) results in damped oscillations, suggesting the onset of the dissipative regime (red solid line). Subsequently, after evolution time $t>3000$, the superfluid regime is reestablished, albeit with some reduction in energy in the minority component due to the excitations created in the majority BEC component.

Numerical simulations were performed using the methods of split-step fast Fourier transform [31,32] and a Runge–Kutta in interaction picture [27]. The integration domain of length $L\in [-6\pi ,6\pi ]$ is employed to accommodate 1024 Fourier modes with a corresponding space step $\Delta x=0.036$. The time step is $\Delta t=0.0005$. The periodic boundary conditions ${\psi }_i(-L/2,t)={\psi }_i(L/2,t)$, $i=1,2$, are adopted to emulate a toroidal trap configuration. The stationary states of the mixture condensate described by GPE (1) are constructed using the Pitaevskii damping procedure [29]. The stability of ground-state solutions is confirmed by observing the long-term evolution of weakly perturbed initial wave profiles. Different initial conditions for the governing GPE are prepared by shifting the position of the minority component in the harmonic trap relative to its minimum at $x=0$.

Now, we estimate the parameter values used for numerical simulations in physical units and compare them with previously reported experimental data. Let us consider a BEC of ⁸⁷Rb atoms prepared in two internal states $|F=1,m_F=-1\rangle$ and $|F=2,m_F=1\rangle$. The s-wave scattering lengths of rubidium atoms in these ground hyperfine states are almost equal $a_1\simeq a_2\simeq 100a_B$, while the inter-component coefficient $a_{12}$ can be independently tuned. The transfer of atoms from one state to the other can be induced by coherent electromagnetic radiation until the desired populations $N_1$, $N_2$ in the corresponding states ${\psi }_1(x,t)$ and ${\psi }_2(x,t)$ are achieved. The mixture condensate is supposed to be held in a toroidal trap with a transverse confinement frequency ${\omega }_{\perp }=2\pi \times 100$ Hz. Then, the units of time and space are given by $\tau =1/{\omega }_{\perp }\simeq 1.6$ ms, $a_{\perp }=\sqrt{\hslash /m{\omega }_{\perp }}\simeq 1\mathsf{\mu }$m, respectively. The length of the integration domain (circumference of the toroidal trap) is equal to $L=12\pi a_{\perp }\simeq 41\mathsf{\mu }$m. The density of the uniformly distributed condensate’s majority component consisting of ${\mathcal{N}}_1={10}^4$ rubidium atoms, required for the calculation of the sound velocity and healing length is $n_1={\mathcal{N}}_1/L{a}_{\perp }^2\simeq 2.1\times {10}^{20}$ m⁻³. Using these data, one can estimate the speed of sound $c_s=(4\pi {\hslash }^2|a_s{|n_1/m^2)}^{1/2}\simeq 2.7$ mm/s and the healing length $\xi =1/\sqrt{8\pi |a_s|n_1}\simeq 0.19\mathsf{\mu }\mathrm{m}$.

It is important to note that the above presented values are associated with uniformly distributed BEC components, taking place for $g_{12}=0$. A localized state of a symbiotic type emerges when there is attraction between the components ($g_{12}=-1.05$), leading to a decrease in the density of the majority component away from the minority one. For a particular case, shown in Figure 1a, the density is reduced by approximately three times. Thus we obtain the corrected values $c_s\simeq 1.6$ mm/s and $\xi \simeq 0.32\mathsf{\mu }\mathrm{m}$.

It is evident from Figure 2b that the transition from a superfluid to a dissipative regime occurs near the critical velocity $\sim 0.1$, which in physical units corresponds to $v_c\simeq 0.07$ mm/s. Therefore the critical velocity leading to suppression of superfluidity in BEC appears to be significantly smaller than the sound velocity, with their ratio being $v_c/c_s\simeq 0.04$. Similar ratios have been reported in experiments with sodium [8] $v_c/c_s\simeq (0.07\pm 0.02)$ and rubidium condensates [7] $v_c/c_s\simeq 0.25$.

The reasons the theoretical prediction overestimates the actual critical velocity have been addressed in many previous works (see e.g., [13]). Among them are the inhomogeneity of the medium, the influence of nonlinearity, the macroscopic size of the moving object, and the dimensionality of the system. All of the quoted factors are inherent to our model. In particular, the size of the localized wave (moving minority component) calculated for the symmetric setting [33] ($N=N_1=N_2,g=g_1=g_2$) $a_0=\sqrt{2\pi }/\left[(g+g_{12})N\right]\simeq 1.2\mathsf{\mu }$m is greater than the healing length $\xi \simeq 0.3\mathsf{\mu }\mathrm{m}$. Thus, the probe object in our model is macroscopic. Currently, uncovering the physical mechanisms behind the onset of the dissipative regime in a BEC superfluid and defining the associated critical velocity is still a challenge.

4. Conclusions

Through numerical simulations, we have studied the superfluid and dissipative regimes in the dynamics of a binary Bose–Einstein condensate where the localized minority component, acting as a single entity (probe object), moves through the majority component, hence, creating excitations. We have designed a configuration in which the minority component is surrounded by the majority one and experiences regular back-and-forth movements under the action of a harmonic trap. Its velocity varies from zero at classical turning points to a maximum value at the lowest point of the harmonic trap. When the velocity of motion exceeds some critical value, strong dissipation occurs, accompanied by a production of quasiparticles in the majority component. After multiple oscillations, the minority component experiences energy losses and slowing down, which eventually cause the return to a non-dissipative mode. It was found that the actual value of the critical velocity leading to the suppression of superfluidity in BEC was smaller than the speed of sound expected from theoretical arguments. The possible reasons for the observed discrepancy were analyzed.