{VERSION 5 0 "IBM INTEL NT" "5.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 
1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 
0 0 1 }{CSTYLE "2D Input" 2 19 "" 0 1 255 0 0 1 0 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 256 "" 1 16 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
257 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 0 0 1 0 0 0 
0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }
{CSTYLE "" -1 261 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 
262 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 
{CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 
0 1 0 1 0 2 2 0 1 }}
{SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT 256 64 "Classical kinetic theory \+
for evaporative cooling in an atom trap" }}{PARA 0 "" 0 "" {TEXT -1 0 
"" }}{PARA 0 "" 0 "" {TEXT -1 1979 "Classical statistical mechanics al
lows one to describe an ensemble of atoms which move in an external po
tential, and which interact with each other by means of elastic collis
ions. During these binary collisions energy is exchanged between the a
toms. The common motion of many non-interacting particles in an extern
al field is described by the classical Liouville equation, which repre
sents a reformulation of the many Newton equations in terms of a singl
e phase-space distribution functions. When an elastic collision integr
al is added to this equation it becomes the Boltzmann equation. The Li
ouville or Boltzmann equations are formulated for continuous distribut
ion functions, but practical realizations are often performed by retur
ning to discrete testparticle equations. These are Newton (Hamilton) e
quations for particles moving in a common potential - which may includ
e a long-distance interaction between the particles - supplemented by \+
a collision term which represents the elastic exchange of energy among
 particle pairs. The Boltzmann equation is interesting insofar as it d
isplays irreversible behaviour caused by the collision term: the Liouv
ille dynamics is based on Newton's equations which are symmetric under
 a reversal of time; the collision integral, however, introduces chang
es to the particle motions such that different initial distributions e
volve towards Maxwell-Boltzmann type equilibrium distributions. Obviou
sly, it cannot work in reverse, i.e., it cannot undo the relaxations a
nd turn an equilibrium distribution (valid at large times) back into t
he original distribution by going backwards in time (the non-uniquenes
s alone would prohibit that). In that sense it contains the arrow of t
ime as observed in nature (2nd law of thermodynamics), but at the same
 time it remains a heuristic equation. For a discussion of this and ot
her features it is worthwhile to consult texts on statistical mechanic
s, a good starting point would be: Gregory H. Wannier: " }{TEXT 258 
19 "Statistical Physics" }{TEXT -1 58 " (paperback reprint by Dover fr
om the 1966 Wiley edition)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 
0 "" 0 "" {TEXT -1 552 "This kinetic motion approach allows one to und
erstand how a sample acquires an equilibrium state through energy exch
ange by elastic collisions. We can start the atoms all with equal velo
city, or we can equidistribute the energies initially over all possibl
e states: when we wait long enough an equilibrium state will develop w
hich is characterized by the average kinetic energy available to the e
nsemble. The steady-state or equilibrium distribution produced by the \+
elastic collisions between the atoms is the well-known Maxwell-Boltzma
nn distribution." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 1038 "Collisions have the tendency to turn a pair of two part
icles of similar-energy into a pair in which one partner is slow, and \+
one is faster. In the present worksheet this is demonstrated by simula
ting evaporative cooling. The idea is to remove the fast particles (in
 evaporative cooling in the coffee or tea cup the hottest molecules le
ave the fluid in the form of steam, thereby leaving an ensemble which \+
is colder on average). In an atom trap one can have a finite barrier h
eight which separates the trap from the environment. The idea is that \+
once an atom has acquired enough energy to go over the barrier it is g
one for good. If we keep the barrier height fixed, the cooling will wo
rk up to some point: once the whole ensemble has cooled off such that \+
the average kinetic energy is well below the barrier height, it is unl
ikely that further eleastic collisions will produce many more atoms of
 sufficient energy to make it over the barrier. Therefore, efficient c
ooling requires a lowering of the trap height as the ensemble cools do
wn." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 640 "O
f course we cannot integrate the full Boltzmann equation in Maple for \+
a three-dimensional system. In fact, a substantial simplification is p
ossible when the system is ergodic. This means that the system visits \+
all energetically accessible points in phase space if given sufficient
 time, and that different sequences of the time history selected at ra
ndom look similar to each other. In such a system it is not important \+
to know (in order to characterize the state) how the convection of the
 particle flow in phase space looks in detail. One can characterize th
e system by its energy dependence alone. The phase-space distribution \+
function " }{TEXT 19 8 "f(r,p,t)" }{TEXT -1 51 " is replaced by the pr
oduct of a density of states " }{TEXT 19 6 "rho(E)" }{TEXT -1 28 " tim
es the level population " }{TEXT 19 4 "f(E)" }{TEXT -1 24 ". The densi
ty of states " }{TEXT 19 6 "rho(E)" }{TEXT -1 242 " is usually a simpl
e function of the potential. Commonly one employes the technique of co
unting the allowed quantum energy levels (taking into account their de
generacy in a multi-dimensional setting such as the spherical harmonic
 oscillator)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 265 "The details of reducing the six-dimensional phase-space \+
Boltzmann equation to a one-dimensional energy-space equation are give
n, e.g., in the paper by J. Walraven et al., Phys. Rev. A 53, 381(1996
). The idea is to discretize the energy range into bins according to \+
" }{TEXT 19 12 "E=(i-1/2)*dE" }{TEXT -1 6 " with " }{TEXT 19 6 "i=1..N
" }{TEXT -1 131 " covering predominantly the energy range from the tra
p bottom up to the barrier height. The time evolution of the level pop
ulation " }{TEXT 19 10 "f_i=f(E_i)" }{TEXT -1 113 " is determined by a
 collision term which involves the product of level populations for tw
o arbitrary states with " }{TEXT 19 3 "E_k" }{TEXT -1 5 " and " }
{TEXT 19 3 "E_l" }{TEXT -1 26 " to connect to the states " }{TEXT 19 
3 "E_i" }{TEXT -1 15 " combined with " }{TEXT 19 3 "E_j" }{TEXT -1 37 
". While determining the evolution of " }{TEXT 19 6 "f(E_i)" }{TEXT 
-1 50 " one sums over all possible k,l while restricting " }{TEXT 19 
15 "E_j=E_k+E_l-E_i" }{TEXT -1 275 ", thus ensuring energy conservatio
n. The collision term is weighted with the probability for collisions \+
to take place (the cross section which is assumed to be energy-indepen
dent at low energies enters). The equation is also weighted by the den
sity of states: on the left by " }{TEXT 19 5 "rho_i" }{TEXT -1 22 ", a
nd on the right by " }{TEXT 19 5 "rho_h" }{TEXT -1 25 ", where the ene
rgy index " }{TEXT 19 14 "h=min(i,j,k,l)" }{TEXT -1 241 ". The selecti
on of the energy index on the right follows from the reduction of the \+
collision integral when the phase-space probability is replaced by the
 energy-level population. Schematically we arrive at a system of diffe
rential equations:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT 19 72 "rho_i diff(f_i,t) = c*add(rho_h*(f_k*f_l - f_i*f_j) , k=1
..N ,l=1..N)   " }{TEXT 2 5 "where" }{TEXT 19 25 "  h=min(i,j,k,l), j=
k+l-i" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 26 "
Here the level population " }{TEXT 19 3 "f_i" }{TEXT -1 74 " is dimens
ionless and the constant c contains the dimensionful quantities " }
{TEXT 19 29 "c=m*sigma*dE^2/(Pi^2*hbar^3) " }{TEXT -1 10 "such that " 
}{TEXT 19 4 "c*dt" }{TEXT -1 83 " represents a measure of the total nu
mber of collisions that occured. Furthermore, " }{TEXT 19 16 "rho_i = \+
rho(E_i)" }{TEXT -1 109 "  is usually a power-law density of states, i
.e., one can show that for the isotropic 3D harmonic oscillator " }
{TEXT 19 5 "rho_i" }{TEXT -1 11 " goes like " }{TEXT 19 5 "E_i^2" }
{TEXT -1 36 ", since more generally it goes like " }{TEXT 19 11 "E_i^(
1/2+d)" }{TEXT -1 24 " for potentials of form " }{TEXT 19 12 "V(r)=a*r
^3/d" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "
" {TEXT -1 58 "This system of differential equations can be solved usi
ng " }{TEXT 19 15 "dsolve[numeric]" }{TEXT -1 289 ", but it is more st
raightforward to apply an extrapolation (Euler) method, i.e., to discr
etize time. To illustrate the method we assume that the energy scale i
s set by the trap depth (in our scaled dimensionless variables we dist
ribute intially all particles equally over the energy range " }{TEXT 
19 6 "E=0..1" }{TEXT -1 77 "). At first we will not remove particles, \+
i.e., we start the simulation with " }{TEXT 19 3 "f_i" }{TEXT -1 41 " \+
representing a step function cut off at " }{TEXT 19 3 "E=1" }{TEXT -1 
60 ", and just watch the evolution of this distribution in time." }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "restart; with(plots):" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 91 "We introduce parameters for the di
mensionless trap height, the total number of energy bins " }{TEXT 19 
1 "N" }{TEXT -1 78 ", and the number of energy bins devoted to the ran
ge above the barrier height " }{TEXT 19 2 "N0" }{TEXT -1 1 ":" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "E_t:=1;" }}}{EXCHG {PARA 0 ">
 " 0 "" {MPLTEXT 1 0 6 "N:=40;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 7 "N0:=10;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "dE:=E_t/(
N-N0);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 162 "The bin energies and t
he density of states are now specified. We begin by choosing the densi
ty of states as appropriate for a harmonic oscillator in 3 dimensions.
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "E:=i->(i-0.5)*dE;" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "rho:=i->E(i)^2;" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "plot([seq([E(i),rho(i)],i=1..N)]);
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 76 "This is the level density of \+
a potential that continues harmonically beyond " }{TEXT 19 3 "E=1" }
{TEXT -1 60 ". We choose all energy levels to be equally populated bel
ow " }{TEXT 19 3 "E=1" }{TEXT -1 24 ", and to be empty above." }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "fv:=array([seq(1,i=1..N-N0),
seq(0,i=N-N0+1..N)]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 173 "The dis
tribution over the energies has to take the density of states into acc
ount. Note that we should really be displaying a histogram - therefore
 the step is not sharp at " }{TEXT 19 3 "E=1" }{TEXT -1 1 "." }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "PL0:=plot([seq([E(i),(rho(i)
*fv[i])],i=1..N)]): display(PL0);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 
15 "The product of " }{TEXT 19 4 "c*dt" }{TEXT -1 116 " has to be smal
l for Euler's extrapolation method to work. The algorithm updates the \+
level population according to: " }{TEXT 19 33 "f_i(t+dt) = f_i(t) + dt
*(df_i/dt)" }{TEXT -1 67 ". Implementation details: the sums on the ri
ght are truncated when " }{TEXT 19 1 "j" }{TEXT -1 176 " becomes negat
ive or exceeds the energy range. The update loop at the end is easily \+
modified to set the distribution to zero for energies above the half-w
ay point, i.e., above " }{TEXT 19 3 "E=1" }{TEXT -1 1 "." }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "cdt:=1/1000;" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 11 "N_col:=100;" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 28 "for icol from 1 to N_col do:" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 21 "for i from 1 to N do:" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 7 "RHS:=0:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "for k from 1 to N
 do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "for l from 1 to N do:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "j:=k+l-i:" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 21 "if j>0 and j<N+1 then" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 66 "h:=min(i,j,k,l): RHS:=RHS+rho(h)/rho(i)*(fv[k]*fv[l]-fv[i]*fv[
j]):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "fi:" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 7 "od: od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "if i<=N \+
then fv[i]:=fv[i]+cdt*RHS: else fv[i]:=0: fi: od:" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 152 "if trunc(icol/10)*10=icol then PL[trunc(icol/10)]:
=plot([seq([E(i),(rho(i)*fv[i])],i=1..N)],title=cat(\"Collision number
 : \",(N/2)*icol*cdt/dE)): fi: od:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 56 "display([PL0,seq(PL[i],i=1..N_col/10)],insequence=tru
e);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 249 "We observe that the atoms
 distribute themselves over the available energy states and reach what
 looks like a Maxwell-Boltzmann distribution that corresponds to some \+
temperature (the distribution peaks at some energy, and has a well-def
ined average)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "E_avg:=add(
E(i)*fv[i]*rho(i),i=1..N)/add(fv[i]*rho(i),i=1..N);" }}}{EXCHG {PARA 
0 "" 0 "" {TEXT -1 112 "Clearly the above result does not provide the \+
correct average, because the distribution has been chopped off at " }
{TEXT 19 10 "(N-1/2)*dE" }{TEXT -1 146 ". In fact, the truncation of t
he distribution at the tail is the key idea behind evaporative cooling
. We will introduce now a removal of atoms at " }{TEXT 19 6 "E(i)>1" }
{TEXT -1 82 ", and look for long-term behaviour. We start with the sam
e distribution as before:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
49 "fv:=array([seq(1,i=1..N-N0),seq(0,i=N-N0+1..N)]);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "PL0:=plot([seq([E(i),(rho(i)*fv[i])
],i=1..N)]):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 78 "We increase the l
ength of the loop and in the update of the energy population " }{TEXT 
19 5 "fv[i]" }{TEXT -1 88 " we set the distribution function to zero f
or atoms above the barrier height, i.e., for " }{TEXT 19 3 "E>1" }
{TEXT -1 213 ". This corresponds to removal of atoms which have acquir
ed such energies. This idea is valid on physical grounds provided the \+
time between collision is long enough for evaporating atoms to disappe
ar from the trap." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "N_col:
=100;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "for icol from 1 to
 N_col do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "for i from 1 to N do:
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "RHS:=0:" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 21 "for k from 1 to N do:" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 21 "for l from 1 to N do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "j:=
k+l-i:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "if j>0 and j<N+1 then" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "h:=min(i,j,k,l): RHS:=RHS+rho(h)/rh
o(i)*(fv[k]*fv[l]-fv[i]*fv[j]):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "f
i:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "od: od:" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 60 "if i<=N-N0 then fv[i]:=fv[i]+cdt*RHS: else fv[i]:=0: \+
fi: od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 152 "if trunc(icol/10)*10=ic
ol then PL[trunc(icol/10)]:=plot([seq([E(i),(rho(i)*fv[i])],i=1..N)],t
itle=cat(\"Collision number : \",(N/2)*icol*cdt/dE)): fi: od:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 75 "display([PL0,seq(PL[i],i=1..
N_col/10)],insequence=true,view=[0..1,0..0.8]);" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 62 "E_avg:=add(E(i)*fv[i]*rho(i),i=1..N)/add(fv[i]
*rho(i),i=1..N);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 454 "One observes
 that the distribution changes rapidly for low collision numbers into \+
a characteristic shape with a maximum level population inside the trap
. With increasing collision number the maximum shifts towards lower en
ergies, but it becomes a rather slow process in the long term, because
 it becomes less likely that particles be evaporated over the high bar
rier. A solution to this problem is to reduce the barrier height as th
e ensemble cools down." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 
0 "" {TEXT -1 225 "The calculated daverage energy is now a more realis
tic measure given that the distribution is bounded over the interval. \+
Before we proceed with the idea of lowering the barrier height we impr
ove the algorithm in two respects." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT 257 11 "Exercise 1:" }}{PARA 0 "" 0 "" {TEXT -1 
174 "Increase the number of collisions by running the loop over longer
 times and observe the modest change in the final distribution. Also p
erform small variations of the product " }{TEXT 19 4 "c*dt" }{TEXT -1 
75 " and verify whether your results are consistent for fixed collisio
n number." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 58 "A good ques
tion to ask is: How much would we gain from an " }{TEXT 19 7 "O(dt^2)
" }{TEXT -1 107 " method? Twice the computational effort per dt-step m
ight be compensated by 4-fold improvement in accuracy." }}{PARA 0 "" 
0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "restar
t; with(plots):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "E_t:=1;" 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "N:=40;" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 7 "N0:=10;" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 15 "dE:=E_t/(N-N0);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 17 "E:=i->(i-0.5)*dE;" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 15 "rho:=i->E(i)^2;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 
55 "We need to farm out the RHS-calculation to a procedure:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "RHSpr:=proc(fv,i) local j,k,
l,h,RHS; global N,cdt,rho;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "RHS:=0
:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "for k from 1 to N do:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "for l from 1 to N do:" }}{PARA 0 ">
 " 0 "" {MPLTEXT 1 0 9 "j:=k+l-i:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
21 "if j>0 and j<N+1 then" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "h:=min
(i,j,k,l): RHS:=RHS+rho(h)/rho(i)*(fv[k]*fv[l]-fv[i]*fv[j]):" }}{PARA 
0 "> " 0 "" {MPLTEXT 1 0 3 "fi:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "
od: od: RHS; end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "N_col:
=100;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "cdt:=1/400;" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "fv:=array([seq(1,i=1..N-N0),
seq(0,i=N-N0+1..N)]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "PL
0:=plot([seq([E(i),(rho(i)*fv[i])],i=1..N)]):" }}}{EXCHG {PARA 0 "" 0 
"" {TEXT -1 168 "The collision loop is modified by introducing a simpl
e predictor-corrector method: extrapolation is used as before in the f
irst step to predict the level population at " }{TEXT 19 4 "t+dt" }
{TEXT -1 33 ", and then the rate of change in " }{TEXT 19 5 "fv[i]" }
{TEXT -1 134 " is recomputed from this predicted population. The real \+
time-step is carried out using an average of both rate of change calcu
lations." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "for icol from 1
 to N_col do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "for i from 1 to N \+
do: RHSp[i]:=RHSpr(fv,i):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "if i<=
N-N0 then gv[i]:=fv[i]+cdt*RHSp[i]: else gv[i]:=0: fi: od:" }}{PARA 0 
"> " 0 "" {MPLTEXT 1 0 53 "for i from 1 to N do: RHS:=0.5*(RHSpr(gv,i)
+RHSp[i]):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "if i<=N-N0 then fv[i]
:=fv[i]+cdt*RHS: else fv[i]:=0: fi: od:" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 159 "if trunc(icol/10)*10=icol then PL[trunc(icol/10)]:=plot([seq(
[E(i),(rho(i)*fv[i])],i=1..N)],title=cat(\"Collision number : \",trunc
((N/2)*icol*cdt/dE))): fi: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 75 "display([PL0,seq(PL[i],i=1..N_col/10)],insequence=true,view=[0
..1,0..0.5]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "E_avg:=add
(E(i)*fv[i]*rho(i),i=1..N)/add(fv[i]*rho(i),i=1..N);" }}}{EXCHG {PARA 
0 "" 0 "" {TEXT 259 11 "Exercise 2:" }}{PARA 0 "" 0 "" {TEXT -1 72 "Ex
plore the stability regime for the improved algorithm, i.e., increase \+
" }{TEXT 19 3 "cdt" }{TEXT -1 68 " and observe at what collision numbe
r the solution becomes unstable." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT 260 11 "Exercise 3:" }}{PARA 0 "" 0 "" {TEXT -1 
211 "Explore the dependence of the results on the number of energy bin
s kept above the trap barrier. In order to make a quantitative compari
son track the maximum of the distribution as a function of collision n
umber." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 254 "To accomplish cooling to low temperatures (i.e., energie
s in the trap) the method becomes inefficient. The obvious strategy is
 to lower the trap height, i.e., to allow particles of lower energy to
 evaporate as time goes on. We lower the trap height from " }{TEXT 19 
7 "E(N-10)" }{TEXT -1 4 " to " }{TEXT 19 6 "E(N/4)" }{TEXT -1 69 " ove
r the time of the run. This is accomplished by defining an index " }
{TEXT 19 6 "N_trap" }{TEXT -1 150 " which varies with time (collision \+
number), and which is used in the update algorithm as the cut-off para
meter beyond which the distribution function " }{TEXT 19 5 "fv[i]" }
{TEXT -1 49 " is set to zero at the end of a propagation step." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 199 "
First we improve the algorithm by avoiding bookkeeping for the over-th
e-barrier states. We discretize the energy range for the (initally) tr
apped states and propagate only those. The double sum over " }{TEXT 
19 3 "k,l" }{TEXT -1 64 " in the calculation of the right hand side is
 extended to cover " }{TEXT 19 2 "N0" }{TEXT -1 26 " additional energy
 levels." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "restart; with(p
lots):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "E_t:=1;" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "N:=40;" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 12 "dE:=E_t/(N);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 17 "E:=i->(i-0.5)*dE;" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 15 "rho:=i->E(i)^2;" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 11 "cdt:=1/400;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 10 "N_col:=50;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "RHSpr:=
proc(fv,i) local j,k,l,h,RHS,fvk,fvl,fvj,N0; global N,cdt,rho;" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "RHS:=0: N0:=15;" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 66 "for k from 1 to N+N0 do: if k<=N then fvk:=fv[k]: e
lse fvk:=0: fi:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "for l from 1 to \+
N+N0 do: if l<=N then fvl:=fv[l]: else fvl:=0: fi:" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 9 "j:=k+l-i:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "if j
>0 and j<=N+N0 then if j<=N then fvj:=fv[j]: else fvj:=0: fi:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "h:=min(i,j,k,l): RHS:=RHS+rho(h)/rh
o(i)*(fvk*fvl-fv[i]*fvj):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "fi:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "od: od: RHS; end:" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 27 "fv:=array([seq(1,i=1..N)]);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "PL0:=plot([seq([E(i),(rho(i)*fv[i])
],i=1..N)]):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "for icol fr
om 1 to N_col do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "for i from 1 t
o N do: RHSp[i]:=RHSpr(fv,i):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "gv
[i]:=fv[i]+cdt*RHSp[i]: od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "for \+
i from 1 to N do: RHS:=0.5*(RHSpr(gv,i)+RHSp[i]):" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 25 "fv[i]:=fv[i]+cdt*RHS: od:" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 159 "if trunc(icol/10)*10=icol then PL[trunc(icol/10)]:=p
lot([seq([E(i),(rho(i)*fv[i])],i=1..N)],title=cat(\"Collision number :
 \",trunc((N/2)*icol*cdt/dE))): fi: od:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 75 "display([PL0,seq(PL[i],i=1..N_col/10)],insequence=tru
e,view=[0..1,0..0.5]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "E
_avg:=add(E(i)*fv[i]*rho(i),i=1..N)/add(fv[i]*rho(i),i=1..N);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "N_TP:=add(fv[i]*rho(i),i=1..
N);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 240 "It is fairly straightforw
ard now to explore the following questions: how does the average energ
y after M collisions depend on the number states kept in the calculati
on above the trap barrier? How sensitive are the results to the resolu
tion " }{TEXT 19 2 "dE" }{TEXT -1 70 "? How do the results differ for \+
other powers in the density of states?" }}{PARA 0 "" 0 "" {TEXT -1 0 "
" }}{PARA 0 "" 0 "" {TEXT -1 238 "Now we come to the main point: we wi
ll lower the barrier height for the trap as a function of time (actual
ly as a function of the number of collisions between the atoms). Over \+
the total number of collisions we reduce the trap height from " }
{TEXT 19 6 "E(N)=1" }{TEXT -1 4 " to " }{TEXT 19 5 "E(N1)" }{TEXT -1 
61 ". The index which is used as a criterion as to how to update " }
{TEXT 19 5 "fv[i]" }{TEXT -1 11 " is called " }{TEXT 19 6 "N_trap" }
{TEXT -1 131 " and is calculated in the collision loop to change gradu
ally. For a start we choose to lower the barrier height by one half, i
.e., " }{TEXT 19 6 "N1=N/2" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 8 "N1:=N/2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
11 "N_col:=100;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "cdt:=1/8
00;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "fv:=array([seq(1,i=1
..N)]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "PL0:=plot([seq([
E(i),(rho(i)*fv[i])],i=1..N)]):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 66 "for icol from 1 to N_col do: N_trap:=N-trunc((N-N1)*(icol/N_co
l));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "for i from 1 to N do: RHSp[
i]:=RHSpr(fv,i):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "if i<=N_trap th
en gv[i]:=fv[i]+cdt*RHSp[i]: else gv[i]:=0: fi: od:" }}{PARA 0 "> " 0 
"" {MPLTEXT 1 0 53 "for i from 1 to N do: RHS:=0.5*(RHSpr(gv,i)+RHSp[i
]):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "if i<=N_trap then fv[i]:=fv[
i]+cdt*RHS: else fv[i]:=0: fi: od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
152 "if trunc(icol/10)*10=icol then PL[trunc(icol/10)]:=plot([seq([E(i
),(rho(i)*fv[i])],i=1..N)],title=cat(\"Collision number : \",(N/2)*ico
l*cdt/dE)): fi: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 77 "disp
lay([PL0,seq(PL[i],i=1..N_col/10)],insequence=true,view=[0..1.0,0..1.0
]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "E_avg:=add(E(i)*fv[i
]*rho(i),i=1..N)/add(fv[i]*rho(i),i=1..N);" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 31 "N_TP:=add(fv[i]*rho(i),i=1..N);" }}}{EXCHG {PARA 0 
"" 0 "" {TEXT -1 146 "Evidently the distribution peaks now at a lower \+
energy value, but for the same number of collisions a smaller fraction
 of bound particles remains." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 
0 "" 0 "" {TEXT 261 11 "Exercise 5:" }}{PARA 0 "" 0 "" {TEXT -1 177 "C
arry out Boltzmann simulations in which the reduction in barrier heigh
t happens faster (or more slowly) by changing the total number of coll
isions (controlled by the parameter " }{TEXT 19 5 "N_col" }{TEXT -1 
226 ", which is proportional to the collision number) in the above loo
p. Is the final result simply a function of collision number, or are t
here differences when the the barrier is lowered at a different rate w
ith collision number?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 
"" {TEXT 262 11 "Exercise 6:" }}{PARA 0 "" 0 "" {TEXT -1 244 "Change t
he power law for the density of states from the harmonic oscillator ca
se and observe the efficiency of evaporative cooling with collision nu
mber. How do the results depend on the power when it is changed in ste
ps of 1/2 up or down from " }{TEXT 19 3 "n=2" }{TEXT -1 1 "?" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 0 "" }}}}{MARK "64 2 0" 20 }{VIEWOPTS 1 1 0 1 1 1803 1 
1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }