Time-dependent variational approach for Bose–Einstein condensates with nonlocal interaction

The nonlinear dynamics of self-bound Bose–Einstein matter waves under the action of an attractive nonlocal and a repulsive local interaction is analyzed by means of a time-dependent variational formalism. The mean-field model described by the Gross–Pitaevskii equation (GPE) is reduced to a single second-order conservative ordinary differential equation admitting a potential function. The stable fixed points and the linear oscillation frequencies of the breather modes are determined. The chosen nonlocal interactions are a Gaussian-shaped potential and a van der Waals-like potential. The variational solutions are compared with direct numerical simulations of the GPE, for each of these nonlocalities. The spatio-temporal frequency spectra of linear waves are also determined.


Introduction
In many situations, nonlocal spatial interactions play a central role. In particular, the experimental observations of dipolar Bose-Einstein condensates (BECs) in systems composed of atoms with a large magnetic moment provided a breakthrough in ultra-cold gas physics [1][2][3]. The creation of self-bound quantum matter waves, or quantum balls [4], has become a concrete possibility to be realized in the near future. Proposals include the self-trapping of mesoscopic atomic clouds by a collective excitation of Rydberg atom pairs [5], additional interactions in dilute BECs of alkali atoms besides the dipolar interaction, creating 2D BEC solitons [6], the role of the orbital angular momentum on the arrest of collapse of vortexfree elliptical clouds of BECs described by Ermakov systems [7,8], and the manipulation of the atomic scattering length through Feshbach resonance to produce stable solitons [9]. In a more general sense, nonlocality is present e.g.in nematic liquid crystals with long-range molecular forces [10], plasmas with a nonlocal response due to heating and ionization processes [11], the arrested collapse of ultra-cold magnetic gases due to quantum fluctuations [12], the crossover from a dilute BEC to a macrodroplet in a dipolar quantum fluid [13], the Rosensweig instability of quantum ferrofluids [14], ballistic transport of atoms [15], or the propagation of light beams in anisotropic nonlinear media [16]. Nonlocality is a decisive influence on the modulational instability of Kerr media [17], as well on the soliton stabilization and collapse arrest in nonlocal nonlinear systems [18].
In the present work, we analyze the nonlinear, timedependent dynamics in a BEC with local and nonlocal interactions by means of a variational treatment, without external confinement. Previous studies in the literature on long-range interaction BECs include, for instance, descriptions restricted to the linear regime [6], specific interactions such as between dark solitons in dipolar BECs [19] by means of the direct simulation of the mean-field equation [20], including an external trapping potential instead of a local coupling [21,22], higher-order dissipative nonlinearities [23], application of the Thomas-Fermi approximation [24], focusing on stationary solutions [25], among other distinctive features. Our time-dependent variational formalism may also be adapted for bosonic systems in a cavity, including experiment [26,27] and theory for cavity back-action [28] and measurement back-action [29].
The dipole-dipole coupling is anisotropic and longrange, making analytical progress more difficult. Therefore, we follow previous works [5,21,[30][31][32][33] proposing a Gaussian kernel, which is particularly well suited for detailed calculations. In addition, a nonlocal coupling reminiscent of the van der Waals potential [5,24,34] is also studied. The specific kernel does not qualitatively change our results, as seen in the comparison between the dynamics for Gaussian and van der Waals-like potentials (section 4). Other possible nonlocal interactions include a Lorentzian shape [22] and singular 1/r n couplings, with n1 [25,35,36].
This work is organized as follows. In section 2, we introduce our basic mean-field model in terms of a nonlocal nonlinear Schrödinger equation and the associated variational formalism, with a Gaussian kernel. In section 3, a suitable time-dependent Ansatz is proposed for the wave function, which reduces the problem to a nonlinear second-order ordinary equation. Being conservative, the resulting dynamical system is associated with a potential function, which allows a sensible interpretation of the linear stability analysis results. Section 4 provides a numerical stability analysis of the variational solutions, by means of direct simulation of the original mean-field, nonlocal nonlinear Schrödinger equation. Section 5 discuss the spectrum of linear excitations. Finally, section 6 contains our conclusions.

Variational approach
The Gross-Pitaevskii equation (GPE) with a nonlocal interaction [37] term can be written as where ÿ is the reduced Planck constant, M is the boson mass and is the coupling constant in terms of the boson-boson s-wave scattering length a. For simplicity we omit writing the explicit time-dependence of the wave function ψ(r, t) in equation (1) and in some equations below. We first consider a Gaussian kernel specified by where the constants W 0 and R are, respectively, the strength and range of the nonlocal interaction. Henceforth, we assume g0 and W 0 <0, so that the local (contact) and nonlocal couplings are respectively repulsive and attractive. Although the sign and strength of g could be tuned by Fesbach's resonance, we exclude attractive contact interactions (g<0) since these are always unstable against collapse in the absence of an external confinement [1,38]. See also the comments below equations (9) and (13) in the following. Notice that our mean-field treatment is valid for dilute condensates where the nonlocal GPE correctly describes the ground state behavior, and so that strongly collisional regimes are avoided. Although the Gaussian kernel is not directly related to physical potentials, it is frequently used [5,21,[30][31][32][33] as an approximate model for nonlocal nonlinearities, since it gives rise to tractable analytical expressions.
The following normalizations to dimensionless variables will henceforth be applied: using for simplicity the same symbols for the original and dimensionless quantities, together with the normalization y y ) Na R 2 1 2 . Some numerical factors were calibrated so that α=γ represents a system where the local and nonlocal influences have the same order of magnitude, as found from inspection of equation (4) taking into account  The functional derivative * d dy = S 0 gives equation (4), while δS/δψ=0 provides the complex conjugate of the GPE.

