2 It is important to note that in contrast to numerical calculations with imaginary time evolution, which only work for stable solutions, this procedure will produce both stable and unstable solutions, and thus uncover yet unexplored parts of the space of solutions of the Gross-Pitaevskii equation. quantities. For dipolar condensates improvements on the simple Gaussian form were made by multiplying it by second-order Hermite polynomials [12]. As an alternative to numerical quantum simulations on multidimensional grids we will extend the variational calculations in such a way that numerically converged results are obtained with significantly reduced computa- tional effort compared to the exact quantum simulations but with similar accuracy. The method is that of cou- pled Gaussian wave packets . It was originally proposed by Heller [13, 14] to describe quantum dynamics in the semiclassical regime, and was successfully applied to the dynamics 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 To analyze the stability of the stationary solutions the dynamical 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 fixed point with α, β = x, y, z ; k, l = 1 , ..., N , determine the stability properties of the solution. If all eigenvalues λ$_{j}$ of the sys- tem are purely imaginary, the motion is confined to the vicinity of the fixed point and quasi-periodic. If one real part or several real parts of the eigenvalues are non-zero, small variations from the fixed point lead to exponential growth of the perturbation. where both the width parameters a k and the scalars γ k are complex quantities, with the latter determining the weight and the phase of the individual Gaussian. In- serting the ansatz (4) into the time-dependent Gross- Pitaevskii equation and applying the time-dependent variational principle where || i Φ( t ) − Hψ ( t ) || 2 is mini- mized by varying Φ, and afterwards Φ is set equal to Φ = ˙ ψ , yields a set of ordinary differential equations for the width parameters a k and the scalars γ$_{k}$ (cf. [15]) The quantities ( v k $_{0}$, V k 2 ) with k = 1 , . . . , N constitute the solution vector to the set of linear equations As a first application we demonstrate the efficiency of the coupled Gaussian wave packet method for conden- sates with attractive 1 /r long-range interaction, for the case of self-trapping ( γ$_{x,y,z}$ = 0). Figure 1 shows the results for the chemical potential as a function of the scaled scattering length N $^{2}$a for superpositions of 1 to 5 Gaussians in comparison with the results of the exact numerical solution. It is evident that all forms reproduce the bifurcation behavior discussed in [6–8]: at a critical point two solutions of the Gross-Pitaevskii equation are born in a tangent bifurcation, one stable (upper branch) and one unstable (lower branch). The numerically ac- curate bifurcation point lies at a$_{cr}$ ≈ − 1 . 025147. It can also be seen that, while the variational calculation with one Gaussian, with a N=1 cr = − 1 . 178, still lies far off the with l = 1 , . . . , N ; m + n = 0 , 2; and x$_{1}$ = x , x$_{2}$ = y , x$_{3}$ = z . Here, V ( x ) = V$_{c}$ + V$_{lr}$ + V$_{t}$ denotes the sum of the contact, the long-range and of external trap potentials. The important and appealing point of this procedure is that all necessary integrals with the trial wave functions g $^{l}$, gk from (4) can be calculated analytically. Stationary variational solutions to the extended Gross- Pitaevskii equation (1) are found by searching for the fixed points of (5), i.e. solving ˙ a k = 0; ˙ γ k = 0 for each k = 1 , ..., N via a 4 N dimensional highly nonlinear root search. The resulting stationary width and weight/phase parameters can then be used to calculate the mean field energy E$_{mf}$ = $^{〈}$Ψ | − Δ + V$_{t}$ + 1 $_{2}$( V$_{c}$ + V$_{lr}$ ) | Ψ 〉 and the chemical potential µ = 〈 Ψ | − Δ + V$_{t}$ + V$_{c}$ + V$_{lr}$ | Ψ 〉 . N ∑ k =1 $^{〈}$g $^{l}$| x m α x n $_{α}$v k $_{0}$| g k $^{〉}$+ 1 2 N ∑ k =1 $^{〈}$g $^{l}$| x m α x n $_{α}$x V k 2 x | g k 〉 = N ∑ k =1 $^{〈}$g $^{l}$| x m α x n $_{α}$V ( x ) | g k 〉 (6) ˙ a k β = − 4( a k $_{β}$) 2 − 1 2 V k 2 $_{,β}$; β = x, y, z ; (5a) ˙ γ k = 2 i ( a k x + a k y + a k $_{z}$) − v k $_{0}$. (5b) ψ ( r , t ) = N ∑ k =1 e i ( a k $_{x}$x $^{2}$+ a k $_{y}$y $^{2}$+ a k $_{z}$z $^{2}$+ γ $^{k}$) ≡ N ∑ k =1 g $^{k}$( a $^{k}$, γ$^{k}$; r ) (4) J = ∂ $^{(}$˙ a k,R α , ˙ a k,I α , ˙ γ $^{k,R}$, ˙ γ k,I ) ∂ $^{(}$a l,R β , a l,I β , γ $^{l,R}$, γl,I ) , (7) -2 -1.6 -1.2 -0.8 -0.4 -1.1 -1 -0.9 -0.8 N$^{2}$a N = 1 N = 2 N = 3 N = 4 N = 5 exact FIG. 1: Chemical potential µ for self-trapped condensates with attractive 1 /r interaction as a function of the scaled scattering length N $^{2}$a obtained by using up to 5 Gaussian wave 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. The inclusion of three Gaussians already well reproduces the exact numerical result.