- #1

svletana

- 21

- 1

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

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]
# schrodinger'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! :)