# Bosonic statistics and the Bose-Einstein condensation¶

Doruk Efe Gökmen -- 30/08/2018 -- Ankara

## Non-interacting ideal bosons¶

Non-interacting bosons is the only system in physics that can undergo a phase transition without mutual interactions between its components.

Let us enumerate the energy eigenstates of a single 3D boson in an harmonic trap by the following program.

In [1]:
Emax = 30
States = []
for E_x in range(Emax):
for E_y in range(Emax):
for E_z in range(Emax):
States.append(((E_x + E_y + E_z), (E_x, E_y, E_z)))
States.sort()
for k in range(Emax):
print '%3d' % k, States[k][0], States[k][1]

  0 0 (0, 0, 0)
1 1 (0, 0, 1)
2 1 (0, 1, 0)
3 1 (1, 0, 0)
4 2 (0, 0, 2)
5 2 (0, 1, 1)
6 2 (0, 2, 0)
7 2 (1, 0, 1)
8 2 (1, 1, 0)
9 2 (2, 0, 0)
10 3 (0, 0, 3)
11 3 (0, 1, 2)
12 3 (0, 2, 1)
13 3 (0, 3, 0)
14 3 (1, 0, 2)
15 3 (1, 1, 1)
16 3 (1, 2, 0)
17 3 (2, 0, 1)
18 3 (2, 1, 0)
19 3 (3, 0, 0)
20 4 (0, 0, 4)
21 4 (0, 1, 3)
22 4 (0, 2, 2)
23 4 (0, 3, 1)
24 4 (0, 4, 0)
25 4 (1, 0, 3)
26 4 (1, 1, 2)
27 4 (1, 2, 1)
28 4 (1, 3, 0)
29 4 (2, 0, 2)


Here it can be perceived that the degeneracy at an energy level $E_n$, which we denote by $\mathcal{N}(E_n)$, is $\frac{(n+1)(n+2)}{2}$. Alternatively, we may use a more systematic approach. We can calculate the number of states at the $n$th energy level as $\mathcal{N}(E_n)=\sum_{E_x=0}^{E_n}\sum_{E_y=0}^{E_n}\sum_{E_z=0}^{E_n}\delta_{(E_x+E_y+E_z),E_n}$, where $\delta_{j,k}$ is the Kronecker delta. In the continuous limit we have the Dirac delta function

$\delta_{j,k}\rightarrow\delta(j-k) =\int_{-\pi}^\pi \frac{\text{d}\lambda}{2\pi}e^{i(j-k)\lambda}$. (1)

If we insert this function into above expression, we get $\mathcal{N}(E_n)=\int_{-\pi}^\pi \frac{\text{d}\lambda}{2\pi}e^{-iE_n\lambda}\left(\sum_{E_x=0}^{E_n}e^{iE_x\lambda}\right)^3$. The geometric sum can be evaluated, hence we have the integral $\mathcal{N}(E_n)=\int_{-\pi}^\pi \frac{\text{d}\lambda}{2\pi}e^{-iE_n\lambda}\left[\frac{1-e^{i\lambda (n+1)}}{1-e^{i\lambda}}\right]^3$. The integration range corresponds to a circular contour $\mathcal{C}$ of radius 1 centered at 0 at the complex plane. If we define $z=e^{i\lambda}$, the integral transforms into $\mathcal{N}(E_n)=\frac{1}{2\pi i}\oint_{\mathcal{C}}\frac{\text{d}z}{z^{n+1}}\left[\frac{1-z^{n+1}}{1-z}\right]^3$. Using the residue theorem, this integral can be evaluated by determining the coefficient of the $z^{-1}$ term in the Laurent series of $\frac{1}{z^{n+1}}\left[\frac{1-z^{n+1}}{1-z}\right]^3$, which is $(n+1)(n+1)/2$. Hence we recover the previous result.

##### Five boson bounded trap model¶

Consider 5 bosons in the harmonic trap, but with a cutoff on the single-particle energies: $E_\sigma\leq 4$. There are $34$ possible single-particles energy states. For this model, the above naive enumeration of these energy states still works. We can label the state of each 5 particle by $\sigma_i$, so that $\{\text{5-particle state}\}=\{\sigma_1,\cdots,\sigma_5\}$. The partititon function for this system is given by $Z(\beta)=\sum_{0\leq\sigma_1\leq\cdots\leq\sigma_5\leq 34}e^{-\beta E(\sigma_1,\cdots,\sigma_5)}$. In the following program, the average occupation number of the ground state per particle is calculated at different temperatures (corresponds to the condensate). However, due to the nested for loops, this method is very inconvenient for higher number of particles.

In [10]:
%pylab inline
import math, numpy as np, pylab as plt

#calculate the partition function for 5 bosons by stacking the bosons in one of the N_states
#number of possible states and counting only a specific order of them (they are indistinguishable)
def bosons_bounded_harmonic(beta, N):
Energy = [] #initialise the vector that the energy values are saved with enumeration
n_states_1p = 0 #initialise the total number of single trapped boson states
for n in range(N + 1):
degeneracy = (n + 1) * (n + 2) / 2  #degeneracy in the 3D harmonic oscillator
Energy += [float(n)] * degeneracy
n_states_1p += degeneracy

n_states_5p = 0 #initialise the total number states of 5 trapped bosons
Z = 0.0 #initialise the partition function
N0_mean = 0.0
E_mean = 0.0
for s_0 in range(n_states_1p):
for s_1 in range(s_0, n_states_1p): #consider the order s_0<s_1... to avoid overcounting
for s_2 in range(s_1, n_states_1p):
for s_3 in range(s_2, n_states_1p):
for s_4 in range(s_3, n_states_1p):
n_states_5p += 1
state = [s_0, s_1, s_2, s_3, s_4] #construct the state of each 5 boson
E = sum(Energy[s] for s in state) #calculate the total energy by above enumeration
Z += math.exp(-beta * E) #canonical partition function
E_mean += E * math.exp(-beta * E) #avg. total energy
N0_mean += state.count(0) * math.exp(-beta * E) #avg. ground level occupation number
return n_states_5p, Z, E_mean, N0_mean

N = 4 #the energy cutoff for each boson
beta = 1.0 #inverse temperature

n_states_5p, Z, E_mean, N0_mean = bosons_bounded_harmonic(beta, N)

print 'Temperature:', 1 / beta, 'Total number of possible states:', n_states_5p, '| Partition function:', Z,\
'| Average energy per particle:', E_mean / Z / 5.0,\
'| Condensate fraction (ground state occupation per particle):', N0_mean / Z / 5.0

cond_frac = []
temperature = []
for T in np.linspace(0.1, 1.0, 10):
n_states_5p, Z, E_mean, N0_mean = bosons_bounded_harmonic(1.0 / T, N)
cond_frac.append(N0_mean / Z / 5.0)
temperature.append(T)

plt.plot(temperature, cond_frac)
plt.title('Condensate? fraction for the $N=5$ bosons bounded trap model ($N_{bound}=%i$)' % N, fontsize = 14)
plt.xlabel('$T$', fontsize = 14)
plt.ylabel('$\\langle N_0 \\rangle$ / N', fontsize = 14)
plt.grid()