Radially symmetric states
A reasonable radially symmetric Ansatz is given by the Gaussian form which is properly normalized, where σ=σ(t) (the width) and β=β(t) are time-dependent variational functions to be determined. The function β is known as the chirp function [39] and is necessary to obtain meaningful differential equations at the end of the time-dependent variational procedure. In the stationary case (∂/∂t=0) one could chose β≡0 from the beginning. In a quantum fluid model (see section 5), the chirp gives a radial contribution to the velocity field, defined via the gradient of the complex phase (action) of the wave function.
Using equation (6) in (5) and integrating over space In passing, the integral of the nonlocal term was readily evaluated using the convolution theorem combined with the spatial Fourier transform (appendix A).
obtained after eliminating β. The pseudo-potential is The contributions can be clearly identified, namely: the first (centrifugal) term in equation (9) is related to the kinetic energy of the condensate and would be neglected in the Thomas-Fermi approximation; the second term, proportional to γ, is due to the local, repulsive coupling; the third term, proportional to α, is due to the nonlocal, attractive coupling. We remark that in the case of an attractive contact interaction (γ<0) the total potential would be not bounded from below and no stable oscillations would be possible. This is due to the dominant contribution of the contact interaction as s  0, as apparent in equations (8) and (9).
The total potential is illustrated in figure 1, where α=20. Notice that increasing the interaction ratio decreases the depth of the potential well. Therefore a strong nonlocal attractive potential leads to a more stably bound state, as expected.
In the general case, it is easy to find η<1 as a necessary condition for obtaining stable minima of V. Analytical results can be derived in specific cases. For very large σ ? 1, equation ( 15 , which is consistently very small for a strong attractive nonlocal interaction. On the other hand, for negligible local interaction and σ=1, one derives s a » - Figure 2 shows the location of the stable fixed point σ (0) numerically obtained, as a function of the nonlocal interaction strength α, for a few choices of η. Notice that there is a minimum α allowing for periodic oscillations of the characteristic length σ. For instance, for η=0 one needs α3.5. Moreover, a smaller η corresponds to a smaller condensate width, due to an enhanced attractive coupling. Figure 3 shows the corresponding angular frequency ω of the linear oscillations, given by w In particular, we note a higher frequency for larger α and smaller η.

Numerical simulations
We here present numerical results showing equilibrium solutions as well as small-amplitude oscillations around the equilibrium for the Gaussian and van der Waals nonlocal potentials. For the Gaussian interaction, the fully 3D GPE (4) is solved numerically, and is compared with numerical solution of equations (8) with the pseudo-potential given by equation (9) as obtained via the variational approach. The nonlocal potential in equation (4) is in the form of a convolution integral, which is calculated efficiently in Fourier space by means of the convolution theorem, using the fast Fourier transform algorithm to numerically calculate the discrete Fourier transform and its inverse. A pseudo-spectral method is used to calculate the spatial derivatives. The pseudo-spectral method is very accurate as long as the solution is well represented on the numerical grid. (Using equation (6), the width of the wave function in the Fourier transformed space can be estimated to be approximately 1/σ, giving a width of the order 1-2; as seen e.g. in figure 2 at large α. We are using a grid spacing of Δx=5/16 giving the Nyqvist wavenumber π/Δx≈10 which is much larger than the spectral width and therefore the solution is resolved well on the numerical grid.) A standard 4th-order Runge-Kutta scheme is used to advance the solution in time. Figures 4 and 5 show the steady-state solutions of the GPE and variational systems, as well as the dynamics when the steadystate solutions are perturbed by 2%. The steady-state of y | | is found by solving the GPE and by iteratively doing an averaging procedure y y y y , which enforces the normalization y y á ñ = | 1. The new solution ψ new is then used as a starting point for a new iteration, until y | | converges to a time-stationary solution. The time difference t N −t 1 is typically chosen to cover one oscillation period of the system to efficiently damp out oscillations. Another averaging scheme is used for the variational solution, s = new s s s , to find the steady-state variational solution. As seen in the top and middle figures 4-5, the comparisons between the 3D GPE solutions and the variational solutions are in general very good, with only minor deviations from the Gaussian Ansatz. To compare the dynamics between the two approaches, the steady-state solutions were perturbed in a manner that the magnitude of ψ is initially 2% larger than the equilibrium solution at the origin x=y=z=0. For the variational approach this is straightforward by choosing a suitable initial condition for σ. For the GPE the initial condition was set by acoordinate stretching procedure of the steady-state solution, again enforcing y y á ñ = | 1 in the perturbed initial condition. Bottom panels of figures 4-5 show the time-evolution of y | | at the origin. Apart from a small difference in the steady-state solutions, the oscillation frequencies also agree to large degree, especially for the lower values of η. For η=0, one can also see oscillations with lower frequency in the 3D GPE solution, seen first as a drift toward lower y | | and then toward larger y | |, indicating that the total (local plus nonlocal) nonlinear potential of the deeply trapped condensate allows multiple oscillation modes.
For comparison, we also carried out numerical simulations using a van der Waals-like interaction [5,24,34], whereC 6 can be related to a strong van der Waals force between Rydberg atoms, while R corresponds to the saturation of the effective interaction, due to the van der Waals shift (see [5] for details). The 3D GPE in normalized units is  An attractive contact interaction (γ<0) would prevent stable nonlinear oscillations, due to the dominant role of the σ −4 term in equation (14) when σ approaches zero. While U(σ) can be expressed in terms of hypergeometric functions [5], a numerically amenable form of the nonlocal interaction potential is also (see appendix B) ò s a s = --- The numerical results for the van der Waals-like nonlocal interaction are presented in figures 6 and 7. The GPE and variational steady-state solutions show excellent agreement for η=0 and α=30 in figure 6, and also good comparison for η=0.15 as seen in figure 7. Also the dynamics show oscillation periods with high degree of agreements between the 3D GPE and variational solutions. Similar as for the Gaussian nonlocal interaction potential, one can also see a lower frequency oscillation in the 3D GPE solution for the case η=0, indicating that the deeply trapped condensate has a more rich dynamics with more than one trapped oscillatory mode.  A set of simulations, similar to the ones in figures 4-7, were carried out for a range of different η and α to produce phase maps shown in figures 8 and 9 for the van der Waals interaction, respectively, showing the central amplitude of ψ and the oscillation frequency for the variational and 3D GPE models. In general the agreement between the variational and 3D GPE solutions are good, with a typical deviation of about 5% or less.

Spectrum of elementary excitations
At this point and below it is more transparent to again employ physical variables. As a complementary study besides the search for localized nonlinear solutions, it is interesting to look for the linear wave spectrum around homogeneous equilibria of the GPE (1). This could be obtained using the Bogoliubov method [38]. Here we use an alternative method by employing the Madelung transformation [41] where n=n(r, t) is the condensate number density and the action S=S(r, t) is related to the condensate velocity field =  S M v , where both n and S are real quantities. Separating the real and imaginary parts of the GPE yields the continuity equation ¶ +  = · ( ) ( ) n nv 0 1 8 t and the momentum transport equation where ¢ = ¶ ¶ ¢ r and where in the last step an integration by parts was performed, assuming the nonlocal interaction to vanish at infinity. Assuming where n 0 is a uniform equilibrium number density and δn and δv are small, first-order perturbations. Linearizing the system and letting the first-order perturbations be plane wave solu- , the result is the dispersion relation in agreement with [33], where 3 denotes the Fourier transform of the nonlocal interaction. In the absence of W or for short wavelengths one regains the Bogoliubov spectrum.

Gaussian kernel
The hydrodynamic equations (18)- (19) and the linear dispersion relation (21) are valid irrespective of the form of the nonlocal interaction. In the case of the model (2), the result from equation (21) is    which shows in a conclusive way an instability (Ω 2 <0) of long wavelengths ((k R) 2 =1), provided η<1. In conclusion, a strong attractive nonlocal interaction destabilizes homogeneous equilibria but allows the existence of stable nonlinear, self-trapped wave packets.
In terms of dimensionless variables t W = W and = k k R, the wave spectrum is which is shown in figure 10. It can be seen that a larger η reduces the unstable (W < 0 2 ) region in wavenumber space.

Van der Waals-like kernel
The Fourier transform of the kernel shown in equation (11) is Comparing equations (24) and (27), we conclude that the long wavelength dispersion relations for the two different kernels are basically the same, up a small discrepancy of order unity in the ∼k 4 term, provided the coupling parameter η is correctly identified and the same characteristic length R is used.
In terms of dimensionless variables t W = W and = k k R, the wave spectrum is p a h f W = -+¯( (¯))¯( ) n R k N k k 2 3 4 , 2 9 2 2 0 3 2 4 which is shown in figure 11, with similar results as for the Gaussian kernel in figure 10.

Conclusions
In this work, the nonlinear dynamics of a completely selfbound (without external confinement) BEC with nonlocal coupling has been analyzed by means of a time-dependent variational formalism. For definiteness, the nonlocal and contact interactions were assumed to be respectively attractive and repulsive, but the proposed method can be easily adapted to alternative situations. It is found that the problem is always reducible to a nonlinear oscillator with a linearly stable equilibrium, whose location and breather frequency can be determined in great detail, specially for a nonlocal Gaussianshaped kernel, which is more amenable to analytical calculations. A van der Waals-like interaction was also considered, which can be related to an attractive effective force in mesoscopic atom clouds due to off-resonant dressing to Rydberg nD states [5]. We have made a comparison between the approximate, variational solutions and the output from the direct simulation of the nonlocal GPE. As a complement to the nonlinear analysis, the linear spectrum of spatio-temporal waves was also addressed and shown to be determined by the Fourier transform of the nonlocal potential. In this context, the Gaussian and van der Waals-like couplings were shown to yield essentially similar dispersion relations. The present results provide insight for the engineering of self-trapped matter waves. Topics to be addressed in future works include the analysis of the existence and stability of BECs with nonradial, rotational states, twisted phonons carrying a finite angular momentum [42], vortex-phonon interactions [43] and twisted nonlinear waves [44], in BECs with long-range interactions.