4 below a p $_{cr}$, ( a p cr − ϵ ) < a < a p $_{cr}$, there exist two additional branches, besides the unstable solution. This proves that the bifurcation is of pitchfork type. Note that the clas- sification of the condensate as unstable for a < a p cr nev- ertheless remains true in physical terms due to the nu- merically small value of ϵ . We also note that for a > a p cr the phase portrait possesses only one elliptic fixed point cooresponding to the stable stationary ground state. ergetically higher branch has purely real stability eigen- values ± λ $^{R}$, corresponding to an unstable excited state e $^{N=1}$, the lower branch possesses purely imaginary eigen- values ± λ I and corresponds to the stable ground state g $^{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 the scattering length is decreased from positive values to- wards the tangent bifurcation, the branch corresponding to the ground state g coupled turns into an unstable state u coupled at a scattering length of a p cr = − 0 . 00359. This is evident from the stability analysis shown in Fig. 3 (b) where the stability eigenvalues for the ground state, cal- culated using 6 Gaussians, are plotted in a small interval of the scattering length around a p $_{cr}$. Above a p cr the eigen- values are purely imaginary, below they are purely real. [Note that in a Bogoliubov analysis this instability should appear as a dynamical instability.] The ground state re- mains unstable down to the tangent bifurcation point at a t cr = − 0 . 01224, where it joins the branch of the unstable excited state. The quality of the calculation using 5 Gaussian wave packets is also demonstrated in Fig. 3 (a) where the re- sults of a numerically grid calculation by imaginary time evolution are shown by crosses. Evidently the numerical results and the results obtained using 5 coupled Gaus- sians excellently agree. The imaginary time calculation, however, can only trace the stable branch of the solution and fails for the unstable branch. Thus it is demonstrated that the Gaussian wave packet method is not only numer- ically accurate but also capable of giving access to regions of the space of solutions of the Gross-Pitaevskii equation with dipolar interaction that are difficult to investigate by conventional numerical full quantum calculations. Is there a chance of observing dipolar BECs on the stable fork arms? The answer probably is no, in the same way as it is in the case of the question of observing the transition to structured ground states, possibly as- sociated with a roton instability, shortly before collapse. The reason is the difficulty of adjusting trap frequencies and the scattering length to the necessary precision in a real experiment. Nevertheless theoretical investigations of this type close to the threshold of instability of dipo- lar condensates are valuable in their own right since they help to understand the nature of the collapse, and thus of “what’s really going on”. The phenomenon of one smooth branch of solutions becoming unstable as a function of a control parame- ter is reminiscent of a pitchfork bifurcation. The two stable solutions on the fork arms which should also be born in a pitchfork bifurcation, and exist in a tiny neigh- borhood ( a p cr − ϵ ) < a < a p $_{cr}$, are numerically hard to trace 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 = − 0 . 036 slightly below a p $_{cr}$. Figure 4 shows con- tours of equal deviation of the mean field energy from that of the ground state in the plane spanned by the two eigenvectors whose eigenvalues are involved in the sta- bility change in Fig 3 (b). The coordinate axes δ$_{1}$, δ$_{2}$ correspond to small variations of the Gaussian parame- ters in the eigenvector directions around the hyperbolic fixed point solution located at the origin. The portrait clearly reveals the existence of two nearby elliptic fixed points corresponding to two additional stable solutions. Therefore, in a small interval ϵ of the scattering length 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 ( δ$_{1}$, δ$_{2}$ in arbitrary units). The figure shows a = − 0 . 0036 close below the pitchfork bi- furcation point, showing three fixed points: Two stable and one hyperbolic. -1 -0.5  0  0.5  1 $_{2}$  -1 -0.5  0  0.5  1 1 [1] A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and T. Pfau, Phys. Rev. Lett. 94 , 160401 (2005). [2] J. Stuhler, A. Griesmaier, T. Koch, M. Fattori, T. Pfau, S. Giovanazzi, P. Pedri, and L. Santos, Phys. Rev. Lett. 95 , 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 77 , 061601(R) (2008). [4] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau, Rep. Prog. Phys. 72 , 126401 (2009). [5] D. O’Dell, S. Giovanazzi, G. Kurizki, and V. M. Akulin, Phys. Rev. Lett. 84 , 5687 (2000).