Characterization of the energy level-structure of a trapped dipolar Bose gas via mean-field parametric resonances
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] J a n Characterization of the energy level-structure of a trapped dipolar Bose gas viamean-field parametric resonances
Asaad R. Sakhel and Roger R. Sakhel Department of Physics, Faculty of Science, Balqa Applied University, Salt 19117, Jordan Department of Physics, Faculty of Science, Isra University, Amman 11622, Jordan (Dated: January 26, 2021)We report parametric resonances (PRs) in the mean-field dynamics of a one-dimensional dipolarBose-Einstein condensate (DBEC) in widely varying trapping geometries. The chief goal is tocharacterize the energy levels of this system by analytical methods and the significance of thisstudy arises from the commonly known fact that in the presence of interactions the energy levelsof a trapped BEC are hard to calculate analytically. The latter characterization is achieved bya matching of the PR energies to energy levels of the confining trap using perturbative methods.Further, this work reveals the role of the interplay between dipole-dipole interactions (DDI) andtrapping geometry in defining the energies and amplitudes of the PRs. The PRs are induced bya negative Gaussian potential whose depth oscillates with time. Moreover, the DDI play a role inthis induction. The dynamics of this system is modeled by the time-dependent Gross- Pitaevskiiequation (TDGPE) that is numerically solved by the Crank-Nicolson method. The PRs are discussedbasing on analytical methods: first, it is shown that it is possible to reproduce PRs by the Lagrangianvariational method that are similar to the ones obtained from TDGPE. Second, the energies at whichthe PRs arise are closely matched with the energy levels of the corresponding trap calculated bytime-independent perturbation theory. Third, the most probable transitions between the trap energylevels yielding PRs are determined by time-dependent perturbation theory. The most significantresult of this work is that we have been able to characterize the above mentioned energy levels of aDBEC in a complex trapping potential.
Keywords: Dipolar Bose-Einstein condensates, variational methods, Crank-Nicolson method, parametricresonances
I. INTRODUCTION
The phenomenon of parametric resonances (PRs) isubiquitious in nature and is a widely examined funda-mental physical property. Today, PRs are one of the out-standing features observed in Bose-Einstein condensates(BECs) [1–18] where resonances are usually the resultof modulating one of the system parameters such as thescattering length [11, 19–22]. They are also generated byexternal means such as a time-dependent trapping geom-etry [11], laser stirring [2], and laser-intensity modulation[1]. Moreover, PRs have been shown to occur in classicalsystems [23, 24] as well as quantum devices such as bire-fringent optical fibers [25], magnetometers [26], super-conducting wave guides [27], and quantum dots [28, 29].Their importance has also been demonstrated in transat-lantic telecommunication fiber optics [30] in connectionto a modulational instability.In this work, we explore PRs in a dipolar BEC(DBEC) within a setting of large-sized traps that wouldallow an examination of the effects of long-range interac-tions. It is known that dipole-dipole interactions (DDI)are capable of shifting the frequency of PRs [1, 31, 32].This shift is sensitive to the trapping geometry [31]. As aresult, one concludes that the interplay between trappinggeometry and DDI determines the energies at which PRsoccur. These energies are equivalent to the energies ofthe trap levels; the trap in which the DBEC is confined.Within this context then, the chief goal is to characterizethe energy-level structure (ELS) of a DBEC in a trap of large size by matching the PR energies to energy levelsof this trap using perturbative methods. Herein, the im-portance of PRs is revealed in characterizing the shape ofthe BEC trap and its energy levels [1]. The significanceof the present investigation arises from the fact that itfacilitates an evaluation of the above-named structure incomplex potentials. This is because it is known that,in the presence of interactions, it is rather hard to an-alytically calculate the energy levels of a trapped BEC.Moreover, via the present research, we hope to motivatefuture experiments that could determine the ELS of morecomplicated traps that may be engineered in the future.That said, the usefulness of the present examination isrevealed from the latter statements. Another goal is toreveal the interplay between trapping geometry and DDI(cf. [20, 33–35]) in defining the frequencies at which PRsoccur. Effects of trapping geometry have already beenstudied earlier such as the influence of the trap aspectratio on the oscillation frequencies [36] and stability of aDBEC [37, 38], except that it hasn’t been related to thestructure of the trap energy levels.At present, we consider a one-dimensional (1D)trapped DBEC that is driven by a negative Gaussianpotential (NGP) whose depth is periodically modulatedwith time. The latter system is simulated numericallyin different traps using the mean-field time-dependentGross-Pitaveskii equation (TDGPE). The 1D DBEC isscanned over a long range of DDI strengths in an at-tempt to detect PRs and to study the above mentionedinterplay. In passing, it should be noted that the 1D Bosegas is a general, important, and well known system thatcan refer to many physical systems, such as optical fibers[25, 39–41]. It is studied to reveal the physics in 1Dwhich is strikingly different than in higher dimensions.The addition of an NGP is for the purpose of causingthe particles to condense into lower energy levels, so tospeak to “catch particles”, and then throw (excite) themto higher energy levels by the oscillating NGP depth.The oscillatory NGP has been earlier shown to work likea modulated contact interaction [7] and could thereforebe viewed analogous to it. It should also be noted thatan NGP has been used to induce a BEC and to studyits growth dynamics [42]. With the above NGP, a laser-light source is modelled whose intensity oscillates withtime [1, 43] and provides a softer stirring of the BEC inenergy space than the stirring by time-dependent spatialmodulations. This method has rarely been used or men-tioned in the BEC literature, and here it is demonstratedthat it leads to significant excitations.The PRs are measured by a quantity that resemblesthe time average of the square of the kinetic energy called“signal energy” [44]. The signal energy is a term bor-rowed from engineering topics for the processing of anoscillating electrical signal. The signal energy has alsobeen found very effective in revealing PRs in one of ourearlier publications [2].One motivation for the present examination arisesfrom the work of Balik et al. [1] who applied a CO -laser generated optical-dipole trap to confine a sampleof Rb atoms. By a modulation of the laser intensity,the authors were able to excite PRs whose frequencieswere found to shift by a change of the laser’s modula-tion depth. In this regard, we make an analogy by ex-plaining the frequency shift as resulting from a changein the depth of the effective time-averaged potential ofthe trapped DBEC. The latter depth is controlled by theDDI, and by scanning the DBEC over a broad range ofthe latter it has been possible to observe PRs at valuesof DDI corresponding to the resonance frequencies. Thisproject has been heavily computational as it required alarge number of runs to locate the DDI regimes wherePRs occur.The most significant result of this work is that we havebeen able to characterize the energy levels of a DBEC ina complex trapping potential, i.e., a trap to which anNGP is added. Other key results are detailed as follows.(i) By a variation of the DDI strength we detect mean-field PRs in the dynamics of a DBEC induced by anNGP whose depth is periodically modulated withtime. The PRs obtained are an inherent featuredeeply embedded in the physics of the DBEC asthey could also be generated by other tools such asthe Lagrangian variational method (LVM) [7, 45–51].(ii) The positions of the PRs are determined by theELS of the trapping geometry. It is shown that theenergies at which PRs occur can be characterized by second-order perturbation theory thus allowingus to determine the above mentioned structure.(iii) It is shown that the depth of the time-averaged ef-fective potential is reduced by an increase in therepulsive DDI strength as a result of which the PRfrequencies are shifted to larger values in the caseof a harmonic oscillator (HO)+NGP. However, inthe case of a BOX+NGP, the PR frequencies dropinstead with rising DDI strength. Thus the inter-play between trapping geometry and DDI controlsthe values of the DDI strength at which PRs arisein the signal energy [44]. Moreover, the amplitudeof the PRs decline with increasing DDI because thedepth of the time-averaged mean-field effective po-tential is reduced with it.(iv) It is found that the occurrence of a PR requires twoconditions: (1) its energy should closely match oneof the trap levels and (2) the transition probabilityto this level should be quite large.The organization of this paper is as follows. In Sec.II the method is presented. In Sec. III the results arepresented and discussed and in Sec. IV the origins of thePRs are explained. Finally in Sec. V we summarize andconclude. In Appendix A technical details for the simu-lations are outlined and in Appendix B a measurementunit is derived.
II. METHOD
In this section, only the rudiments of the method areoutlined. The reader is therefore referred to Ref.[52] formore details on the method and the numerics involved.
A. Basic Units
In the present work, lengths and energies are in unitsof the trap a ho = p ¯ h/ ( m ¯ ω ) and ¯ h ¯ ω , respectively, where¯ ω = ( ω x ω y ω z ) / is the geometric average of the trap-ping frequencies along the coordinate axes, and m is themass of the atom. It should be however noted, that themodes in the y - and z -directions are frozen since thepresent system is described by the TDGPE that is re-duced from 3D to 1D by integrating out the transversedirections. Time is unitless where t = ¯ ωτ , τ being thetime in seconds. B. Systems
The systems considered are 1D strongly repulsiveDBECs that are confined in a few different trapping ge-ometries. The confining potentials are power-law trapsthat have the form V tr ( x ) = 12 (cid:12)(cid:12)(cid:12)(cid:12) xL x (cid:12)(cid:12)(cid:12)(cid:12) p , (1)where p is the trapping exponent, and L x a length scalethat shapes V tr ( x ) so that there is some flatness in theneighborhood of x = 0 when p is much larger than 2. V tr ( x ) is in units of ¯ h ¯ ω whereas x and L x are in a ho .It should be emphasized that power-law traps have beenrealized in the United Kingdom using spatial light mod-ulation [53].The DBECs are excited by an NGP whose depth os-cillates with time. Experimentally, the NGP is generatedby the application of a focusing red-laser beam [42, 54–61]. The NGP has been considered in theoretical investi-gations [62–67] as well. Its importance has been demon-strated in an architecture of quantum computation [66],the transportation of BECs [59], and the extraction ofatoms from a BEC [63, 67]. Solitons [68] in an NGP andproperties of a BEC in a harmonic trap plus an eccen-tric NGP [65] have also been examined. The red-detunedlaser beam interacts with the BEC in such a way so asto introduce an NGP into it by phase-imprinting [69].Adding to this the requirement of an oscillating inten-sity, the NGP is modelled by V DT ( x, t ) = [ A + δA cos(Ω t )] exp( − βx ) , (2)where Ω = 2 πf is the driving frequency, A the princi-pal depth, δA the modulation amplitude, and as usual1 / √ β a measure of the NGP width. A and δA are inunits of ¯ h ¯ ω , β in a − ho , and Ω is in units of ¯ ω . As alreadynoted in the introduction, the analog of the present exci-tation method is the periodic modulation of a scatteringlength such as a ( t ) = a bg + δa sin( ωt ) [19–21], where δa is the modulation amplitude, and a bg an unperturbedbackground. The analog of a bg is the A and of δa the δA .In one of our earlier publications [7], it has been verifiedthat the effects on the BEC arising from the NGP withoscillating depth [Eq.(2)] are indeed similar to those from a ( t ) above if a bg < δV DT ( x, t ) = δA cos(Ω t ) e − βx , generates an oscillatingforce along the length of the DBEC that is given by thepotential gradient δF x = − ∂δV DT ( x, t ) ∂x = 2 βxδA cos(Ω t ) e − βx . (3)This force transfers a momentum from the oscillatingNGP to the BEC that is symmetric about x = 0 andreads ∆ p x = Z t δF x dt = 2 βxe − βx δA sin(Ω t )Ω . (4)The latter is maximal at x = ± / √ β and minimal at x = 0 and the edges of the trap. C. Dipolar Interactions and their Control
A note on the manipulation of the DDI is in orderhere. Consider the DDI potential given by [52, 70] U dd ( R ) = C dd (1 − θ ) | R | , (5)where R = r − r ′ is the relative position vector of twodipoles at r and r ′ , θ the angle between R and the ori-entation of the dipoles, C dd = E α / (¯ h ¯ ωǫ ) [70, 71] (intrap units) with E being an external electric field, ǫ thepermittivity of free space, and α the static polarizabil-ity. Thus, the DDI can then be induced and tuned byan external field E [72] to various orders of magnitude[73] and in polar molecules the dipole moment can be setup to 10 times larger than in atomic systems. The DDIalso occur naturally [74–78] if the atoms possess a mag-netic dipole moment ¯ µ in which case C dd = µ ¯ µ / (4 π ¯ h ¯ ω )[52, 70, 76] (in trap units) with µ the permittivity of freespace. The control of DDI has also been further explainedin the article by Lahaye et al. [70]. Moreover, by usingthe linear Stark effect, it is possible to excite Rydberg-dressed atoms to very high principle quantum numbersto achieve large dipole moments like p ∼ D , where D is the ratio between the dipolar and s -wave scatteringlength [79]. In this regard, we justify the use of largevalues of the DDI parameter G dd in the present work.Before we continue, it should be emphasized that theDDI in the present work is repulsive. D. Mean-Field Gross-Pitaevskii Equation
The trapped DBEC is described by the time-dependent Gross-Pitaevskii equation (TDGPE). It is re-duced from the 3D to the 1D form by integrating out thecontributions in the transverse direction (see e.g. Refs.[52, 80]). We thus consider a cigar-shaped DBEC thatis elongated along the x –axis with strong radial confine-ment in the transverse direction. The dynamics in thetransverse direction is frozen in the radial ground state φ ( ρ ) = 1 d ρ √ π e − ρ / (2 d ρ ) , (6)where d ρ is the width of the Gaussian. The reduced 1DTDGPE reads then as in Ref.[52] i ∂ψ ( x, t ) ∂t = (cid:20) − ∂ ∂x + V tr ( x ) + V DT ( x, t ) + G D | ψ ( x, t ) | +4 π G dd Z + ∞−∞ dk x π e − ik x x | e ψ ( k x , t ) | j D ( τ x ) (cid:21) ψ ( x, t ) , (7)where e ψ ( k x , t ) is the Fourier transform of ψ ( x, t ). Thelatter is the longitudinal wave function which is normal-ized to one vis-´ a -vis R + ∞−∞ | ψ ( x, t ) | dx = 1 with ψ ( x, t ) inunits of a − / ho . j D ( τ x ) is a function given by j D ( τ x ) = √ πd ρ Z + ∞−∞ dτ y e − τ y h D ( τ ) , (8)with τ y = d ρ k y / √ τ = q τ x + τ y , and h D ( τ ) = 1 √ πd ρ h − √ πe τ τ { − erf ( τ ) } i . (9) G D and G dd are respectively the 1D s-wave and dipo-lar interaction parameters (see Sec. II F below). Thefourth term on the right-hand-side of (7) introduces theusual mean-field s-wave interaction nonlinearity. The lastterm introduces the mean-field dipolar nonlinearity de-rived from the long-range DDI. It is a special integralthat is designed to eliminate the singularity in the DDIpotential.Equation (7) is solved numerically using the famoussplit-step Crank-Nicolson (CN) method [52, 80] in realtime. It is first solved in imaginary time to initializethe BEC in the trapping geometry (1) superimposed onwhich is the NGP (2) without the modulated part, thatis V (0) DT ( x ) = A exp( − βx ) . (10)In the second step, the DBEC is driven by the NGP (2)in real time to examine its ensuing dynamics. For fur-ther technical details, the reader is referred to AppendixA. The codes used for solving the TDGPE have beenwritten by the group of Antun Balaz in Belgrade andare fully explained in Ref.[52] for recent versions treatingDBECs and also Ref. [80] for earlier versions on ordinaryBECs. Numerous other TDGPE codes are available bythis group [81–87] and have been used extensively. E. Gross-Pitaevskii Energies
The Gross-Pitaevskii (GP) energies are as usual eval-uated via E GP ( t ) = Z + ∞−∞ dx " (cid:12)(cid:12)(cid:12)(cid:12) ∂ψ ( x, t ) ∂x (cid:12)(cid:12)(cid:12)(cid:12) +[ V tr ( x ) + V DT ( x, t )] | ψ ( x, t ) | + 12 G D | ψ ( x, t ) | +2 π G dd Z + ∞−∞ dk x π e − ik x x | e ψ ( x, t ) | j D ( τ x ) | ψ ( x, t ) | (cid:21) . (11) The time average is then computed using h E GP i t = 1 T Z T E GP ( t ) dt. (12) F. Parameters
The 3D s -wave and dipolar interaction parameters aredefined as G = 4 πN a and G dd = 3 N a dd , where a and a dd are the s-wave and dipolar scattering lengths, and N thenumber of particles. The interaction parameters acting in1D here, G D and G dd , are obtained after the reductionof the 3D TDGPE to the 1D form (7). As a result of thelatter, the G and G dd are divided by a factor 2 πd ρ so that G D = G πd ρ and G dd = G dd πd ρ , (13)where d ρ is the width of the integrated-out wave functionin the transverse direction. The G D and G dd are inputdirectly into the code without explicit evaluation via N , a , a dd , and d ρ . These parameters define the strength ofthe s-wave and dipolar-interaction nonlinearities of theTDGPE. The systems are simulated from x = − L to x = + L with L = 30 for a harmonic oscillator (HO)trap and quartic trap (QT), and L = 51 . A = − β = 4, and Ω = 2 πf with f = 10. Thevalues of δA applied are 5, 10, and 20. The DBECs arescanned along a range of G dd ranging from 0 to 400 insteps of 2 and for each G dd a run was performed. Valuesof G D = 50 , ,
150 were used for each set of the runs. G D , G dd L , a , a dd , and d ρ are all in units of a ho . G. Signal Energy
For a time-dependent physical observable f ( t ), thesignal energy D is formally defined by the integral D = Z + ∞−∞ | f ( t ) | dt, (14)where f ( t ) is taken to be the mean-field kinetic energy E kin ( t ) = 12 Z + ∞−∞ (cid:12)(cid:12)(cid:12)(cid:12) ∂ψ ( x, t ) ∂x (cid:12)(cid:12)(cid:12)(cid:12) dx, (15) ψ ( x, t ) being the time-dependent wavefunction describ-ing the DBEC. We refer to the signal energy by D kin because it is derived from E kin . E kin is in units of ¯ h ¯ ω and D kin in (¯ h ¯ ω ) ¯ ω . However, since we cannot numeri-cally integrate to infinity, neither timely nor spatially, welimit the integral in Eq.(14) from t = 0 to some time T that is sufficiently long to reveal enough of the dynamicalproperties, and in Eq.(15) to the length of the simulationgrid from x = − L to L . The usefulness of the signalenergy in displaying important properties about excita-tions in a driven BEC has already been demonstrated inour recent work [2]. There, it has been argued that thesignal energy is tantamount to the time-average of thesquared amplitude of an oscillating signal describing adynamic variable. It can be likened to the root-mean-squared value of an alternating voltage or current. Since E kin ( t ) is found here again to be oscillating with time,it is quite convenient to apply the latter concept to itsmeasurement. The reason for using the kinetic energystems additonally from the fact that it is an importantproperty of the condensate [88]. Other quantities, suchas the potential energy, zero-point energy, and the radialsize could also have been used here because they revealthe same PRs with the same properties as obtained for E kin . Further, E kin has been used in a number of previ-ous examinations [7, 46, 89–92] that further demonstrateits importance.As far as the experimental measurement of D kin isconcerned, this can be performed as follows. Whilethe BEC is excited by the focusing red-laser beam ofan oscillating intensity, an in-situ recording of densityprofiles n ( x, t ) is performed as a function of time bya CCD camera using a nondestructive method as thatdemonstrated in the experiment of Onofrio et al. [93].The latter performed a repeated in-situ non-destructivephase-contrast imaging that has been used by Andrews et al. [94] and Ketterle et al. [95]. After the n ( x, t ) arerecorded, they are Fourier-transformed to obtain the cor-responding momentum-density distributions n ( p , t ) thatcan be used to compute E kin ( t ) via E kin ( t ) = 12 m Z dp p n ( p , t ) . (16)Having obtained E kin , D kin is then evaluated easily. III. RESULTSA. Resonances in a Box Potential
Figure 1 displays D kin / D (0) kin vs G dd for a few differentvalues of G D and δA in a BOX+NGP trap. D (0) kin is themaximum signal energy for G D = 50 at δA = 20 in therange of G dd considered and is used for normalizationof the signals in all frames. Several PRs are discoveredwhose amplitudes decline with increasing values of G dd .This is in line with Refs.[75, 96], where it has been foundthat the DDI reduce the amplitude of DBEC dynam-ics. An increase in the DDI causes therefore a weakerresponse to the oscillating NGP as it turns out that therepulsive DDI reduce the depth of the effective mean-fieldpotential (see Sec. III D below) and with it the occupancyof the NGP so that lesser particles are excited. . . mag. × δA = 5 . mag. × G dd D k i n / D ( ) k i n ( A ) [ × f = 10, β = 4 p = 100, L x = 51 . G D = 50, A = −
30 4003503002502001501005001 . G dd D k i n ( × ) ( B ) G D = 100 4003503002502001501005001 . S i g n a l E n e r g y D k i n / D ( ) k i n G dd D k i n ( × ) ( C ) G D = 150 4003503002502001501005001 . . G dd FIG. 1: Signal energy D kin of a driven one-dimensional DBECin a box as a function of the dipolar interaction parameter G dd [Eq.(13)]. D kin is normalized by D (0) kin , the maximumof D kin at δA = 20 in the G dd -range considered. G D is the s − wave interaction parameter [Eq.(13)]. The box is generatedfrom Eq.(1) using L x = 51 . p = 100. The DBEC isexcited by an NGP whose principal depth A is periodicallymodulated by an amplitude δA via Eq.(2) with A = − f = 10, and β = 4. Various values are used for δA and G D .Frame ( A ) and thick-solid line: D kin at G D = 50 with δA =5. The signal is magnified (mag.) by 60 times ( × ) to enablecomparison with other cases; thin-dashed line: δA = 10 . × B, C ): as in ( A ); but for G D = 100 and 150, respectively. A and δA are in units of ¯ h ¯ ω ; L x , G D , G dd in a ho , whereas β in a − ho . D kin and D (0) kin are in (¯ h ¯ ω ) ¯ ω . Qualitatively, the same features are observed for allvalues of G D in Fig. 1. Notably, an increase in G D shiftsthe whole spectrum backwards (to the left) keeping thesame distance ∆ G dd between each pair of peaks, demon-strating that ∆ G dd is not influenced by G D . For ex-ample, the third peak from the left in frame ( A ) at G dd ∼
200 is shifted backwards by an amount of 50 in frame( B ) from its position in ( A ), and by an amount 100 in × δA = 5 (mag × G dd D k i n / D ( ) k i n ( A ) f = 10, β = 4 p = 2, L x = 1 G D = 50, A = − . G dd D k i n / D ( ) k i n ( B ) G D = 100 4003503002502001501005001 . S i g n a l E n e r g y D k i n / D ( ) k i n G dd D k i n / D ( ) k i n ( C ) G D = 150 4003503002502001501005000 . G dd FIG. 2: As in Fig. 1; but for an HO trap with p = 2, L x = 1,magnifications of ×
30 for δA = 5 and × δA = 10, and D (0) kin is the maximum of D kin at δA = 20 in frame ( A ). A and δA are in units of ¯ h ¯ ω ; L x , G D , G dd in a ho , whereas β in a − ho . D kin and D (0) kin are in (¯ h ¯ ω ) ¯ ω . frame ( C ). That is, an increase of G D by 50 causes all theresonance positions to shift backwards by the same valueof 50, (similarly for an increase by 100) being a ratherunprecedented and remarkable feature. Thus there seemsto be no effect for the interplay between s -wave interac-tions and DDI on the principle features of the spectrum of D kin in a box potential. This demonstrates that s -waveinteractions and DDI work similarly in determining thepositions of PRs. The above results reveal informationabout the energy-level structure of the DBEC and shapeof the trapping potential. The equidistance of peaks mir-rors the confinement homogeneity, i.e. the flatness of thebox. This can be particularly concluded from a compar-ison with the HO potential in the next section.The reason for the PRs and the shift of their po-sitions with G D are discussed in Sec. IV basing onwell-known theoretical methods: the Wenzel-Kramer-Brillouin (WKB) approximation, time-independent andtime-dependent perturbation theory, and LVM. Most im- portantly, it is shown that the PRs arise whenever thetime-averaged energy of the DBEC closely matches oneof the energy levels of the trap+NGP. B. Resonances in a Harmonic Trap
In Fig. 2, PRs are also encountered in an HO+NGPtrap, except that the separation ∆ G dd is not uniform asin the box but rather increases with G dd . Henceforth,this mirrors confinement inhomogeneity. Further, thisreveals the role of the interplay between trapping geome-try and DDI in determining the PR energies. A uniformtrap leads to equidistant whereas a nonuniform one tonon-equidistant PRs along the G dd axis. What is pecu-liar though, is that the shift of the whole spectrum as aresult of changing G D is also observed here in the sameuniform manner as reported in Fig. 1. Increasing G D by50 causes the PR peaks to shift backwards by 50 alongthe G dd axis. C. Resonances in a Quartic Trap
In Fig. 3 for a QT+NGP, there is only one well-defined PR at δA = 5 in the same range of G dd consid-ered. At the larger δA , a disordered excitation patternarises. Compared to Figs. 1 and 2, this is rather surpris-ing since a pattern similar to that for the HO+NGP traphad been anticipated. Thus for a QT at stronger drivingquite a larger number of modes is excited than in the HOtrap and box. As such, there exist trapping geometriesthat under certain conditions do not support ordered PRpatterns. Once again it can be seen that the geometry ofthe trap strongly influences the PR phenomenon and itspattern. Therefore, it can be used to control the PRs. Asin the previous figures, an increase of G D by 50 causesthe whole PR spectrum to shift by 50 along the G dd axis. D. Effective Potential
To this end, it is useful to state the reason for thedecline of the amplitude of a PR with G dd as observedin Figs. 1 and 2. This decline is attributed to a change inthe depth of the time-averaged effective potential givenby h V eff ( x ) i t = * V DT ( x, t ) + G D | ψ ( x, t ) | +4 π G dd Z + ∞−∞ dk x π e − ik x x j D ( τ x ) | e ψ ( k x , t ) | + + V tr ( x ) , (17)with × δA = 5 (mag × G dd D k i n / D ( ) k i n ( A ) p = 4, L x = 1 δA = +10, β = 4 G D = 50, A = −
30 4003503002502001501005001 . . . . . G dd D k i n / D ( ) k i n ( B ) G D = 100 4003503002502001501005001 . . . . . S i g n a l E n e r g y D k i n / D ( ) k i n G dd D k i n / D ( ) k i n ( C ) G D = 150 4003503002502001501005001 . . . . . G dd FIG. 3: As in Fig. 1; but for a QT with p = 4, L x = 1, magni-fication of × δA = 5 and 10, and D (0) kin is a normalizationfactor. A and δA are in units of ¯ h ¯ ω ; L x , G D , G dd in a ho ,whereas β in a − ho . D kin and D (0) kin are in (¯ h ¯ ω ) ¯ ω h· · · i = 1 T Z T ( · · · ) dt, (18) T being the total simulation time. Fig. 4( A ) shows as anexample h V eff ( x ) i t for the HO+NGP of Fig. 2( A ). Fig.4( B ) displays the minimum of h V eff ( x ) i t at x = 0 forprevious systems in Figs. 1 and 2 as a function of the PRvalues of G dd . h V eff ( x ) i t becomes shallower as G dd risescausing lesser particles to occupy the NGP. Hence, lesserparticles contribute to the strength of the excitations asthey are thrown out of the NGP by its oscillating depththereby causing the drop in the PR amplitudes D kin . No-tably, h V eff ( x ) i t closely follows the spatial shape of theNGP and it can be therefore argued that the influence ofthe NGP on the DBEC is reduced as h V eff ( x ) i t becomesshallower. Further the depth h V eff ( x ) i t being controlledby G dd plays also a decisive role in determining the PRenergies. . . . . G dd = 29 . x h V e ff ( x ) i t δA = 10(HO+NGP)( A ) − − − − − − − − − − G D = 100, BOX+NGP G dd V e ff , m i n δA = 10( B ) − − − − − FIG. 4: Frame ( A ): Time-averaged effective potential h V eff ( x ) i t [Eq.(17)] for the HO+NGP in Fig. 2( A ) at δA = 10and the indicated resonant values of G dd . Solid line: 29.9297;dotted line: 67.8894; triple-dotted line: 120.449; dashed-dotted line: 191.988; dashed-double-dotted line: 292.727.Frame ( B ): Minima V eff,min = h V eff (0) i t at the PR valuesof G dd for some of the systems presented in Figs. 1 and 2at δA = 10: Solid circles: BOX+NGP at G D = 100; opencircles: same at 50; open squares: HO+NGP at G D = 50;solid squares: same at 100. h V eff ( x ) i t , V eff,min , A , and δA arein units of ¯ h ¯ ω , whereas L x , x , G D , and G dd in units of a ho ,and β in a − ho . IV. ORIGINS OF THE PARAMETRICRESONANCES
The reasons for the appearance of the PRs at cer-tain values of G dd are deeply encrypted in the numericalsolutions of the TDGPE and hard to decipher. There-fore, there is no other way to circumvent this prob-lem but to seek qualitative explanations basing on othermodels and methods, such as LVM [7, 45–51, 97], theWKB approximation, as well as time-dependent andtime-independent perturbation theory [98, 99]. It shouldbe noted that LVM has been applied earlier to treat thesame HO+NGP system as the present one, but withoutDDI [7]. The latter investigation also demonstrated thepresence of PRs although in a different framework. As-tonishingly, by using LVM in the present work it has beenpossible to generate a few resonances in a manner sim-ilar to those obtained in Figs. 1– 3. The LVM analysisof these PRs can be used in a qualitative manner to castmore light on them and to reach a better understand-ing of their origin. The time-averaged GP energies ofsome of our systems at PR were matched to energy levelscalculated by using the perturbation theory correspond-ing to either the HO+NGP, BOX+NGP, or QT+NGPtraps. This enables us to characterize their level struc-ture. By applying time-dependent perturbation theory,it has been demonstrated that each PR corresponds toa certain transition between an energy level m of thetrap and the energy level n to which the PR has beenmatched, identified by the maximum transitional prob-ability between a set of other probabilities to this state n . A. LVM analysis
1. Euler-Lagrange Differential Equation
The Euler-Lagrange equation for the dynamics of thewidth w ( t ) of the DBEC in a 1D harmonic trap can befound in Ref.[45]. We thus add to this LVM equationthe contribution coming from the NGP with oscillatingdepth [Eq.(2)]. Further, using our definitions for G dd and G D [Eq.(13)], this LVM equation becomes¨ w + w = 1 w + G D √ πw − β [ A + δA cos(Ω t )] w (1 + βw ) / − √ π G dd w c ( κ ) (19)with κ = d ρ /w , c ( κ ) = 1 + 10 κ − κ − κ d ( κ )(1 − κ ) , (20)and d ( κ ) = arctanh √ − κ √ − κ . (21)The width at time t = 0 is used as an initial conditionfor solving the LVM equation (19) and is taken to be thevalue w = w ( t = 0) at which the DBEC is in equilib-rium. As in Refs.[7, 46], w is obtained by solving theequation w − w − G D √ πw + Aβw (1 + βw ) / + 2 √ π G dd w c ( κ ) = 0 , (22)at κ = d ρ /w that is gotten from Eq.(19 by setting δA = 0 and ¨ w = 0.
2. Frequency of Breathing Mode
In this section, we revisit the frequency of the breath-ing mode explored earlier in some of our publications[7, 46], except that this time the effects of the DDI areadded. As before, the breathing-mode frequency canbe obtained by a linearization of Eq.(19) via w ( t ) = w + δw ( t ). However, to avoid mathematical complex-ities, it is easier to equivalently take the differential ofboth sides of (19) with respect to w at w . Consequently,one gets δ ¨ w + Ω B δw = 0 (23)with Ω B given byΩ B ( t ) = 1 + 3 w + r π G D w +[ A + δA cos(Ω t )] β (1 − βw )(1 + βw ) / +2 √ π G dd w " ∂c ( κ ) ∂w (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) w − w c ( κ ) . (24)The term between brackets yields ∂c ( κ ) ∂w (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) κ − w c ( κ ) = κ d ρ (1 − κ ) − − κ + 12 κ − κ +9 κ ( κ + 4) p − κ arctanh q − κ ! . (25)Again, it can be concluded that an imaginary value ofΩ B is obtained if1 + 3 w + r π G D w < − [ A + δA cos(Ω t )] β (1 − βw )(1 + βw ) / − √ π G dd w " ∂c ( κ ) ∂w (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) w − w c ( κ ) , (26)that leads to a damping of the DBEC oscillations.
3. LVM kinetic energy
The LVM kinetic energy is given by [7, 45, 46] e kin ( t ) = 12 [ ˙ w ( t )] + 12[ w ( t )] , (27)and is used now for f ( t ) in Eq.(14) for evaluating D kin and checking the possibility of generating PRs by theLVM equation (19).
4. LVM Results
The LVM equation (19) is solved numerically withMathematica over the same range of G dd as in Sec. IIIusing similar parameters for G D , A , δA , L x , and β , ex-cept for Ω which is set to other values for the best re-sponse of the system to changes in G dd . For example,an integer value for Ω yields resonant behavior in w ( t )and D kin contrary to a real value [7]. In this regard, itshould be emphasized that we principally aim at a qual-itative comparison with TDGPE that would help us inexplaining the PRs of Sec. III.In Fig. 5, the signal energy D kin obtained from e kin ( t ), where the latter is given by (27), is graphed as afunction of G dd and demonstrates that LVM surprisinglygenerates PRs along the same lines as the ones displayedin Sec. III. Another astonishing result is that these reso-nances are shifted in their positions along the G dd axiswhen G D is changed. This is similar to what has beenreported earlier in Figs. 1– 3 of the TDGPE results, ex-cept that in LVM they shift in the opposite direction andtheir amplitudes and shapes in Fig. 5 change somewhatbecause of this shift. This opposite behavior is an arte-fact of the model arising from the fact that LVM relieson the variational Gaussian Ansatz [45–47] ψ ( x ) = 1 π / w e − x w + iβx (28)to evaluate the mean-field Lagrangian [7, 45–51]. ThisAnsatz is much less flexible than the numerical solutionto the TDGPE resulting in this different behavior. Itshould further be iterated that it is mostly suitable forBECs in an HO trap and thus needs to be developed forother types of traps. Nevertheless the discovery of theLVM PRs as presented in Fig. 5 strongly substantiatesour results on the TDGPE PRs and supports our claimsthat they are an inherent feature in the DBECs, andnot just a result of some other influence such as noise ornumerical chaos.Figure 6 displays the time-average of Ω B ( t ) [Eq.(24)],that is h Ω B i t graphed against G dd for various G D . Fora larger G D , the height of the curve drops. This can beconnected to the PR shifting to larger resonant valuesof G dd in Fig. 5. Inspecting Eq.(24), one can see thatthe term proportional to G D is linearly added to thatproportional to G dd . Thus when h Ω B i t is changed as aresult of varying G D , the position of the correspondingPR is updated.
50, 450, 210, 4 G D = 10, Ω = 2 G dd ( a ho ) L V M D k i n / D ( ) k i n (HO+NGP) . . . . . . . . FIG. 5: Parametric resonances generated by the Euler-Lagrange differential equation (19) for a 1D DBEC in thesame HO+NGP trap of Fig. 2. The graph displays the sig-nal energy D kin of the LVM kinetic energy e kin ( t ) [Eq.(27)]versus the DDI interaction parameter G dd . D kin is normal-ized by D (0) kin = 1 × . Different values for G D and Ω areconsidered: 10 and 2, respectively (solid line); 10, 4 (dashedline); 50, 2 (dotted line); and 50, 4 (dashed-dotted line). Theinitial conditions used in solving the Euler-Lagrange equation(19) at time t = 0 are the initial “speed” ˙ w = ˙ w (0) = 0 andthe values of the equilibrium w = w (0) satisfying Eq.(22) ateach G dd with the value of d ρ in κ = d ρ /w set to 0.6. G D , G dd , and L x are in units of a ho whereas A in ¯ h ¯ ω and β in a − ho . D kin and D (0) kin are in units of (¯ h ¯ ω ) ¯ ω . G D = 10 G dd L V M h Ω B i t Ω = 2 (HO+NGP) . . . . . . FIG. 6: As in Fig. 5; but for the time average h Ω B i t of theLVM breathing-mode frequency Ω B ( t ) given by [Eq.(24)] atvarious values of G D . Solid line: G D = 10; dashed-double-dotted line: 20; triple-dotted line: 30; double-dotted line: 40;dotted line: 50. Ω B ( t ) is in units of ¯ ω . B. Energy Levels of the HO+NGP Trap viaPerturbation Theory
The goal of this section is to demonstrate that the GPtime-averaged total energies h E GP i t [Eq.(12)] at whichthe driven DBEC resonates can be generated from firstand second-order perturbation theory thereby revealinga certain ELS via the quantum numbers n correspond-ing to h E GP i t . The HO+NGP is considered here first.Within this context, the NGP is taken to be a small per-0 TABLE I: Match of the time averaged GP energies h E GP i t [Eq.(12)] with the energy levels E n obtained from second-orderperturbation theory. The system considered here is the HO+NGP trap in Fig. 2( A ) for δA = 10. The E n are obtained fromEqs.(29–33). From left to right the table lists G dd , h E GP i t , E n at quantum state n , the number of states M found for thesecond-order correction (31), and the difference between the matched energies. h E GP i t and E n are in units of ¯ h ¯ ω , and G dd in a ho . G dd h E GP i t E n n M |h E GP i t − E n | turbation within a large HO trap. This is in view of thefact that the size of our simulated system is 60 ( a ho ) cast-ing a very a high potential energy at the trap edges ofthe order of ∼
450 (¯ hω ). Therefore, it is reasonable toapply time-independent perturbation theory to calculatethe energy levels of the DBEC. The first order correctionto the energy is thus E (1) n = h φ n ( x ) | V DT ( x ) | φ n ( x ) i = A √ π n n ! Z + ∞−∞ dx H n ( x ) e − ( β +1) x , (29)where φ n ( x ) are the HO functions being the solutionsto the non-interacting Schr¨odinger equation with an HOpotential and are given by φ n ( x ) = 1 p n n ! √ π H n ( x ) e − x / . ( n = 1 , , , · · · )(30)The second-order correction arising from a number oflevels M is thus E (2) n = M X kk = n |h φ n ( x ) | V DT ( x ) | φ k ( x ) i| E (0) k − E (0) n , = A π M X kk = n n + k n ! k !( k − n ) × (cid:12)(cid:12)(cid:12)(cid:12)Z + ∞−∞ dx H n ( x ) H k ( x ) e − ( β +1) x (cid:12)(cid:12)(cid:12)(cid:12) . (31)The above Eqs.(29) and (31) can be evaluated numeri-cally by Mathematica. The energy level E n of the DBECin the HO+NGP is then the sum of the latter correctionsplus the unperturbed energy level E (0) n = (cid:18) n + 12 (cid:19) (32) such that E n = E (0) n + E (1) n + E (2) n . (33)Table I lists for example the h E GP i t of the DBEC atthe resonant values of G dd for the system in Fig. 2(A)at δA = 10 with corresponding values E n obtained fromEqs.(29–33). The values of h E GP i t and E n agree verywell after matching and consequently an analytical ELScan be deduced.It makes sense to state that when h E GP i t matcheswith an E n the PR occurs. Therefore, a Green’s functionof the form G ( E ) ∼ E − E n + i Γcould be proposed that accounts for these PRs, where Γis a width and E = h E GP i t . The latter h E GP i t increasewith G dd while |h V eff (0) i| [see Sec. III D and Eq.(17)]drops with it. The quantum number n rises as well indi-cating that the PRs shift to higher frequencies. In fact,this can be related to the results of the experiment ofBalik et al. [1] who excited PRs via a CO laser with anintensity modulated via a depth h . They showed that areduction in h shifts the PRs to higher frequencies. Ourfindings above are analogous to theirs for this case. C. Energy Levels in a BOX+NGP trap viaPerturbation Theory
The goal of this section is the same as the prior one,but for the box potential. The wave function φ n ( x ) for a1D box of length 2 L x is appropriately φ n ( x ) = 1 √ L x cos (cid:18) nπx L x (cid:19) , ( n = 1 , , , · · · ) (34)again being the solution to the box-potential Schr¨odingerequation in the noninteracting case. Similar to Eqs.(29)and (31) this yields1 TABLE II: As in Table I; but for the BOX+NGP in Fig. 1( A ) with the corrections to the energies obtained via φ n ( x ) given byEq.(34). h E GP i t and E n are in units of ¯ h ¯ ω , and G dd in a ho . G dd h E GP i t E n n M |h E GP i t − E n | E (1) n = AL x Z + L x − L x dx (cid:20) cos (cid:18) nπx L x (cid:19)(cid:21) e − βx (35)and E (2) n = 8 A π M X kk = n k − n ) × (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)Z + L x − L x dx cos (cid:18) nπx L x (cid:19) cos (cid:18) kπx L x (cid:19) e − βx (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (36)where the unperturbed energies are given by E (0) n = n π L x (37)(which are in trap units). Accordingly, Table II lists thesame quantities as in Table I; but for the BOX+NGP. Itis found again that perturbation theory is able to repro-duce the GP energies at which the PRs occur. It makessense also here to state that when h E GP i t matches E n a PR arises. Notably, the PRs correspond to very highenergy levels n because the size of the system 2 L x isquite large ( L x = 51 .
2) yielding small differences in the E n and thereby a large density of states. However, thebehavior of h E GP i t with G dd is opposite to that of theHO+NGP. The reduction in the depth of the effectivepotential with G dd causes the PR energies to drop shift-ing them to lower frequencies instead. Thus the DDI inthe case of a BOX+NGP cannot prevent the bosons fromfalling deeper into the NGP. But why does this happenfor the box and not the HO? Two reasons come to mind:(1) the density in the BOX+NGP is lower than in theHO+NGP as a result of which the dipolar nonlinearity isweaker in the box; (2) the quantum pressure in the HOtrap is larger than in the box. This finding, in conjunc-tion with the one in the previous section, demonstratesthe role of an interplay between the trapping geometryand interactions in defining the positions of the PRs inthe spectrum of D kin . D. Energy Levels in a QT+NGP via PerturbationTheory
The energy levels of a general power-law trap V tr ( x ) = α | x | ν (38)can be obtained from WKB theory as described in ad-vanced quantum mechanics textbooks [99] and yields (bysetting ¯ h = m = 1) E (0) n = α "(cid:18) n + 12 (cid:19) r π α Γ (cid:0) ν + (cid:1) Γ( ν + 1) νν + 2 . (39)The left and right turning points a and b , respectively,are obtained from a = − b = − E (0) n α ! /ν , (40)where in the present case α = 1 / ν = 4. The WKBwave functions are determined by the usual matching ofthe solutions at the turning points and read φ n ( x ) = R p κ n ( x ) e − R ax κ n ( y ) dy ; x ≪ a R p k n ( x ) sin "Z bx k n ( y ) dy + π ; a < x < bR p κ n ( x ) e − R xb κ n ( y ) dy ; x ≫ b (41)with k n ( x ) = q E (0) n − V tr ( x )] ,κ n ( x ) = ik ( x ) , (42)the wave vectors. R is a constant determined by normal-izing the ground-state wave function at n = 0, that is R + ∞−∞ | φ ( x ) | = 1 so that2 TABLE III: As in Table I; but for the QT+NGP in Fig. 3(A) with the corrections to the energies obtained via φ n ( x ) given byEq.(41) with the wave vectors (42). h E GP i t and E n are in units of ¯ h ¯ ω , and G dd in a ho . G dd h E GP i t E n n M |h E GP i t − E n | Z a −∞ dx (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) R p κ ( x ) e − R ax κ ( y ) dy (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + Z ba dx (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) R p k ( x ) sin "Z bx k ( y ) dy + π + Z + ∞ b dx (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) R p κ ( x ) e − R ax κ ( y ) dy (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 1 . (43)The first-order correction to the unperturbed energy E (0) n [Eq.(39)] via the NGP perturbation becomes E (1) n = 4 R A Z ba dx k n ( x ) sin "Z bx k n ( x ) dx + π e − βx . (44)It should be noted, that the wave functions in the regimesoutside the QT have not been considered in evaluating E (1) n since they do not overlap with the NGP and there-fore do not make any contribution. Along the same lines,the second-order correction yields E (2) n = X mm = n E (0) m − E (0) n × (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) R A Z ba dx p k n ( x ) k m ( x ) e − βx × sin "Z bx k n ( y ) dy + π sin "Z bx k m ( y ) dy + π , (45)where E (0) n is given by (39).Table III shows again that h E GP i t at the PRs in theQT can be closely matched to the energies of second-order perturbation theory. It should be emphasized thatwell-behaved, exact analytic solutions to the Schr¨odingerequation with a quartic oscillator are hitherto unknown.Although, the Heun function [100] may provide a solutionwithin a restricted range in the neighborhood of the trapcenter, outside the latter it diverges to very large values.We have found that even a numerical solution by Math-ematica yields divergent solutions away from the center(not shown here). E. Transition Probabilities
1. Harmonic Trap
In what follows, the transition probabilities at whichPRs occur in an HO+NGP trap are evaluated. Basingon the assumption that the NGP and DDI are pertur-bations inducing transitions between the different statesin the traps, the probabilities for these can be computedaccording to time-dependent perturbation theory. In thisregard, the probability for a transfer from some state m to state n is from standard quantum mechanics textbooksgiven by P ( n, m, t ) = | c n,m ( t ) | , (46)where as usual c n,m ( t ) = λi ¯ h Z t dt ′ e i ( E (0) m − E (0) n ) t ′ / ¯ h h φ m ( x ) | V ( t ′ ) | φ n ( x ) i (47)with V ( t ′ ) being a time-dependent perturbation. Con-sidering first φ n ( x ) to be the HO states given by (30)and V ( t ′ ) the time-dependent NGP (2), then it is easyto show that the probability introduced by the NGP is(setting ¯ h = 1) P DT ( n, m, t ) = (cid:12)(cid:12)(cid:12)(cid:12) λi √ π n + m n ! m ! × (cid:26)Z + ∞−∞ H m ( x ) H n ( x ) e − ( β +1) x dx (cid:27) × " AF ( m, n, t ) + δA ( F +Ω ( m, n, t ) + F − Ω ( m, n, t ) , (48)where F ( m, n, t ) = e i ( m − n ) t − i ( m − n ) ,F +Ω ( m, n, t ) = e i ( m − n +Ω) t − i ( m − n + Ω) ,F − Ω ( m, n, t ) = e i ( m − n − Ω) t − i ( m − n − Ω) . (49)On the other hand, to compute the transition probabilitybecause of the DDI, the Fourier transform (FT) of the 1D3dipolar potential V D ( x − x ′ ), such as the one presentedby Sinha and Santos [101], needs to be applied to avoidthe singularity in the DDI. This FT reads for a quasi 1DBose gas ˜ V D ( k ) = 4 αd ℓ " − σe σ Γ(0 , σ ) , (50)were Γ(0 , σ ) is the incomplete gamma function, and σ = k ℓ /
4. Here ℓ = p ¯ h/ ( µ ¯ ω ) is a length scale with µ = m/ α has values be-tween − / φ = π/
2) and 1 ( φ = 0) where φ is the anglebetween the dipoles and the longitudinal axis of the trap.In the present case ℓ = d ρ is the width of the wave func-tion in the transverse direction. Appendix B shows howto cast (50) in trap units, that is ˜ V D → ˜ V D / (¯ hω ).Following Muruganandam [45], the spatial integral re-quired in Eq.(47) for the transition between states n and m , namely h φ n ( x ) | V ( t ′ ) | φ m ( x ) i → N Z + ∞−∞ dx Z + ∞−∞ dx ′ φ ⋆n ( x ) φ m ( x ′ ) V D ( x − x ′ ) , (51)can be rewritten using the FT [Eq.(50)] with its strengthgiven by (B5) and the FT of the HO functions e φ n ( k ) [7]where e φ n ( k ) = FT[ φ n ( x )] = 1 p √ π n n ! ( − i ) n H n ( k ) e − k / . (52)This then yields that h φ n ( x ) | V ( t ′ ) | φ m ( x ) i → N π Z + ∞−∞ dk e φ ⋆n ( k ) e φ m ( − k ) ˜ V D ( k ) . (53)As such, the probability for the DDI reads then P DDI ( n, m, t ) = λ (cid:12)(cid:12)(cid:12)(cid:12) G dd ( − i ) n + m √ π n + m n ! m ! × Z + ∞−∞ H n ( k ) H m ( − k ) e − k × " − k d ρ e k d ρ / Γ(0 , k d ρ / F ( m, n, t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (54)The above two probabilities P DT ( n, m, t ) [Eq.(48)] and P DDI ( n, m, t ) can be easily evaluated numerically byMathematica for the integer values of m and n . These TABLE IV: Most probable transitions to the states n at whichthe GP energies are matched in Table I for the HO+NGPtrap. The transition probability is computed by Eqs.(48),(54), and (55). G dd is in units of a ho . G dd m n P ( n, m, t ) /λ two probabilities are dependent on each other and there-fore the total probability is P ( n, m, t ) = P DT ( n, m, t ) P DDI ( n, m, t ) . (55)Table IV displays P ( n, m, t ) (in units of arbitrary λ )for the most probable transitions to the states n matchedin Table I for a length of one driving cycle t = 0 .
1. Inthis, it is assumed that λ ≪
1. It should be noted, thata significant number of m ↔ n transitions turned out tobe prohibited with identically zero values for P ( n, m, t ).This explains the absence of PRs for values of G dd witha time-averaged GP that nevertheless can be matchedto one of the energy levels given by perturbation theoryin Sec. IV B. Thus there are two conditions for PRs tooccur: (1) its energy should match one of the trap levelsand (2) the transition to this level should have a highprobability.
2. Box Potential
As in the previous section, the goal is to find themost probable transitions at which the PRs arise in theBOX+NGP trap. In this case, P DT ( n, m, t ) = (cid:12)(cid:12)(cid:12)(cid:12) λi L x × (Z + L x − L x cos (cid:18) mπx L x (cid:19) cos (cid:18) nπx L x (cid:19) e − βx dx ) × " AF ( m, n, t ) + δA ( F +Ω ( m, n, t ) + F − Ω ( m, n, t )) (56)and4 TABLE V: As in Table IV; but for the BOX+NGP trap. G dd is in units of a ho . G dd m n P ( n, m, t ) /λ P DDI ( n, m, t ) = (cid:12)(cid:12)(cid:12)(cid:12) λi L x G dd × Z + K − K sin (cid:16) kL x − nπ (cid:17) kL x − nπ + sin (cid:16) kL x + nπ (cid:17) kL x + nπ × sin (cid:16) kL x − mπ (cid:17) kL x − mπ + sin (cid:16) kL x + mπ (cid:17) kL x + mπ × " − k d ρ e k d ρ / Γ[0 , k d ρ / F ( m, n, t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (57)where the FT of cos ( nπx/ (2 L x )) inside the box is F T [cos ( nπx/ (2 L x ))] =2 L x sin (cid:16) kL x − nπ (cid:17) kL x − nπ + sin (cid:16) kL x + nπ (cid:17) kL x + nπ . (58)In (57), K is a cutoff momentum introduced to simplifythe integration, but which would still give the same resultas going to infinite values. Table V is as in Table IV, butfor the box. The most probable transitions are found tocome from relatively high states m towards n matched inTable II.Finally, it must be noted that because of the diver-gent wavefunction of the QT+NGP already mentioned inSec. IVE, an analytical evaluation of the P DT and P DDI is currently extremely difficult if not impossible, and istherefore left for the future.
V. SUMMARY AND CONCLUSIONS
In summary then, we have reported mean-field para-metric resonances (PRs) in a one-dimensional (1D) dipo-lar Bose-Einstein condensate (DBEC) excited by a nega-tive Gaussian potential (NGP) with a periodically oscil-lating depth. The PRs have been detected by the signalenergy [44], a quantity that closely resembles the time-average of the square of the energy (in this work the kinetic energy). It is similar to the root-mean-squared(RMS) value of an oscillating electrical signal arising froman alternating voltage or current source. The latter NGPwas for modelling a red laser light source with oscillatingintensity in a manner similar to Refs.[1, 43].By a matching of the PR energies to energy levelscomputed by time-independent perturbation theory wewere able to characterize the energy-level structure (ELS)of a DBEC in a complex trap, i.e., a trap to which anNGP has been added. Within the purpose of the lattercharacterization, a few different traps have been appliedto test the effect of confining geometry and its inter-play with dipole-dipole interactions (DDI) on the PRs.It turns out the DDI play an important role in definingthe amplitude of the PRs and their energies. The DDIreduces the depth of the effectve mean-field potential,thereby controlling the occupancy of the NGP, and inturn the positions and strengths of the PRs.The key feature of these PRs is that their propertiesare sensitive to the ELS of the confining potential, i.e.,external trap+NGP. The trapping geometry determinesthe wavefunction that describes the system, and in turnthe wavefunction in conjunction with the NGP and DDIdetermine the transition probabilities between differentquantum states of the confinement. The PRs correspondto the most probable transitions determined from time-dependent perturbation theory. Nevertheless, it is rathersurprising that these PRs correspond only to specific val-ues of the principal quantum numbers of high transitionalprobabilities, say m and n , and that they do not arise atother values ( m , n ) of comparable, or slightly lower, prob-abilities. One might even begin to assume the presenceof magic numbers for n and m , although this assertionneeds to be proven.Most importantly, DBEC PRs could be astonishinglyproduced by the Lagrangian variational method (LVM)[46] that approximates the wavefunctions by a simpleGaussian Ansatz. However, it was not possible to useLVM for the box and quartic trap because the GaussianAnsatz had been originally designed for the harmonictrapping case. The LVM PRs are a further manifesta-tion of the fact that these resonances are an inherent fea-ture present in the DBECs awakened by external drivingagents.The present work can also be somewhat related toa previous study on resonances and dynamical fragmen-tation in a stirred BEC [102] in which a series of PRsin the total energy arises as the rotational frequency isincreased. In this regard, the authors conclude that frag-mentation of the gas appears simultaneously with theresonant absorption of energy and angular momentumfrom the external excitation agent. Whereas their PRsare associated with fragmentation of the BEC, in thepresent mean-field GP formulation there is no fragmen-tation. The PRs of the current work arise because of res-onant absorption of momentum from the dynamic NGP.This momentum arises from the time-varying force Eq.(3)which is the gradient of the NGP that acts along the5length of the BEC.The previous results in Secs. III A - C are also a mani-festation for the role of a mean-field dipolar nonlinearityin defining the properties of the PRs and − via its in-terplay with the trapping geometry − also their spacings∆ G dd . One question that arises is if similar results couldbe obtained by simulations of the same systems using themany-body Schr¨odinger equation with DDI. In this case,the DDI do not arise in the form of a nonlinearity andit would be interesting to make comparisons between themany-body effects of the DDI and that of the DDI non-linearity. Acknowledgement
The calculations were performed on the PARADOX-IV supercomputing facility at the Scientic ComputingLaboratory of the Institute of Physics Belgrade, sup-ported in part by the Ministry of Education, Science andTechnological Development of the Republic of Serbia un-der project No. ON171017. We thank William J Mullin(UMASS, Amherst USA) for insightful comments on anearlier version of this manuscript that helped us to im-prove it substantially. The authors declare that there areno conflicts of interest.
Appendix A: Numerics
The TDGPE [Eq.(7)] is solved numerically via thesplit-step Crank-Nicolson (CN) method [52, 80] in realtime for the HO ( p = 2), QT ( p = 4), and box potentialtraps. For the HO and QT, it is propagated along a gridof N x = 6000 pixels of size ∆ x = 0 .
01 with L x = 1in Eq.(1). The time step chosen is ∆ t = 1 × − ,the number of time steps in the transient run is set to N pas = 10 , and the same in the final run N run = 10 .For the box potential N x = 2048, ∆ x = 0 . L x = 51 . t = 0 . N pas = 200000, and N run = 200000. Inall cases, the system is initialized with a stationary NGPof depth A = −
30 via the imaginary-time CN method fora number of N stp = 2 × time steps, and then takenthrough N pas = 20000 steps in the transient, and endingwith N run = 20000 steps in the final run. The time step∆ t is the same as in the corresponding real-time simu-lations. In the CN codes that have been used, the finalrun (NRUN) just complements the transient one (NPAS).The latter could be the stage during which the BEC is al-lowed to evolve and relax to a stable state. The final runcould then be the stage where the BEC is optionally ex-cited by an external agent in order to examine its ensuingdynamics. The authors of the code originally separatedthe transient and final runs for organizational purposes. This is in order to have data files in the relaxation stageof the BEC towards a stable state separate from the onewhere this stable BEC is suddenly excited by an externaldriving force. In contrast, the present work considers aBEC that is continuously excited in both transient andfinal runs, so there is no difference between both of them.Sets of runs were performed each at a fixed value of G D . For each G D , the system was scanned along arange of G dd from 0 to 400 in steps of 2. It should beemphasized that for each G dd there was a run. This largenumber of runs has been conducted in the form of parallelarray jobs on the high performance computational clusterof the Scientific Computing Laboratory of the Institute ofPhysics in Belgrade, Serbia. Each simulation took aboutfive days to finish. In essence, this has been a heavilycomputational project. Appendix B: Units of V D ( k ) Beginning with the dipole moment d where [52] d = µ ˜ µ π (B1)and the dipolar scattering length a dd = µ ˜ µ m π ¯ h , (B2)˜ µ being the magnetic dipole moment and µ the magneticpermeability, we substitute (B2) into (B1) to eliminate˜ µ so that d = 3¯ h a dd /m . Then the rescaling d → d / (¯ hω ) yields˜ d = d ¯ hω = 3¯ h a dd m ¯ hω = 3˜ a dd a ho , (B3)where ˜ a dd = a dd /a ho . Hence, the coefficient of Eq.(50)becomes in trap units4 α ˜ d d ρ = 4 α ˜ d ρ a dd a ho , (B4)with ˜ d ρ = d ρ /a ho . Using Eq.(13) for G dd , one eventu-ally gets that the coefficient with α = − / αd ℓ → π G dd a ho N . (B5) [1] S. Balik, A. L. Win, and M. D. Havey,Phys. Rev. A , 023404 (2009), URL https://link.aps.org/doi/10.1103/PhysRevA.80.023404 . [2] R. R. Sakhel and A. R. Sakhel, J. Phys:Condens. Matter , 315401 (2020), URL https://iopscience.iop.org/article/10.1088/1361-648X/ab7f06 .[3] J. H. V. Nguyen, M. C. Tsatsos, D. Luo, A. U. J.Lode, G. D. Telles, V. S. Bagnato, and R. G.Hulet, Phys. Rev. X , 011052 (2019), URL https://link.aps.org/doi/10.1103/PhysRevX.9.011052 .[4] T. Chen, K. Shibata, Y. Eto, T. Hirano, andH. Saito, Phys. Rev. A , 063610 (2019), URL https://link.aps.org/doi/10.1103/PhysRevA.100.063610 .[5] C.-X. Zhu, W. Yi, G.-C. Guo, and Z.-W.Zhou, Phys. Rev. A , 023619 (2019), URL https://link.aps.org/doi/10.1103/PhysRevA.99.023619 .[6] Z.-C. Li, Q.-H. Jiang, Z. Lan, W. Zhang, andL. Zhou, Phys. Rev. A , 033617 (2019), URL https://link.aps.org/doi/10.1103/PhysRevA.100.033617 .[7] A. R. Sakhel and R. R. Sakhel, Journal of Low Temper-ature Physics , 120 (2018), ISSN 1573-7357, URL https://doi.org/10.1007/s10909-017-1826-7 .[8] P. Molignini, L. Papariello, A. U. J. Lode, andR. Chitra, Phys. Rev. A , 053620 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.98.053620 .[9] S. Robertson, F. Michel, and R. Parentani,Phys. Rev. D , 056003 (2018), URL https://link.aps.org/doi/10.1103/PhysRevD.98.056003 .[10] S. Lellouch, M. Bukov, E. Demler, and N. Goldman,Phys. Rev. X , 021015 (2017).[11] I. Vidanovi´c, H. Al-Jibbouri, A. Balaˇz, and A. Pel-ster, Physica Scripta T149 , 014003 (2012), URL https://doi.org/10.1088%2F0031-8949%2F2012%2Ft149%2F014003 .[12] A. Posazhennikova, M. Trujillo-Martinez, andJ. Kroha, Phys. Rev. Lett. , 225304 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.116.225304 .[13] D. Kobyakov, V. Bychkov, E. Lundh, A. Bezett, andM. Marklund, Phys. Rev. A , 023614 (2012), URL https://link.aps.org/doi/10.1103/PhysRevA.86.023614 .[14] J.-K. Xue, G.-Q. Li, A.-X. Zhang, andP. Peng, Phys. Rev. E , 016606 (2008), URL https://link.aps.org/doi/10.1103/PhysRevE.77.016606 .[15] P. Engels, C. Atherton, and M. A. Hoe-fer, Phys. Rev. Lett. , 095301 (2007), URL https://link.aps.org/doi/10.1103/PhysRevLett.98.095301 .[16] A. I. Nicolin, R. Carretero-Gonz´alez, and P. G.Kevrekidis, Phys. Rev. A , 063609 (2007), URL https://link.aps.org/doi/10.1103/PhysRevA.76.063609 .[17] M. Kr¨amer, C. Tozzo, and F. Dalfovo,Phys. Rev. A , 061602 (2005), URL https://link.aps.org/doi/10.1103/PhysRevA.71.061602 .[18] C. Tozzo, M. Kr¨amer, and F. Dalfovo,Phys. Rev. A , 023613 (2005), URL https://link.aps.org/doi/10.1103/PhysRevA.72.023613 .[19] T. Chen and B. Yan, Phys.Rev. A , 063615 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.98.063615 .[20] Sabari, Subramaniyan and Kumar, R. Kishor,Eur. Phys. J. D , 48 (2018), URL https://doi.org/10.1140/epjd/e2018-80354-2 .[21] Cairncross, William and Pelster, Axel,Eur. Phys. J. D , 106 (2014), URL https://doi.org/10.1140/epjd/e2014-40835-x .[22] I. Vidanovi´c, A. Balaˇz, H. Al-Jibbouri, andA. Pelster, Phys. Rev. A , 013618 (2011), URL https://link.aps.org/doi/10.1103/PhysRevA.84.013618 .[23] V. I. Arnold, A. Weinstein, and K. Vogtmann, Mathe- matical Methods of Classical Mechanics, Graduate Textsin Mathematics , vol. 60 (Springer,Berlin, 1989).[24] L. D. Landau and E. Lifshitz,
Mechanics, Course ofTheoretical Physics , vol. 1 (Butterworth-Heinemann,Oxford, 1976).[25] A. Armaroli and F. Biancalana, Phys.Rev. A , 063848 (2013), URL https://link.aps.org/doi/10.1103/PhysRevA.87.063848 .[26] F. Beato and A. Palacios-Laloy, EPJ Quantum Tech-nology , 1 (2020).[27] N. V. Fomin, O. L. Shalaev, and D. V. Shant-sev, Journal of Applied Physics , 8091(1997), https://doi.org/10.1063/1.365417, URL https://doi.org/10.1063/1.365417 .[28] M. Calvo, Phys. Rev. B , 10953 (1999), URL https://link.aps.org/doi/10.1103/PhysRevB.60.10953 .[29] G. Hackenbroich, B. Rosenow, and H. A. Wei-denm¨uller, Phys. Rev. Lett. , 5896 (1998), URL https://link.aps.org/doi/10.1103/PhysRevLett.81.5896 .[30] F. Matera, A. Mecozzi, M. Romagnoli, andM. Settembre, Opt. Lett. , 1499 (1993), URL http://ol.osa.org/abstract.cfm?URI=ol-18-18-1499 .[31] G. Bismut, B. Pasquiou, E. Mar´echal, P. Pedri,L. Vernac, O. Gorceix, and B. Laburthe-Tolra,Phys. Rev. Lett. , 040404 (2010), URL https://link.aps.org/doi/10.1103/PhysRevLett.105.040404 .[32] K. G´oral and L. Santos, Phys.Rev. A , 023613 (2002), URL https://link.aps.org/doi/10.1103/PhysRevA.66.023613 .[33] C. Mishra and R. Nath, Phys.Rev. A , 033633 (2016), URL https://link.aps.org/doi/10.1103/PhysRevA.94.033633 .[34] F. W¨achtler and L. Santos, Phys.Rev. A , 043618 (2016), URL https://link.aps.org/doi/10.1103/PhysRevA.94.043618 .[35] B. Schulz, S. Sala, and A. Saenz, New Jour-nal of Physics , 065002 (2015), URL https://doi.org/10.1088%2F1367-2630%2F17%2F6%2F065002 .[36] A. R. P. Lima and A. Pelster, Phys.Rev. A , 041604 (2011), URL https://link.aps.org/doi/10.1103/PhysRevA.84.041604 .[37] F. Ancilotto and F. Toigo, Phys.Rev. A , 023617 (2014), URL https://link.aps.org/doi/10.1103/PhysRevA.89.023617 .[38] R. M. Wilson and J. L. Bohn, Phys.Rev. A , 023623 (2011), URL https://link.aps.org/doi/10.1103/PhysRevA.83.023623 .[39] G. P. Agrawal, Phys. Rev. Lett. , 880 (1987), URL https://link.aps.org/doi/10.1103/PhysRevLett.59.880 .[40] F. Abdullaev, S. Darmanyan, A. Kobyakov,and F. Lederer, Physics Letters A ,213 (1996), ISSN 0375-9601, URL .[41] S. Ambomo, C. M. Ngabireng, P. T. Dinda,A. Labruy`ere, K. Porsezian, and B. Kalithasan,J. Opt. Soc. Am. B , 425 (2008), URL http://josab.osa.org/abstract.cfm?URI=josab-25-3-425 .[42] Michael C. Garrett, Adrian Ratnapala, Eikbert D. vanOoijen, Christopher J. Vale, Kristian Weegink, Se-bastian K. Schnelle, Otto Vainio, Norman R. Heck-enberg, Halina Rubinsztein-Dunlop, and MatthewJ. Davis, Phys. Rev. A , 013630 (2011), URL https://doi.org/10.1103/PhysRevA.83.013630 .[43] L. W. Clark, L.-C. Ha, C.-Y. Xu, and C. Chin, Phys. Rev. Lett. , 155301 (2015).[44] A. V. Oppenheim and G. C. Verghese,
Signals, Systems,and Inference (Prentice Hall Signal Processing Series,Pearson, 2015), 1st ed.[45] P. Muruganandam and S. K. Adhikari, Las. Phys. ,813 (2012).[46] Roger Sakhel and Asaad Sakhel, J. Phys. B: At. Mol.Opt. Phys. , 105301 (2017).[47] V. M. Perez-Garcia, H. Michinel, J. I. Cirac, M. Lewen-stein, and P. Z¨oller, Phys. Rev. Lett. , 5320 (1996).[48] V. M. Perez-Garcia, H. Michinel, J. I. Cirac, M. Lewen-stein, and P. Z¨oller, Phys. Rev. A , 1424 (1997).[49] Hamid Al-Jibbouri, Ivana Vidanovic, Antun Balaz, andAxel Pelster, J. Phys. B: At. Mol. Opt. Phys. , 065303(2013).[50] G. M. Falco, J. Phys. B: At. Mol. Opt. Phys. , 215303(2009).[51] Qi Wei, Liang Zhao-Xin, and Zhang Zhi-Dong, ChinesePhysics Letters , 060303 (2013).[52] R. Kishor Kumar, Luis E. Young-S., DuˇssanVudragovi´c, Antun Balaˇz, Paulsamy Muru-ganandam, S.K.Adhikari, Computer PhysicsCommunications , 117 (2015), URL .[53] G. D. Bruce, S. L. Bromley, G. Smirne,L. Torralbo-Campo, and D. Cassettari,Phys. Rev. A , 053410 (2011), URL https://link.aps.org/doi/10.1103/PhysRevA.84.053410 .[54] M. Hammes, D. Rychtarik, H.-C. N¨agerl, andR. Grimm, Phys. Rev. A , 051401 (2002), URL https://link.aps.org/doi/10.1103/PhysRevA.66.051401 .[55] C. Tuchendler, A. M. Lance, A. Browaeys,Y. R. P. Sortais, and P. Grangier,Phys. Rev. A , 033425 (2008), URL https://link.aps.org/doi/10.1103/PhysRevA.78.033425 .[56] D. M. Stamper-Kurn, H.-J. Miesner, A. P.Chikkatur, S. Inouye, J. Stenger, and W. Ket-terle, Phys. Rev. Lett. , 2194 (1998), URL https://link.aps.org/doi/10.1103/PhysRevLett.81.2194 .[57] D. Comparat, A. Fioretti, G. Stern,E. Dimova, B. L. Tolra, and P. Pillet,Phys. Rev. A , 043410 (2006), URL https://link.aps.org/doi/10.1103/PhysRevA.73.043410 .[58] D. Jacob, E. Mimoun, L. D. Sarlo,M. Weitz, J. Dalibard, and F. Gerbier,New J. Phys. , 065022 (2011), URL https://doi.org/10.1088%2F1367-2630%2F13%2F6%2F065022 .[59] T. L. Gustavson, A. P. Chikkatur, A. E. Lean-hardt, A. G¨orlitz, S. Gupta, D. E. Pritchard, andW. Ketterle, Phys. Rev. Lett. , 020401 (2001), URL https://link.aps.org/doi/10.1103/PhysRevLett.88.020401 .[60] M. D. Barrett, J. A. Sauer, and M. S. Chap-man, Phys. Rev. Lett. , 010404 (2001), URL https://link.aps.org/doi/10.1103/PhysRevLett.87.010404 .[61] M. Schulz, H. Crepaz, F. Schmidt-Kaler, J. Es-chner, and R. Blatt, J. Mod. Opt. , 1619 (2007),https://doi.org/10.1080/09500340600861740, URL https://doi.org/10.1080/09500340600861740 .[62] N. P. Proukakis, J. Schmiedmayer, and H. T. C.Stoof, Phys. Rev. A , 053603 (2006), URL https://link.aps.org/doi/10.1103/PhysRevA.73.053603 .[63] R. B. Diener, Biao Wu, M. G. Raizen, andQ. Niu, Phys. Rev. Lett. , 070401 (2002), URL https://link.aps.org/doi/10.1103/PhysRevLett.89.070401 .[64] T. Aioi, T. Kadokura, T. Kishimoto, andH. Saito, Phys. Rev. X , 021003 (2011), URL https://link.aps.org/doi/10.1103/PhysRevX.1.021003 .[65] H. Uncu, D. Tarhan, E. Demiralp, and O. E.M¨ustecaplioglu, Las. Phys. , 331 (2008), URL https://doi.org/10.1134/S1054660X08030237 .[66] C. Weitenberg, S. Kuhr, K. Mølmer, and J. F.Sherson, Phys. Rev. A , 032322 (2011), URL https://link.aps.org/doi/10.1103/PhysRevA.84.032322 .[67] A. V. Carpentier, J. Belmonte-Beitia, H. Michinel,and M. I. Rodas-Verde, J. of Mod. Opt. , 2819(2008), https://doi.org/10.1080/09500340802209763,URL https://doi.org/10.1080/09500340802209763 .[68] N. G. Parker, N. P. Proukakis, M. Leadbeater, andC. S. Adams, Phys. Rev. Lett. , 220401 (2003), URL https://link.aps.org/doi/10.1103/PhysRevLett.90.220401 .[69] C. J. Pethick and S. H., Bose-Einstein Condensation inDilute Gases (Cambridge University Press, Cambridge,2008), 2nd ed.[70] T. Lahaye, C. Menotti, L. Santos, M. Lewen-stein, and T. Pfau, Reports on Progressin Physics , 126401 (2009), URL https://doi.org/10.1088%2F0034-4885%2F72%2F12%2F126401 .[71] D. H. J. O’Dell, S. Giovanazzi, and C. Eber-lein, Phys. Rev. Lett. , 250401 (2004), URL https://link.aps.org/doi/10.1103/PhysRevLett.92.250401 .[72] S. Yi and L. You, Phys. Rev. A , 041604 (2000), URL https://link.aps.org/doi/10.1103/PhysRevA.61.041604 .[73] J. L. Bohn, M. Cavagnero, and C. Ticknor,New Journal of Physics , 055039 (2009), URL https://doi.org/10.1088%2F1367-2630%2F11%2F5%2F055039 .[74] K. G´oral, K. Rza¸˙zewski, and T. Pfau,Phys. Rev. A , 051601 (2000), URL https://link.aps.org/doi/10.1103/PhysRevA.61.051601 .[75] A. J. Olson, D. L. Whitenack, and Y. P.Chen, Phys. Rev. A , 043609 (2013), URL https://link.aps.org/doi/10.1103/PhysRevA.88.043609 .[76] T. Koch, T. Lahaye, J. Metz, B. Fr¨ohlich, A. Gries-maier, and T. Pfau, Nat. Phys. , 218 (2008), URL https://doi.org/10.1038/nphys887 .[77] M. Lu, N. Q. Burdick, S. H. Youn, and B. L.Lev, Phys. Rev. Lett. , 190401 (2011), URL https://link.aps.org/doi/10.1103/PhysRevLett.107.190401 .[78] S. H. Youn, M. Lu, U. Ray, and B. L.Lev, Phys. Rev. A , 043425 (2010), URL https://link.aps.org/doi/10.1103/PhysRevA.82.043425 .[79] A. Filinov and M. Bonitz, Phys.Rev. A , 043628 (2012), URL https://link.aps.org/doi/10.1103/PhysRevA.86.043628 .[80] P. Muruganandam and S. K. Adhikari, ComputerPhysics Communications , 1888 (2009), URL .[81] D. Vudragovi´c, I. Vidanovi´c, A. Balaˇz, P. Mu-ruganandam, and S. K. Adhikari, ComputerPhysics Communications , 2021 (2012), URL .[82] V. Lonˇcar, A. Balaˇz, A. Bogojevi´c, S. ˇSkrbi´c, P.Muruganandam, and S. K. Adhikari, ComputerPhysics Communications , 406 (2016), URL .[83] V. Lonˇcar, L. E. Young-S., S. ˇSkrbi´c,P. Muruganandam, S. K. Adhikari, andA. Balaˇz, Computer Physics Communica-tions , 190 (2016), ISSN 0010-4655, URL .[84] L. Young-S., D. Vudragovi´c, P. Muruganan-dam, S. K. Adhikari, and A. Balaˇz, ComputerPhysics Communications , 209 (2016), URL .[85] B. Satari´c, V. Slavni´c, A. Belic, A. Balaˇz, P.Muruganandam, and S. K. Adhikari, ComputerPhysics Communications , 411 (2016), URL .[86] L. E. Young-S., P. Muruganandam, S. K.Adhikari, V. Lonˇcar, D. Vudragovi´c, andA. Balaˇz, Computer Physics Communica-tions , 503 (2017), ISSN 0010-4655, URL .[87] R. Kishor Kumar, V. Lonˇcar, P. Muruganandam,S. K. Adhikari, and A. Balaz, Computer Physics Com-munications , 74 (2019), ISSN 0010-4655, URL .[88] E. E. Nikitin and L. P. Pitaevskii, arxiv:cond-mat/050868v41 (2005).[89] Roger R. Sakhel, Asaad R. Sakhel, Humam B. Ghassib,J. Low. Temp. Phys. , 177 (2013).[90] R. R. Sakhel, A. R. Sakhel, H. B. Ghassib,and A. Balaˇz, The European Physical Jour-nal D , 66 (2016), ISSN 1434-6079, URL https://doi.org/10.1140/epjd/e2016-60085-2 .[91] A. R. Sakhel, Physica B: Condensed Mat-ter , 72 (2016), ISSN 0921-4526, URL .[92] R. R. Sakhel and A. R. Sakhel, Journal of Low Temper-ature Physics , 106 (2019), ISSN 1573-7357, URL https://doi.org/10.1007/s10909-018-2068-z .[93] R. Onofrio, C. Raman, J. M. Vogels, J. R. Abo-Shaeer,A. P. Chikkatur, and W. Ketterle, Phys. Rev. Lett. ,2228 (2000).[94] M. R. Andrews, M.-O. Mewes, N. J. van Druten, D. S.Durfee, D. M. Kurn, and W. Ketterle, Science , 84(1996).[95] W. Ketterle, D. S. Durfee, and D. M. Stamper-Kurn,in Proceedings of the Enrico Fermi International Schoolof Physics , edited by M. Inguscio, S. Stringari, andC. E. Wieman (IOS Press, Amsterdam, 1999), vol. CXL,p. 67.[96] J. T. Mendon¸ca, H. Ter¸cas, and A. Gam-mal, Phys. Rev. A , 063610 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.97.063610 .[97] Z.-Ch. Wang, H.-W. Yang, and S. Yin, Eur. Phys. J. D , 117 (2002).[98] Eugene Merzbacher, Quantum Mechanics (John Wileyand Sons, Inc., New York, 1998), Third ed.[99] D. J. Griffiths and D. F. Schr¨oter,
Introduction to Quan-tum Mechanics (Cambridge University Press, 2018), 3rded.[100] A. Ronveaux, ed.,
Heun’s Differential Equations (Ox-ford, New York, 1995), First ed.[101] S. Sinha and L. Santos, Phys.Rev. Lett. , 140406 (2007), URL https://link.aps.org/doi/10.1103/PhysRevLett.99.140406 .[102] M. C. Tsatsos and A. U. J. Lode, Journal of Low Tem-perature Physics , 171 (2015), ISSN 1573-7357, URL https://doi.org/10.1007/s10909-015-1335-5https://doi.org/10.1007/s10909-015-1335-5