Populating the interactive namespace from numpy and matplotlib
Temperature: 1.0 Total number of possible states: 575757 | Partition function: 17.3732972183 | Average energy per particle: 1.03133265311 | Condensate fraction (ground state occupation per particle): 0.446969501933


Here we see that all particles are in the ground states at very low temperatures this is a simple consequence of Boltzmann statistics. At zero temperature all the particles populate the ground state. Bose-Einstein condensation is something else, it means that a finite fraction of the system is in the ground-state for temperatures which are much higher than the gap between the gap between the ground-state and the first excited state, which is one, in our system. Bose-Einstein condensation occurs when all of a sudden a finite fraction of particles populate the single-particle ground state. In a trap, this happens at higher and higher temperatures as we increase the particle number.

Alternatively, we can characterise any single particle state $\sigma=0,\cdots,34$ by an occupation number $n_\sigma$. Using this occupation number representation, the energy is given by $E=n_0E_0+\cdots + n_{34}E_{34}$, and the partition function is $Z(\beta)=\sum^{N=5}_{n_0=0}\cdots\sum^{N=5}_{n_{34}=0}e^{-\beta(n_0E_0+\cdots + n_{34}E_{34})}\delta_{(n_0+\cdots + n_{34}),N=5}$. Using the integral representation of the Kronecker delta given in (1), and evaluating the resulting sums, we have

$Z(\beta)=\int_{-\pi}^\pi\frac{\text{d}\lambda}{2\pi}e^{-iN\lambda}\Pi_{E=0}^{E_\text{max}}[f_E(\beta,\lambda)]^{\mathcal{N}(E)}$. (2)

### The bosonic density matrix¶

Distinguishable particles: The partition function of $N$ distinguishable particles is given by $Z^D(\beta)=\int \text{d}\mathbf{x}\rho(\mathbf{x},\mathbf{x},\beta)$, where $\mathbf{x}=\{0,\cdots,N-1\}$, i.e. the positions of the $i$th particle; and $\rho$ is the $N$ distinguishable particle density matrix. If the particles are non-interacting (ideal), then the density matrix can simply be decomposed into $N$ single particle density matrices as

$\rho^{D,\text{ideal}}(\mathbf{x},\mathbf{x}',\beta)=\Pi_{i=0}^{N-1}\rho(x_i,x_i',\beta)$, (3)

with the single particle density matrix $\rho(x_i,x_i',\beta)=\sum_{\lambda_i=0}^{\infty}\psi_{\lambda_i}(x_i)\psi_{\lambda_i}^{*}(x'_i)e^{-\beta E_{\lambda_i}}$, where $\lambda_i$ is the energy eigenstate of the $i$th particle. That means that the quantum statistical paths of the two particles are independent. More generally, the interacting many distinguishable particle density matrix is

$\rho^{D}(\mathbf{x},\mathbf{x}',\beta)=\sum_{\sigma}\Psi_{\sigma}(\mathbf{x})\Psi_{\sigma}^{*}(\mathbf{x}')e^{-\beta E_{\sigma}}$, (4)

where the sum is done over the all possible $N$ particle states $\sigma=\{\lambda_0,\cdots,\lambda_{N-1}\}$. The interacting paths are described by the paths whose weight are modified through Trotter decomposition, which correlates those paths.

Indistinguishable particles: The particles $\{0,\cdots,N-1\}$ are indistinguishable if and only if

$\Psi_{\sigma_\text{id}}(\mathbf{x})=\xi^\mathcal{P}\Psi_{\sigma_\text{id}}(\mathcal{P}\mathbf{x})$ $\forall \sigma$, (5)

where they are in an indistinguishable state ${\sigma_\text{id}}$, $\mathcal{P}$ is any $N$ particle permutation and the species factor $\xi$ is $-1$ (antisymmetric) for fermions, and $1$ (symmetric) for bosons. Here we focus on the bosonic case. Since there are $N!$ such permutations, if the particles are indistinguishable bosons, using (5) we get $\frac{1}{N!}\sum_{\mathcal{P}}\Psi_\sigma(\mathcal{P}x)=\Psi_\sigma(\mathbf{x})$, i.e. $\Psi_\sigma(x)=\Psi_{\sigma_\text{id}}(x)$. Furthermore, from a group theory argument it follows that $\frac{1}{N!}\sum_{\mathcal{P}}\Psi_\sigma(\mathcal{P}x)=0$ otherwise (fermionic or distinguishable). This can be expressed in a more compact form as

$\frac{1}{N!}\sum_{\mathcal{P}}\Psi_\sigma(\mathcal{P}x)=\delta_{{\sigma_\text{id}},\sigma}\Psi_\sigma(x)$. (6)

