Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Simulating 1D time-independent Bose Einstein Condensation

  1. Jun 26, 2017 #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.

    Code (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
        print('Please insert only the iteration number')

    # 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])
    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! :)
  2. jcsd
  3. Jun 29, 2017 #2


    User Avatar

    Staff: Mentor

    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.
  4. Jul 3, 2017 #3
    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!
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted

Similar Discussions: Simulating 1D time-independent Bose Einstein Condensation
  1. 1D simulator (Ruby) (Replies: 0)