seen in Fig. 2 (b), reproduces the biconcave shape of the condensate as does the numerical solution. Thus the method of coupled Gaussians is a viable and full-fledged alternative to direct numerical solutions of the Gross- Pitaevskii equation for dipolar condensates. correct result, the inclusion of only one more Gaussian brings the chemical potential curve already close to the numerical result, and practically no improvement is vis- ible in Fig. 1 when 3 or more Gaussians are included. Using 5 coupled Gaussians the exact bifurcation point is reproduced with an accuracy of 10 − $^{6}$. Similar results are obtained 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-Gaussian biconcave “blood-cell-like” shape. To demonstrate the power of the coupled Gaussian wave packet method, we choose a set of such parameters. We consider an axisym- metric trap with (particle number scaled) trap frequen- cies N $^{2}$γ$_{z}$ = 25200 along the polarization direction of the dipoles and N $^{2}$γ$_{ρ}$ = 3600 in the plane perpendicular to it (corresponding to an aspect ratio of λ = γ$_{z}$/γ$_{ρ}$ = 7). For this set of parameters we show in Fig. 2 (a) the 3  56000  57000  58000  59000  60000  61000  62000 -0.04 -0.03 -0.02 -0.01  0  0.01 E mf 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   58700  58725  58750  58775 2 3 4 5 6 E mf Number of Gaussians n (a) v  0  20  40  60  80  100  120  0  0.02  0.04  0.06  0.08  0.1 = 0.0 = 0.008  = 0.016  (b) v n FIG. 3: (a) Mean field energy of a dipolar condensate for (particle number scaled) trap frequencies N $^{2}$γ$_{z}$ = 25200 and N $^{2}$γ$_{ρ}$ = 3600 as a function of the scattering length. In the variational calculation with one Gaussian a stable ground state (g $^{N=1}$) and an unstable excited state (e$^{N=1}$) emerge in a tangent bifurcation. Using coupled Gaussians two unstable states emerge (labeled u $^{coupled}$), of which the lower one turns into a stable ground state (g $^{coupled}$) in a pitchfork bifurca- tion. (b) Stability eigenvalues λ of the pitchfork bifurcation point for calculations with 6 coupled Gaussians, scattering length in rectangle marked in (a). Real and imaginary parts of two selected eigenvalues of the Jacobian (7) as a function of the scattering length. For a < a p cr = − 0 . 00359 the solu- tion is unstable with one pair of real eigenvalues. At a p cr the real eigenvalues vanish in a pitchfork bifurcation and a stable ground state forms with purely imaginary eigenvalues. Only those eigenvalues involved in the stability change are shown. 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 the numerical one (triangles) at different z coordinates. Both so- lutions show a biconcave shaped condensate. The figures are for (particle number scaled) trap frequencies N $^{2}$γ$_{z}$ = 25200 and N $^{2}$γ$_{ρ}$ = 3600, and scattering length a = 0. 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 6 Gaussians would be indistinguishable in the figure from those obtained using 5 Gaussians, and the results for 2– 4 Gaussians are not shown for the sake of clarity of the figure. 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 branches of solution are born in a tangent bifurcation. The en- convergence behavior of the mean field energy. We com- pare the variational solution as the number of Gaussian wave packets is increased from 2 to 6 with the mean field energy value of a numerical lattice calculation (imagi- nary time evolution combined with FFT) with a grid size of 128 × 512, at scattering length a = 0 as an example. The mean field energy for one Gaussian is E$_{mf}$ = 60361 E$_{d}$ and lies far outside the vertical energy scale. Evidently the numerical value is more than excel- lently reproduced by 5 and 6 coupled Gaussians. The behavior for other scattering lengths is similar. Also the wave function nicely converges, and moreover, as can be