Decay mechanisms of superflow of Bose-Einstein condensates in ring traps
Masaya Kunimi Ippei Danshita
cond-mat.quant-gas























E-mail: masaya.kunimi@yukawa.kyoto-u.ac.jp

Present address : Department of Physics, Kindai University, Higashi-Osaka, Osaka 577-8502, Japan
Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan


















We study the supercurrent decay of a Bose-Einstein condensate in a ring trap combined with a repulsive barrier potential. In recent experiments, Kumar et al. [Phys. Rev. A 95, 021602(R) (2017)] have measured the dependence of the decay rate on the temperature and the barrier strength. However, the origin of the decay observed in the experiment remains unclear. We calculate the rate of supercurrent decay due to thermally activated phase slips (TAPS) by using the Kramers formula based on the Gross-Pitaevskii mean-field theory. The resulting decay rate is astronomically small compared to that measured in the experiment, thus excluding the possibility of TAPS as the decay mechanism. Alternatively, we argue that three-body losses can be relevant to the observed decay and predict that one can observe supercurrent decay via TAPS by decreasing the number of atoms.










Decay mechanisms of superflow of Bose-Einstein condensates in ring traps
    Ippei Danshita
    December 30, 2023
========================================================================






§ INTRODUCTION



Systems of ultracold gases confined in ring-shaped traps have served as an ideal platform for studies of superfluidity <cit.>. Their exquisite controllability has led to experimental observations of various fundamental properties of superfluids, including persistent currents <cit.>, critical velocities <cit.>, and hysteresis loops <cit.>. Such systems attract growing interest also as an atomic analog of superconducting quantum interference device <cit.>, which constitutes a basic element for atomtronic circuits <cit.>. 

A recent experiment performed by the NIST group has raised a puzzling question concerning the superfluidity in a ring-shaped Bose-Einstein condensate (BEC) <cit.>. In this experiment, they prepared a current-carrying BEC with winding number 1 as an initial state and measured the lifetime of this state. To induce the decay of the persistent current, they exposed the BEC to a repulsive potential barrier, which is generated by a blue detuned laser beam. They found that the observed decay rate is significantly dependent on the barrier strength U_0 and the temperature T. One might naively expect that this result could be interpreted as a decay of the persistent current via thermally activated phase slips (TAPS) <cit.>, whose decay rate Γ obeys the Arrhenius law, i.e., Γ∝ e^-E_ B^-/(k_ BT). Here E_ B^- denotes the energy barrier separating the state with winding number 1 from that with winding number zero (see Fig. <ref>).  However, the NIST group presented a rough estimation of E_ B^- to argue that the observed decay rate is inconsistent with the Arrhenius law. Thus the decay mechanism of the persistent currents remains unclear.



The above-mentioned observation of Ref. <cit.> challenges our common understanding of the superfluidity in the sense that TAPS has been established as a universal decay mechanism applicable to various superfluid systems at finite temperatures, such as superfluid ^4He <cit.>, superconductors <cit.>, spin superfluids <cit.>, and one-dimensional Bose gases in optical lattices <cit.>. Resolving this puzzle is important also for engineering and controlling the atomtronic circuits, which require quantitative understanding of the persistent current in the presence of repulsive potential barriers <cit.>.

In this paper, we first reexamine the possibility of TAPS as a decay mechanism by quantitatively computing the nucleation rate of TAPS without any fitting parameters. Our quantitative result shows that the energy barrier in Ref. <cit.> was considerably underestimated. Nevertheless, the computed decay rate completely disagrees with the value measured in Ref. <cit.>, thus confirming that TAPS is irrelevant to the observed decay. Alternatively, we find that three-body losses can induce the decay of superflow in the experimental setup and discuss the possibility of this decay mechanism in the experiment. This indicates that one may observe decay via TAPS by optimizing some experimental parameters in such a way that TAPS is enhanced. We explicitly identify parameter regions in which TAPS dominates over three-body losses and predict the decay rate via TAPS as a reference that can be directly compared with future experiments.

This paper is organized as follows. In Sec. <ref>, we explain the system considered in this paper. In Secs. <ref> and <ref>, we show our results for the energy barrier, lifetime, and real-time dynamics within the three-body loss in the quasi-2D and 3D systems. In Sec. <ref>, we summarize our results. In the Appendix, we explain the details of the numerical method used in this paper.





§ MODEL AND METHODS


We consider a BEC in a combined potential of a ring trap and a repulsive barrier to mimic the situation in the experiment of Ref. <cit.>. The interatomic interaction is sufficiently weak so that the dynamics of the superfluid order parameter ψ(x, y, z, t) is quantitatively described within the mean-field approximation. Specifically, the three-dimensional (3D) Gross-Pitaevskii (GP) equation with a phenomenological three-body loss term <cit.> is given by

    iħ∂/∂ tψ(x, y, z, t)
       =[-ħ^2/2m∇^2+U(x,y,z)+g|ψ(x, y, z, t)|^2]ψ(x, y, z, t)
          -iħ/2L_3|ψ(x, y, z, t)|^4ψ(x, y, z, t),

