Pitchfork bifurcations in blood-cell shaped dipolar Bose-Einstein condensates Stefan Rau, J¨org Main, Patrick K¨oberle, and Gu¨nter Wunner Institut fu¨r Theoretische Physik 1, Universita¨t Stuttgart, 70550 Stuttgart, Germany (Dated: January 8, 2010) WedemonstratethatthemethodofcoupledGaussianwavepacketsisafull-fledgedalternativeto directnumericalsolutions ofthe Gross-Pitaevskiiequationof condensateswith electromagnetically induced attractive 1/r interaction, or with dipole-dipole interaction. Moreover, Gaussian wave packets are superior in that they are capable of producing both stable and unstable stationary solutions,andthusofgivingaccesstoyetunexploredregionsofthespaceofsolutionsoftheGross- Pitaevskiiequation. Weapplythemethodtoclarifythetheoreticalnatureofthecollapsemechanism 0 of blood-cell shaped dipolar condensates: On the route to collapse the condensate passes through 1 a pitchfork bifurcation, where the ground state itself turns unstable, before it finally vanishes in a 0 tangent bifurcation. 2 PACSnumbers: 67.85.-d,03.75.Hh,05.30.Jp,05.45.-a n a J Bose-Einstein condensates (BECs) with dipole-dipole withϑ(cid:48) theanglebetweenr−r(cid:48) andtheaxisofanexter- 8 interactionhavebecomeanactiveandexcitingfieldofre- nalmagneticfield. Forcompletenesswewillalsoconsider ] searchbecausetheyofferthepossibilityoftuningtherel- thecaseofanisotropic“gravity-like”attractive1/rlong- s ativestrengthsoftheshort-rangeisotropiccontactinter- range interaction, a g action and the anisotropic long-range dipole interaction - by manipulating the s-wave scattering length via Fesh- (cid:90) |ψ(r(cid:48),t)|2 t V (r)=−2N d3r(cid:48) . (3) n bach resonances, and thus of studying a wealth of new lr |r−r(cid:48)| a phenomena that occur as one crosses the whole range u from dominance of the contact interaction to that of According to O’Dell et al. [5] this interaction could be q . the dipole interaction. The experimental realization of electromagnetically induced by exposing the condensate at a BEC of chromium atoms [1–3], which possess a strong atoms to an appropriately arranged set of triads of laser m magnetic dipole moment, has given additional impetus beams. The appealing feature of such “monopolar” - to the field (for a comprehensive list of references see condensates is that they can be self-trapping, i.e. ex- d the recent review by Lahaye et al. [4]). In the dilute ist without an external trapping potential. The equa- n limit,thetheoreticaldescriptionofthesecondensatescan tions above have been brought into dimensionless form o c be done in the framework of the Gross-Pitaevskii equa- by introducing natural units, which for monopolar in- [ tion (GPE). This nonlinear Schr¨odinger equation has teraction (Vmono = −u/r) are [6–8] the “Bohr radius” 1 been solved in the literature so far by simple variational au = (cid:126)2/(mu) for lengths, the “Rydberg energy” Eu = v ansatzes, where the mean-field energy is minimized, e.g., (cid:126)2/(2ma2u) for energies and ωu = Eu/(cid:126) for frequen- 3 withtheconjugategradientmethodorbyimaginarytime cies. Natural units for dipolar atoms with magnetic mo- 5 evolution. In this Letter we will show that the method ment µ are [9] the dipole length ad =µ0µm/(2π(cid:126)2), the 2 ofcoupledGaussianwavepacketsisanadequatealterna- dipole energy Ed =(cid:126)2/(2ma2d) and the dipole frequency .1 tivetosolvingtheGPEofBECswithlong-rangeinterac- ωd =Ed/(cid:126). The quantities γx,y,z in (1) denote the trap- 1 tions. Moreover,wewillshowthatthemethodissuperior ping frequencies in the three spatial directions measured 0 in that it also yields unstable stationary solutions, and in the respective frequency units, N is the number of 0 1 thus opens access to regions of the space of solutions of bosons,andathescatteringlengthinunitsofau andad, : theGPEunexploredheretofore. Asanapplicationofthe respectively. v method we will analyze in detail the theoretical nature The most obvious way of solving the Gross-Pitaevskii i X of the collapse mechanism of dipolar BECs. equation (1) is its direct numerical integration on multi- r The GPE for ultracold gases with long-range interac- dimensional grids using, e.g., fast-Fourier techniques. a tions, described by the interatomic potential V (r), has The stationary ground state can be obtained by imagi- lr the form nary time evolution. These calculations, however, may turn out laborious, and physical insight can often be idψ(r,t)=(cid:2)−∆+γ2x2+γ2y2+γ2z2 gainedusingapproximate,inparticular,variationalsolu- dt x y z tions. Acommonapproachemployedfordeterminingthe +8πNa|ψ(r,t)|2+Vlr(r)(cid:3)ψ(r,t), (1) dynamics and stability of condensates both with contact interaction only [10, 11] and with additional long-range where for dipolar interaction we have interaction [8] is to assume a simple Gaussian form of (cid:90) 1−3cos2ϑ(cid:48) the wave function, with time-dependent width parame- V (r)=N d3r(cid:48) |ψ(r(cid:48),t)|2. (2) lr |r−r(cid:48)|3 ters and phase, and to investigate the dynamics of these 2 quantities. Fordipolarcondensatesimprovementsonthe It is important to note that in contrast to numerical simple Gaussian form were made by multiplying it by calculations with imaginary time evolution, which only second-order Hermite polynomials [12]. work for stable solutions, this procedure will produce As an alternative to numerical quantum simulations both stable and unstable solutions, and thus uncover on multidimensional grids we will extend the variational yet unexplored parts of the space of solutions of the calculations in such a way that numerically converged Gross-Pitaevskii equation. results are obtained with significantly reduced computa- tional effort compared to the exact quantum simulations Toanalyzethestabilityofthestationarysolutionsthe but with similar accuracy. The method is that of cou- dynamicalequations(5)aresplitintoreal(R)andimag- pled Gaussian wave packets. It was originally proposed inary (I) parts and linearized around the fixed points. by Heller [13, 14] to describe quantum dynamics in the The eigenvalues of the Jacobian matrix J at the fixed semiclassical regime, and was successfully applied to the point dynamics of molecules and atoms in external fields [15]. ∂(cid:0)a˙k,R,a˙k,I,γ˙k,R,γ˙k,I(cid:1) Theideaistochoosetrialwavefunctionswhicharesuper- J = α α , (7) positions of N different Gaussians centered at the origin ∂(cid:0)al,R,al,I,γl,R,γl,I(cid:1) β β N N withα,β =x,y,z; k,l=1,...,N,determinethestability ψ(r,t)=(cid:88)ei(akxx2+akyy2+akzz2+γk) ≡(cid:88)gk(ak,γk;r) propertiesofthesolution. Ifalleigenvaluesλ ofthesys- j k=1 k=1 tem are purely imaginary, the motion is confined to the (4) vicinity of the fixed point and quasi-periodic. If one real where both the width parameters ak and the scalars γk partorseveralrealpartsoftheeigenvaluesarenon-zero, are complex quantities, with the latter determining the small variations from the fixed point lead to exponential weight and the phase of the individual Gaussian. In- growth of the perturbation. serting the ansatz (4) into the time-dependent Gross- Pitaevskii equation and applying the time-dependent --00..44 variational principle where ||iΦ(t) − Hψ(t)||2 is mini- --00..88 mized by varying Φ, and afterwards Φ is set equal to Φ = ψ˙, yields a set of ordinary differential equations for µµ --11..22 N = 1 N = 4 the width parameters ak and the scalars γ (cf. [15]) --11..66 N = 2 N = 5 k N = 3 exact --22 1 a˙k =−4(ak)2− Vk ; β =x,y,z; (5a) --11..11 --11 --00..99 --00..88 β β 2 2,β NN22aa γ˙k =2i(ak+ak+ak)−vk. (5b) x y z 0 FIG. 1: Chemical potential µ for self-trapped condensates The quantities (vk,Vk) with k =1,...,N constitute the with attractive 1/r interaction as a function of the scaled 0 2 solution vector to the set of linear equations scattering length N2a obtained by using up to 5 Gaussian wave packets in comparison with the result of the exact nu- N N merical solution of the stationary Gross-Pitaevskii equation. (cid:88)(cid:10)gl|xmxnvk|gk(cid:11)+ 1(cid:88)(cid:10)gl|xmxnxVkx|gk(cid:11) Note that all forms yield a tangential bifurcation diagram, α α 0 2 α α 2 with a stable (upper) and an unstable (lower) branch. The k=1 k=1 inclusionofthreeGaussiansalreadywellreproducestheexact N =(cid:88)(cid:10)gl|xmxnV(x)|gk(cid:11) (6) numerical result. α α k=1 As a first application we demonstrate the efficiency of with l = 1,...,N; m+n = 0,2; and x = x, x = y, the coupled Gaussian wave packet method for conden- 1 2 x =z. Here,V(x)=V +V +V denotesthesumofthe sates with attractive 1/r long-range interaction, for the 3 c lr t contact, the long-range and of external trap potentials. case of self-trapping (γ = 0). Figure 1 shows the x,y,z The important and appealing point of this procedure is results for the chemical potential as a function of the that all necessary integrals with the trial wave functions scaled scattering length N2a for superpositions of 1 to gl,gk from (4) can be calculated analytically. 5 Gaussians in comparison with the results of the exact Stationary variational solutions to the extended Gross- numericalsolution. Itisevidentthatallformsreproduce Pitaevskii equation (1) are found by searching for the the bifurcation behavior discussed in [6–8]: at a critical fixed points of (5), i.e. solving a˙k = 0;γ˙k = 0 for each point two solutions of the Gross-Pitaevskii equation are k = 1,...,N via a 4N dimensional highly nonlinear root born in a tangent bifurcation, one stable (upper branch) search. Theresultingstationarywidthandweight/phase and one unstable (lower branch). The numerically ac- parameters can then be used to calculate the mean field curate bifurcation point lies at a ≈ −1.025147. It can cr energy E = (cid:10)Ψ|−∆+V + 1(V +V )|Ψ(cid:11) and the also be seen that, while the variational calculation with mf t 2 c lr chemical potential µ = (cid:104)Ψ|−∆+V +V +V |Ψ(cid:105). one Gaussian, with aN=1 = −1.178, still lies far off the t c lr cr 3 correct result, the inclusion of only one more Gaussian seen in Fig. 2 (b), reproduces the biconcave shape of brings the chemical potential curve already close to the the condensate as does the numerical solution. Thus the numerical result, and practically no improvement is vis- method of coupled Gaussians is a viable and full-fledged ible in Fig. 1 when 3 or more Gaussians are included. alternative to direct numerical solutions of the Gross- Using 5 coupled Gaussians the exact bifurcation point is Pitaevskii equation for dipolar condensates. reproducedwithanaccuracyof10−6. Similarresultsare obtained in the presence of a trapping potential. (a) 62000 We now turn to dipolar condensates. Previous studies [12, 16] have shown that in certain regions of the param- 61000 eter space dipolar condensates assume a non-Gaussian 60000 biconcave “blood-cell-like” shape. To demonstrate the power of the coupled Gaussian wave packet method, we mf 59000 chooseasetofsuchparameters. Weconsideranaxisym- E metric trap with (particle number scaled) trap frequen- 58000 gN = 1 cies N2γz = 25200 along the polarization direction of eN = 1 the dipoles and N2γρ =3600 in the plane perpendicular 57000 guccoouupplleedd to it (corresponding to an aspect ratio of λ = γ /γ = numerical z ρ 7). For this set of parameters we show in Fig. 2 (a) the 56000 -0.04 -0.03 -0.02 -0.01 0 0.01 a (a) 58775 (b) n 6000 Emf 5588772550 v λ 240000 000 RIme λλ -2000 -4000 58700 -6000 2 3 4 5 6 -0.005 -0.0045 -0.004 -0.0035 -0.003 -0.0025 -0.002 Number of Gaussians a (b) 120 1 8000 nv FIG. 3: (a) Mean field energy of a dipolar condensate for z = 0.0 (particle number scaled) trap frequencies N2γ =25200 and Ψ 60 z 40 zz == 00..000186 N2γρ = 3600 as a function of the scattering length. In the 20 variational calculation with one Gaussian a stable ground 0 state (gN=1) and an unstable excited state (eN=1) emerge in 0 0.02 0.04 0.06 0.08 0.1 ρ atangentbifurcation. UsingcoupledGaussianstwounstable statesemerge(labeleducoupled),ofwhichtheloweroneturns FIG. 2: (a) Convergence of the mean-field energy with in- into a stable ground state (gcoupled) in a pitchfork bifurca- creasing number of coupled Gaussian wave packets (squares) tion. (b) Stability eigenvalues λ of the pitchfork bifurcation and comparison with the value obtained by a lattice calcula- point for calculations with 6 coupled Gaussians, scattering tion with grid size 128×512 (dashed line), which lies ener- length in rectangle marked in (a). Real and imaginary parts geticallyhigherthantheexactconvergedvariationalsolution of two selected eigenvalues of the Jacobian (7) as a function (solid line). (b) Comparison of the variational wave func- of the scattering length. For a < apcr = −0.00359 the solu- tionfor6coupledGaussians(solidcurves)withvaluesofthe tion is unstable with one pair of real eigenvalues. At apcr the numericalone(triangles)atdifferentz coordinates. Bothso- realeigenvaluesvanishinapitchforkbifurcationandastable lutions show a biconcave shaped condensate. The figures are ground state forms with purely imaginary eigenvalues. Only for (particle number scaled) trap frequencies N2γ = 25200 those eigenvalues involved in the stability change are shown. z and N2γ =3600, and scattering length a=0. ρ Figure 3 (a) shows, for the same set of trap frequen- convergence behavior of the mean field energy. We com- cies, the results for the mean field energy of the conden- pare the variational solution as the number of Gaussian sate as a function of the scattering length a (in units a ) d wavepacketsisincreasedfrom2to6withthemeanfield for a wave function with one Gaussian, and for 5 cou- energy value of a numerical lattice calculation (imagi- pled Gaussian wave packets. Results obtained using 6 nary time evolution combined with FFT) with a grid Gaussians would be indistinguishable in the figure from size of 128 × 512, at scattering length a = 0 as an those obtained using 5 Gaussians, and the results for 2– example. The mean field energy for one Gaussian is 4 Gaussians are not shown for the sake of clarity of the E = 60361E and lies far outside the vertical energy figure. mf d scale. Evidently the numerical value is more than excel- Similar to the above findings for monopolar conden- lently reproduced by 5 and 6 coupled Gaussians. The sates, and as is known from previous variational calcula- behavior for other scattering lengths is similar. Also the tions[9]fordipolarcondensates, forN =1twobranches wave function nicely converges, and moreover, as can be of solution are born in a tangent bifurcation. The en- 4 ergetically higher branch has purely real stability eigen- below ap, (ap −(cid:15))<a<ap, there exist two additional cr cr cr values ±λR, corresponding to an unstable excited state branches,besidestheunstablesolution. Thisprovesthat eN=1,thelowerbranchpossessespurelyimaginaryeigen- the bifurcation is of pitchfork type. Note that the clas- values ±λI and corresponds to the stable ground state sification of the condensate as unstable for a < ap nev- cr gN=1. At the bifurcation point the branches of the sta- ertheless remains true in physical terms due to the nu- bility eigenvalues merge and vanish. merically small value of (cid:15). We also note that for a>ap cr The situation is different if the condensate wave func- the phase portrait possesses only one elliptic fixed point tion is described by more than one Gaussian. As the cooresponding to the stable stationary ground state. scattering length is decreased from positive values to- wards the tangent bifurcation, the branch corresponding to the ground state gcoupled turns into an unstable state 1 ucoupled at a scattering length of ap = −0.00359. This cr is evident from the stability analysis shown in Fig. 3 (b) 0.5 where the stability eigenvalues for the ground state, cal- culatedusing6Gaussians, areplottedinasmallinterval of the scattering length around ap. Above ap the eigen- δ 1 0 cr cr values are purely imaginary, below they are purely real. [NotethatinaBogoliubovanalysisthisinstabilityshould -0.5 appear as a dynamical instability.] The ground state re- mains unstable down to the tangent bifurcation point at at =−0.01224,whereitjoinsthebranchoftheunstable -1 cr -1 -0.5 0 0.5 1 excited state. δ 2 The quality of the calculation using 5 Gaussian wave packets is also demonstrated in Fig. 3 (a) where the re- FIG.4: Contourplotofthemeanfieldenergywiththeeigen- sults ofa numerically grid calculationbyimaginary time vectorscorrespondingtotheeigenvaluesofFig.3(b)lineariz- evolution are shown by crosses. Evidently the numerical ing the vicinity of the fixed point (δ1,δ2 in arbitrary units). The figure shows a = −0.0036 close below the pitchfork bi- results and the results obtained using 5 coupled Gaus- furcation point, showing three fixed points: Two stable and sians excellently agree. The imaginary time calculation, one hyperbolic. however, canonlytrace thestable branchof thesolution andfailsfortheunstablebranch. Thusitisdemonstrated thattheGaussianwavepacketmethodisnotonlynumer- Is there a chance of observing dipolar BECs on the icallyaccuratebutalsocapableofgivingaccesstoregions stable fork arms? The answer probably is no, in the ofthespaceofsolutionsoftheGross-Pitaevskiiequation same way as it is in the case of the question of observing with dipolar interaction that are difficult to investigate the transition to structured ground states, possibly as- by conventional numerical full quantum calculations. sociated with a roton instability, shortly before collapse. The phenomenon of one smooth branch of solutions The reason is the difficulty of adjusting trap frequencies becoming unstable as a function of a control parame- and the scattering length to the necessary precision in a ter is reminiscent of a pitchfork bifurcation. The two real experiment. Nevertheless theoretical investigations stable solutions on the fork arms which should also be of this type close to the threshold of instability of dipo- borninapitchforkbifurcation,andexistinatinyneigh- larcondensatesarevaluableintheirownrightsincethey borhood (ap − (cid:15)) < a < ap, are numerically hard to help to understand the nature of the collapse, and thus cr cr trace and therefore not plotted in the figure. Their ex- of “what’s really going on”. istence, and the pitchfork type of the bifurcation, how- ever, can be proven by looking at the “phase portrait” plotted in Fig. 4 at a value of the scattering length a = −0.036 slightly below ap. Figure 4 shows con- cr [1] A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and tours of equal deviation of the mean field energy from T. Pfau, Phys. Rev. Lett. 94, 160401 (2005). thatofthegroundstateintheplanespannedbythetwo [2] J. Stuhler, A. Griesmaier, T. Koch, M. Fattori, T. Pfau, eigenvectors whose eigenvalues are involved in the sta- S. Giovanazzi, P. Pedri, and L. Santos, Phys. Rev. Lett. bility change in Fig 3 (b). The coordinate axes δ ,δ 95, 150406 (2005). 1 2 correspond to small variations of the Gaussian parame- [3] Q.Beaufils,R.Chicireanu,T.Zanon,B.Laburthe-Tolra, E. Mar´echal, L. Vernac, J.-C. Keller, and O. Gorceix, ters in the eigenvector directions around the hyperbolic Phys. Rev. A 77, 061601(R) (2008). fixed point solution located at the origin. The portrait [4] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and clearly reveals the existence of two nearby elliptic fixed T. Pfau, Rep. Prog. Phys. 72, 126401 (2009). points corresponding to two additional stable solutions. [5] D. O’Dell, S. Giovanazzi, G. Kurizki, and V. M. Akulin, Therefore, in a small interval (cid:15) of the scattering length Phys. Rev. Lett. 84, 5687 (2000). 5 [6] I. Papadopoulos, P. Wagner, G. Wunner, and J. Main, stein, and P. Zoller, Phys. Rev. A 56, 1424 (1997). Phys. Rev. A 76, 053604 (2007). [12] S.Ronen,D.C.E.Bortolotti,andJ.L.Bohn,Phys.Rev. [7] H.Cartarius,J.Main,andG.Wunner,Phys.Rev.A77, Lett. 98, 030406 (2007). 013618 (2008). [13] E. J. Heller, J. Chem. Phys. 65, 4979 (1976). [8] H. Cartarius, T. Fabˇciˇc, J. Main, and G. Wunner, Phys. [14] E. J. Heller, J. Chem. Phys. 75, 2923 (1981). Rev. A 78, 013615 (2008). [15] T. Fabˇciˇc, J. Main, and G. Wunner, Phys. Rev. A 79, [9] P. Ko¨berle, H. Cartarius, T. Fabˇciˇc, J. Main, and 043416 and 043417 (2009). G. Wunner, New Journal of Physics 11, 023017 (2009). [16] O. Dutta and P. Meystre, Phys. Rev. A 75, 053604 [10] V. M. P´erez-Garc´ıa, H. Michinel, J. I. Cirac, M. Lewen- (2007). stein, and P. Zoller, Phys. Rev. Lett. 77, 5320 (1996). [11] V. M. P´erez-Garc´ıa, H. Michinel, J. I. Cirac, M. Lewen-