Prelude

First, let's load some packages and set plot defaults.

using Plots, LaTeXStrings
gr(legend=false,titlefontsize=12,size=(500,300),colorbar=false,grid=false)
using FourierGPE
import VortexDistributions:unwrap

Gross-Pitaevskii equation

We are going to solve the Gross-Pitaevskii equation

\[ i\hbar\frac{\partial \psi(x,t)}{\partial t}=\left(-\frac{\hbar^2\partial_x^2}{2m}+V(x,t)+g|\psi|^2\right)\psi \]

for particular initial and boundary conditions.

Bright soliton

The bright soliton provides a good test of any numerical simulation of the Gross-Pitaevskii equation as it involves a delicate balance between kinetic dispersion and the attractive nonlinearity. An initial state with finite momentum also tests the periodicity of the fft method since the soliton will eventually wrap around the domain.

The bright soliton wavefunction with wavenumber $k$ describing its collective motion is

\[ \psi_s(x)= \sqrt{\frac{N_s}{2\xi_s}}\textrm{sech}(x/\xi_s)e^{ikx} \]

where the soliton scale $\xi_s$ for $N_s$ particles is given by

\[ \xi_s \equiv \frac{2}{|g|N_s} \]

Simulation

Potential function

import FourierGPE.V
V(x,t) = 0.0 |> complex

V (generic function with 3 methods)

Units

In any numerical calculations we should have a clear understanding of our choice of physical units.

In length unit $\xi_s$, and time unit

\[ t_s\equiv \frac{m \xi_s^2}{\hbar}, \]

and rescaled wavefunction $\bar\psi = \psi\sqrt{\xi_s}$, our dimensionless form of the equation of motion is

\[ i\frac{\partial \bar\psi(\bar x,\bar t)}{\partial \bar t}=\left(-\frac{\bar\partial_x^2}{2}+\bar g|\bar\psi|^2\right)\bar\psi \]

where the dimensionless interaction parameter is

\[ \bar g \equiv \frac{m\xi_s}{\hbar^2}g < 0 \]

Initialize the simulation

Create the sim struct holding all parameters, with predefined grids.

L = (60.0,)
N = (512,)
sim = Sim(L,N)
@unpack_Sim sim;

Parameters

μ = 25.0
g = -0.01
γ = 0.0
Ns = 200
ξs = 2/abs(g)/Ns
us = 20
tf = 1π |> Float64
Nt = 150
t = LinRange(0.,tf,Nt);

Initial condition

We initialize the bright soliton with dimensionless velocity $u_s$ as

x = X[1]
ψs(x) = sqrt(Ns/2ξs)*sech(x/ξs)*exp(im*us*x)
ψi = ψs.(x)
ϕi = kspace(ψi,sim)

@pack_Sim! sim; #finally, pack everything up for simulation.

Evolve in k-space

Now we have everything we need to evolve in k-space

@time sol = runsim(sim);

4.515066 seconds (19.36 M allocations: 1.127 GiB, 5.77% gc time, 87.82% c
ompilation time)

Plot the solution

We plot the density of atoms, and the phase after removing the Galilean boost.

ψp = xspace(sol[end],sim).*exp.(-im*us*x) # drop translation factor for plot
p1 = plot(x,abs2.(ψp))
xlabel!(L"x");ylabel!(L"|\psi|^2")
p2 = plot(x,unwrap(angle.(ψp)))
xlabel!(L"x");ylabel!(L"\textrm{phase} (\psi)")
p = plot(p1,p2,layout=(2,1),size=(600,400))

The phase is constant over the soliton, as it would be in the lab frame. To visualize the motion, we can make an animation, saved to the media folder.

anim = @animate for i=1:Nt-8
    ψ = xspace(sol[i],sim)
    y = abs2.(ψ)
    plot(x,y,fill=(0, 0.2),size=(600,150),legend=false,xticks=false,yticks=false,axis=false)
end

filename = "brightsoliton.gif"
gif(anim,joinpath(@__DIR__,"../media",filename),fps = 25);

---