Simulating 1D time-independent Bose Einstein Condensation

In summary: It's a power, rather than a simple number like you would see in most physics formulas. Second, you're not taking into account the fact that the particles are interacting with each other. This is important because it means the wavefunction will be more complex and won't converge to a simple form like a Gaussian.In summary, your code doesn't seem to be working as it should.
  • #1
svletana
21
1
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.

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! :)
 
Technology news on Phys.org
  • #2
I'm not very fluent in Python, but your code looks all right. I can however comment on the approach.

svletana said:
# 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.

svletana said:
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.
svletana said:
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.
 
  • Like
Likes jim mcnamara and QuantumQuest
  • #3
DrClaude said:
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!
 

Related to Simulating 1D time-independent Bose Einstein Condensation

1. What is Bose Einstein Condensation (BEC)?

Bose Einstein Condensation is a quantum phenomenon that occurs when a group of bosons, which are particles with integer spin, occupy the same quantum state at low temperatures. This results in the formation of a macroscopic wavefunction, also known as a Bose Einstein condensate.

2. What is the significance of studying 1D time-independent BEC?

Studying 1D time-independent BEC allows us to understand the properties of BEC in a simplified and controlled setting. It also allows us to investigate the effects of external factors such as trapping potentials and particle interactions on the condensate's behavior.

3. How is 1D time-independent BEC simulated?

1D time-independent BEC can be simulated using numerical techniques such as the Gross-Pitaevskii equation, which describes the dynamics of the wavefunction of a BEC. Other methods include mean-field theories and Monte Carlo simulations.

4. What are the applications of studying 1D time-independent BEC?

The study of 1D time-independent BEC has various applications in fields such as quantum computing, precision measurements, and creating new states of matter. It also has implications for understanding the behavior of other quantum systems.

5. What are the challenges in simulating 1D time-independent BEC?

One of the main challenges in simulating 1D time-independent BEC is accurately modeling the interactions between particles, which can be complex and difficult to predict. Another challenge is dealing with the effects of external factors such as temperature and external forces on the condensate's behavior. Additionally, numerical simulations can be computationally intensive and require advanced techniques and resources.

Similar threads

  • Programming and Computer Science
Replies
1
Views
995
Replies
1
Views
2K
  • Programming and Computer Science
Replies
15
Views
1K
  • Programming and Computer Science
Replies
6
Views
2K
  • Thermodynamics
Replies
3
Views
1K
  • Programming and Computer Science
Replies
5
Views
1K
  • Programming and Computer Science
Replies
2
Views
2K
  • Programming and Computer Science
Replies
17
Views
2K
  • Programming and Computer Science
3
Replies
97
Views
7K
  • Programming and Computer Science
Replies
1
Views
1K
Back
Top