Boltzmann equation simulation for a trapped Fermi gas of atoms O. Goulkoa,b,, F. Chevyc, C. Lobod aDepartment of Applied Mathematics and Theoretical Physics, Universityof Cambridge, Centre for Mathematical Sciences, Cambridge, CB3 0WA, UnitedKingdom bPresent address: Physics Department, Arnold Sommerfeld Centerfor Theoretical Physics, and Centerfor NanoScience, Ludwig-Maximilians-Universit¨at, Theresienstraße 37, 80333 Munich, Germany 2 cLaboratoire KastlerBrossel, Ecole Normale Sup´erieure, CNRS, UPMC, 24 rue Lhomond, 75231 Paris Cedex 05, France 1 dSchool of Mathematics, Universityof Southampton, Highfield, Southampton, SO17 1BJ, UnitedKingdom 0 2 n a J Abstract 0 3 ThedynamicsofaninteractingFermigasofatomsatsufficientlyhightemperaturescanbeefficientlystudied via a numericalsimulation of the Boltzmann equation. In this workwe describe in detail the setup we used ] recently to study the oscillations of two spin-polarised fermionic clouds in a trap. We focus here on the s a evaluation of interparticle interactions. We compare different ways of choosing the phase space coordinates g of a pair of atoms after a successful collision and demonstrate that the exact microscopic setup has no - t influence on the macroscopic outcome. n a Keywords: Boltzmann equation, cold atoms, interacting Fermi gas u PACS: 03.75.Ss, 03.75.Hh q . t a m - 1. Introduction d n The Boltzmann simulation of the semiclassical regime of trapped Fermi gases [1] has become a crucial o tool to understand recent experiments [2, 3, 4, 5] where spin-polarised atomic clouds undergo strong col- c lisions. The simulation is one of the very few methods which allow us to calculate accurately e.g. spin [ transport coefficients in ways that can then be quantitatively compared with experiment. In a previous 1 work [6] we studied the collision of two spin-polarised fermionic clouds and obtained excellent agreement v with experimental results. The goal of the present paper is to describe in detail the numerical setup that 5 3 was used there and to present tests of the simulationand other informationwhich is potentially very useful 2 for those wishing to implement this important method. 6 Our setup is related to that of [7] but we introduce several new modifications [8]. For instance, we . 1 place an artificial grid over the continuous coordinate space that allows us to evaluate particle collisions 0 efficiently. However, the main difference between this work and [7] is the method used to choose the new 2 positions and momenta of two colliding particles after the collision. A quantum mechanical collision is a 1 fundamentally random process with constraints given by the symmetries of the system. There are several : v ways to implement these constraints within the framework of a molecular dynamics simulation. We will i X discussthe physicsofseveralcollisionalsetups andalsocomparethe numericalresultsassociatedwith them indetail. Wefindthat,althoughtheyexhibitverysimilarmacroscopicbehaviour,therearedistincttechnical r a advantages which make some of the methods preferable over others. In Sec. 2 of this paper we introduce the general numerical setup. We describe the method of test particles and explain what conditions need to be fulfilled for a semiclassical collision to take place and how to efficiently search the system for suitable collision partners. The different ways of choosing the new positionsandmomentaafteracollisionarediscussedindetailinSec.3. InSec.4weshowhowtodetermine Email addresses: [email protected] (O.Goulko), [email protected] (F.Chevy), [email protected] (C.Lobo) Preprint submitted toCNSNS January 31, 2012 the optimal simulation parameters and present several tests of the numerical method, in particular also of the different collisional setups. A summary of our results and an outlook on future work can be found in Sec. 5. 2. Numerical setup 2.1. Definitions Weconsiderasystemoftwo-componentfermionswithequalmassm,labelledbythespinindexs= , and work in units in which ~ = k = 1. Fermions of opposite spin can interact via s-wave scattering{↑an↓d} B the cross section is given by 4πa2 σ = , (1) 1+a2p2 /4 rel where a is the scatteringlength and p =p p is the relative momentum of the two atoms. We assume rel ↑ ↓ − that the system is in the normal phase and that the temperature is sufficiently high, so that the two spin distributions can be described semiclassically in terms of functions f (r,p,t). Within the local density s approximationthe equilibrium distribution for fermions is given by the Fermi-Dirac distribution 1 f(FD)(r,p)= , (2) s e(p2/2m+V(r)−µs)/T +1 where T is the temperature and µ the chemical potential defined by the normalisation condition s d3p N =N/2= d3rn (r,t)= d3r f (r,p,t). (3) s Z s Z Z (2π)3 s The atom number for each species is assumed to be equal, N = N = N/2, and will be held fixed during ↑ ↓ the simulation. The trapping potential is assumed to be harmonic, 1 V(r)= m(ω2x2+ω2y2+ω2z2), (4) 2 x y z with the three trapping frequencies ω , ω and ω . In the following we will only consider the isotropic case x y z ω = ω = ω or a cigar shaped trap with ω < ω = ω ω . We denote the geometric average of the x y z z x y ⊥ ≡ trap frequencies as ω =(ω ω ω )1/3. 0 x y z We define the Fermi energyfromthe energy distributionof the non-interactinggasat zerotemperature, ε˜ = T˜ = (3N)1/3ω . This energy scale marks the frontier between the classical and quantum regimes. F F 0 For T & T˜ the fermions behave essentially like classical particles and can be described by the Maxwell- F Boltzmann distribution, f(MB)(r,p)=N ω03e−(p2/2m+V(r))/T. (5) s sT3 The Fermi energy can also be used to determine the typical scales of the cloud, ε˜ = k˜2/2m = mω2R2/2, F F i i where k˜ is the Fermi momentum and R = k˜ /(mω ) are the Thomas-Fermi radii in the three spatial F i F i directions. These quantities give the widths of the zero-temperature momentum and density distributions respectivelyandarethereforeusefultodescribethe extentofthe cloudinmomentumandcoordinatespace. The time-evolution of the distribution function f (r,p,t) is given by the Boltzmann equation, s ∂ f +(p/m) ∇ f ∇ V ∇ f = I[f ,f ], (6) t s r s r p s s s · − · − where the left-hand side represents the propagation of the atoms in the potential and the right-hand side stands for the collision integral, which depends on the particle statistics. Here, we consider only collisions betweenatomscarryingoppositespins. Indeed,forultracoldfermions,thePauliprincipleforbidsinteractions 2 between particles of the same spin [1]. At temperatures T &T˜ when the Maxwell-Boltzmann distribution F is applicable, the collision integral takes the form d3p dσ p p I [f ,f ]= s dΩ | s− s|[f f f′f′], (7) class s s Z (2π)3 Z dΩ m s s− s s where the indices s and s label the two colliding atoms of opposite spin, the primed variables refer to quantities after the collision and Ω is the angle between the incoming and outgoing relative momenta. The Pauli exclusion principle only allows fermions to scatter into a previously unoccupied quantum state. This reduces the scattering probability by the Pauli blocking term proportional to (1 f′)(1 f′). Taking this − s − s into account, the collision integral reads d3p dσ p p I[f ,f ]= s dΩ | s− s|[f f (1 f′)(1 f′) f′f′(1 f )(1 f )]. (8) s s Z (2π)3 Z dΩ m s s − s − s − s s − s − s At temperatures T & T˜ the Pauli terms are close to one, and therefore the collision integral (8) tends to F the classical form (7). 2.2. Test particles Toefficientlysimulatethe evolutionofthecontinuousdistributionfunctionf (r,p,t)weusethe method s of test particles. These are point-like particles which form a discrete approximation to f (r,p,t) through s δ-functions. In order for this approximation to be accurate we will represent each fermion by several test particles [7, 9, 10]. The higher the ratio N˜/N of test particles to atoms, the more precisely the continuous distribution function will be approximated, N N˜s f (r,p,t) s (2π)3δ(r r (t))δ(p p (t)). (9) s → N˜ − i − i s Xi=1 Physical observables need to be rescaled, for instance the test particle cross section becomes σ˜ = σ(N/N˜). A generic thermal expectation value of a single-particle observable X(r,p) can be easily calculated within the test particle picture, as the integration reduces to a sum over all test particles, 1 d3p X = d3r f (r,p,t)X(r,p) (10) h i N Z (2π)3 s Xs s N˜ 1 = X(r ,p ). (11) N˜ i i Xi=1 We introduce a discrete time step ∆t, such that during each time step the test particles propagate withoutcolliding,followingtheir classicaltrajectories. Atthe endofeachtime stepcollisionsbetweenthem are evaluated. In a harmonic potential the trajectories are given by r (t ) = r (t )cos(ω ∆t)+(p (t )/mω )sin(ω ∆t), (12) i n+1 i n i i n i i p (t ) = p (t )cos(ω ∆t) r (t )mω sin(ω ∆t). (13) i n+1 i n i i n i i − Note that since the time step is fixed, the trigonometric functions only need to be evaluated once during the entire simulation, so that using the exact solution is more efficient than using the Verlet algorithm, as for instance in Refs. [7, 9], in which case the accelerations need to be recalculated for every time step. The Verlet algorithm is more general as it is applicable for any potential. But since in this work we will only consider harmonic potentials, we will use the exact solution instead. We evaluate the probabilityofatwo-particlecollisioninthe samewayasdescribedin[7]. Firstwe must test whether a given pair of test particles reaches the point of closest approach during the present time 3 step. Thisconditionis importantto preventparticlesfromattempting to collidewitheachother repeatedly over several consecutive time steps, an issue which will be further addressed below. If the closest approach condition is true we check if the minimal distance d of the test particles fulfils the classicalcondition for min scattering: πd2 <σ˜. Ifthisconditionisalsosatisfiedweproposeacollisionatthetimeofclosestapproach. min However due to Pauli statistics, even if the classical conditions for scattering are fulfilled, a collision can only take place if the new state of the particles was previously unoccupied. To take this into account we calculatethe quantummechanicalscatteringprobabilitygivenbythe Pauliterm(1 f′)(1 f′)andaccept − s − s or reject the collision according to this probability. Clearly, the point-like particle picture is unsuitable for the calculation of this probability. To return to a continuous distribution we therefore have to smear out the δ-functions representing the test particles, e.g. by Gaussians in position and momentum space: e−(p−pi)2/wp2 δ(p p ) (14a) − i → (√πw )3 p e−(x−xi)2/wx2 e−(y−yi)2/wy2 e−(z−zi)2/wz2 δ(r r ) . (14b) i − → √πw √πw √πw x y z The widths of these Gaussians, w , w , w and w , need to be tuned so that on the one hand fluctuations x y z p due to the discrete nature of the test particle picture are smoothed out, but on the other hand the physical structure of the distribution function f remains preserved [7]. The first condition is equivalent to s w w (N/N˜)1/3, (15) p r ≫ where we introduced w =(w w w )1/3, the geometric averageof the spatial widths. The secondcondition r x y z implies w R and w k˜ . i i p F ≪ ≪ WealsorequirethatthesmearingbytheGaussianfunctionspreservetemperature-dependentdegeneracy effects of the Fermi distribution, particularly close to the Fermi surface. This implies w k˜ (T/T˜ ) and w R (T/T˜ ). (16) p F F i i F ≪ ≪ Notethatthesmearingwidthinmomentumspaceisisotropic,whileinpositionspacethesmearingwidthcan be different depending onthe spatialdirection, if the correspondingtrap frequencies are unequal. Since the Thomas-Fermiradii are inversely proportionalto the correspondingtrap frequencies it is sensible to choose w = w ω /ω for the spatial widths. Furthermore Eq. (16) together with the definition of the Thomas- i r 0 i Fermi radiiimply w =mω w . Hence all four smearing widths canbe reduced to only one free parameter. p 0 r At very low temperatures the margin given by the conditions (15) and (16) becomes so narrow that it is impossible to find smearing widths satisfying both conditions, without having to significantly increase the number of test particles. This limits the applicability of this setup to temperatures above approximately 0.2T˜ for N˜/N = 10. Moreover, at very low temperatures the system undergoes a phase transition into a F superfluid state. This method does not include the relevant degrees of freedom of the superfluid Fermi gas and is only applicable in the normal phase. 2.3. Auxiliary grid The main numerical challenge is to efficiently evaluate collisions between the test particles. The total number of pairs of opposite spin is N˜2/4, which is an unfavourable scaling given that we want to use large particle numbers and a high test particle to particle ratio. In this work we develop a more efficient method thantocheckallpossibletestparticlepairs. Thekeyobservationisthatsincethecrosssectionisdecreasing with increasing relative momentum it can never exceed σ˜ = 4πa2N/N˜ and consequently the maximal max distancetwocollidingtestparticlescanhaveisd =2a N/N˜. Havingthisinmindwesuperposeacubic max q grid with cell size d on the continuous space. The grid has finite extent which can be set by demanding max that a certain proportion, for instance at least 95%, of the test particles are within the grid. For a cigar shaped trap the grid must have larger extent in the axial direction. At the end of each time step we move systematically through all grid cells starting in one of the corners and note all possible collision partners 4 that fulfil the classical scattering conditions. To make sure that we do not miss collisions due to boundary effects, for each particle we check not only all opposite spin particles in the same grid cell, but also in all neighbouring cells (the ones sharing a face, an edge or a vertex with the given cell). This ensures that all particles in a sphere of radius d around a given particle are definitely accounted for. This makes a total max of 33 = 27 cells for each particle, however to avoid double counting we only need to evaluate cells in the positive direction, which means on average 14 cells per particle. A small systematic error source remains with this setup. If the relative velocity of two particles is large they canbe innon-neighbouringcells atthe beginningandatthe endofatime step, althoughinthe course of the time step they come within each other’s allowed collision range. Such a possible collision will then not be accounted for. However this systematic error can be minimised by choosing the time step to be sufficiently small and also by choosing the cell size to be larger than d . Also note that for large relative max velocities the cross section is small and collisions between very fast particles are rare events. After having searched the entire grid for classically allowed collision pairs we proceed to choose which collisionswillindeedtakeplace. Todosowe consecutivelyselectrandompairsfromthe list. We thenprop- agateboth particlestothe pointoftheirclosestapproach,letthem scatter(the exactsetupfordetermining thenewpositionsandmomentaafterscatteringisdescribedinSec. 3)andthenpropagatethembacktothe originaltime. To account for quantum statistics we then calculate the Pauli blocking factors using the new distributions f′ = f (r′,p′). The continuous distributions are obtained by replacing the δ-functions in (9) s s according to (14). If we obtain a value f′ >1 from the summation we set f′ =1. The probability that the s s collisionis acceptedis thengivenby (1 f′)(1 f′). Regardlessifa collisionis acceptedorrejectedneither − s − s of the particles concerned is allowed to collide again with another particle during the present time step. If the collision is accepted we keep the new positions and momenta. If the collision is rejected we return to the values before the collision. This procedure is repeated until all possible pairs have been evaluated. 3. Collisional setup Since the picture of colliding point-like particles with well-defined positions and momenta is a classical interpretation of a quantum mechanical scattering process, it is unsurprising that there must be some ambiguity in the implementation of the collisional setup. In the semiclassical particle picture each collision has12degreesoffreedom: three positionandthree momentumcomponents foreachofthe twoparticles,or equivalently three components of the total and relative positions and momenta. In quantum mechanics we consider wave packets rather than particles and concepts like particle position or momentum are not well- defined. Instead the system is described by operators which, depending on the symmetries of the system, commute with the Hamiltonian and define quantum numbers corresponding to conserved quantities. The concept of a trajectory does not exist, as a particle is not localised in phase space but rather smeared out over a certain phase space volume in accordance with Heisenberg’s uncertainty principle. Therefore, in the collisional setup, it is not necessary to preserve for instance the exact positions of the two atoms during a collision. In the presence of an axially symmetric external potential for instance, the only conserved quantities are the total energy and the angular momentum component L . z 3.1. Angular momentum preserving setup The method used in this work and in [6] can be motivated as follows. As collisions are local, we can disregard for the moment the external potential and go to the centre of mass frame of the two particles. Motivatedby the analogueofclassicalscatteringwe wishto conservethe totalmomentum, the totalenergy and the total angular momentum of the system during the collision. Conservation of the total momentum p=p +p and the total energy E =E =(p2+p2)/2m together imply the conservationof the modulus s s kin s s ofthe relative momentum p = p p , since p2+p2 =2(p2+p2)=4mE. The direction ofthe relative rel | s− s| rel s s momentum vector can change during the collision. Finally we must also conserve the angular momentum L = r p . We can satisfy all these constraints by rotating both the relative momentum and relative rel rel × positionvectorsbythesamearbitraryanglearoundtheL-axis. Theangleofthisrotationisthe onlydegree offreedomofthecollisionandisdeterminedatrandom. Fromthenewvaluesp′ andr′ the newvaluesof rel rel 5 the momenta and positions of the individual particles can be recoveredusing total momentum conservation and the centre of mass coordinates respectively. So far we have ignored the external potential, which would be justified if the collisions took place when the atoms were infinitesimally close to each other and therefore experiencing the same potential. However, thetwocollidingparticlesarenotexactlyatthesamepositionbeforethecollisionandtheirrelativeposition changes after a successful collision. Thus potential energy is not conserved and total energy conservation is not exact during a collision in general. Nevertheless, we will show below that such changes in the total energy are negligibly small for almost any collision. Furthermore, these small changes cancel each other in a many-particle system which experiences many collisions. 3.2. Energy preserving setup It is possible to preserve energy conservation exactly by employing a different setup, for instance as in [7]. Inthis reference,therelativepositionstaysfixedduringacollisionandthe relativemomentumvectoris rotated by a random angle in space (such a rotation has two degrees of freedom). As a direct consequence of the unrestricted rotation this setup violates angular momentum conservation. While it is true that this violation would occur naturally in non-axially symmetric potentials even for non-interacting atoms, this collisionalsetupprovidesanadditional,unphysicalchangeof angularmomentum. Atany rate,since in this work we consider either isotropic or axially symmetric potentials, either the total angular momentum or its axial component L are conserved. z 4. Tests and optimal parameters How accurately the numerical setup represents the physical picture depends crucially on the values of the simulation parameters, in particular N˜/N, the time step ∆t and the smearing widths w and w . For r p all tests described below we use N = 10000 atoms and N˜/N = 10, which is sufficient for temperatures T 0.2T˜ . F ≥ The optimal value of ∆t depends on the physical parameters of the system. Obvious requirements are thatthetimestepmustbesmallerthanthetypicaltimebetweentwocollisionsandthattheaveragedistance travelledduringatimestepmustbemuchsmallerthanthediameterofthecrosssection, v ∆t σ˜ /π. rel h i ≪ h i Another constraint is that the time step must be smaller than the half-period with respect to thpe largest trap frequency, ∆t < π/ω , but unless the aspect ratio is extremely large this condition is much weaker max than the other ones. There is no lower bound on the time step, however the simulation slows down with decreasing ∆t. We first describe the tuning of the parameters and the corresponding tests of the simulation with the angular momentum preserving collisional setup 3.1. We perform the same tests as in [7] to demonstrate that even with the different collisional setup we obtain very good agreement. Then we explicitly compare the different collisional setups in Sec. 4.4. 4.1. Equilibrium collision rates To obtainthe correctdynamical properties,for instance the damping time of excitationmodes, we need toensurethattheequilibriumcollisionrateobservedinthesimulationcorrespondstothecorrecttheoretical value. The number of collisions in a given time interval can be very easily obtained from the simulation, simply by counting all test particle collisions and then multiplying by the ratio (N/N˜). The theoretical value for the collision rate γ in the presence of Pauli blocking is given by the following integral, d3p d3p dσ γ = d3r s s dΩ p p f f (1 f′)(1 f′). (17) block Z Z (2π)3(2π)3 Z dΩ| s− s| s s − s − s After inserting the Fermi-Dirac distribution this integral can be calculated numerically [7]. The numerical setup also provides a powerful tool to artificially switch off Pauli blocking. This allows us to separately check for errors related to the general numerical setup and the calculation of the Pauli blocking factors. 6 Without Pauli blocking the integral is simpler and can be solved analytically for the Maxwell-Boltzmann distribution, d3p d3p dσ γ = d3r s s dΩ p p f f (18) noblock Z Z (2π)3(2π)3 Z dΩ| s− s| s s = N↑N↓π2ωT032 (cid:18)1+ ma12Tema12TEi(cid:18)−ma12T(cid:19)(cid:19), (19) whereEi(x)= x (et/t)dtistheexponentialintegral. Furthermorewecanobtainthetheoreticalprediction −∞ for the Pauli blRocking probability by solving the integral (18) for the Fermi-Dirac distribution. The proba- bility p that two classically colliding fermions will indeed scatter is then given by the ratio of γ to Pauli block this integral. To find the optimal value for ∆t for each system we measure the collision rate in the absence of Pauli blockingfordecreasingvaluesofthetimestepandcompareittothetheoreticalprediction. Thetimestepis small enough when we reach good agreement. It is important to switch off Pauli blocking for the tuning of thetime step, asinthe presenceofPauliblockingthe collisionrateissensitivetothe valuesofthesmearing widths, which canobscure inaccuraciesdue to a time step that is too large. After having found the optimal time step we check that repeated (unphysical) collisions of the same particle pair are rare. This has always been found to be the case with our collisional setup. We then use this optimal value for ∆t to establish the optimal values of the smearing widths. We first identify the allowed interval for w and w given by the p r conditions (15) and (16) and perform the same collision rate matching as described above, this time with Pauli blocking switched on. We find that the optimal widths lie between w = w /(mω ) = 1.0l for the r p 0 ho lowestandw =w /(mω )=2.0l forthe highesttemperaturesusedinouranalysis,wherel =1/√mω r p 0 ho ho 0 is the natural harmonic oscillator length unit. The measured collision rates for the optimal choice of parameters with and without Pauli blocking to- getherwiththetheoreticalpredictionsareshowninFig.1. Forsufficientlyhightemperaturestheagreement is very good. At very low temperatures it becomes increasingly difficult to resolve the conditions (15) and (16)fortheGaussiansmearingwidthssimultaneously. Forlargervaluesofthescatteringlengththisproblem gets worse since test particles with a large relative distance can scatter on each other [7] and therefore the continuous distribution function needs to be resolved accurately over larger scales. 5 10 γ /(N ω ) γ /(N ω ) block s 0 block s 0 γ /(N ω ) γ /(N ω ) 4 noblock s 0 8 noblock s 0 p p Pauli Pauli 3 6 2 4 1 2 0 0 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 ~ ~ T/T T/T F F Figure1: (Colouronline)Theequilibriumcollisionratesperparticle(inunitsofω0)withandwithoutPauliblocking,aswell as the Pauli probability for a successful scattering versus temperature for k˜Fa = 1 (left) and k˜Fa = 2 (right). The lines | | | | correspondtothetheoretical predictionandthesymbolstothevaluesobtainedwiththesimulation. 4.2. Equilibrium energy distributions Another important test is to check that the system thermalises to the correct equilibrium energy distri- bution, independently of the initial distribution. In the presence of Pauli blocking the energy is distributed 7 according to the Fermi-Dirac distribution, whereas without Pauli blocking the particles will be distributed according to the Maxwell-Boltzmann distribution. Figures 2 and 3 show the results of this test for a low temperature system with k˜ a = 1 and isotropic trap frequencies. The parameter values for the time step F | | and the smearing widths are optimal. We performed two tests of the thermalisation. First we initialised the systemaccordingto the Fermi-Diracdistribution for T =0.2T˜ (Fig.2)and ranthe simulationwithout F Pauliblocking. Afterashorttime thesystemthermalisedaccordingtothe Maxwell-Boltzmanndistribution for T = 0.31T˜ . Note that the temperatures of the two distributions are not necessarily equal, since the F equipartitiontheoremdoes notholdfor the Fermi-Diracdistribution. Whenchangingfromonedistribution to the other the average energy of the system remains conserved and hence the temperature of the new equilibrium state is different. Figure 3 shows the corresponding results for the reverse situation: the initial distribution is Maxwell-Boltzmann with T =0.31T˜ and Pauli blocking is switched on. It is clearly visible F from both figures that the correct equilibrium distribution is always attained at the end of the simulation. This agreement improves further at higher temperatures. 0.12 0.12 start distribution end distribution 0.1 M.-B. distribution 0.1 M.-B. distribution F.-D. distribution F.-D. distribution 0.08 0.08 s s N N E)/ 0.06 E)/ 0.06 (s (s f f 0.04 0.04 0.02 0.02 0 0 0 20 40 60 80 100 0 20 40 60 80 100 E/ω E/ω 0 0 0.01 0.01 start distribution end distribution M.-B. distribution M.-B. distribution 0.001 0.001 F.-D. distribution F.-D. distribution s s N N E)/ 0.0001 E)/ 0.0001 (s (s 2)f 2)f E 1e-05 E 1e-05 3/ 3/ 0 0 ω ω ( ( 1e-06 1e-06 1e-07 1e-07 0 20 40 60 80 100 0 20 40 60 80 100 E/ω E/ω 0 0 Figure 2: (Colour online) The equilibrium energy distributions without Pauli blocking. The top panel shows the energy distributionson alinearscale, the bottom panel shows the energydistributions scaled byω3/E2 onalogarithmicscale. The 0 start distribution (left) is Fermi-Dirac (F.-D.) at T = 0.2T˜F and the end distribution (right) is Maxwell-Boltzmann (M.-B.) withthesameaverageenergybutatT =0.31T˜F,asexpected. 4.3. Collective excitations Collective excitations emerge when a many-particle system is perturbed away from equilibrium. Here we will discuss three different excitation modes, the sloshing mode (also known as dipole or Kohn mode), the breathing mode (monopole mode) and the quadrupole mode. We will confirm that the simulation gives 8 0.12 0.12 start distribution end distribution 0.1 M.-B. distribution 0.1 M.-B. distribution F.-D. distribution F.-D. distribution 0.08 0.08 s s N N E)/ 0.06 E)/ 0.06 (s (s f f 0.04 0.04 0.02 0.02 0 0 0 20 40 60 80 100 0 20 40 60 80 100 E/ω E/ω 0 0 0.01 0.01 start distribution end distribution M.-B. distribution M.-B. distribution 0.001 0.001 F.-D. distribution F.-D. distribution s s N N E)/ 0.0001 E)/ 0.0001 (s (s 2)f 2)f E 1e-05 E 1e-05 3/ 3/ 0 0 ω ω ( ( 1e-06 1e-06 1e-07 1e-07 0 20 40 60 80 100 0 20 40 60 80 100 E/ω E/ω 0 0 Figure3: (Colouronline)TheequilibriumenergydistributionswithPauliblocking. Thetoppanelshowstheenergydistributions onalinearscale,thebottompanelshowstheenergydistributionsscaledbyω3/E2onalogarithmicscale. Thestartdistribution 0 (left)isMaxwell-Boltzmann(M.-B.)atT =0.31T˜Fandtheenddistribution(right)isFermi-Dirac(F.-D.)withthesameaverage energybutatT =0.2T˜F,asexpected. 9 0.3 0.2 0.1 z >/l 0 z < -0.1 -0.2 -0.3 0 2π 4π 6π 8π ω t z Figure4: Simulationoftheequilibriumsloshingmode z /lz,wherelz =1/√mωz. Themodeisundampedandtheoscillation h i frequencyequalsωz. the correct frequencies and damping properties of these modes. The following tests were performed for a spherical trap ω =ω =ω =ω . x y z 0 Thesloshingmodeisexcitedbyasmalldisplacementofthe centreofmassfromitsequilibriumposition, or equivalently by a short-lived force represented by an additional linear term in the potential. The time evolution of each of the three centre of mass coordinates r is well-known to be an undamped oscillation i h i with the correspondingharmonic oscillatorfrequency ω . Figure 4 shows suchan oscillationfor a systemat i k˜ a =1 and T =0.2T˜ . F F | | Thebreathingmodecanbeexcitedbycompressingorexpandingthecloud. Inasphericaltrapthisyields an undamped oscillation of r2 with frequency 2ω around its equilibrium value [11]. Figure 5 illustrates 0 this mode for a system withhk˜ia =1 and T =0.2T˜ . F F | | By perturbing the system via a short-lived small increase in one or several of the trap frequencies we can excite the quadrupole mode. We excite the quadrupole mode Q(t) = x2 y2 by applying the h i − h i perturbation ∆p = cx and ∆p = cy with c = 0.2mω in the same way as in [7]. Unlike the sloshing x y 0 − andthe breathingmode this mode isdamped. The frequencyofthe quadrupolemode athightemperatures approaches the ideal gas value 2ω , while at low temperatures it is closer to the hydrodynamic frequency 0 √2ω . Figure 6 shows the quadrupole mode for k˜ a = 1 in the high and in the low temperature regime. 0 F | | In both cases the damping of the mode is clearly visible. The frequency of the oscillation can be extracted from the corresponding Fourier transform Q(ω) = ∞dtQ(t)eiωt and is in agreement with the theoretical 0 predictionwhilethedampingcanbeextractedfromRitsimaginarypart. Weseethatthedampingisstronger at the lower temperature T =0.4T˜ when collisions are more frequent (see Fig. 1 for a plot of the collision F rate). 4.4. Comparison of different collisional setups To study the impact of the collisional setup on the outcome of the simulation we have implemented the setups 3.1 (angular momentum preserving) and 3.2 (energy preserving) and confirmed that they generate the same results for both equilibrium and non-equilibrium systems. First we analyse the magnitude of the changesinthetotalenergywiththeangularmomentumconservingsetup. Figure7showsatypicalhistogram of the relative change of the energy of a particle pair ∆E/E = (E E )/E after a successful init final init init − collision. The histogram is sharply peaked around zero. In this typical example the majority of collisions (approximately84%)conservethe energyofthe collidingpairwithanaccuracyofupto ∆E/E 10−4. init | |≤ Since collisions are frequent and the energy changes have random sign we also observe large cancellation 10