where m is the mass of the atom, U(x, y, z)≡ (1/2)mω_ r^2(r-R)^2+(1/2)mω_z^2z^2+U_ ext(r) represents the trap potential and the barrier potential, r = (x,y), ω_ r and ω_z are the trap frequency of the radial direction and z direction, respectively, R is the mean radius of the ring trap, g≡ 4πħ^2a_ s/m is the coupling constant, a_ s is the s-wave scattering length, and L_3 is the three-body loss rate. We note that the three-body loss term is set to zero except for the calculations of Fig. <ref> in Sec. <ref>. The barrier potential is created by dithering the Gaussian laser beam. The time-averaged potential U_ ext(r) is given by <cit.>

    U_ ext(r)   =U_0/2{ erf[√(2)/w(x-R+l_ d/2)].
       .- erf[√(2)/w(x-R-l_ d/2)]}e^-2y^2/w^2,

where U_0 is the potential strength, w is the 1/e^2 width of the laser beam, l_ d is the length of the dither, and erf(·) is the error function. 

In Sec. <ref>, we show the quasi-two-dimensional (2D) results. The quasi-2D GP equation is derived by substituting ψ(x, y, z, t)=Ψ(r, t)[1/√(√(π) a_z)]e^-(z/a_z)^2/2 into Eq. (<ref>):

    iħ∂/∂ tΨ(r, t)
       =[-ħ^2/2m∇_ 2D^2+U_ 2D(r)+g_ 2D|Ψ(r, t)|^2]Ψ(r, t)
         -iħ/2L_ 3, 2D|Ψ(r, t)|^4Ψ(r, t),

where ∇^2_ 2D≡∂^2/∂ x^2+∂^2/∂ y^2 is the two-dimensional Laplacian, U_ 2D(r)≡ (1/2)mω_ r^2(r-R)^2+U_ ext(r) is the two-dimensional external potential, g_ 2D≡√(8π)ħ^2a_ s/(m a_z) is the two-dimensional coupling constant <cit.>,  a_z≡√(ħ/(mω_z)) is the harmonic-oscillator length in the z direction, and L_3, 2D≡ L_3/(√(3)π a_z^2) is the quasi-2D three-body loss rate. We also note that the L_ 3,2D term is set to zero except for the calculations of Fig. <ref>. in Sec. <ref>.

To obtain excitation spectra, we linearize the GP equation around a stationary solution Ψ(r, t)=e^-i μ t/ħΦ(r), where μ is the chemical potential. The linearized GP equation, i.e., the Bogoliubov equation, is given by

    H_ B[ u_i(r); v_i(r) ]   =ϵ_i
    [ u_i(r); v_i(r) ]
    ,
    
    H_ B   ≡[              ℒ    g_ 2DΦ(r)^2; -g_ 2DΦ^∗(r)^2           -ℒ^∗ ]
    ,
    ℒ   ≡ -ħ^2/2m∇^2_ 2D+U_ 2D(r)-μ+2g_ 2D|Φ(r)|^2,

where ϵ_i is the excitation energy and u_i(r) and v_i(r) are eigenfunctions of the excited state labeled by i.

In addition to the Bogoliubov equation, we need to diagonalize the following matrix: H_ E≡σ_3H_ B, where σ_3= diag(+1,-1). Let λ_i be an eigenvalue of the matrix H_ E, which is a real value in contrast to the eigenvalue of the Bogoliubov equation ϵ_i∈ℂ.

In the experiment of Ref. <cit.>, they prepare a state with the winding number W=1 as the initial state. Here, we calculate the decay from W=1 state (metastable state) to W=0 state (ground state). According to the Kramers formula <cit.>, the nucleation rate can be written as

    Γ   ≡Γ_–Γ_+,
    Γ_±   ≡|ϵ_ DI|/ħ/2π∏_n^'√(λ_n^ s±/|λ_n^ u|)e^-E_ B^±/(k_ BT),

where Γ_- (Γ_+) is a nucleation rate of the decay process (the acceleration process) (see Fig. <ref>), ϵ_ DI is the frequency of the unstable mode and λ_n^ u and λ_n^ s± are the eigenvalues of the matrix H_ E, E_ B^± is the energy barrier of the acceleration and the decay process, k_ B is the Boltzmann constant, and T is the temperature of the system. In the product of Eq. (<ref>), we omit the zero modes (λ_n=0). We define the nucleation rate (inverse of the lifetime) as the difference between Γ_- and Γ_+. We note that the Kramers formula is valid for E_ B^±≫ k_ BT.


Here, we explain the numerical method used in this work. In the quasi-2D calculations of the GP equation and the Bogoliubov equation, we use the almost same methods as those used in our previous work  <cit.>, i.e., the space discretization is performed by the discrete variable representation method <cit.> (see also the Appendix <ref>) and seeking the unstable stationary solution of the GP equation is based on the pseudoarclength continuation method <cit.> and the Newton method. For the real-time dynamics, we use the pseudospectral method. The typical numerical meshes taken in this work are 129× 129. 

In the 3D calculation, we use the standard discretization method for the calculations of the energy barrier (see Sec. <ref>) and pseudospectral method for the calculations of the real-time dynamics. The typical numerical meshes taken in this work are 129× 129× 65.

It is worth emphasizing that although the GP equation has been extensively used for studying superfluidity of ring-shaped BEC <cit.>, decay rates via TAPS at 2D, to our knowledge, have never been quantitatively calculated because of the difficulty in finding the unstable solutions.





§ RESULTS






 §.§ Quasi-2D results






Here, we show our theoretical results for quasi-2D systems, which are compared to the experimental ones. In the following calculations, we consider ^23Na atom and set the system parameters as the case I in Refs. <cit.>, which corresponds to the lowest temperature and the strongest two dimensionality in the experiment. The specific values are as follows: m=3.82× 10^-26  kg, (ω_ r, ω_z)/(2π)=(258  Hz, 974  Hz), R=22.4 μ m, a_ s=2.75  nm, w=6 μ m, and l_ d=21.8 μ m. We use the particle number N≃ 3.75× 10^5, which comes from the ground-state particle number at the chemical potential μ=μ_ exp≡ h× 2.91  kHz. Here, μ_ exp is the chemical potential measured in the experiment.




  §.§.§ Density and phase profiles



First, we show the density and phase profiles of the order parameter for the metastable and unstable states in Fig. <ref>. We see a solitonic vortex (SV) at the low-density region under the barrier potential in the unstable solution [Figs. <ref> (b) and <ref> (d)]. The NIST group assumed the presence of such a SV in the unstable state when they estimated the energy barrier <cit.>. Our results support the validity of their assumption.




  §.§.§ Energy barrier and lifetime



Next, we show the energy barrier as a function of the strength of the external potential in Fig. <ref>. There we compare our results with those estimated in Ref. <cit.> and see that the latter is considerably underestimated. This happens because the energy barrier was estimated by the energy of the SV [Eq. (2) in Ref. <cit.>] for a uniform system. The present system is nonuniform due to the trap potential and the external potential, which significantly modify the quantitative size of the energy barrier.







We then show the lifetime of the superflow due to the TAPS in Fig. <ref>. We see that the lifetime due to TAPS is astronomically longer than that observed in Ref. <cit.> (1-10  s). The origin of such long lifetime is the large energy barrier compared to the temperature. In the present case, the energy barrier is typically E_ B/h∼ 200  kHz (see Fig. <ref>) and the temperature is given by k_ BT/h≃ 0.6  kHz (T=30  nK). Consequently, we obtain the extremely long lifetime Γ^-1∼ e^E_ B/(k_ BT)∼ e^300. From the above results, we conclude that the TAPS is irrelevant to the decay of superflow observed in Ref. <cit.>. 





  §.§.§ Three-body loss





Instead of the TAPS, we investigate the effects of three-body losses. The reason why we consider three-body losses is that the time scale of the NIST experiment (∼ 5 s) is longer than that of the typical cold atomic gas experiments (≲ 1 s). Hence the three-body loss, which is negligible in the usual cases, may become relevant to the supercurrent decay. Specifically, the superflow state is a metastable state, which may be broken by three-body losses.

In order to investigate the effects of three-body losses in quasi-2D systems, we numerically solve Eq. (<ref>) with L_3,2D≠0. Here, we set L_3=1.1× 10^-30  cm^-6/s, which is the three-body loss rate of the ^23Na atom <cit.>. The initial condition is the ground-state solution of the nondissipative GP equation (L_3=0) for U_0=0 and N≃ 4.40× 10^5, which comes from the ground-state particle number at μ=(10/9)μ_ exp. The reason why we put a factor 10/9 to μ_ exp is that μ_ exp is measured at the end of the experiment <cit.>. According to Ref. <cit.>, the atom loss reduces the chemical potential by 10%. Therefore, the initial particle number is evaluated by (10/9)μ_ exp. 

