Assignment 3, Question 3

nmax:=50;
hbar:=1;
ω:=1;
ϵ:=(10^(-3));
mass:=1;
NN:=(10^5);
kbTcrit:=0.94 hbar ω (NN^(1/3)); (*Eq. (3.14) *)
kbTemp:= 0.5 kbTcrit;
Ncond:=NN(1-((((kbTemp/kbTcrit)))^3));
tildeμ:=0;
σ:=Sqrt[(hbar/(mass ω))];
EE[n_]:=hbar ω n;
calN:=(1/(((π (σ^2)))^(1/4)));
Nsamples:=32;
rmax:=8σ;
dr := rmax/Nsamples;
mbar[n_]:=(1/(Exp[ (EE[n] -tildeμ)/kbTemp]-1));
u[n_,l_,r_]:=Sqrt[(((2^(n+2))(n-l)!)/((π^(1/2))(2n+1)!!))] (1/(σ^(l+3))) Exp[-((r^2)/(2 (σ^2)))](r^l)LaguerreL[(n-l),l+(1/2),((r^2)/( (σ^2)))];
DensThermal:=Table[{dr mm,Sum[Sum[If[EvenQ[nn-ll],mbar[nn] ((u[nn,ll,dr mm])^2)(((2 ll + 1))/((4π))),0],{ll,0,nn}],{nn,1,nmax}]},{mm,1,Nsamples}]
DensBEC[r_]:=Ncond (calN^6) Exp[-((r^2)/( (σ^2)))]

In this part we check if we added enough excited states to accomodate the uncondensed atoms (nmax) . However we allow ourselves to be stingy since it takes too long for the density otherwise:

NN-Ncond

Sum of all excited state populations :

Sum[Sum[If[EvenQ[nn-ll],mbar[nn](2 ll + 1),0],{ll,0,nn}],{nn,1,nmax}]

Plot[{DensBEC[r]},{r,0,4 σ},AxesLabel ->{Style["r",FontSize->18,Black], Style["ρ(r)", FontSize->18,Black]}]

plotthermal=ListPlot[DensThermal,PlotRange-> All,AxesLabel ->{Style["r",FontSize->18,Black], Style["ρ(r)", FontSize->18,Black]}]

If you compare this with Fig . 2.2  in Pethik and Smith . Our solution is somehow cut off in the cusp near r -> 0 due to the limited number of modes taken .

ampl=20;

FitPlot=Plot[ampl Exp[-((r^2)/( (2 kbTemp/(mass (ω^2))) ))],{r,0, rmax}]