Boltzmann equation simulation for a trapped Fermi gas of atoms

The dynamics of an interacting Fermi gas of atoms at sufficiently high temperatures can be efficiently studied via a numerical simulation of the Boltzmann equation. In this work we 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 evaluation of interparticle interactions. We compare different ways of choosing the phase space coordinates of a pair of atoms after a successful collision and demonstrate that the exact microscopic setup has no influence on the macroscopic outcome.


Introduction
The Boltzmann simulation of the semiclassical regime of trapped Fermi gases [1] has become a crucial tool to understand recent experiments [2,3,4,5] where spin-polarised atomic clouds undergo strong collisions. 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 work [6] we studied the collision of two spin-polarised fermionic clouds and obtained excellent agreement with experimental results. The goal of the present paper is to describe in detail the numerical setup that was used there and to present tests of the simulation and other information which is potentially very useful for those wishing to implement this important method.
Our setup is related to that of [7] but we introduce several new modifications [8]. For instance, we place an artificial grid over the continuous coordinate space that allows us to evaluate particle collisions efficiently. However, the main difference between this work and [7] is the method used to choose the new positions and momenta of two colliding particles after the collision. A quantum mechanical collision is a fundamentally random process with constraints given by the symmetries of the system. There are several ways to implement these constraints within the framework of a molecular dynamics simulation. We will discuss the physics of several collisional setups and also compare the numerical results associated with them in detail. We find that, although they exhibit very similar macroscopic behaviour, there are distinct technical 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 positions and momenta after a collision are discussed in detail in Sec. 3. In Sec. 4 we show how to determine 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.

Definitions
We consider a system of two-component fermions with equal mass m, labelled by the spin index s = {↑, ↓} and work in units in which = k B = 1. Fermions of opposite spin can interact via s-wave scattering and the cross section is given by where a is the scattering length and p rel = p ↑ − p ↓ is the relative momentum of the two atoms. We assume 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 s (r, p, t). Within the local density approximation the equilibrium distribution for fermions is given by the Fermi-Dirac distribution where T is the temperature and µ s the chemical potential defined by the normalisation condition 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, with the three trapping frequencies ω x , ω y and ω z . In the following we will only consider the isotropic case ω x = ω y = ω z or a cigar shaped trap with ω z < ω x = ω y ≡ ω ⊥ . We denote the geometric average of the trap frequencies as ω 0 = (ω x ω y ω z ) 1/3 . We define the Fermi energy from the energy distribution of the non-interacting gas at zero temperature, ε F =T F = (3N ) 1/3 ω 0 . This energy scale marks the frontier between the classical and quantum regimes. For T T F the fermions behave essentially like classical particles and can be described by the Maxwell-Boltzmann distribution, The Fermi energy can also be used to determine the typical scales of the cloud, wherek F is the Fermi momentum and R i =k F /(mω i ) are the Thomas-Fermi radii in the three spatial directions. These quantities give the widths of the zero-temperature momentum and density distributions respectively and are therefore useful to describe the extent of the cloud in momentum and coordinate space.
The time-evolution of the distribution function f s (r, p, t) is given by the Boltzmann equation, 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 between atoms carrying opposite spins. Indeed, for ultracold fermions, the Pauli principle forbids interactions between particles of the same spin [1]. At temperatures T T F when the Maxwell-Boltzmann distribution is applicable, the collision integral takes the form 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 ( Taking this into account, the collision integral reads At temperatures T T F the Pauli terms are close to one, and therefore the collision integral (8) tends to the classical form (7).

Test particles
To efficiently simulate the evolution of the continuous distribution function f s (r, p, t) we use the method of test particles. These are point-like particles which form a discrete approximation to f s (r, p, t) through δ-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 of test particles to atoms, the more precisely the continuous distribution function will be approximated, Physical observables need to be rescaled, for instance the test particle cross section becomesσ = σ(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, We introduce a discrete time step ∆t, such that during each time step the test particles propagate without colliding, following their classical trajectories. At the end of each time step collisions between them are evaluated. In a harmonic potential the trajectories are given by 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 probability of a two-particle collision in the same way as described in [7]. First we must test whether a given pair of test particles reaches the point of closest approach during the present time step. This condition is important to prevent particles from attempting to collide with each other 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 min of the test particles fulfils the classical condition for scattering: πd 2 min <σ. If this condition is also satisfied we propose a collision at the time of closest approach. 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 calculate the quantum mechanical scattering probability given by the Pauli term (1 − f ′ s )(1 − f ′ s ) and accept 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: The widths of these Gaussians, w x , w y , w z and w p , need to be tuned so that on the one hand fluctuations 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 s remains preserved [7]. The first condition is equivalent to where we introduced w r = (w x w y w z ) 1/3 , the geometric average of the spatial widths. The second condition implies w i ≪ R i and w p ≪k F . We also require that the smearing by the Gaussian functions preserve temperature-dependent degeneracy effects of the Fermi distribution, particularly close to the Fermi surface. This implies Note that the smearing width in momentum space is isotropic, while in position space the smearing width can be different depending on the spatial direction, if the corresponding trap frequencies are unequal. Since the Thomas-Fermi radii are inversely proportional to the corresponding trap frequencies it is sensible to choose w i = w r ω 0 /ω i for the spatial widths. Furthermore Eq. (16) together with the definition of the Thomas-Fermi radii imply w p = mω 0 w r . Hence all four smearing widths can be reduced to only one free parameter. 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 F forÑ /N = 10. Moreover, at very low temperatures the system undergoes a phase transition into a 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.

Auxiliary grid
The main numerical challenge is to efficiently evaluate collisions between the test particles. The total number of pairs of opposite spin isÑ 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 than to check all possible test particle pairs. The key observation is that since the cross section is decreasing with increasing relative momentum it can never exceedσ max = 4πa 2 N/Ñ and consequently the maximal distance two colliding test particles can have is d max = 2a N/Ñ . Having this in mind we superpose a cubic grid with cell size d max on the continuous space. The grid has finite extent which can be set by demanding 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 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 max around a given particle are definitely accounted for. This makes a total of 3 3 = 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 can be in non-neighbouring cells at the beginning and at the end of a time step, although in the 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 max . Also note that for large relative 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 collisions will indeed take place. To do so we consecutively select random pairs from the list. We then propagate both particles to the point of their closest approach, let them scatter (the exact setup for determining the new positions and momenta after scattering is described in Sec. 3) and then propagate them back to the original time. To account for quantum statistics we then calculate the Pauli blocking factors using the new The continuous distributions are obtained by replacing the δ-functions in (9) according to (14). If we obtain a value f ′ s > 1 from the summation we set f ′ s = 1. The probability that the collision is accepted is then given by ( Regardless if a collision is accepted or rejected neither 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.

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 has 12 degrees of freedom: three position and three momentum components for each of the two particles, 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 welldefined. 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 .

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. Motivated by the analogue of classical scattering we wish to conserve the total momentum, the total energy and the total angular momentum of the system during the collision. Conservation of the total momentum p = p s + p s and the total energy E = E kin = (p 2 s + p 2 s )/2m together imply the conservation of the modulus of the relative momentum p rel = |p s − p s |, since p 2 + p 2 rel = 2(p 2 s + p 2 s ) = 4mE. The direction of the relative momentum vector can change during the collision. Finally we must also conserve the angular momentum L = r rel × p rel . We can satisfy all these constraints by rotating both the relative momentum and relative position vectors by the same arbitrary angle around the L-axis. The angle of this rotation is the only degree of freedom of the collision and is determined at random. From the new values p ′ rel and r ′ rel the new values of the momenta and positions of the individual particles can be recovered using 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, the two colliding particles are not exactly at the same position before the collision and their relative position 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.

Energy preserving setup
It is possible to preserve energy conservation exactly by employing a different setup, for instance as in [7]. In this reference, the relative position stays fixed during a collision and the relative momentum vector is 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 collisional setup provides an additional, unphysical change of angular momentum. At any rate, since in this work we consider either isotropic or axially symmetric potentials, either the total angular momentum or its axial component L z are conserved.

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 , the time step ∆t and the smearing widths w r and w p . For all tests described below we use N = 10000 atoms andÑ /N = 10, which is sufficient for temperatures The optimal value of ∆t depends on the physical parameters of the system. Obvious requirements are that the time step must be smaller than the typical time between two collisions and that the average distance travelled during a time step must be much smaller than the diameter of the cross section, v rel ∆t ≪ σ /π. Another constraint is that the time step must be smaller than the half-period with respect to the largest trap frequency, ∆t < π/ω max , but unless the aspect ratio is extremely large this condition is much weaker 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.

Equilibrium collision rates
To obtain the correct dynamical properties, for instance the damping time of excitation modes, we need to ensure that the equilibrium collision rate observed in the simulation corresponds to the correct theoretical 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/Ñ ). The theoretical value for the collision rate γ in the presence of Pauli blocking is given by the following integral, 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.
Without Pauli blocking the integral is simpler and can be solved analytically for the Maxwell-Boltzmann distribution, where Ei(x) = x −∞ (e t /t)dt is the exponential integral. Furthermore we can obtain the theoretical prediction for the Pauli blocking probability by solving the integral (18) for the Fermi-Dirac distribution. The probability p Pauli that two classically colliding fermions will indeed scatter is then given by the ratio of γ block to this integral.
To find the optimal value for ∆t for each system we measure the collision rate in the absence of Pauli blocking for decreasing values of the time step and compare it to the theoretical prediction. The time step is small enough when we reach good agreement. It is important to switch off Pauli blocking for the tuning of the time step, as in the presence of Pauli blocking the collision rate is sensitive to the values of the smearing widths, which can obscure inaccuracies due 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 p and w r given by the 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 r = w p /(mω 0 ) = 1.0l ho for the lowest and w r = w p /(mω 0 ) = 2.0l ho for the highest temperatures used in our analysis, where l ho = 1/ √ mω 0 is the natural harmonic oscillator length unit.
The measured collision rates for the optimal choice of parameters with and without Pauli blocking together with the theoretical predictions are shown in Fig. 1. For sufficiently high temperatures the agreement is very good. At very low temperatures it becomes increasingly difficult to resolve the conditions (15) and (16) for the Gaussian smearing widths simultaneously. For larger values of the scattering length this problem 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.

Equilibrium energy distributions
Another important test is to check that the system thermalises to the correct equilibrium energy distribution, independently of the initial distribution. In the presence of Pauli blocking the energy is distributed 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 F a| = 1 and isotropic trap frequencies. The parameter values for the time step and the smearing widths are optimal. We performed two tests of the thermalisation. First we initialised the system according to the Fermi-Dirac distribution for T = 0.2T F (Fig. 2) and ran the simulation without Pauli blocking. After a short time the system thermalised according to the Maxwell-Boltzmann distribution for T = 0.31T F . Note that the temperatures of the two distributions are not necessarily equal, since the equipartition theorem does not hold for the Fermi-Dirac distribution. When changing from one distribution 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 F and Pauli blocking is switched on. It is clearly visible from both figures that the correct equilibrium distribution is always attained at the end of the simulation. This agreement improves further at higher temperatures.

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  the correct frequencies and damping properties of these modes. The following tests were performed for a spherical trap ω x = ω y = ω z = ω 0 . The sloshing mode is excited by a small displacement of the centre of mass from its equilibrium position, 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 i is well-known to be an undamped oscillation with the corresponding harmonic oscillator frequency ω i . Figure 4 shows such an oscillation for a system at |k F a| = 1 and T = 0.2T F .
The breathing mode can be excited by compressing or expanding the cloud. In a spherical trap this yields an undamped oscillation of r 2 with frequency 2ω 0 around its equilibrium value [11]. Figure 5 illustrates this mode for a system with |k F a| = 1 and T = 0.2T 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) = x 2 − y 2 by applying the perturbation ∆p x = −cx and ∆p y = cy with c = 0.2mω 0 in the same way as in [7]. Unlike the sloshing and the breathing mode this mode is damped. The frequency of the quadrupole mode at high temperatures approaches the ideal gas value 2ω 0 , while at low temperatures it is closer to the hydrodynamic frequency √ 2ω 0 . Figure 6 shows the quadrupole mode for |k F a| = 1 in the high and in the low temperature regime. 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(ω) = ∞ 0 dtQ(t)e iωt and is in agreement with the theoretical prediction while the damping can be extracted from its imaginary part. We see that the damping is stronger at the lower temperature T = 0.4T F when collisions are more frequent (see Fig. 1 for a plot of the collision rate).

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 changes in the total energy with the angular momentum conserving setup. Figure 7 shows a typical histogram of the relative change of the energy of a particle pair ∆E/E init = (E final − E init )/E init after a successful collision. The histogram is sharply peaked around zero. In this typical example the majority of collisions (approximately 84%) conserve the energy of the colliding pair with an accuracy of up to |∆E/E init | ≤ 10 −4 . Since collisions are frequent and the energy changes have random sign we also observe large cancellation effects, such that the total energy of the system is conserved with an accuracy of the order of 10 −7 at any given time. This is smaller than for instance the energy deviations due to Verlet's algorithm quoted in [7].  To understand the good accuracy of energy conservation we also make a qualitative theoretical estimate. During the collision the kinetic energy is conserved, such that it is sufficient to consider the potential energy V of the two particles, where (x i , y i , z i ) are the coordinates of the particle i and ω ⊥ = ω x = ω y is the radial trap frequency. We can recast this expression in terms of the centre of mass and the relative coordinates (X, Y, Z) and (x, y, z) respectively, The position of the centre of mass is not affected by the collision. The variation of energy is thus only due to the contribution of the relative motion V rel defined by We now define a second frame (x ′ , y ′ , z ′ ) such that the z ′ -axis is parallel to the relative angular momentum L of the pair. In addition, since the trap is rotationally invariant around z, we can choose the initial coordinate system such that L is in the (x, z) plane, in which case y ′ = y. If θ is the angle between L and the z axis, the two sets of coordinates are related by the following relations By conservation of angular momentum, the relative coordinates evolve in the (x ′ , y ′ ) plane. We thus have z ′ = 0, x ′ = d cos ϕ and y ′ = d sin ϕ, where d is the relative distance between the two particles. In terms of d and ϕ, the relative potential energy is During the collision d and θ are conserved. The variation of energy is thus only due to the rotation in the (x ′ , y ′ ) and we obtain where ϕ i and ϕ f correspond to the initial and final values of ϕ respectively. As an approximation, we normalise to the average total energy of the pair,Ē ≈ 6T , and obtain As shown above the distance between the particles can be at most d max = 2a N/Ñ . Hence for the range of parameters considered in this study the constant prefactor is of order 10 −4 . The trigonometric functions further reduce the variation of the energy and are responsible for the pointed shape of the histogram. Note that the method of test particles is important for a small energy variation. Due toÑ being larger than N the cross section and hence the average distance between two scattering partners is reduced, which leads to a smaller ∆E/E init (in our case by a factor ofÑ /N = 10). To demonstrate this effect Fig. 8 shows the variance of the energy change as a function of the ratioÑ /N . The data is well-fitted by a line with slope −2 on a double-logarithmic scale, which confirms the relation (25). Analogous scaling laws can be shown also for the dependence on the temperature and the scattering length. We also compared the time evolutions of the centre of mass coordinates in response to a perturbation. The differences between the two setups were found to be smaller than the statistical fluctuations. Taking these results into account we conclude that the two setups 3.1 and 3.2 are equivalent on a macroscopical scale and we have the freedom to choose whichever option appears more convenient.
We adopt the angular momentum conserving setup 3.1, not only because of this property, but also since the setup from [7] has a small technical disadvantage. Since the new direction of the relative momentum vector is chosen uniformly on a sphere, in many cases the particles are found to approach each other again after a successful collision. This implies that in the following time step they are likely to undergo another collision. The proportion of such events compared to the total number of collisions can reach up to 40% at the lowest temperature considered in this study and amounts to around 10% at T =T F . To avoid this overcounting of collisions one needs to implement an additional routine that prevents particles from colliding with each other repeatedly within short time intervals. In other words the molecular chaos assumption does not hold with the setup from [7] and needs to be enforced artificially. In our setup repeated scattering is so rare that this effect can be neglected. This is because the relative position and the relative momentum vectors of the particle pair are rotated by the same random angle and hence their relative angle is preserved. After a collision takes place the particles continue to move on trajectories that diverge from each other.

Outlook
In addition to the two collisional setups described above we propose a third very general collisional setup which can be explored in future work.
To preserve macroscopical averages we keep the centre of mass coordinates and momentum fixed and work with the six degrees of freedom for the relative momentum and position. Based on the idea of a delocalised particle pair, we assume that the relative phase space coordinates of the particles are distributed according to a Gaussian. We then choose the new relative position and momentum randomly according to the Gaussian probability distribution with the additional constraints given by the symmetries of the problem. Energy and L z conservation for instance reduce the problem by two degrees of freedom. If energy and total angular momentum are conserved the problem is reduced by four degrees of freedom and so on.
It is important to stress however, that the task of defining a probability measure on a complicated hypersurface in the phase space (in our case defined by the constraints of energy and angular momentum conservation) is a mathematically highly nontrivial problem and will pose a challenge for the numerical implementation.

Conclusion
We have implemented a Boltzmann equation approach to the simulation of the semiclassical regime of two-component Fermi gases. We described how to implement the method using test particles and stressed the advantages of the auxiliary grid. We also presented an extensive range of tests for the simulation. The main focus of this work was a discussion of the different possibilities for implementing collisions, namely the angular momentum conserving and the energy conserving setups. We studied in detail the error in the energy of the former scheme and showed that it is indeed extremely small, so that it becomes a valid alternative way of calculating the collisional integral and has distinct advantages since it minimises repeated scattering events. In future work we intend to apply the method to systems with components of unequal mass which bring a greater variety of dynamical regimes and which can provide a more stringent test of the Boltzmann equation approach to these systems.