The time dependence of the potential strength U_0 is determined by mimicking the experimental sequences (see Fig. 1 in Ref. <cit.>). Specifically, we start with U_0=0 at t = 0. The first 1 s is the preparation stage for the W=1 state. The strength is linearly ramped up in 300  ms, kept constant (U_0=1.1μ_ exp) in 400  ms, and linearly ramped down in 300  ms. During this ramp-up process, we rotate the barrier potential at the angular velocity Ω_0≡ħ/(m R^2)≃ 2π× 0.88  Hz. For 1  s≤ t ≤ 3  s, the potential strength is kept zero. Finally, U_0 is linearly ramped up to the desired value in 70  ms and kept constant after that.

We use the fourth-order Runge-Kutta method for the time evolution and the pseudospectral method for space discretization. The numerical meshes are used for 128× 128 and the time step Δ t for Δ t=0.5 μ s.

The time dependence of the angular momentum of the z component per particle is shown in Fig. <ref>, where the expression of the angular momentum is given by 

    L_z(t)≡∫ drΨ^∗(r, t)l_zΨ(r,t),
    
    l_z≡ -iħ(x∂/∂ y-y∂/∂ x).

We can see a sudden jump of the angular momentum immediately after t=3 s. This is due to the ramp-up of the potential strength U_0 during 3 s≤ t≤ 3.07 s, which does not change the winding number. After this sudden decrease, the angular momentum decays at a certain time for U_0≥ 0.7μ_ exp. This decay can be attributed to the fact that the energy barrier E_ B^± diminishes as the total particle number decreases via three-body losses. This result means that three-body losses can induce the decay of superflow within a few tens of seconds at zero temperature.

Here, we comment on the fluctuations of the angular momentum around zero after the phase slip. After the vortex passes through the ring trap, various oscillation modes are excited. Hence the system has weak dissipation (three-body loss), these excitations do not decay. Consequently, the angular momentum strongly fluctuates around L_z=0.




  §.§.§ Discussion





In Sec. <ref>, we found that three-body losses can induce the decay of superflow at zero temperature. The remaining problem is whether this mechanism is the origin of the decay observed in the NIST experiment or not. To investigate this, we compare the decay time scale of the experiment with the theoretical one. Before discussing our results, we briefly explain how to extract the lifetime in the NIST experiment <cit.>. The experimentalists measured the winding number of the BEC as a function of the temperature, the barrier strength, and the time by using the interference technique <cit.>. The ensemble average of the winding number can be interpreted as the angular momentum per particle. They fitted the results to the exponential function exp(-Γ t) [see Fig. 2 (a) in Ref. <cit.>]. The experimental lifetime is given by 1/Γ.


The experimental results show that the lifetime is about 6.7  s for U_0=0.6884μ_ exp. In the theoretical calculations, we obtain that the time for the instability to set in is about 20.4  s for U_0=0.6884μ_ exp. This result shows that the theoretical decay time scale is about three times longer than the experimental one. This means that the experimental results cannot be attributed solely to the three-body losses. In addition to this, our calculations cannot explain the temperature dependence of the lifetime.

An obvious reason for this discrepancy is the lack of finite-temperature effects in our GP calculations. These effects should enhance the decay of superflow. Moreover, other mechanisms, which accelerate the decay of superflow, may exist in the actual experimental setup. For example, one-body loss due to the inelastic collision between the atoms and the background gases, inelastic collision of the light scattering, and inevitable experimental noise, and so on. Hence our results should be interpreted as the upper bound of the decay time scale. The hybrid effects of the three-body loss and the finite temperature effects are a possible scenario of the decay of superflow in the NIST experiment.


To confirm whether or not the hybrid effects appear in the experiment, we need to perform the dynamical simulations including the finite temperature effects and the three-body loss and compare them with the experiments. For example, the truncated-Wigner approximation (TWA) <cit.>, the Zaremba-Nikuni-Griffin (ZNG) formalism <cit.>, and the stochastic projected GP (SPGP) method <cit.> have been used for simulating dynamics of BEC at finite temperatures in previous literature. Among them, the ZNG formalism cannot describe the nucleation process of the vortices due to thermal fluctuations. The TWA at finite temperature and the SPGP method are possible ways to attack this problem. To investigate the decay problem by using such kind of methods will be our future work.

Nevertheless, we can at least predict that, if the decay is indeed due to three-body losses, one can experimentally observe the decay via TAPS by suppressing the energy barrier in the following way: we seek the parameter region where the lifetime due to the TAPS becomes 1-100  ms by tuning the particle number in the trap. In this time scale, we can neglect the effects of three-body losses. Figure <ref> shows the lifetime of the superflow for N≃ 2.7× 10^4 [note that the particle number in the experiment is O(10^5)], which corresponds to the ground-state particle number for μ=h× 0.5  kHz. Here, other parameters such as trap frequencies are fixed. If these results agree with future experiments, we can conclude that the superflow decay via TAPS occurs without being affected by three-body losses.




 §.§ 3D results




In Sec. <ref>, we showed the quasi-2D results. However, the confinement in the z direction is not sufficiently strong for quantitatively justifying the quasi-2D approximation. In fact, the experimental parameters are given by ω_z=2π× 974 Hz, ω_ r=2π× 256 Hz, k_ BT=k_ B× 30 nK≃ h× 600 Hz, and μ=h× 2.91 kHz. From these values, we can obtain the following inequality; μ> ħω_z>k_ BT>ħω_ r, which is inconsistent with a quasi-2D condition ħω_z≫max(k_ BT, μ, ħω_ r). In the experimental setup, the thermal excitations along the z direction are suppressed and the spatial dependence of the order parameter on the z direction is not of the Gaussian form. This means that the conclusions in the previous section should be treated carefully. 

In this section, we show the 3D calculation results in order to clarify the validity of the quasi-2D calculations performed in the previous section. We show the energy barrier and the real-time dynamics with the same parameters as those of the quasi-2D systems. We will see that the quasi-2D and 3D results are qualitatively similar. 




  §.§.§ Density and phase profile





First, we show the density and phase profiles in the 3D systems for the metastable and unstable stationary states in Fig. <ref>. Here, the particle number is N_ 3D≃ 4.31× 10^5, which is evaluated as the particle number of the ground state at the chemical potential μ_ exp by solving the GP equation. We see the SV at the low-density region in the unstable solution [Figs. <ref> (b) and <ref> (f)]. The existence of the SV reflects the z dependence of the density [Fig. <ref> (d)]. Comparing Fig. <ref> (c) with Fig. <ref> (d), we can see that the density around x∼ R of the unstable state is lower than that of the metastable state. This is due to the existence of the SV. These results are consistent with the quasi-2D ones (see Fig. <ref>).




  §.§.§ Energy barrier







Next, we show the energy barrier in Fig. <ref>. The 3D results are shown by the green triangles. We can find that the values of the energy barrier for 3D systems are close to those of the experimental one compared to the 2D ones. In the following, we explain the reason why the energy barriers in 3D are so different from that in 2D (the difference is roughly 200 kHz).

The difference between 2D and 3D systems can be attributed to the discrepancy in the critical strength of the external potential U_ c, where the energy barrier for the W=1 states vanishes. The critical values are given by U_ c^ 2D≃ 0.960μ_ exp and U_ c^ 3D≃ 0.839μ_ exp, which are determined by the time evolution of the GP equation as performed in Sec. <ref> and Sec. <ref>. This result indicates that one should compare the energy barriers in 2D and 3D systems as a function of U_0/U_ c instead of U_0/μ_ exp.  In Fig. <ref>, we show the energy barrier as a function of rescaled external potential strength U_0/U_ c. This result shows that the difference between the 2D and 3D energy barrier at the same U_0/U_ c is about 50 kHz, which is small compared to the difference that is extracted from the same U_0/μ_ exp.

Our results show that the values of the 3D energy barrier are small compared to those of the 2D ones. Nevertheless, our conclusion that the TAPS is irrelevant to the supercurrent decay observed in the NIST experiment still remains. The reasons are as follows. The typical energy barrier for the 3D system is given by E_ B∼ h× 100 kHz. This means that the ratio between the energy barrier and the temperature is given by E_ B/(k_ BT)∼ 100 for T=30 nK. Although we do not calculate the prefactor of the Kramers formula in 3D systems due to the high computational costs, we can estimate the lifetime by using the 2D prefactor. The lifetime for T=30 nK and U_0=0.6884μ_ exp becomes Γ^-1∼ 10^57 s, which is much longer than the experimental timescales. From this 3D calculation, we can confirm that the TAPS is irrelevant in the NIST experiment.




  §.§.§ Three-body loss





Finally, we show the results for the effects of three-body losses in the 3D systems. We solve the 3D dissipative GP equation (<ref>). The numerical procedures are the same as the quasi-2D case. The initial condition is the ground-state solution of the nondissipative GP equation (L_3=0) for U_0=0 and N≃ 5.32× 10^5, which comes from the ground-state particle number at μ=(10/9)μ_ exp. We use the split-step method for the time evolution and the pseudospectral method for space discretization. The numerical meshes and the time step used here are 128× 128× 64 and Δ t=0.5 μ s, respectively.

Figure <ref> shows the time evolution of the z component of the angular momentum for various potential strengths. Comparing these results with those shown in Fig. <ref>, the qualitative behavior is rather similar to the quasi-2D one. We also find that the time for the instability to set in is about 20.3 s for U_0=0.6884μ_ exp. This value is almost the same as  the quasi-2D case (see Sec. <ref>). From these results, we conclude that the effects of the z direction for the three-body loss induced decay are not significant. 