By definition, the bosonic density matrix should be $\rho^\text{bose}(\mathbf{x},\mathbf{x}',\beta)=\sum_{\sigma=\{\sigma_\text{id}\}}\Psi_\sigma(\mathbf{x})\Psi^{*}_\sigma(\mathbf{x}')e^{-\beta E_\sigma}=\sum_{\sigma}\delta_{{\sigma_\text{id}},\sigma}\Psi_\sigma(\mathbf{x})\Psi^{*}_\sigma(\mathbf{x}')e^{-\beta E_\sigma}$, i.e. a sum over all $N$ particle states which are symmetric. If we insert Eqn. (6) here in the latter equality, we get $\rho^\text{bose}(\mathbf{x},\mathbf{x}',\beta)=\frac{1}{N!}\sum_\sigma\Psi_\sigma(\mathbf{x})\sum_\mathcal{P}\Psi^{*}_\sigma(\mathcal{P}\mathbf{x}')e^{-\beta E_\sigma}$. Exchanging the sums, we get $\rho^\text{bose}(\mathbf{x},\mathbf{x}',\beta)=\frac{1}{N!}\sum_\mathcal{P}\sum_\sigma\Psi_\sigma(\mathbf{x})\Psi^{*}_\sigma(\mathcal{P}\mathbf{x}')e^{-\beta E_\sigma}$. In other words, we simply have

$\boxed{\rho^\text{bose}(\mathbf{x},\mathbf{x}',\beta)=\frac{1}{N!}\sum_\mathcal{P}\rho^D(\mathbf{x},\mathcal{P}\mathbf{x}',\beta)}$, (7)

that is the average of the distinguishable density matrices over all permutations of $N$ particles.

For ideal bosons, we have $\boxed{\rho^\text{bose, ideal}(\mathbf{x},\mathbf{x}',\beta)=\frac{1}{N!}\sum_\mathcal{P}\rho(x_0,\mathcal{P}x_0',\beta)\rho(x_1,\mathcal{P}x_1',\beta)\cdots\rho(x_{N-1},\mathcal{P}x_{N-1}',\beta)}$. (8)

The partition function is therefore

$Z^\text{bose}(\beta)=\frac{1}{N!}\int \text{d}x_0\cdots\text{d}x_{N-1}\sum_\mathcal{P}\rho^D(\mathbf{x},\mathcal{P}\mathbf{x},\beta)=\frac{1}{N!}\sum_\mathcal{P}Z_\mathcal{P}$, (9)

i.e. an integral over paths and an average over all permutations. We should therefore sample both positions and permutations.

For fermions, the sum over permutations $\mathcal{P}$ involve a weighting with factor $(-1)^{\mathcal{P}}$:

$\rho^\text{fermi}(\mathbf{x},\mathbf{x}',\beta)=\frac{1}{N!}\sum_\mathcal{P}(-1)^\mathcal{P}\rho^D(\mathbf{x},\mathcal{P}\mathbf{x}',\beta)$

Therefore for fermions corresponding path integrals are nontrivial, and they involve Grassmann variables (see e.g. Negele, Orland https://www.amazon.com/Quantum-Many-particle-Systems-Advanced-Classics/dp/0738200522 ).

#### Sampling permutations¶

The following Markov-chain algorithm samples permutations of $N$ elements on a list $L$. The permutation function for the uniformly distributed $\mathcal{P}$ is $Y_N=\sum_\mathcal{P}1=N!$.

In [ ]:
import random

N = 3 #length of the list
statistics = {}
L = range(N) #initialise the list
nsteps = 10
for step in range(nsteps):
i = random.randint(0, N - 1) #pick two random indices i and j from the list L
j = random.randint(0, N - 1)
L[i], L[j] = L[j], L[i] #exchange the i'th and j'th elements
if tuple(L) in statistics:
statistics[tuple(L)] += 1 #if a certain configuration appears again, add 1 to its count
else:
statistics[tuple(L)] = 1 #if a certain configuration for the first time, give it a count of 1
print L
print range(N)
print

for item in statistics:
print item, statistics[item]


Let us look at the permutation cycles and their frequency of occurrence:

In [26]:
import random

N = 20 #length of the list
stats = [0] * (N + 1) #initialise the "stats" vector
L = range(N) #initialise the list
nsteps = 1000000 #number of steps
for step in range(nsteps):
i = random.randint(0, N - 1) #pick two random indices i and j from the list L
j = random.randint(0, N - 1)
L[i], L[j] = L[j], L[i] #exchange the i'th and j'th elements in the list L
#Calculate the lengths of the permutation cycles in list L
if step % 100 == 0: #i.e. at each 100 steps
cycle_dict = {} #initialise the permutation cycle dictionary
for k in range(N): #loop over the list length,where keys (k) represent the particles
cycle_dict[k] = L[k] #and the values (L) are for the successors of the particles in the perm. cycle
while cycle_dict != {}: #i.e. when the cycle dictionary is not empty?
starting_element = cycle_dict.keys()[0] #save the first (0th) element in the cycle as the starting element
cycle_length = 0 #initialise the cycle length
old_element = starting_element #ancillary variable
while True:
cycle_length += 1 #increase the cycle length while...
new_element = cycle_dict.pop(old_element) #get the successor of the old element in the perm. cycle
if new_element == starting_element: break #the new element is the same as the first one (cycle complete)
else: old_element = new_element #move on to the next successor in the perm. cycle
stats[cycle_length] += 1 #increase the number of occurrences of a cycle of that length by 1
for k in range(1, N + 1): #print the cycle lengths and their number of occurrences
print k, stats[k]

1 10130
2 5008
3 3395
4 2438
5 1969
6 1659
7 1403
8 1260
9 1118
10 949
11 943
12 833
13 778
14 745
15 642
16 618
17 610
18 553
19 530
20 492


The partition function of permutations $\mathcal{P}$ on a list of lentgth $N$ is $Y_N=\sum_\mathcal{P}\text{weight}(\mathcal{P})$. Let $z_n$ be the weight of a permutation cycle of length $n$. Then, the permutation $[0,1,2,3]\rightarrow[0,1,2,3]$, which can be represented as $(0)(1)(2)(3)$, has the weight $z_1^4$; similarly, $(0)(12)(3)$ would have $z_1^2z_2$, etc.

Generally, the cycle $\{n_1,\cdots,n_{k-1},\text{last element}\}$, i.e. the cycle containing the last element, has a length $k$, with the weight $z_k$. The remaining $N-k$ elements have the partition function $Y_{(N-k)}$. Hence, the total partition function is given by $Y_N=\sum_{k=1}^Nz_k\{\text{# of choices for} \{n_1,\cdots,n_{k-1}\}\}\{\text{# of cycles with} \{n_1,\cdots,n_{k}\}\}Y_{N-k}$

$\implies Y_N=\sum_{k=1}^N z_k{{N-1}\choose{k-1}}(k-1)!Y_{N-k}$ which leads to the following recursion formula

$\boxed{Y_N=\frac{1}{N}\sum_{k=1}^N z_k\frac{N!}{(N-k)!}Y_{N-k}, (\text{with }Y_0=1)}$. (10)

Using the convolution property, we can regard the $l+1$ bosons in a permutation cycle of length $l$ at temperatyre $1/\beta$ as a single boson at a temperature $1/(l\beta)$.

Example 1: Consider the permutation $[0,3,1,2]\rightarrow[0,1,2,3]$ consists of the following permutation cycle $1\rightarrow 2 \rightarrow 3 \rightarrow 1$ of length 3 ($\mathcal{P}=(132)$). This corresponds to the partition function $Z^\text{bose}_{(0)(132)}(\beta)=\int \text{d}x_0\rho(x_0,x_0,\beta)\int\text{d}x_1\int\text{d}x_2\int\text{d}x_3\rho(x_1,x_3,\beta)\rho(x_3,x_2,\beta)\rho(x_2,x_1,\beta)$. Using the convolution property, we have: $\int\text{d}x_3\rho(x_1,x_3,\beta)\rho(x_3,x_2,\beta)=\rho(x_1,x_2,2\beta)\implies\int\text{d}x_2\rho(x_1,x_2,2\beta)\rho(x_2,x_1,\beta)=\rho(x_1,x_1,3\beta)$. The single particle partition function is defined as $z(\beta)=\int\text{d}\mathbf{x}\rho(\mathbf{x},\mathbf{x},\beta) =\left[ \int\text{d}x\rho(x,x,\beta)\right]^3$.

$\implies Z^\text{bose}_{(0)(132)}(\beta)=\int \text{d}x_0\rho(x_0,x_0,\beta)\int\text{d}x_1\rho(x_1,x_1,3\beta)=z(\beta)z(3\beta)$.

Example 2: $Z^\text{bose}_{(0)(1)(2)(3)}(\beta)=z(\beta)^4$.

Simulation of bosons in a harmonic trap: (Carefully note that here are no intermediate slices in the sampled paths, since the paths are sampled from the exact distribution.)

In [ ]:
import random, math, pylab, mpl_toolkits.mplot3d

#3 dimensional Levy algorithm, used for resampling the positions of entire permutation cycles of bosons
#to sample positions
def levy_harmonic_path(k, beta):
#direct sample (rejection-free) three coordinate values, use diagonal density matrix
#k corresponds to the length of the permutation cycle
xk = tuple([random.gauss(0.0, 1.0 / math.sqrt(2.0 *
math.tanh(k * beta / 2.0))) for d in range(3)])
x = [xk] #save the 3 coordinate values xk into a 3d vector x (final point)
for j in range(1, k): #loop runs through the permutation cycle
#Levy sampling (sample a point given the latest sample and the final point)
Upsilon_1 = (1.0 / math.tanh(beta) +
1.0 / math.tanh((k - j) * beta))
Upsilon_2 = [x[j - 1][d] / math.sinh(beta) + xk[d] /
math.sinh((k - j) * beta) for d in range(3)]
x_mean = [Upsilon_2[d] / Upsilon_1 for d in range(3)]
sigma = 1.0 / math.sqrt(Upsilon_1)
dummy = [random.gauss(x_mean[d], sigma) for d in range(3)] #direct sample the j'th point
x.append(tuple(dummy)) #construct the 3d path (permutation cycle) by appending tuples
return x

#(Non-diagonal) harmonic oscillator density matrix, used for organising the exchange of two elements
#to sample permutations
def rho_harm(x, xp, beta):
Upsilon_1 = sum((x[d] + xp[d]) ** 2 / 4.0 *
math.tanh(beta / 2.0) for d in range(3))
Upsilon_2 = sum((x[d] - xp[d]) ** 2 / 4.0 /
math.tanh(beta / 2.0) for d in range(3))
return math.exp(- Upsilon_1 - Upsilon_2)

N = 256 #number of bosons
T_star = 0.3
beta = 1.0 / (T_star * N ** (1.0 / 3.0)) #??
nsteps = 1000000
positions = {} #initial position dictionary
for j in range(N): #loop over all particles, initial permutation is identity (k=1)
a = levy_harmonic_path(1, beta) #initial positions (outputs a single 3d point)
positions[a[0]] = a[0] #positions of particles are keys for themselves in the initial position dict.
for step in range(nsteps):
boson_a = random.choice(positions.keys()) #randomly pick the position of boson "a" from the dict.
perm_cycle = [] #initialise the permutation cycle
while True: #compute the permutation cycle of the boson "a":
perm_cycle.append(boson_a) #construct the permutation cycle by appending the updated position of boson "a"
boson_b = positions.pop(boson_a) #remove and return (pop) the position of "a", save it as a temp. var.
if boson_b == perm_cycle[0]: break #if the cycle is completed, break the while loop
else: boson_a = boson_b #move boson "a" to position of "b" and continue permuting
k = len(perm_cycle) #length of the permutation cycle
#SAMPLE POSITIONS:
perm_cycle = levy_harmonic_path(k, beta) #resample the particle positions in the current permutation cycle
positions[perm_cycle[-1]] = perm_cycle[0] #assures that the new path is a "cycle" (last term maps to the first term)
for j in range(len(perm_cycle) - 1): #update the positions of bosons
positions[perm_cycle[j]] = perm_cycle[j + 1] #construct the "cycle": j -> j+1
#SAMPLE PERMUTATION CYCLES by exchanges:
#Pick two particles and attempt an exchange to sample permutations (with Metropolis acceptance rate):
a_1 = random.choice(positions.keys()) #pick the first random particle
b_1 = positions.pop(a_1) #save the random particle to a temporary variable
a_2 = random.choice(positions.keys()) #pick the second random particle
b_2 = positions.pop(a_2) #save the random particle to a temporary variable
weight_new = rho_harm(a_1, b_2, beta) * rho_harm(a_2, b_1, beta) #the new Metropolis acceptance rate
weight_old = rho_harm(a_1, b_1, beta) * rho_harm(a_2, b_2, beta) #the old Metropolis acceptance rate
if random.uniform(0.0, 1.0) < weight_new / weight_old:
positions[a_1] = b_2 #accept
positions[a_2] = b_1
else:
positions[a_1] = b_1 #reject
positions[a_2] = b_2

#Figure output:
fig = pylab.figure()
ax = mpl_toolkits.mplot3d.axes3d.Axes3D(fig)
ax.set_aspect('equal')
list_colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k']
n_colors = len(list_colors)
dict_colors = {}
i_color = 0
# find and plot permutation cycles:
while positions:
x, y, z = [], [], []
starting_boson = positions.keys()[0]
boson_old = starting_boson
while True:
x.append(boson_old[0])
y.append(boson_old[1])
z.append(boson_old[2])
boson_new = positions.pop(boson_old)
if boson_new == starting_boson: break
else: boson_old = boson_new
len_cycle = len(x)
if len_cycle > 2:
x.append(x[0])
y.append(y[0])
z.append(z[0])
if len_cycle in dict_colors:
color = dict_colors[len_cycle]
ax.plot(x, y, z, color + '+-', lw=0.75)
else:
color = list_colors[i_color]
i_color = (i_color + 1) % n_colors
dict_colors[len_cycle] = color
ax.plot(x, y, z, color + '+-', label='k=%i' % len_cycle, lw=0.75)
# finalize plot
pylab.title('$N=%i$, $T^*=%s$' % (N, T_star))
pylab.legend()
ax.set_xlabel('$x$', fontsize=16)
ax.set_ylabel('$y$', fontsize=16)
ax.set_zlabel('$z$', fontsize=16)
ax.set_xlim3d([-8, 8])
ax.set_ylim3d([-8, 8])
ax.set_zlim3d([-8, 8])
pylab.savefig('snapshot_bosons_3d_N%04i_Tstar%04.2f.png' % (N, T_star))
pylab.show()


But we do know that for the harmonic trap, the single 3-dimensional particle partition function is given by $z(\beta)=\left(\frac{1}{1-e^{-\beta}}\right)^3$. The permutation cycle of length $k$ corresponds to $z_k=z(k\beta)=\left(\frac{1}{1-e^{-k\beta}}\right)^3$. Hence, using (9) and (10), we have that

$Z^\text{bose}_N=Y_N/{N!}=\frac{1}{N}\sum_{k=1}^N z_k Z^\text{bose}_{N-k}, (\text{with }Z^\text{bose}_0=1)$. (11)

(Due to Landsberg, 1961 http://store.doverpublications.com/0486664937.html) This recursion relation relates the partition function of a system of $N$ ideal bosons to the partition function of a single particle and the partition functions of systems with fewer particles.

In [36]:
import math, pylab

def z(k, beta):
return 1.0 / (1.0 - math.exp(- k * beta)) ** 3 #partition function of a single particle in a harmonic trap

def canonic_recursion(N, beta): #Landsberg recursion relations for the partition function of N bosons
Z = [1.0] #Z_0 = 1
for M in range(1, N + 1):
Z.append(sum(Z[k] * z(M - k, beta) \
for k in range(M)) / M)
return Z #list of partition functions for boson numbers up to N

N = 256 #number of bosons
T_star = 0.5 #temperature
beta = 1.0 / N ** (1.0 / 3.0) / T_star
Z = canonic_recursion(N, beta) #partition function
pi_k = [(z(k, beta) * Z[N - k] / Z[-1]) / float(N) for k in range(1, N + 1)] #probability of a cycle of length k
# graphics output
pylab.plot(range(1, N + 1), pi_k, 'b-', lw=2.5)
pylab.ylim(0.0, 0.01)
pylab.xlabel('cycle length $k$', fontsize=16)
pylab.ylabel('cycle probability $\pi_k$', fontsize=16)
pylab.title('Cycle length distribution ($N=%i$, $T^*=%s$)' % (N, T_star), fontsize=16)
pylab.savefig('plot-prob_cycle_length.png')

phase = [pi[k+1] - pi[k] for k in range(1, N+1)]

# graphics output
pylab.plot(range(1, N + 1), pi_k, 'b-', lw=2.5)
pylab.ylim(0.0, 0.01)
pylab.xlabel('cycle length $k$', fontsize=16)
pylab.ylabel('cycle probability $\pi_k$', fontsize=16)
pylab.title('Cycle length distribution ($N=%i$, $T^*=%s$)' % (N, T_star), fontsize=16)
pylab.savefig('plot-prob_cycle_length.png')


Since we have an analytical solution to the problem, we can now implement a rejection-free direct sampling algorithm for the permutations.

In [40]:
import math, random

def z(k, beta): #partition function of a single particle in a harmonic trap
return (1.0 - math.exp(- k * beta)) ** (-3)

def canonic_recursion(N, beta): #Landsberg recursion relation for the partition function of N bosons in a harmonic trap
Z = [1.0]
for M in range(1, N + 1):
Z.append(sum(Z[k] * z(M - k, beta) for k in range(M)) / M)
return Z

def make_pi_list(Z, M): #the probability for a boson to be in a permutation length of length up to M?
pi_list = [0.0] + [z(k, beta) * Z[M - k] / Z[M] / M for k  in range(1, M + 1)]
pi_cumulative = [0.0]
for k in range(1, M + 1):
pi_cumulative.append(pi_cumulative[k - 1] + pi_list[k])
return pi_cumulative

def naive_tower_sample(pi_cumulative):
eta = random.uniform(0.0, 1.0)
for k in range(len(pi_cumulative)):
if eta < pi_cumulative[k]: break
return k

def levy_harmonic_path(dtau, N): #path sampling (to sample permutation positions)
beta = N * dtau
x_N = random.gauss(0.0, 1.0 / math.sqrt(2.0 * math.tanh(beta / 2.0)))
x = [x_N]
for k in range(1, N):
dtau_prime = (N - k) * dtau
Upsilon_1 = 1.0 / math.tanh(dtau) + 1.0 / math.tanh(dtau_prime)
Upsilon_2 = x[k - 1] / math.sinh(dtau) +  x_N / math.sinh(dtau_prime)
x_mean = Upsilon_2 / Upsilon_1
sigma = 1.0 / math.sqrt(Upsilon_1)
x.append(random.gauss(x_mean, sigma))
return x

### main program starts here ###
N = 8 #number of bosons
T_star = 0.1 #temperature
beta = 1.0 / N ** (1.0 / 3.0) / T_star
n_steps = 1000
Z = canonic_recursion(N, beta) #{N} boson partition function
for step in range(n_steps):
N_tmp = N #ancillary
x_config, y_config, z_config = [], [], [] #initialise the configurations in each 3 directions
while N_tmp > 0: #iterate through all particles
pi_sum = make_pi_list(Z, N_tmp)
k = naive_tower_sample(pi_sum)
x_config += levy_harmonic_path(beta, k)
y_config += levy_harmonic_path(beta, k)
z_config += levy_harmonic_path(beta, k)
N_tmp -= k #reduce the number of particles that are in the permutation cycle of length k


### Physical properties of the 1-dimensional classical and bosonic systems¶

• Consider 2 non-interacting distinguishable particles in a 1-dimensional harmonic trap:
In [5]:
import random, math, pylab

#There are only two possible cases: For k=1, we sample a single position (cycle of length 1),
#for k=2, we sample two positions (a cycle of length two).
def levy_harmonic_path(k):
x = [random.gauss(0.0, 1.0 / math.sqrt(2.0 * math.tanh(k * beta / 2.0)))] #direct-sample the first position
if k == 2:
Ups1 = 2.0 / math.tanh(beta)
Ups2 = 2.0 * x[0] / math.sinh(beta)
x.append(random.gauss(Ups2 / Ups1, 1.0 / math.sqrt(Ups1)))
return x[:]

def pi_x(x, beta):
sigma = 1.0 / math.sqrt(2.0 * math.tanh(beta / 2.0))
return math.exp(-x ** 2 / (2.0 * sigma ** 2)) / math.sqrt(2.0 * math.pi) / sigma

beta = 2.0
nsteps = 1000000
#initial sample has identity permutation
low = levy_harmonic_path(2) #tau=0
high = low[:] #tau=beta
data = []
for step in xrange(nsteps):
k = random.choice([0, 1])
low[k] = levy_harmonic_path(1)[0]
high[k] = low[k]
data.append(high[k])

list_x = [0.1 * a for a in range (-30, 31)]
y = [pi_x(a, beta) for a in list_x]
pylab.plot(list_x, y, linewidth=2.0, label='Exact distribution')
pylab.hist(data, normed=True, bins=80, label='QMC', alpha=0.5, color='green')
pylab.legend()
pylab.xlabel('$x$',fontsize=14)
pylab.ylabel('$\\pi(x)$',fontsize=14)
pylab.title('2 non-interacting distinguishable 1-d particles',fontsize=14)
pylab.xlim(-3, 3)
pylab.savefig('plot_A1_beta%s.png' % beta)

• Consider two non-interacting indistinguishable bosonic quantum particles in a one-dimensional harmonic trap:
In [6]:
import math, random, pylab, numpy as np

def z(beta):
return 1.0 / (1.0 - math.exp(- beta))

def pi_two_bosons(x, beta): #exact two boson position distribution
pi_x_1 = math.sqrt(math.tanh(beta / 2.0)) / math.sqrt(math.pi) * math.exp(-x ** 2 * math.tanh(beta / 2.0))
pi_x_2 = math.sqrt(math.tanh(beta)) / math.sqrt(math.pi) * math.exp(-x ** 2 * math.tanh(beta))
weight_1 = z(beta) ** 2 / (z(beta) ** 2 + z(2.0 * beta))
weight_2 = z(2.0 * beta) / (z(beta) ** 2 + z(2.0 * beta))
pi_x = pi_x_1 * weight_1 + pi_x_2 * weight_2
return pi_x

def levy_harmonic_path(k):
x = [random.gauss(0.0, 1.0 / math.sqrt(2.0 * math.tanh(k * beta / 2.0)))]
if k == 2:
Ups1 = 2.0 / math.tanh(beta)
Ups2 = 2.0 * x[0] / math.sinh(beta)
x.append(random.gauss(Ups2 / Ups1, 1.0 / math.sqrt(Ups1)))
return x[:]

def rho_harm_1d(x, xp, beta):
Upsilon_1 = (x + xp) ** 2 / 4.0 * math.tanh(beta / 2.0)
Upsilon_2 = (x - xp) ** 2 / 4.0 / math.tanh(beta / 2.0)
return math.exp(- Upsilon_1 - Upsilon_2)

beta = 2.0
list_beta = np.linspace(0.1, 5.0)
nsteps = 10000
low = levy_harmonic_path(2)
high = low[:]
fract_one_cycle_dat, fract_two_cycles_dat = [], []

for beta in list_beta:
one_cycle_dat = 0.0 #initialise the permutation fractions for each temperature
data = []
for step in xrange(nsteps):
# move 1 (direct-sample the positions)
if low[0] == high[0]: #if the cycle is of length 1
k = random.choice([0, 1])
low[k] = levy_harmonic_path(1)[0]
high[k] = low[k] #assures the cycle
else: #if the cycle is of length 2s
low[0], low[1] = levy_harmonic_path(2)
high[1] = low[0] #assures the cycle
high[0] = low[1]
one_cycle_dat += 1.0 / float(nsteps) #calculate the fraction of the single cycle cases
data += low[:] #save the position histogram data
# move 2 (Metropolis for sampling the permutations)
weight_old = (rho_harm_1d(low[0], high[0], beta) * rho_harm_1d(low[1], high[1], beta))
weight_new = (rho_harm_1d(low[0], high[1], beta) * rho_harm_1d(low[1], high[0], beta))
if random.uniform(0.0, 1.0) < weight_new / weight_old:
high[0], high[1] = high[1], high[0]

fract_one_cycle_dat.append(one_cycle_dat)
fract_two_cycles_dat.append(1.0 - one_cycle_dat) #save the fraction of the two cycles cases

#Exact permutation distributions for all temperatures
fract_two_cycles = [z(beta) ** 2 / (z(beta) ** 2 + z(2.0 * beta)) for beta in list_beta]
fract_one_cycle = [z(2.0 * beta) / (z(beta) ** 2 + z(2.0 * beta)) for beta in list_beta]

#Graphics output:
list_x = [0.1 * a for a in range (-30, 31)]
y = [pi_two_bosons(a, beta) for a in list_x]
pylab.plot(list_x, y, linewidth=2.0, label='Exact distribution')
pylab.hist(data, normed=True, bins=80, label='QMC', alpha=0.5, color='green')
pylab.legend()
pylab.xlabel('$x$',fontsize=14)
pylab.ylabel('$\\pi(x)$',fontsize=14)
pylab.title('2 non-interacting bosonic 1-d particles',fontsize=14)
pylab.xlim(-3, 3)
pylab.savefig('plot_A2_beta%s.png' % beta)
pylab.show()
pylab.clf()

fig = pylab.figure(figsize=(10, 5))

ax = fig.add_subplot(1, 2, 1)
ax.plot(list_beta, fract_one_cycle_dat, linewidth=4, label='QMC')
ax.plot(list_beta, fract_one_cycle, linewidth=2, label='exact')
ax.legend()
ax.set_xlabel('$\\beta$',fontsize=14)
ax.set_ylabel('$\\pi_2(\\beta)$',fontsize=14)
ax.set_title('Fraction of cycles of length 2',fontsize=14)

ax = fig.add_subplot(1, 2, 2)
ax.plot(list_beta, fract_two_cycles_dat, linewidth=4, label='QMC')
ax.plot(list_beta, fract_two_cycles, linewidth=2,label='exact')
ax.legend()
ax.set_xlabel('$\\beta$',fontsize=14)
ax.set_ylabel('$\\pi_1(\\beta)$',fontsize=14)
ax.set_title('Fraction of cycles of length 1',fontsize=14)

pylab.savefig('plot_A2.png')
pylab.show()
pylab.clf()

<Figure size 432x288 with 0 Axes>
<Figure size 432x288 with 0 Axes>

We can use dictionaries instead of lists. The implementation is in the following program.

Here we also calculate the correlation between the two particles, i.e. sample of the absolute distance $r$ between the two bosons. The comparison between the resulting distribution and the distribution for the distinguishable case corresponds to boson bunching (high weight for small distances between the bosons).

In [7]:
import math, random, pylab

def prob_r_distinguishable(r, beta): #the exact correlation function for two particles
sigma = math.sqrt(2.0) / math.sqrt(2.0 * math.tanh(beta / 2.0))
prob = (math.sqrt(2.0 / math.pi) / sigma) * math.exp(- r ** 2 / 2.0 / sigma ** 2)
return prob

def levy_harmonic_path(k):
x = [random.gauss(0.0, 1.0 / math.sqrt(2.0 * math.tanh(k * beta / 2.0)))]
if k == 2:
Ups1 = 2.0 / math.tanh(beta)
Ups2 = 2.0 * x[0] / math.sinh(beta)
x.append(random.gauss(Ups2 / Ups1, 1.0 / math.sqrt(Ups1)))
return x[:]

def rho_harm_1d(x, xp, beta):
Upsilon_1 = (x + xp) ** 2 / 4.0 * math.tanh(beta / 2.0)
Upsilon_2 = (x - xp) ** 2 / 4.0 / math.tanh(beta / 2.0)
return math.exp(- Upsilon_1 - Upsilon_2)

beta = 0.1
nsteps = 1000000
low_1, low_2 = levy_harmonic_path(2)
x = {low_1:low_1, low_2:low_2}
data_corr = []
for step in xrange(nsteps):
# move 1
a = random.choice(x.keys())
if a == x[a]:
dummy = x.pop(a)
a_new = levy_harmonic_path(1)[0]
x[a_new] = a_new
else:
a_new, b_new = levy_harmonic_path(2)
x = {a_new:b_new, b_new:a_new}
r = abs(x.keys()[1] - x.keys()[0])
data_corr.append(r)
# move 2
(low1, high1), (low2, high2) = x.items()
weight_old = rho_harm_1d(low1, high1, beta) * rho_harm_1d(low2, high2, beta)
weight_new = rho_harm_1d(low1, high2, beta) * rho_harm_1d(low2, high1, beta)
if random.uniform(0.0, 1.0) < weight_new / weight_old:
x = {low1:high2, low2:high1}

#Graphics output:
list_x = [0.1 * a for a in range (0, 100)]
y = [prob_r_distinguishable(a, beta) for a in list_x]
pylab.plot(list_x, y, linewidth=2.0, label='Exact distinguishable distribution')
pylab.hist(data_corr, normed=True, bins=120, label='Indistinguishable QMC', alpha=0.5, color='green')
pylab.legend()
pylab.xlabel('$r$',fontsize=14)
pylab.ylabel('$\\pi_{corr}(r)$',fontsize=14)
pylab.title('Correlation function of non-interacting 1-d bosons',fontsize=14)
pylab.xlim(0, 10)
pylab.savefig('plot_A3_beta%s.png' % beta)
pylab.show()
pylab.clf()

<Figure size 432x288 with 0 Axes>

### 3-dimensional bosons¶

#### Isotropic trap¶

In [264]:
import random, math, numpy, sys, os
import matplotlib.pyplot as plt

def harmonic_ground_state(x):
return math.exp(-x ** 2)/math.sqrt(math.pi)

def levy_harmonic_path_3d(k):
x0 = tuple([random.gauss(0.0, 1.0 / math.sqrt(2.0 *
math.tanh(k * beta / 2.0))) for d in range(3)])
x = [x0]
for j in range(1, k):
Upsilon_1 = 1.0 / math.tanh(beta) + 1.0 / \
math.tanh((k - j) * beta)
Upsilon_2 = [x[j - 1][d] / math.sinh(beta) + x[0][d] /
math.sinh((k - j) * beta) for d in range(3)]
x_mean = [Upsilon_2[d] / Upsilon_1 for d in range(3)]
sigma = 1.0 / math.sqrt(Upsilon_1)
dummy = [random.gauss(x_mean[d], sigma) for d in range(3)]
x.append(tuple(dummy))
return x

def rho_harm_3d(x, xp):
Upsilon_1 = sum((x[d] + xp[d]) ** 2 / 4.0 *
math.tanh(beta / 2.0) for d in range(3))
Upsilon_2 = sum((x[d] - xp[d]) ** 2 / 4.0 /
math.tanh(beta / 2.0) for d in range(3))
return math.exp(- Upsilon_1 - Upsilon_2)

N = 512
T_star = 0.8
list_T = numpy.linspace(0.8,0.1,5)
beta = 1.0 / (T_star * N ** (1.0 / 3.0))
cycle_min = 10
nsteps = 50000
data_x, data_y, data_x_l, data_y_l = [], [], [], []

for T_star in list_T:
# Initial condition
filename = 'data_boson_configuration_N%i_T%.1f.txt' % (N,T_star)
positions = {}
if os.path.isfile(filename):
f = open(filename, 'r')
for line in f:
a = line.split()
positions[tuple([float(a[0]), float(a[1]), float(a[2])])] = \
tuple([float(a[3]), float(a[4]), float(a[5])])
f.close()
if len(positions) != N:
sys.exit('ERROR in the input file.')
print 'starting from file', filename
else:
for k in range(N):
a = levy_harmonic_path_3d_anisotropic(1)
positions[a[0]] = a[0]
print 'Starting from a new configuration'

# Monte Carlo loop
for step in range(nsteps):
# move 1: resample one permutation cycle
boson_a = random.choice(positions.keys())
perm_cycle = []
while True:
perm_cycle.append(boson_a)
boson_b = positions.pop(boson_a)
if boson_b == perm_cycle[0]:
break
else:
boson_a = boson_b
k = len(perm_cycle)
data_x.append(boson_a[0])
data_y.append(boson_a[1])

if k > cycle_min:
data_x_l.append(boson_a[0])
data_y_l.append(boson_a[1])
perm_cycle = levy_harmonic_path_3d(k)
positions[perm_cycle[-1]] = perm_cycle[0]
for k in range(len(perm_cycle) - 1):
positions[perm_cycle[k]] = perm_cycle[k + 1]

# move 2: exchange
a_1 = random.choice(positions.keys())
b_1 = positions.pop(a_1)
a_2 = random.choice(positions.keys())
b_2 = positions.pop(a_2)
weight_new = rho_harm_3d(a_1, b_2) * rho_harm_3d(a_2, b_1)
weight_old = rho_harm_3d(a_1, b_1) * rho_harm_3d(a_2, b_2)
if random.uniform(0.0, 1.0) < weight_new / weight_old:
positions[a_1] = b_2
positions[a_2] = b_1
else:
positions[a_1] = b_1
positions[a_2] = b_2

f = open(filename, 'w')
for a in positions:
b = positions[a]
f.write(str(a[0]) + ' ' + str(a[1]) + ' ' + str(a[2]) + ' ' +
str(b[0]) + ' ' + str(b[1]) + ' ' + str(b[2]) + '\n')
f.close()

# Analyze cycles, do 3d plot
import pylab, mpl_toolkits.mplot3d

fig = pylab.figure()
ax = mpl_toolkits.mplot3d.axes3d.Axes3D(fig)
ax.set_aspect('equal')
n_colors = 10
list_colors = pylab.cm.rainbow(numpy.linspace(0, 1, n_colors))[::-1]
dict_colors = {}
i_color = 0
positions_copy = positions.copy()
while positions_copy:
x, y, z = [], [], []
starting_boson = positions_copy.keys()[0]
boson_old = starting_boson
while True:
x.append(boson_old[0])
y.append(boson_old[1])
z.append(boson_old[2])
boson_new = positions_copy.pop(boson_old)
if boson_new == starting_boson: break
else: boson_old = boson_new
len_cycle = len(x)
if len_cycle > 2:
x.append(x[0])
y.append(y[0])
z.append(z[0])
if len_cycle in dict_colors:
color = dict_colors[len_cycle]
ax.plot(x, y, z, '+-', c=color, lw=0.75)
else:
color = list_colors[i_color]
i_color = (i_color + 1) % n_colors
dict_colors[len_cycle] = color
ax.plot(x, y, z, '+-', c=color, label='k=%i' % len_cycle, lw=0.75)
pylab.title(str(N) + ' bosons at T* = ' + str(T_star))
pylab.legend()
ax.set_xlabel('$x$', fontsize=16)
ax.set_ylabel('$y$', fontsize=16)
ax.set_zlabel('$z$', fontsize=16)
xmax = 6.0
ax.set_xlim3d([-xmax, xmax])
ax.set_ylim3d([-xmax, xmax])
ax.set_zlim3d([-xmax, xmax])
pylab.savefig('plot_boson_configuration_N%i_T%.1f.png' %(N,T_star))
pylab.show()
pylab.clf()

#Plot the histograms
list_x = [0.1 * a for a in range (-50, 51)]
y = [harmonic_ground_state(a) for a in list_x]
pylab.plot(list_x, y, linewidth=2.0, label='Ground state')
pylab.hist(data_x, normed=True, bins=120, alpha = 0.5, label='All bosons')
pylab.hist(data_x_l, normed=True, bins=120, alpha = 0.5, label='Bosons in longer cycle')
pylab.xlim(-3.0, 3.0)
pylab.xlabel('$x$',fontsize=14)
pylab.ylabel('$\pi(x)$',fontsize=14)
pylab.title('3-d non-interacting bosons $x$ distribution $N= %i$, $T= %.1f$' %(N,T_star))
pylab.legend()
pylab.savefig('position_distribution_N%i_T%.1f.png' %(N,T_star))
pylab.show()
pylab.clf()

plt.hist2d(data_x_l, data_y_l, bins=40, normed=True)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.title('The distribution of the $x$ and $y$ positions')
plt.colorbar()
plt.xlim(-3.0, 3.0)
plt.ylim(-3.0, 3.0)
plt.show()

starting from file data_boson_configuration_N512_T0.8.txt

starting from file data_boson_configuration_N512_T0.6.txt

<Figure size 432x288 with 0 Axes>
starting from file data_boson_configuration_N512_T0.5.txt

<Figure size 432x288 with 0 Axes>
starting from file data_boson_configuration_N512_T0.3.txt

<Figure size 432x288 with 0 Axes>
starting from file data_boson_configuration_N512_T0.1.txt

<Figure size 432x288 with 0 Axes>

#### Anisotropic trap¶

We can imitate the experiments that imitate 1-d bosons in cigar shaped anisotropic harmonic traps, and 2-d bosons in pancake shaped anisotropic harmonic traps.

In [8]:
%pylab inline
import random, math, numpy, os, sys

def levy_harmonic_path_3d_anisotropic(k, omega):
sigma = [1.0 / math.sqrt(2.0 * omega[d] *
math.tanh(0.5 * k * beta * omega[d])) for d in xrange(3)]
xk = tuple([random.gauss(0.0, sigma[d]) for d in xrange(3)])
x = [xk]
for j in range(1, k):
Upsilon_1 = [1.0 / math.tanh(beta * omega[d]) +
1.0 / math.tanh((k - j) * beta * omega[d]) for d in range(3)]
Upsilon_2 = [x[j - 1][d] / math.sinh(beta * omega[d]) + \
xk[d] / math.sinh((k - j) * beta * omega[d]) for d in range(3)]
x_mean = [Upsilon_2[d] / Upsilon_1[d] for d in range(3)]
sigma = [1.0 / math.sqrt(Upsilon_1[d] * omega[d]) for d in range(3)]
dummy = [random.gauss(x_mean[d], sigma[d]) for d in range(3)]
x.append(tuple(dummy))
return x

def rho_harm_3d_anisotropic(x, xp, beta, omega):
Upsilon_1 = sum(omega[d] * (x[d] + xp[d]) ** 2 / 4.0 *
math.tanh(beta * omega[d] / 2.0) for d in range(3))
Upsilon_2 = sum(omega[d] * (x[d] - xp[d]) ** 2 / 4.0 /
math.tanh(beta * omega[d] / 2.0) for d in range(3))
return math.exp(- Upsilon_1 - Upsilon_2)

omegas = numpy.array([[4.0, 4.0, 1.0], [1.0, 5.0, 1.0]])

for i in range(len(omegas[:,1])):
N = 512
nsteps = 100000
omega_harm = 1.0
omega = omegas[i,:]
for d in range(3):
omega_harm *= omega[d] ** (1.0 / 3.0)
T_star = 0.5
T = T_star * omega_harm * N ** (1.0 / 3.0)
beta = 1.0 / T
print 'omega: ', omega
# Initial condition
if i == 0:
filename = 'data_boson_configuration_anisotropic_N%i_T%.1f_cigar.txt' % (N,T_star)
elif i == 1:
filename = 'data_boson_configuration_anisotropic_N%i_T%.1f_pancake.txt' % (N,T_star)
positions = {}
if os.path.isfile(filename):
f = open(filename, 'r')
for line in f:
a = line.split()
positions[tuple([float(a[0]), float(a[1]), float(a[2])])] = \
tuple([float(a[3]), float(a[4]), float(a[5])])
f.close()
if len(positions) != N:
sys.exit('ERROR in the input file.')
print 'starting from file', filename
else:
for k in range(N):
a = levy_harmonic_path_3d_anisotropic(1,omega)
positions[a[0]] = a[0]
print 'Starting from a new configuration'
for step in range(nsteps):
boson_a = random.choice(positions.keys())
perm_cycle = []
while True:
perm_cycle.append(boson_a)
boson_b = positions.pop(boson_a)
if boson_b == perm_cycle[0]: break
else: boson_a = boson_b
k = len(perm_cycle)
perm_cycle = levy_harmonic_path_3d_anisotropic(k,omega)
positions[perm_cycle[-1]] = perm_cycle[0]
for j in range(len(perm_cycle) - 1):
positions[perm_cycle[j]] = perm_cycle[j + 1]
a_1 = random.choice(positions.keys())
b_1 = positions.pop(a_1)
a_2 = random.choice(positions.keys())
b_2 = positions.pop(a_2)
weight_new = (rho_harm_3d_anisotropic(a_1, b_2, beta, omega) *
rho_harm_3d_anisotropic(a_2, b_1, beta, omega))
weight_old = (rho_harm_3d_anisotropic(a_1, b_1, beta, omega) *
rho_harm_3d_anisotropic(a_2, b_2, beta, omega))
if random.uniform(0.0, 1.0) < weight_new / weight_old:
positions[a_1], positions[a_2] = b_2, b_1
else:
positions[a_1], positions[a_2] = b_1, b_2

f = open(filename, 'w')
for a in positions:
b = positions[a]
f.write(str(a[0]) + ' ' + str(a[1]) + ' ' + str(a[2]) + ' ' +
str(b[0]) + ' ' + str(b[1]) + ' ' + str(b[2]) + '\n')
f.close()

import pylab, mpl_toolkits.mplot3d
fig = pylab.figure()
ax = mpl_toolkits.mplot3d.axes3d.Axes3D(fig)
ax.set_aspect('equal')
n_colors = 10
list_colors = pylab.cm.rainbow(numpy.linspace(0, 1, n_colors))[::-1]
dict_colors = {}
i_color = 0
positions_copy = positions.copy()
while positions_copy:
x, y, z = [], [], []
starting_boson = positions_copy.keys()[0]
boson_old = starting_boson
while True:
x.append(boson_old[0])
y.append(boson_old[1])
z.append(boson_old[2])
boson_new = positions_copy.pop(boson_old)
if boson_new == starting_boson: break
else: boson_old = boson_new
len_cycle = len(x)
if len_cycle > 2:
x.append(x[0])
y.append(y[0])
z.append(z[0])
if len_cycle in dict_colors:
color = dict_colors[len_cycle]
ax.plot(x, y, z, '+-', c=color, lw=0.75)
else:
color = list_colors[i_color]
i_color = (i_color + 1) % n_colors
dict_colors[len_cycle] = color
ax.plot(x, y, z, '+-', c=color, label='k=%i' % len_cycle, lw=0.75)
pylab.legend()
ax.set_xlabel('$x$', fontsize=16)
ax.set_ylabel('$y$', fontsize=16)
ax.set_zlabel('$z$', fontsize=16)
xmax = 8.0
ax.set_xlim3d([-xmax, xmax])
ax.set_ylim3d([-xmax, xmax])
ax.set_zlim3d([-xmax, xmax])
if i == 0:
pylab.title(str(N) + ' bosons at T* = ' + str(T_star) + ' cigar potential')
pylab.savefig('position_distribution_N%i_T%.1f_cigar.png' %(N,T_star))
elif i == 1:
pylab.title(str(N) + ' bosons at T* = ' + str(T_star) + ' pancake potential')
pylab.savefig('position_distribution_N%i_T%.1f_pancake.png' %(N,T_star))
pylab.show()

Populating the interactive namespace from numpy and matplotlib
omega:  [4. 4. 1.]
starting from file data_boson_configuration_anisotropic_N512_T0.5_cigar.txt

omega:  [1. 5. 1.]
starting from file data_boson_configuration_anisotropic_N512_T0.5_pancake.txt


There it is found that the critical temperature for Bose-Einstein condensation is around $T^*\sim 0.9$.

## To do:¶

• Calculate the pair correlation function
In [ ]: