Simulating 1D time-independent Bose Einstein Condensation

[svletana]

Hello! I'm trying to simulate a one dimensional time independent BEC, I hope this is the right place to ask for help.

First of all, here's my code in Python.

import sys
import numpy as np
import matplotlib.pyplot as plt

if len(sys.argv) == 1:
    niter = 100
elif len(sys.argv) == 2:
    niter = int(sys.argv[1]) # number of iterations
else:
    print('Please insert only the iteration number')
    exit(1)

# useful functions

# Laplacian operator for the second derivative
def Laplacian(x):
    h = x[1]-x[0] # assume uniformly spaced points
    n = len(x)
    M = -2 * np.diag([1]*n)
    for i in range(1,n):
        M[i,i-1] = 1
        M[i-1,i] = 1
    return M/h ** 2# for normalizing eigenstates
def Normalizate(U, x):
    h = x[1] - x[0] # assume uniformly spaced points
    n = len(x)
    for j in range(0, n):
        suma = 0.0
        for i in range(1, n):
            suma = suma + U[i, j]**2

        suma = suma * h
        rnorm = 1 / np.sqrt(suma)
        # print j,’ integral (sin normalizar) =’,rnorm
        # Normalization
        rsign = 1
        if U[1, j] < 0:
            rsign = -1

        rnorm = rnorm * rsign
        for i in range(0, n):
            U[i, j] = U[i, j] * rnorm

# parameters for Rubidium 87 (a.u.)
m = 87 # mass
a = 110 # scattering length
N = 500 # particle number
g = 4 * np.pi * a * N / m # coefficient for gross pitaevskii equation
w = 1 # frequancy of trapping potential
L = 1.5
dx = 600

# initialize wavefunction
x = np.linspace(-L, L, dx)
phi = [20 * np.exp(- (i ** 2)) for i in x]

W = 0.5 * m * (w ** 2) * (x ** 2)

# main loop
for i in range(niter):

    # calculate density
    density = [abs(p * p) for p in phi]

    # Schrödinger's eq in matrix form
    T = -0.5 / m * Laplacian(x) # kinetic energy
    P = g * np.diag(density)  #GPE non linear term

    H = T + np.diag(W) + P #hamiltonian
    E, U = np.linalg.eigh(H) #eigenvalues and eigenstates

    Normalizate(U, x)

    phi = U[:, 0]print('chemical potential: ' + str(E[0]))
### Theorical value for the chemical potential, only works for the Thomas Fermi approximation!
ff = (3*g*N*np.sqrt(m/2)*w)**(2/3)
print('theoretical value (TF): ' + str(ff))

# plt.plot(x,W,'--k')
plt.plot(x,U[:, 0])
plt.show()

It's not very complex up to this point. What my program does is start the simulation with some random wavefunction. I calculate the density and build the Gross-Pitaevskii Equation (GPE) and solve by calculating eigenvalues and eigenstates. Afterwards, I normalize the eigenstates and start over by taking the new eigenstate, calculating its density and solving the GPE once again (I'm working only with the ground state). It's a fixed point iteration.

I'm hoping that after a good number of iterations, the wavefunction will look like a Gaussian, which should be the ground state for a 1D BEC.

However, if you run the program you will find the wavefunction doesn't seem to converge to any particular waveform. The chemical potential values also don't match what I predict (around 0.5ω for a small number of particles, where ω is the trapping potential's frequency, and the theoretical value for the TF approximation for a large number of particles).

Also, I don't see how the particle number is involved in this equation other than in the g coefficient. How is that supposed to change my result aside from a different effective potential?

I'll gladly take any other suggestion or advice! Thanks in advance! :)

 

[DrClaude]

I'm not very fluent in Python, but your code looks all right. I can however comment on the approach.

# parameters for Rubidium 87 (a.u.)
m = 87 # mass
a = 110 # scattering length
N = 500 # particle number
g = 4 * np.pi * a * N / m # coefficient for gross pitaevskii equation
w = 1 # frequancy of trapping potential
L = 1.5
dx = 600

Two problems here. First, the potential strength g is not the same for a 1D condensate than for a 3D, see M. Olshanii, Phys. Rev. Lett. 81, 938 (1998). Second, the atomic unit of mass is not the same as the atomic mass unit. In atomic units, the mass of rubidium is not 87.

It's not very complex up to this point. What my program does is start the simulation with some random wavefunction. I calculate the density and build the Gross-Pitaevskii Equation (GPE) and solve by calculating eigenvalues and eigenstates. Afterwards, I normalize the eigenstates and start over by taking the new eigenstate, calculating its density and solving the GPE once again (I'm working only with the ground state). It's a fixed point iteration.

This should work for small values of the nonlinearity. The value of g you are using is about 7944, so I am not surprised that it is not working. You should start with $g \sim 1$ and make sure it works then. Also, are you normalizing the initial wave function? For it to converge, you should start as close as possible to the expected final state, so I would take the ground state of the corresponding potential without the non-linear term.

Also, I don't see how the particle number is involved in this equation other than in the g coefficient. How is that supposed to change my result aside from a different effective potential?

If you normalize $\psi$ (which I wouldn't call the wave function, but that's another story), to 1, the particle number only enters through $g$. At the level of approximation of the GRoss-Pitaevskii equation, different BECs will behave the same provided that $g$ is the same (if you change species, hence the mass and the scattering length, you simply need to change the number of atoms to get the same result). Remember however that $|\psi|^2$ is then related to the total number of condensed atoms.

 

[svletana]

I'm not very fluent in Python, but your code looks all right. I can however comment on the approach.Two problems here. First, the potential strength g is not the same for a 1D condensate than for a 3D, see M. Olshanii, Phys. Rev. Lett. 81, 938 (1998). Second, the atomic unit of mass is not the same as the atomic mass unit. In atomic units, the mass of rubidium is not 87.This should work for small values of the nonlinearity. The value of g you are using is about 7944, so I am not surprised that it is not working. You should start with $g \sim 1$ and make sure it works then. Also, are you normalizing the initial wave function? For it to converge, you should start as close as possible to the expected final state, so I would take the ground state of the corresponding potential without the non-linear term.
If you normalize $\psi$ (which I wouldn't call the wave function, but that's another story), to 1, the particle number only enters through $g$. At the level of approximation of the GRoss-Pitaevskii equation, different BECs will behave the same provided that $g$ is the same (if you change species, hence the mass and the scattering length, you simply need to change the number of atoms to get the same result). Remember however that $|\psi|^2$ is then related to the total number of condensed atoms.

Thanks for your help! :)
Instead of setting the real g value my teacher told me to try any value, so I did that. There was a problem with normalization as well, that was messing up the whole calculations. Thanks again!