§ SUMMARY



We have investigated the supercurrent decay of BEC in ring traps to reveal the origin of the decay observed in the experiment of Ref. <cit.>. First, we reconsidered the decay via the TAPS based on the Kramers formula in the quasi-2D systems. Our results show that the lifetime via TAPS is astronomically long. This means that the TAPS is not relevant in the experiment of Ref. <cit.>. The same conclusion has also been obtained from the 3D calculations.

Next, we performed the numerical simulations of the GP equation with the three-body loss term in the quasi-2D and 3D systems and found that three-body losses can induce the decay of superflow at zero temperature. Comparing the decay time scale of the theory and experiment, we found the decay time scale is about three times longer than the experimental one. These results indicate that the experimental results cannot be attributed solely to the three-body loss. We proposed a possible scenario that the hybrid effects of the three-body loss and the finite temperature is the origin of the decay observed in the NIST experiment.  We also proposed that one can enhance the TAPS by decreasing the energy barrier in order to observe the decay via the TAPS in future experiments. 

We finally emphasize that the effects of particle loss, to our knowledge, have not been considered in the context of superflow decay before the present work. In this sense, our results may open up a new possibility in the study of superfluidity. Specifically, consideration of such effects is expected to be relevant also to advanced superfluid systems, including exciton polariton BEC <cit.> and ultracold atomic gases with controllable particle losses <cit.>.






We thank G. K. Campbell for providing us the experimental data and useful comments, A. Kumar for providing us with the experimental data, and P. Naidon for critical reading of the manuscript. M.K. thanks S. Goto for his advice regarding numerical calculations. I.D. thanks L. Mathey for useful discussions. M.K. was supported by Grant-in-Aid for JSPS Research Fellow Grant No. JP16J07240. I.D. was supported by KAKENHI grants from JSPS: Grants No. 15H05855, No. 25220711, No. 18K03492, and No. 18H05228, by research grant from CREST, JST, and by Q-LEAP program of MEXT, Japan.

















§ DVR METHOD



In this appendix, we explain the DVR method <cit.>. For numerical calculations, we consider a box region [-L_x'/2, +L_x'/2]×[-L_y'/2, +L_y'/2], where L_x'≡ L_x/a_ r and L_y'≡ L_y/a_ r are the dimensionless system sizes for x and y directions, respectively, and a_ r≡√(ħ/(mω_ r)) is the harmonic-oscillator length of the radial direction. We discretize the space into N_x× N_y lattice points. The condensate wave function is expanded as a series of the DVR functions:

    Ψ'(r', t')   =∑_i=-(N_x-1)/2^(N_x-1)/2 ∑_j=-(N_y-1)/2^(N_y-1)/2Ψ'_i j(t')f_i(x')g_j(y'),

where N_x and N_y are odd natural numbers and the DVR functions f_i(x') and g_j(y') are defined by

    f_i(x')≡1/N_x√(Δ x)sin[π(x'/Δ x -i)]/sin[π(x'/Δ x -i)/N_x],
    
    g_j(y')≡1/N_y√(Δ y)sin[π(y'/Δ y -j)]/sin[π(y'/Δ y -j)/N_y].

Here, the mesh sizes of the x and y directions are given by Δ x≡ L_x'/N_x and Δ y≡ L_y'/N_y. We can show that, on the mesh points, which are x_i'≡Δ x× i and y_j'≡Δ y× j, the DVR functions satisfy f_i(x_k)=δ_i,k/√(Δ x) and g_j(y_k')=δ_j,k/√(Δ y). We can also show that the DVR functions satisfy the orthonormal relations

    ∫^+L_x'/2_-L_x'/2d x'f_i(x')f_j(x')=δ_i j,
       ∫^+L_y'/2_-L_y'/2d y'g_i(y')g_j(y')=δ_i j,

where we used the Gauss quadrature to derive these relations <cit.>. We note that the wave function on the mesh points satisfies Ψ'(x_i', y_j', t')=Ψ'_i j(t')/√(Δ xΔ y).

From the above relations, we can derive the matrix elements of the  first and second derivative terms

    T^(1)_x, i j   ≡∫^+L_x'/2_-L_x'/2d x' f_i(x')d/d x'f_j(x')
       =
    
    π/N_xΔ x(-1)^i-j/sin[π(i-j)/N_x]   (i≠j)
    
    0   (i=j)
    
    ,


    T^(1)_y, i j   ≡∫^+L_y'/2_-L_y'/2d y' g_i(x')d/d y'g_j(y')
       =
    
    π/N_yΔ y(-1)^i-j/sin[π(i-j)/N_y]   (i≠j)
    
    0   (i=j)
    
    ,


    T_x, i j^(2)   ≡ -1/2∫^+L_x'/2_-L_x'/2d x'f_i(x')d^2/d x'^2f_j(x')
       =
    
    (-1)^i-j/Δ x^2π^2/N_x^2cos[π(i-j)/N_x]/sin^2[π(i-j)/N_x]   (i≠j)
    π^2/6Δ x^2(1-1/N_x^2)    (i=j)
    
    ,
    
    T_y, i j^(2)   ≡ -1/2∫^+L_y'/2_-L_y'/2d y'g_i(y')d^2/d y'^2g_j(y')
       =
    
    (-1)^i-j/Δ y^2π^2/N_y^2cos[π(i-j)/N_y]/sin^2[π(i-j)/N_y]   (i≠j)
    π^2/6Δ y^2(1-1/N_y^2)    (i=j)
    
    .


One of the advantages of the DVR method is that the matrix elements of the potential and the interaction terms are diagonal. Although the matrix elements of the kinetic terms have the off-diagonal elements, the matrix is sparse. This is the advantage to the calculation of the Bogoliubov equation. The other advantage is that the DVR method has higher accuracy compared to the standard discretization method (centered difference method). The DVR method can represent up to the polynomials of degree O(N_i) exactly, where N_i (i=x,y) is the number of lattice points for i direction. This is in contrast to the standard discretization method, which can represent up to the quadratic polynomials exactly. This fact comes from the Gauss quadrature in the derivation. 

We note that another type of basis has been used to calculate the GP and the Bogoliubov equation in Ref. <cit.>.




999
Ryu2007
C. Ryu, M. F. Andersen, P. Cladé, V. Natarajan, K. Helmerson, and W. D. Phillips, Phys. Rev. Lett. 99, 260401 (2007).
Ramanathan2011
A. Ramanathan, K. C. Wright, S. R. Muniz, M. Zelan, W. T. Hill, C. J. Lobb, K. Helmerson, W. D. Phillips, and G. K. Campbell, Phys. Rev. Lett. 106, 130401 (2011).
Moulder2012
S. Moulder, S. Beattie, R. P. Smith, N. Tammuz, and Z. Hadzibabic, Phys. Rev. A 86, 013629 (2012).
Beattie2013
S. Beattie, S. Moulder, R. J. Fletcher, and Z. Hadzibabic, Phys. Rev. Lett. 110, 025301 (2013).
Neely2013
T. W. Neely, A. S. Bradley, E. C. Samson, S. J. Rooney, E. M. Wright, K. J. H. Law, R. Carretero-González, P. G. Kevrekidis, M. J. Davis, and B. P. Anderson, Phys. Rev. Lett. 111, 235301 (2013).
Wright2013
K. C. Wright, R. B. Blakestad, C. J. Lobb, W. D. Phillips, and G. K. Campbell, Phys. Rev. Lett. 110, 025302 (2013).
Wright2013_2
K. C. Wright, R. B. Blakestad, C. J. Lobb, W. D. Phillips, and G. K. Campbell, Phys. Rev. A 88, 063633(2013).
Ryu2013
C. Ryu, P. W. Blackburn, A. A. Blinova, and M. G. Boshier, Phys. Rev. Lett. 111, 205301 (2013).
Ryu2014
C. Ryu, K. C. Henderson, and M. G. Boshier, New J. Phys. 16, 013046 (2014).
Eckel2014
S. Eckel, J. G. Lee, F. Jendrzejewski, N. Murray, C. W. Clark, C. J. Lobb, W. D. Phillips, M. Edwards, and G. K. Campbell, Nature (London) 506, 200 (2014).
Jendrzejewski2014
F. Jendrzejewski, S. Eckel, N. Murray, C. Lanier, M. Edwards, C. J. Lobb, and G. K. Campbell, Phys. Rev. Lett. 113, 045305 (2014).
Eckel2014_2
S. Eckel, F. Jendrzejewski, A. Kumar, C. J. Lobb, and G. K. Campbell, Phys. Rev. X 4, 031052 (2014).
Corman2014
L. Corman, L. Chomaz, T. Bienaimé, R. Desbuquois, C. Weitenberg, S. Nascimbène, J. Dalibard, and J. Beugnon, Phys. Rev. Lett. 113, 135302 (2014).
Kumar2017
A. Kumar, S. Eckel, F. Jendrzejewski, and G. K. Campbell, Phys. Rev. A 95, 021602(R) (2017).
Clarke2004
J. Clarke and A. I. Braginski, The SQUID Handbook (Wiley-VCH, Weinheim, 2004).
Seaman2007
B. T. Seaman, M. Krämer, D. Z. Anderson, and M. J. Holland, Phys. Rev. A 75, 023615 (2007).
Pepino2009
R. A. Pepino, J. Cooper, D. Z. Anderson, and M. J. Holland, Phys. Rev. Lett. 103, 140405 (2009).
Langer_Fisher1967
J. S. Langer and M. E. Fisher, Phys. Rev. Lett. 19, 560 (1967).
Langer1967
J. S. Langer and V. Ambegaokar, Phys. Rev. 164, 498 (1967).
McCumber1970
D. E. McCumber and B. I. Halperin, Phys. Rev. B 1, 1054 (1970).
Langer1970
J. S. Langer and J. D. Reppy, in Progress in Low Temperature Physics, edited by C. J. Gorter (North-Holland, Amsterdam, 1970), Vol. 6, Chap. 1.
Reppy1992
J. D. Reppy, J. Low Temp. Phys. 87, 205 (1992).
Halperin2010
B. I. Halperin, G. Refael, and E. Demler, Int. J. Mod. Phys. B 24, 4039 (2010).
Kim2016
S. K. Kim, S. Takei, and Y. Tserkovnyak, Phys. Rev. B 93, 020402(R) (2016).
Tanzi2016
L. Tanzi, S. S. Abbate, F. Cataldini, L. Gori, E. Lucioni, M. Inguscio, G. Modugno, and C. D'Errico, Sci. Rep. 6, 25965 (2016).
Kunimi2017
M. Kunimi and I. Danshita, Phys. Rev. A 95, 033637 (2017).
Kagan1998
Y. Kagan, A. E. Muryshev, and G. V. Shlyapnikov, Phys. Rev. Lett. 81, 933 (1998).
Kumar2017_pb
A. Kumar, private communications.
Petrov2000_2
D. S. Petrov, M. Holzmann, and G. V. Shlyapnikov, Phys. Rev. Lett. 84, 2551 (2000).
Langer1969
J. S. Langer, Ann. Phys. 54, 258 (1969).
Hanggi1990
P. Hänggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 62, 251 (1990).
Baye1986
D. Baye and P.-H. Heenen, J. Phys. A: Math. Gen. 19, 2041 (1986).
Keller1987
H. B. Keller, Lectures on Numerical Methods in Bifurcation Problems (Springer-Verlag, Berlin/Heidelberg/New York/Tokyo, 1987).
Kunimi2015
M. Kunimi and Y. Kato, Phys. Rev. A 91, 053608 (2015).
Brand2001
J. Brand and W. P. Reinhardt, J. Phys. B: At. Mol. Opt. Phys. 34, L113 (2001).
Piazza2009
F. Piazza, L. A. Collins, and A. Smerzi, Phys. Rev. A 80, 021601(R) (2009).
Mathey2014
A. C. Mathey, C. W. Clark, and L. Mathey, Phys. Rev. A 90, 023604 (2014).
Mateo2015
A. Muñoz Mateo, A. Gallemí, M. Guilleumas, and R. Mayol, Phys. Rev. A 91, 063625 (2015).
Snizhko2016
K. Snizhko, K. Isaieva, Y. Kuriatnikov, Y. Bidasyuk, S. Vilchinskii, and A. Yakimenko, Phys. Rev. A 94, 063642 (2016).
Supplemental_Kumar
Supplemental materials of Ref. <cit.>.
Stamper-Kurn1998
D. M. Stamper-Kurn, M. R. Andrews, A. P. Chikkatur, S. Inouye, H.-J. Miesner, J. Stenger, and W. Ketterle, Phys. Rev. Lett. 80, 2027 (1998).
Campbell_pb
G. K. Campbell, private communication.


Blakie2008
P. B. Blakie, A. S. Bradley, M. J. Davis, R. J. Ballagh, and C. W. Gardiner, Adv. Phys.57, 363 (2008).
Polkovnikov2010
A. Polkovnikov, Ann. Phys. 325, 1790 (2010).
Griffin2009
A. Griffin, T. Nikuni, and E. Zaremba, Bose-Condensed Gases at Finite Temperatures, (Cambridge University Press, Cambridge, 2009).
Rooney2013
S. J. Rooney, T. W. Neely, B. P. Anderson, and A. S. Bradley, Phys. Rev. A 88, 063620 (2013).
Deng2010
H. Deng, H. Haug, and Y. Yamamoto, Rev. Mod. Phys. 82, 1489 (2010).
Labouvie2016
R. Labouvie, B. Santra, S. Heun, and H. Ott, Phys. Rev. Lett. 116, 235302 (2016).
Tomita2017
T. Tomita, S. Nakajima, I. Danshita, Y. Takasu, and Y. Takahashi, Sci. Adv. 3, e1701513 (2017).
Numerical_recipes
W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes, 3rd ed. (Cambridge University Press, Cambridge, 2007).
Schneider1999
B. I. Schneider and D. L. Feder, Phys. Rev. A 59, 2232 (1999).





