Pitchfork bifurcations in blood-cell shaped dipolar Bose-Einstein condensates
PPitchfork bifurcations in blood-cell shaped dipolar Bose-Einstein condensates
Stefan Rau, J¨org Main, Patrick K¨oberle, and G¨unter Wunner
Institut f¨ur Theoretische Physik 1, Universit¨at Stuttgart, 70550 Stuttgart, Germany (Dated: November 15, 2018)We demonstrate that the method of coupled Gaussian wave packets is a full-fledged alternative todirect numerical solutions of the Gross-Pitaevskii equation of condensates with electromagneticallyinduced attractive 1 /r interaction, or with dipole-dipole interaction. Moreover, Gaussian wavepackets are superior in that they are capable of producing both stable and unstable stationarysolutions, and thus of giving access to yet unexplored regions of the space of solutions of the Gross-Pitaevskii equation. We apply the method to clarify the theoretical nature of the collapse mechanismof blood-cell shaped dipolar condensates: On the route to collapse the condensate passes througha pitchfork bifurcation, where the ground state itself turns unstable, before it finally vanishes in atangent bifurcation. PACS numbers: 67.85.-d, 03.75.Hh, 05.30.Jp, 05.45.-a
Bose-Einstein condensates (BECs) with dipole-dipoleinteraction have become an active and exciting field of re-search because they offer the possibility of tuning the rel-ative strengths of the short-range isotropic contact inter-action and the anisotropic long-range dipole interactionby manipulating the s -wave scattering length via Fesh-bach resonances, and thus of studying a wealth of newphenomena that occur as one crosses the whole rangefrom dominance of the contact interaction to that ofthe dipole interaction. The experimental realization ofa BEC of chromium atoms [1–3], which possess a strongmagnetic dipole moment, has given additional impetusto the field (for a comprehensive list of references seethe recent review by Lahaye et al. [4]). In the dilutelimit, the theoretical description of these condensates canbe done in the framework of the Gross-Pitaevskii equa-tion (GPE). This nonlinear Schr¨odinger equation hasbeen solved in the literature so far by simple variationalansatzes, where the mean-field energy is minimized, e.g.,with the conjugate gradient method or by imaginary timeevolution. In this Letter we will show that the methodof coupled Gaussian wave packets is an adequate alterna-tive to solving the GPE of BECs with long-range interac-tions. Moreover, we will show that the method is superiorin that it also yields unstable stationary solutions, andthus opens access to regions of the space of solutions ofthe GPE unexplored heretofore. As an application of themethod we will analyze in detail the theoretical natureof the collapse mechanism of dipolar BECs.The GPE for ultracold gases with long-range interac-tions, described by the interatomic potential V lr ( r ), hasthe form i dd t ψ ( r , t ) = (cid:2) − ∆ + γ x x + γ y y + γ z z +8 πN a | ψ ( r , t ) | + V lr ( r ) (cid:3) ψ ( r , t ) , (1)where for dipolar interaction we have V lr ( r ) = N (cid:90) d r (cid:48) − ϑ (cid:48) | r − r (cid:48) | | ψ ( r (cid:48) , t ) | . (2) with ϑ (cid:48) the angle between r − r (cid:48) and the axis of an exter-nal magnetic field. For completeness we will also considerthe case of an isotropic “gravity-like” attractive 1 /r long-range interaction, V lr ( r ) = − N (cid:90) d r (cid:48) | ψ ( r (cid:48) , t ) | | r − r (cid:48) | . (3)According to O’Dell et al. [5] this interaction could beelectromagnetically induced by exposing the condensateatoms to an appropriately arranged set of triads of laserbeams. The appealing feature of such “monopolar”condensates is that they can be self-trapping, i.e. ex-ist without an external trapping potential. The equa-tions above have been brought into dimensionless formby introducing natural units, which for monopolar in-teraction ( V mono = − u/r ) are [6–8] the “Bohr radius” a u = (cid:126) / ( mu ) for lengths, the “Rydberg energy” E u = (cid:126) / (2 ma u ) for energies and ω u = E u / (cid:126) for frequen-cies. Natural units for dipolar atoms with magnetic mo-ment µ are [9] the dipole length a d = µ µm/ (2 π (cid:126) ), thedipole energy E d = (cid:126) / (2 ma ) and the dipole frequency ω d = E d / (cid:126) . The quantities γ x,y,z in (1) denote the trap-ping frequencies in the three spatial directions measuredin the respective frequency units, N is the number ofbosons, and a the scattering length in units of a u and a d ,respectively.The most obvious way of solving the Gross-Pitaevskiiequation (1) is its direct numerical integration on multi-dimensional grids using, e.g., fast-Fourier techniques.The stationary ground state can be obtained by imagi-nary time evolution. These calculations, however, mayturn out laborious, and physical insight can often begained using approximate, in particular, variational solu-tions. A common approach employed for determining thedynamics and stability of condensates both with contactinteraction only [10, 11] and with additional long-rangeinteraction [8] is to assume a simple Gaussian form ofthe wave function, with time-dependent width parame-ters and phase, and to investigate the dynamics of these a r X i v : . [ c ond - m a t . qu a n t - g a s ] J a n quantities. For dipolar condensates improvements on thesimple Gaussian form were made by multiplying it bysecond-order Hermite polynomials [12].As an alternative to numerical quantum simulationson multidimensional grids we will extend the variationalcalculations in such a way that numerically convergedresults are obtained with significantly reduced computa-tional effort compared to the exact quantum simulationsbut with similar accuracy. The method is that of cou-pled Gaussian wave packets . It was originally proposedby Heller [13, 14] to describe quantum dynamics in thesemiclassical regime, and was successfully applied to thedynamics of molecules and atoms in external fields [15].The idea is to choose trial wave functions which are super-positions of N different Gaussians centered at the origin ψ ( r , t ) = N (cid:88) k =1 e i ( a kx x + a ky y + a kz z + γ k ) ≡ N (cid:88) k =1 g k ( a k , γ k ; r )(4)where both the width parameters a k and the scalars γ k are complex quantities, with the latter determining theweight and the phase of the individual Gaussian. In-serting the ansatz (4) into the time-dependent Gross-Pitaevskii equation and applying the time-dependentvariational principle where || i Φ( t ) − Hψ ( t ) || is mini-mized by varying Φ, and afterwards Φ is set equal toΦ = ˙ ψ , yields a set of ordinary differential equations forthe width parameters a k and the scalars γ k (cf. [15])˙ a kβ = − a kβ ) − V k ,β ; β = x, y, z ; (5a)˙ γ k = 2 i ( a kx + a ky + a kz ) − v k . (5b)The quantities ( v k , V k ) with k = 1 , . . . , N constitute thesolution vector to the set of linear equations N (cid:88) k =1 (cid:10) g l | x mα x nα v k | g k (cid:11) + 12 N (cid:88) k =1 (cid:10) g l | x mα x nα x V k x | g k (cid:11) = N (cid:88) k =1 (cid:10) g l | x mα x nα V ( x ) | g k (cid:11) (6)with l = 1 , . . . , N ; m + n = 0 ,
2; and x = x , x = y , x = z . Here, V ( x ) = V c + V lr + V t denotes the sum of thecontact, the long-range and of external trap potentials.The important and appealing point of this procedure isthat all necessary integrals with the trial wave functions g l , g k from (4) can be calculated analytically.Stationary variational solutions to the extended Gross-Pitaevskii equation (1) are found by searching for thefixed points of (5), i.e. solving ˙ a k = 0; ˙ γ k = 0 for each k = 1 , ..., N via a 4 N dimensional highly nonlinear rootsearch. The resulting stationary width and weight/phaseparameters can then be used to calculate the mean fieldenergy E mf = (cid:10) Ψ | − ∆ + V t + ( V c + V lr ) | Ψ (cid:11) and thechemical potential µ = (cid:104) Ψ | − ∆ + V t + V c + V lr | Ψ (cid:105) . It is important to note that in contrast to numericalcalculations with imaginary time evolution, which onlywork for stable solutions, this procedure will produceboth stable and unstable solutions, and thus uncoveryet unexplored parts of the space of solutions of theGross-Pitaevskii equation.To analyze the stability of the stationary solutions thedynamical equations (5) are split into real ( R ) and imag-inary ( I ) parts and linearized around the fixed points.The eigenvalues of the Jacobian matrix J at the fixedpoint J = ∂ (cid:0) ˙ a k,Rα , ˙ a k,Iα , ˙ γ k,R , ˙ γ k,I (cid:1) ∂ (cid:0) a l,Rβ , a l,Iβ , γ l,R , γ l,I (cid:1) , (7)with α, β = x, y, z ; k, l = 1 , ..., N , determine the stabilityproperties of the solution. If all eigenvalues λ j of the sys-tem are purely imaginary, the motion is confined to thevicinity of the fixed point and quasi-periodic. If one realpart or several real parts of the eigenvalues are non-zero,small variations from the fixed point lead to exponentialgrowth of the perturbation. -2-1.6-1.2-0.8-0.4 -1.1 -1 -0.9 -0.8 µ N aN = 1N = 2N = 3-2-1.6-1.2-0.8-0.4 -1.1 -1 -0.9 -0.8 µ N aN = 4N = 5exact FIG. 1: Chemical potential µ for self-trapped condensateswith attractive 1 /r interaction as a function of the scaledscattering length N a obtained by using up to 5 Gaussianwave packets in comparison with the result of the exact nu-merical solution of the stationary Gross-Pitaevskii equation.Note that all forms yield a tangential bifurcation diagram,with a stable (upper) and an unstable (lower) branch. Theinclusion of three Gaussians already well reproduces the exactnumerical result. As a first application we demonstrate the efficiency ofthe coupled Gaussian wave packet method for conden-sates with attractive 1 /r long-range interaction, for thecase of self-trapping ( γ x,y,z = 0). Figure 1 shows theresults for the chemical potential as a function of thescaled scattering length N a for superpositions of 1 to5 Gaussians in comparison with the results of the exactnumerical solution. It is evident that all forms reproducethe bifurcation behavior discussed in [6–8]: at a criticalpoint two solutions of the Gross-Pitaevskii equation areborn in a tangent bifurcation, one stable (upper branch)and one unstable (lower branch). The numerically ac-curate bifurcation point lies at a cr ≈ − . a N=1cr = − . − . Similar results areobtained in the presence of a trapping potential.We now turn to dipolar condensates. Previous studies[12, 16] have shown that in certain regions of the param-eter space dipolar condensates assume a non-Gaussianbiconcave “blood-cell-like” shape. To demonstrate thepower of the coupled Gaussian wave packet method, wechoose a set of such parameters. We consider an axisym-metric trap with (particle number scaled) trap frequen-cies N γ z = 25200 along the polarization direction ofthe dipoles and N γ ρ = 3600 in the plane perpendicularto it (corresponding to an aspect ratio of λ = γ z /γ ρ =7). For this set of parameters we show in Fig. 2 (a) the E m f Number of Gaussians n(a) v 0 20 40 60 80 100 120 0 0.02 0.04 0.06 0.08 0.1 Ψ ρ z = 0.0z = 0.008 z = 0.016 (b) vn FIG. 2: (a) Convergence of the mean-field energy with in-creasing number of coupled Gaussian wave packets (squares)and comparison with the value obtained by a lattice calcula-tion with grid size 128 ×
512 (dashed line), which lies ener-getically higher than the exact converged variational solution(solid line). (b) Comparison of the variational wave func-tion for 6 coupled Gaussians (solid curves) with values of thenumerical one (triangles) at different z coordinates. Both so-lutions show a biconcave shaped condensate. The figures arefor (particle number scaled) trap frequencies N γ z = 25200and N γ ρ = 3600, and scattering length a = 0. convergence behavior of the mean field energy. We com-pare the variational solution as the number of Gaussianwave packets is increased from 2 to 6 with the mean fieldenergy value of a numerical lattice calculation (imagi-nary time evolution combined with FFT) with a gridsize of 128 × a = 0 as anexample. The mean field energy for one Gaussian is E mf = 60361 E d and lies far outside the vertical energyscale. Evidently the numerical value is more than excel-lently reproduced by 5 and 6 coupled Gaussians. Thebehavior for other scattering lengths is similar. Also thewave function nicely converges, and moreover, as can be seen in Fig. 2 (b), reproduces the biconcave shape ofthe condensate as does the numerical solution. Thus themethod of coupled Gaussians is a viable and full-fledgedalternative to direct numerical solutions of the Gross-Pitaevskii equation for dipolar condensates. E m f a(a) g N = 1 e N = 1 g coupled u coupled numerical-6000-4000-2000 0 2000 4000 6000-0.005 -0.0045 -0.004 -0.0035 -0.003 -0.0025 -0.002 λ a(b) Re λ Im λ FIG. 3: (a) Mean field energy of a dipolar condensate for(particle number scaled) trap frequencies N γ z = 25200 and N γ ρ = 3600 as a function of the scattering length. In thevariational calculation with one Gaussian a stable groundstate (g N=1 ) and an unstable excited state (e
N=1 ) emerge ina tangent bifurcation. Using coupled Gaussians two unstablestates emerge (labeled u coupled ), of which the lower one turnsinto a stable ground state (g coupled ) in a pitchfork bifurca-tion. (b) Stability eigenvalues λ of the pitchfork bifurcationpoint for calculations with 6 coupled Gaussians, scatteringlength in rectangle marked in (a). Real and imaginary partsof two selected eigenvalues of the Jacobian (7) as a functionof the scattering length. For a < a pcr = − . a pcr thereal eigenvalues vanish in a pitchfork bifurcation and a stableground state forms with purely imaginary eigenvalues. Onlythose eigenvalues involved in the stability change are shown. Figure 3 (a) shows, for the same set of trap frequen-cies, the results for the mean field energy of the conden-sate as a function of the scattering length a (in units a d )for a wave function with one Gaussian, and for 5 cou-pled Gaussian wave packets. Results obtained using 6Gaussians would be indistinguishable in the figure fromthose obtained using 5 Gaussians, and the results for 2–4 Gaussians are not shown for the sake of clarity of thefigure.Similar to the above findings for monopolar conden-sates, and as is known from previous variational calcula-tions [9] for dipolar condensates, for N = 1 two branchesof solution are born in a tangent bifurcation. The en-ergetically higher branch has purely real stability eigen-values ± λ R , corresponding to an unstable excited statee N=1 , the lower branch possesses purely imaginary eigen-values ± λ I and corresponds to the stable ground stateg N=1 . At the bifurcation point the branches of the sta-bility eigenvalues merge and vanish.The situation is different if the condensate wave func-tion is described by more than one Gaussian. As thescattering length is decreased from positive values to-wards the tangent bifurcation, the branch correspondingto the ground state g coupled turns into an unstable stateu coupled at a scattering length of a pcr = − . a pcr . Above a pcr the eigen-values are purely imaginary, below they are purely real.[Note that in a Bogoliubov analysis this instability shouldappear as a dynamical instability.] The ground state re-mains unstable down to the tangent bifurcation point at a tcr = − . pitchfork bifurcation. The twostable solutions on the fork arms which should also beborn in a pitchfork bifurcation, and exist in a tiny neigh-borhood ( a pcr − (cid:15) ) < a < a pcr , are numerically hard totrace and therefore not plotted in the figure. Their ex-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 = − .
036 slightly below a pcr . Figure 4 shows con-tours of equal deviation of the mean field energy fromthat of the ground state in the plane spanned by the twoeigenvectors whose eigenvalues are involved in the sta-bility change in Fig 3 (b). The coordinate axes δ , δ correspond to small variations of the Gaussian parame-ters in the eigenvector directions around the hyperbolicfixed point solution located at the origin. The portraitclearly reveals the existence of two nearby elliptic fixedpoints corresponding to two additional stable solutions.Therefore, in a small interval (cid:15) of the scattering length below a pcr , ( a pcr − (cid:15) ) < a < a pcr , there exist two additionalbranches, besides the unstable solution. This proves thatthe bifurcation is of pitchfork type. Note that the clas-sification of the condensate as unstable for a < a pcr nev-ertheless remains true in physical terms due to the nu-merically small value of (cid:15) . We also note that for a > a pcr the phase portrait possesses only one elliptic fixed pointcooresponding to the stable stationary ground state. -1 -0.5 0 0.5 1 δ -1-0.5 0 0.5 1 δ FIG. 4: Contour plot of the mean field energy with the eigen-vectors corresponding to the eigenvalues of Fig. 3 (b) lineariz-ing the vicinity of the fixed point ( δ , δ in arbitrary units).The figure shows a = − . Is there a chance of observing dipolar BECs on thestable fork arms? The answer probably is no, in thesame way as it is in the case of the question of observingthe transition to structured ground states, possibly as-sociated with a roton instability, shortly before collapse.The reason is the difficulty of adjusting trap frequenciesand the scattering length to the necessary precision in areal experiment. Nevertheless theoretical investigationsof this type close to the threshold of instability of dipo-lar condensates are valuable in their own right since theyhelp to understand the nature of the collapse, and thusof “what’s really going on”. [1] A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, andT. Pfau, Phys. Rev. Lett. , 160401 (2005).[2] J. Stuhler, A. Griesmaier, T. Koch, M. Fattori, T. Pfau,S. Giovanazzi, P. Pedri, and L. Santos, Phys. Rev. Lett. , 150406 (2005).[3] Q. Beaufils, R. Chicireanu, T. Zanon, B. Laburthe-Tolra,E. Mar´echal, L. Vernac, J.-C. Keller, and O. Gorceix,Phys. Rev. A , 061601(R) (2008).[4] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, andT. Pfau, Rep. Prog. Phys. , 126401 (2009).[5] D. O’Dell, S. Giovanazzi, G. Kurizki, and V. M. Akulin,Phys. Rev. Lett. , 5687 (2000). [6] I. Papadopoulos, P. Wagner, G. Wunner, and J. Main,Phys. Rev. A , 053604 (2007).[7] H. Cartarius, J. Main, and G. Wunner, Phys. Rev. A ,013618 (2008).[8] H. Cartarius, T. Fabˇciˇc, J. Main, and G. Wunner, Phys.Rev. A , 013615 (2008).[9] P. K¨oberle, H. Cartarius, T. Fabˇciˇc, J. Main, andG. Wunner, New Journal of Physics , 023017 (2009).[10] V. M. P´erez-Garc´ıa, H. Michinel, J. I. Cirac, M. Lewen-stein, and P. Zoller, Phys. Rev. Lett. , 5320 (1996).[11] V. M. P´erez-Garc´ıa, H. Michinel, J. I. Cirac, M. Lewen- stein, and P. Zoller, Phys. Rev. A , 1424 (1997).[12] S. Ronen, D. C. E. Bortolotti, and J. L. Bohn, Phys. Rev.Lett. , 030406 (2007).[13] E. J. Heller, J. Chem. Phys. , 4979 (1976).[14] E. J. Heller, J. Chem. Phys. , 2923 (1981).[15] T. Fabˇciˇc, J. Main, and G. Wunner, Phys. Rev. A ,043416 and 043417 (2009).[16] O. Dutta and P. Meystre, Phys. Rev. A75