Numerical solution of the nonlinear Schrödinger equation using smoothed-particle hydrodynamics

We formulate a smoothed-particle hydrodynamics numerical method, traditionally used for the Euler equations for fluid dynamics in the context of astrophysical simulations, to solve the nonlinear Schrödinger equation in the Madelung formulation. The probability density of the wave function is discretized into moving particles, whose properties are smoothed by a kernel function. The traditional fluid pressure is replaced by a quantum pressure tensor, for which a robust discretization is found. We demonstrate our numerical method on a variety of numerical test problems involving the simple harmonic oscillator, soliton-soliton collision, Bose-Einstein condensates, collapsing singularities, and dark matter halos governed by the Gross-Pitaevskii-Poisson equation. Our method is conservative, applicable to unbounded domains, and is automatically adaptive in its resolution, making it well suited to study problems with collapsing solutions.

Figure 1

Convergence properties of the SPH code for obtaining the ground state of the Schrödinger equation with a one-dimensional (1D) simple harmonic oscillator potential. Second-order convergence is achieved by increasing particle number $N$ (and decreasing smoothing length as $1/N$). The obtained ground-state solution for $\rho$ at the various resolutions is shown in the inset (the shades of teal correspond to the shades of the circles in the convergence plot: higher resolution approaches the analytic answer, which is shown with the thick gray line).

Figure 2

Time evolution of oscillating wave function in simple harmonic oscillator potential. The SPH approach captures the dynamics well. The analytic solution is shown in thin gray lines.

Figure 3

Evolution of two solitons in simple harmonic oscillator potential. The SPH code captures the profile well, which consists of two solitons oscillating in opposite directions.

Figure 4

Collision of two bright solitons evolved under the NLSE. Collision occurs at $t=5$, during which fringes are formed. The solitons return to their original profile after interaction.

Figure 5

Collision of two dark solitons evolved under the NLSE. Collision occurs at $t=2.5$.

Figure 6

Ground state of 2D BEC with various values of the parameter $g$, as obtained by the SPH approach and compared to results from high resolution numerical simulations in Ref. [8].

Figure 7

Comparison of the blow-up solution of the focusing NLSE with numerical solutions obtained via our SPH method and a second-order finite difference scheme for comparison. The SPH solution produces a collapse with the right scaling while the finite difference scheme shows focusing-defocusing oscillations early on due to discretization effects. In the insets, we show the location of the SPH particles at three different times, and draw circles around each that have radii proportional to the adaptive smoothing length of the particle. The inset at the third time frame has an aggregation of 12 particles near the singularity.

Figure 8

3D BEC dark matter halo profiles calculated with our SPH method for various values of the parameter $g$. The profiles approach the Thomas-Fermi approximation for increasing values of $g$.