## Abstract

Using the maximum entropy principle, a kinetic model equation is proposed to simplify the intricate collision term in the semi-classical Boltzmann equation for dilute quantum gases in the normal phase. The kinetic model equation keeps the main properties of the Boltzmann equation, including conservation of mass, momentum and energy, the entropy dissipation property, and rotational invariance. It also produces the correct Prandtl numbers for the Fermi gases. To validate the proposed model, the kinetic model equation is numerically solved in the hydrodynamic and kinetic flow regimes using the asymptotic preserving scheme. The results agree well with those of the quantum Euler and Boltzmann equations.

## 1. Introduction

The experimental realization of quantum degeneracy in ultracold atomic gases has attracted intensive research effort to understand the interacting quantum systems (see Dalfovo *et al.* (1999) and Giorgini *et al.* (2008) for reviews). The controllability of the interactions, energy and spin population make these quantum systems ideal to study the crossover from a Bardeen–Cooper–Schrieffer superfluid to Bose–Einstein condensation, which is ubiquitous in high-temperature superconductivity, neutron stars, nuclear matter and quark-gluon plasma. To probe the properties of quantum-degenerate and strongly interacting Fermi gases, elementary oscillations are often excited. Also, to detect the low-lying excitations, one usually switches off the external trapping potential and allows the atomic cloud to expand freely. In these cases, the systems are not in thermodynamic equilibrium. Theoretically, the Euler and Navier–Stokes equations (Schäfer 2010; Cao *et al.* 2011), and the semi-classical Boltzmann equation (Uehling & Uhlenbeck 1933), can model the dynamics of dilute quantum gas in the normal phase. While the Navier–Stokes equations are applicable only in the hydrodynamic flow regime with very small Knudsen numbers (the Knudsen number is defined as *Kn*=*λ*/*L* or *ω*/*ν*, where *λ* and *ν* are the mean free path and mean collision frequency, *L* and *ω* are the characteristic length scale and frequency of the fluid system), the semi-classical Boltzmann equation is applicable for a wider range of conditions (Nikuni & Griffin 1998; Struchtrup 2005). Specifically, through the Chapman–Enskog expansion (Chapman & Cowling 1970), the Euler and Navier–Stokes equations are, respectively, recovered as the zeroth-order and first-order approximations in the power series expansion of the distribution function in terms of the Knudsen number.

The Knudsen number in the experiments of ultracold quantum gases is usually large. For instance, the Bose gases in the normal phase are usually in the transition regime where *Kn*∼1 (Pitaevskii & Stringari 2003; Buggle *et al.* 2005), say, in the experiment of Meppelink *et al.* (2009), about 1.3×10^{9} atoms are confined in harmonic trap with the trap frequency *ω*_{rad}=2*π*×90 Hz in the radial direction and *ω*_{ax}=2*π*×1.46 Hz in the axial direction. The average collision frequency is *ν*≈90 Hz, and the Knudsen number is about 6 and 0.1. Further incensement of atom numbers will reduce the Knudsen number, but this is forbidden by the strong three-body loss. Although the two-component Fermi gases in the normal phase can be in the hydrodynamic regime in the central region of the cloud, in the surface region, the gas is always collisionless where *Kn*≫1 (O’Hara *et al.* 2002; Regal & Jin 2003; Wright *et al.* 2007). This indicates that the Navier–Stokes equations will fail to describe the dynamics of the trapped quantum gas. Alternatively, one turns to the semi-classical Boltzmann equation to study the low-lying collective oscillations and free expansion of the atomic clouds. However, owing to the complicated nature of the nonlinear integro-differential collision integral, the Boltzmann equation is difficult to handle both analytically and numerically. For this reason, it is useful to consider kinetic model equations that greatly simplify the collision integral, but retain its main physical properties, e.g. conservation of mass, momentum and energy, validation of the entropy dissipation property, and recovery of the correct Prandtl number (ratio between the shear viscosity and thermal conductivity). The semi-classical Bhatnager–Gross–Krook (BGK) model based on the relaxation-time approximation is the well-known kinetic model to simplify the intricate collision integral (see Bhatnagar *et al.* (1954) for the BGK model for classical gases). The only difference in the semi-classical version is that the Maxwell–Boltzmann equilibrium distribution function is replaced by the Fermi–Dirac or Bose–Einstein equilibrium distribution function, depending on the statistic types of the particles. Based on this model, analytical expressions for the mode frequency and damping of the low-lying collective oscillations in quantum gases have been obtained by the method of moments (Khawaja *et al.* 2000; Bruun & Smith 2007; Riedl *et al.* 2008) which agree with the experimental data to some extent. Note that these analytical expressions contain a spatially averaged relaxation time which is determined by the shear viscosity. This coincides with most of the experiments where the temperature is spatially homogenous so that one does not have to consider the thermal conductivity.

In some experiments, however, the thermal conductivity has to be taken into account. For instance, the cross-dimensional rethermalization problem where the temperature in one direction is initially increased by parametric heating (Monroe *et al.* 1993; DeMarco *et al.* 1999; Loftus *et al.* 2002; Regal *et al.* 2003). Generally speaking, the shear relaxation time and thermal relaxation time of the quantum gas are different. Since the BGK model has only one relaxation time, it is not suitable for problems where the shear viscosity and thermal conductivity are both important. If the relaxation time is chosen to produce correct shear viscosity, it cannot produce correct thermal conductivity, and vise versa.

In this paper, we propose a kinetic model equation which simplifies the collision integral of the semi-classical Boltzmann equation but keeps its main physical characters. Not only does it retain the desired properties of the BGK model equation, but it also produces the correct Prandtl number for the Fermi gases. The paper is organized as follows. We introduce the semi-classical Boltzmann equation and its equilibrium properties in §2, and develop a kinetic model to simplify the collision integral of the semi-classical Boltzmann equation in §3. In §4, we demonstrate how to recover the correct Prandtl number for the Fermi gases. Finally, to validate the proposed model, we numerically solve the kinetic model equation and compare with the results of other numerical methods in §5.

## 2. Semi-classical Boltzmann equation

For a single-component quantum gas, the semi-classical Boltzmann equation takes the form of (Uehling & Uhlenbeck 1933)
2.1
where *f*(**r**,**p**,*t*) is the one-body distribution function, with **r** being the spatial coordinates, **p** the momentum coordinates, *t* the time, *m* the atom mass and *U*(**r**,*t*) the effective potential. *I*_{coll}[*f*] is the collision integral, which takes the complicated form of
2.2
where *w*(*g*,*ϑ*) d*Ω* is the cross section with the relative velocity between the collision particles *g* and the post-collision angle *ϑ*. The element of solid angle is given by , where *ψ* is the azimuthal angle about *g*. The differential d*ϕ*_{1}∝*d***p**_{1} (Uehling & Uhlenbeck 1933). The subscript 1 denotes functions and variables pertaining to the second molecule in the collision, while the primes denote the functions and variables after the collision. The Boltzmann equation for classical gases is recovered at *θ* = 0, while for *θ* = −1 and 1, the semi-classical Boltzmann equation describes the dynamics of the Fermi and Bose gases, respectively.

Once the distribution function *f*(**r**,**p**,*t*) is known, one can get the macroscopic quantities such as the mass density *ρ*, macroscopic velocity **v**, internal energy density *e*, pressure tensor *P* and heat flux **Q**. For the gas in a *D*-dimensional (*D*=1,2,3) momentum space, these macroscopic quantities are given by
2.3
where *h* is the Planck constant and **u**=**p**/*m*−**v** is the peculiar velocity. The subscripts *α* and *β* represent the spatial coordinates. The Einstein summation convention for repeated indexes will be used throughout this paper.

Introducing the entropy function (Müller 1985)
2.4
to the Boltzmann equation (2.1), one gets the equilibrium distribution function
2.5
where *μ* is the chemical potential, *k*_{B} is the Boltzmann constant and *T* is the temperature.

When the collisions between atoms are so frequent that the local equilibrium can always be established, one can approximate the distribution function *f* by the equilibrium distribution function *f*_{eq}, which yields
2.6
and
2.7
*e*=*P*_{αα}/2, and **Q**=0, where *δ*_{αβ} is the Kronecker’s delta function, and
2.8
is the Bose–Einstein (*θ*=1) or Fermi–Dirac (*θ*=−1) function, with the local fugacity
2.9
and gamma function *Γ*(*n*). For the Bose gases, 0<*z*≤1 (with *z*=1 for the onset of Bose–Einstein condensation), while for the Fermi gases, *z*>0. When , , the quantum gas is in the near-classical limit, where the equilibrium distribution function is very close to the Maxwellian equilibrium distribution function for classical gases,
2.10
and the behaviour of quantum gases is similar to the classical ones. While (or ), the Bose (or Fermi) gas is in the nearly degenerate limit.

## 3. Kinetic model equation

Several considerations should be taken into account when simplifying the collision integral of the semi-classical Boltzmann equation (Struchtrup 2005). First of all, the simplified collision term must guarantee conservation of mass, momentum and energy. Secondly, like the Boltzmann equation, the kinetic model equation must be rotationally invariant. Thirdly, in the equilibrium state, the distribution function must reduce to the Fermi–Dirac (or Bose–Einstein) equilibrium distribution function for the Fermi (or Bose) gases. Fourthly, the entropy dissipation property (H-theorem), which states that the production of the entropy is always positive and vanishes only if the system is in equilibrium, must be satisfied. Finally, the shear viscosity and thermal conductivity derived from the kinetic model equation should coincide with those obtained from the semi-classical Boltzmann equation.

For classical gas dynamics, various kinetic models have been proposed, say, the BGK model (Bhatnagar *et al.* 1954), S model (Shakhov 1968) and Liu model (Liu 1990). None of them satisfies all the above five requirements. Among them, the well-known and mostly used model is the BGK model, which simplifies the collision integral as
3.1
with a single relaxation time *τ*. It satisfies the first four requirements, but fails the fifth. This problem is tackled in the ellipsoidal statistical BGK (ES–BGK) model, where an additional parameter was introduced (Holway 1966; Andries *et al.* 2000). Here, we shall construct a semi-classical kinetic model analogous to the classical ES–BGK model, expecting it retains the main physical properties of the collision integral of the semi-classical Boltzmann equation.

### (a) Ellipsoidal reference distribution function

We approximate the collision integral by −(*f*−*f*_{r})/*τ* so that the semi-classical Boltzmann equation (2.1) becomes
3.2
where *f*_{r} is the reference distribution function. When *f*_{r} reduces to *f*_{eq}, it becomes the BGK model equation. Note, the collision integral must satisfy the conservation conditions for mass, momentum and energy.

Notice that the reference distribution function in the classical ES–BGK model is determined by maximizing the Shannon entropy (Holway 1966), which is equivalent to the maximization of the classical entropy function *S*(*f*_{r}) (see equation (2.4) with *θ*=0) subjected to the given constraints. Such an entropy maximization principle gives us the least biased description about a system for which only partial information (such as mass, momentum, and pressure) is known. Following the entropy maximization principle for indistinguishable bosons and fermions, one can get the reference distribution function under some constraints or information (Kapur 1989). Here, equivalently, we shall get the reference distribution function *f*_{r} by maximizing the corresponding entropy function *S*(*f*_{r}) under some constraints. Before doing so, let us recall that the special reference distribution function—the equilibrium distribution function *f*_{eq}—can be obtained by the maximum entropy principle under the constraints that the mass density, momentum and energy are conserved, , and . As will be shown in §4, the resulting BGK model equation cannot produce the correct Prandtl number. Therefore, we introduce additional information,
3.3
where *W*_{αα}=*P*_{αα} for the conservation of energy.

Let us first assume that *W*_{αβ}, a quantity that is related to the pressure tensor, is known. By means of the Lagrange’s multipliers *z*_{r},*C*_{α} and , we maximize the following function
and get the following anisotropic reference distribution function (the conservation of momentum leads to *C*_{α}=0)
3.4
with the multipliers *z*_{r}(**r**,*t*) and the matrix *λ*_{αβ}(**r**,*t*) being determined through the following equations
3.5a
and
3.5b
where |⋅| denotes the determinant of a matrix.

Although the Lagrange’s method gives an extremum, whether it gives a maximum or minimum has to be further determined. Suppose there exists another distribution function *f*_{r}′ which satisfies the same constraints as *f*_{r} does. Consequently, we have . According to equation (2.4), the entropy difference between the two distribution functions is
3.6
Let . Since it is a convex function, we have for any pairs of positive numbers *x*_{1} and *x*_{2}, where the equality holds when *x*_{1}=*x*_{2}. Consequently, choosing *x*_{1}=*f*_{r} and *x*_{2}=*f*_{r}′, one could find that *S*(*f*_{r})−*S*(*f*_{r}′) is always positive (or zero) when *f*_{r}≠*f*_{r}′ (or *f*_{r}=*f*_{r}′). Therefore, the reference distribution function *f*_{r}, as given by equation (3.4), maximizes the entropy.

We now determine *W*_{αβ} to completely get the reference distribution function. Indeed, there are infinite possible values for *W*_{αβ} as long as *W*_{αα}=*P*_{αα}. A simple choice leading to the rotational invariance of the kinetic model equation is (Holway 1966)
3.7
where the hydrodynamic pressure *p* is defined as the average of the diagonal components of the pressure tensor *P*_{αβ}, i.e. *p*=*D*^{−1}*P*_{αα}, and *b* is a free parameter. However, for *f*_{r} to be of a finite norm, the matrix *λ*_{αβ} should be positive definite, that is, its eigenvalues must be non-negative. As can seen from equation (3.5b), this is equivalent to requiring *W*_{αβ} to be positive definite. From equation (3.7), it follows that *b* should satisfy
3.8
for *D*=2 and *D*=3. Specifically, when *D*=3, the region is the same as that in the ES–BGK model for the classical gases. Note that for *D*=1, *P*_{αβ} vanishes and *b* loses its meaning.

When *b*=0, the kinetic model reduces to the quantum BGK model. So the present kinetic model contains the quantum BGK model as a special case. Also, in the near classical limit , the kinetic model reduces to the ES–BGK model for classical gases (Holway 1966): from equations (3.4), (3.5) and (3.7), we have
3.9
where *λ*_{αβ}=*W*_{αβ}/*ρ*.

Now the kinetic model is completely described by equations (3.2), (3.4), (3.5) and (3.7), except for the relaxation time *τ* and parameter *b*. As it will be shown in §4, the two parameters will be determined by the coefficients of shear viscosity and thermal conductivity derived from the semi-classical Boltzmann equation.

### (b) Road to local equilibrium

From equation (3.2), one can see that collisions drive the system towards the states with anisotropic momentum distributions. This raises a question whether the local equilibrium can be achieved or not. To address this question, we consider a spatially homogeneous problem, where equation (3.2) reduces to ∂*f*/∂*t*=−(*f*−*f*_{r})/*τ*. Multiplying this equation by *v*_{α}*v*_{β}−*v*^{2}*δ*_{αβ}/3 and integrating over the momentum spaces, we have
3.10
where *σ*_{αβ}=*P*_{αβ}−*pδ*_{αβ} is the traceless stress tensor. Since *b*<1, the traceless stress tensor *σ*_{αβ} always approaches zero. This means that, eventually, we have
3.11
Then, from equations (2.6), (2.7) and (3.5), we have , where
3.12
It can be shown that is a monotonically decreasing function of *z* (see the Appendix A). Therefore, leads to *z*_{r}=*z*. Furthermore, with equations (2.6), (2.7) and (3.5), we have *λ*_{αβ}=*δ*_{αβ}*k*_{B}*T*/*m*. This means that in equilibrium, the reference distribution function *f*_{r} reduces to the equilibrium one *f*_{eq}.

The road to local equilibrium can thus be understood as follows: as *f* approaches *f*_{r}, *f*_{r} gradually adjusts itself to the local equilibrium distribution function (since |*σ*_{αβ}| decreases), and eventually we have
3.13

In fact, equation (3.13) can be interpreted in another way. That is, if *f*=*f*_{r}, we must have *f*=*f*_{r}=*f*_{eq}. When *f*=*f*_{r}, all the moments of *f*−*f*_{r} vanish. Since the equality holds when *σ*_{αβ}=0, we have *W*_{αβ}=*pδ*_{αβ}, i.e. equation (3.11). As we have discussed above from equation (3.11) we can obtain equation (3.13).

### (c) Rotational invariance of the kinetic model equation

The semi-classical Boltzmann equation is invariant under the rotation of spatial coordinates and the corresponding rotation of momentum coordinates. Since the reference distribution function *f*_{r} has a direction-dependent tensor *λ*_{αβ}, it may seem that *f*_{r} is not rotationally invariant, and hence equation (3.2) cannot be used to simplify the semi-classical Boltzmann equation. This is not true. In fact, as the momentum coordinates rotate, the pressure tensor *P*_{αβ} and *λ*_{αβ} change in a unique way that makes *f*_{r} rotationally invariant.

For simplicity, we first consider the case where only the *x*,*y* spatial coordinates rotate. That is, the new spatial coordinates **r**′=(*x*′,*y*′,*z*′) are
3.14
Of course, the momentum space should rotate exactly the same way, where the new momentum components **p**=(*p*_{x}′,*p*_{y}′,*p*_{z}′) are
3.15
After the rotation, the reference distribution function becomes
3.16
with
3.17
where *W*_{α′β′}′=(1−*b*)*p*′*δ*_{α′β′}+*bP*_{α′β′}′ and *p*′=*D*^{−1}*P*_{α′α′}′ are defined in terms of the new pressure tensor (because of the rotation of momentum coordinates)
3.18

Under the spatial and momentum rotations, the left side of the kinetic model equation (3.2) remains unchanged, just like the semi-classical Boltzmann equation. Therefore, to demonstrate the rotational invariance of the kinetic model equation, one has to prove *f*_{r}=*f*_{r}′. This can be done as follows. First, with the help of mathematical software such as Maple, it can be easily shown that |*W*_{αβ}|=|*W*_{α′β′}′|. Then the comparison between equations (3.5) and (3.17) leads to , which is equivalent to *z*_{r}=*z*_{r}′. Also, we have |*λ*_{αβ}|=|*λ*_{α′β′}′|. With this result, we find that the proof of is equivalent to the proof of . The latter can be verified by Maple.

Since the kinetic model equation (3.2) is rotational invariant under the above two-dimensional rotations, it is definitely rotational invariant for general three-dimensional rotations. This is because a three-dimensional rotation can be decomposed into two-dimensional rotations, i.e. firstly the *x*,*y* (and correspondingly *p*_{x},*p*_{y}) axes rotate, secondly the *y*,*z* axes rotate, and finally the *x*,*z* axes rotate.

### (d) Proof of the H-theorem

The H-theorem shows that the entropy production due to collision is always non-negative and vanishes only if the system is in equilibrium. According to equation (2.4), the entropy production due to the collision in the kinetic model equation (3.2) is
3.19
For the Bose gases (*θ*=1), is a convex function. Applying the corresponding convexity inequality to equation (3.19), one gets
3.20
where the two equalities hold when *f*=*f*_{r} (the last integral is non-negative and equals zero when *f*=*f*_{r}, since is non-negative). For the Fermi gases, one can use the convexity inequality for . Similarly, one has
3.21
where the equality holds when *f*=*f*_{r}. Therefore, for both Bose and Fermi gases, we have *τS*^{(c)}≥*S*(*f*_{r})−*S*(*f*), where the equality holds only when *f*=*f*_{r}.

In §3*b*, we have already shown that when *f*=*f*_{r}, we will have *f*=*f*_{r}=*f*_{eq}. Thus, to prove the H-theorem, one has to prove *S*(*f*_{r})≥*S*(*f*). It is more convenient to introduce the distribution function *f*_{s} which maximizes the entropy function (2.4) under the constraints of reproducing the actual density, macroscopic velocity and pressure tensor from *f*. Accordingly, *f*_{s} has the same form of *f*_{r}, but with *z*_{r}, *λ*_{αβ} and *W*_{αβ} being, respectively, replaced by *z*_{s}, , and *P*_{αβ}:
3.22
with
3.23
By definition, *S*(*f*)≤*S*(*f*_{s}), where the equality holds when *f*=*f*_{s}.

Now we have to prove *S*(*f*_{r})≥*S*(*f*_{s}). From equations (2.4) and (3.4), we obtain *S*(*f*_{r})−*S*(*f*_{s}) = *ρ*[*y*(*z*_{r})−*y*(*z*_{s})]/*m*, where . Numerical calculation shows that *y*(*x*) is a decreasing function of *x*. Therefore, to prove the H-theorem, it is equivalent to prove *z*_{r}≤*z*_{s}.

For *W*_{αβ} and *P*_{αβ} which satisfy equation (3.7) with *b*≥−1/(*D*−1), we can use the Brunn–Minkowsky inequality |*W*_{αβ}|≥|*P*_{αβ}| (Andries *et al.* 2000), so that equations (3.5) and (3.23) give . Consequently, we have *z*_{r}≤*z*_{s}. This means that the entropy production is always positive and the entropy reaches its maximum only when *f*=*f*_{r}=*f*_{eq}. Therefore, the kinetic model equation (3.2) satisfies the H-theorem.

## 4. Recovery of the correct Prandtl number

The relaxation time *τ* and parameter *b* can be determined by equating the shear viscosity and thermal conductively derived from the kinetic model equation (3.2) with those obtained from the semi-classical Boltzmann equation (2.1).

### (a) Shear viscosity and thermal conductivity from the kinetic model equation

We first calculate the shear viscosity and thermal conductivity of the kinetic model equation by means of the Chapman–Enskog expansion (Chapman & Cowling 1970). The basic idea of this expansion is to expand the distribution function around the local equilibrium in terms of the Knudsen number, which gives the Euler equations at the zeroth-order. For the first-order approximation (Navier–Stokes level), we assume a solution of the kinetic model equation (3.2) of the form:
4.1
where *f*^{(1)}(**r**,**p**,*t*) denotes a small deviation from the local equilibrium.

To calculate *f*^{(1)}(**r**,**p**,*t*), we need to express the reference distribution function *f*_{r} in terms of the local macroscopic quantities and the Knudsen number. For a system not far away from local equilibrium, the traceless stress tensor *σ*_{αβ} is of order . It can be shown from equation (3.5) that , , and
4.2
which allows the power-series expansion of the reference distribution function *f*_{r} in terms of the Knudsen number as
4.3
Substituting equations (4.1) and (4.3) into equation (3.2), we get the expression for *f*^{(1)} as
where the last term can be expressed as , with being given by (see Nikuni & Griffin (1998) for detailed calculations with *D*=3)

On substituting equation (4.1) into expressions for the heat flux *Q* and pressure tensor *P*_{αβ}, we recover the Fourier’s law of heat conduction *Q*=−*κ*∇_{r}*T*, with the thermal conductivity given by
4.4
and the Newton’s friction law for the stress tensor *σ*_{αβ}=−*η*(∂*v*_{α}/∂*β*+∂*v*_{β}/∂*α*) (*D*=2 or 3, *α*≠*β*), with the shear viscosity given by
4.5

Note that in the BGK model, the shear viscosity is *η*=*pτ*. Also, in the near-classical limit (), the shear viscosity and thermal conductivity depend on the product of relaxation time *τ* and pressure *p* only, i.e. *κ*=(*D*+2)*τpk*_{B}/2*m* and *η*=*pτ*/(1−*b*). So the Prandtl number *Pr*=(*D*+2)*k*_{B}*η*/2*mκ*=2/3 for classical monatomic gases can be reproduced when *b* takes the minimum allowed value of (Holway 1966). When the quantum effects are taken into account, the expression for the shear viscosity remains unchanged but the thermal conductivity depend also on the local fugacity *z*. If one continues to use the Prandtl number *Pr*=(*D*+2)*k*_{B}*η*/2*mκ* for quantum gases, following equations (4.4) and (4.5), the Prandtl number of the semi-classical kinetic model equation (3.2) is
4.6

Since *b* must be larger than to ensure that the reference distribution function does not diverge in the momentum space and to guarantee the H-theorem, the kinetic model equation can be used to reproduce a Prandtl number of the semi-classical Boltzmann equation if it is not smaller than .

### (b) Shear viscosity and thermal conductivity of the Fermi gases

We consider two-component balanced Fermi gases. In most experiments, the two components move together and one needs only to use one distribution function (Riedl *et al.* 2008). Owing to Pauli’s exclusion principle, the *s*-wave scattering happens between atoms with different spins. The collision integral between different components is given by equation (2.2) with d*ϕ*_{1}=*d***p**_{1}/*h*^{3} and the differential cross section *w*(*g*,*ϑ*)=*a*^{2}/[1+(*πmga*/*h*)^{2}], where *a* is the *s*-wave scattering length.

Applying the Chapman–Enskog expansion (Chapman & Cowling 1970) to the semi-classical Boltzmann equation (Kavoulakis *et al.* 2000; Watabe *et al.* 2010), one obtains the thermal conductivity as
4.7
and the shear viscosity as
4.8
where
4.9
4.10
and
4.11
Here, , , , and *T*_{a}=*h*^{2}/4*π*^{2}*mk*_{B}*a*^{2}. Note that when the *s*-wave scattering length *a* is very small (corresponding to ), the results are the same with those reported in Watabe *et al.* (2010).

Derived from the semi-classical Boltzmann equation, the numerical calculation of equations (4.7) and (4.8) shows that the Prandtl number . This is larger than the minimum Prandtl number given by the kinetic model equation (the reason is for the Fermi gases), e.g. see the solid and dashed lines in figure 1*b*. Consequently, when *b* takes the values of
4.12
the correct Prandtl numbers can be recovered. Note that *b* increases from the minimum value at the near-classical limit to its maximum value 1 at the degenerate limit (see the solid line in figure 1*c* or figure 1*d*). It is interesting to note that when and the fugacity *z*≈22, we have *b*≈0 (figure 1*c*), where the kinetic model equation reduces to the quantum BGK model equation.

When the value of *b* is determined, the relaxation time *τ* can be determined by equalling equation (4.5) with equation (4.8), which yields
4.13

### (c) Prandtl number of the Bose gases

For the dilute Bose gases under the *s*-wave interaction, the Prandtl number from the Chapman–Enskog expansion of the semi-classical Boltzmann equation is (Nikuni & Griffin 1998; Kavoulakis *et al.* 2000)
4.14
where *I*_{A} and *I*_{B} are still given by equations (4.9) and (4.10), respectively, but *z* in equation (4.11) is replaced by −*z*.

Numerical calculation shows that . Since for the Bose gases, in general, the kinetic model equation could not yield the correct Prandtl number (figure 2). Figure 2*b* shows that our kinetic model only recovers the true Prandtl number in the classical limit where the fugacity *z* approaches zero. As *z* increases, the difference in the Prandtl number between our model and that of the full Boltzmann equation grows. However, when , the Prandtl number given by our model is closer to the realistic one than that of the quantum BGK model, where *b*=0.

## 5. Numerical simulations

In this section, we consider a one-dimensional shock wave problem to numerically validate the proposed kinetic model equation for the Fermi gases in the normal phase. We consider the case where the *s*-wave scattering length *a* is very small, so that *T*/*T*_{a} approaches zero in equations (4.9) and (4.10). For simplicity, we ignore the effective potential. Note that some numerical solutions of the quantum Euler equation (Yang *et al.* 2007; Hu & Jin 2011), the semi-classical BGK model equation (Shi & Yang 2008; Yang & Hung 2009; Muljadi & Yang 2011), and the semi-classical Boltzmann equation (Filbet *et al.* 2011) have already been obtained.

### (a) Normalization

We first introduce the following dimensionless variables
5.1
where *L* is the characteristic length of the flow field, *T*_{0} is the reference temperature, *z*_{0} is the reference fugacity and *ρ*_{0} is the reference mass density given by equation (2.6) when *T*=*T*_{0} and *z*=*z*_{0}. With these dimensionless variables, the mass density, pressure and internal energy density are given by
5.2
The *hat* symbol will hereafter be omitted. Then equation (3.2) becomes
5.3
where *f*_{r} can be obtained from equation (3.4), divided by , with *z*_{r} and *λ*_{αβ} being uniquely determined by
5.4
and the parameter *ϵ* is related to the Knudsen number as
5.5
which reduces to in the classical limit when *z*=0.

### (b) Asymptotic preserving scheme

Several numerical schemes have been proposed to solve the BGK kinetic model equation (5.3) when *b*=0 (Shi & Yang 2008; Yang & Hung 2009; Muljadi & Yang 2011). These methods can be extended to solve the general case where *b*≠0. Also, various numerical schemes developed for dynamics of classical kinetic model equations can be used to solve equation (5.3). Here, we adopt the asymptotic preserving scheme (Filbet & Jin 2011), which allows the consistent treatment of problems with a broad range of Knudsen numbers. This property is important because the whole flow field could be in both hydrodynamic and kinetic flow regimes. In the hydrodynamic limit, this scheme can capture macroscopic gas dynamics even if the small scale determined by the relaxation time *τ* is not numerically resolved. The computational accuracy is guaranteed by the fact that, using the Chapman–Enskog expansion, this numerical scheme yields the correct Euler equations when we allow *τ* to approach zero at the fixed spatial and time steps. Therefore, the computation of hydrodynamic flow can be as fast and accurate as transition and collisionless flows.

In the temporal discretization of equation (5.3), the stream part is treated explicitly but the collision part is treated implicitly, resulting in
5.6
where Δ*t* is the time step and the variables with superscript *n* denote the values of these variables at the *n*th time step. The spatial discretization of *ξ*_{x}∇_{x}*f*^{n} is given by a second-order finite volume scheme with minmod slope limiter. Given the value of *f*^{n}, *f*^{n+1} is calculated according to
5.7
where can be solved explicitly, because the macroscopic quantities at the (*n*+1)th time step can be solved explicitly.

Firstly, according to the conservation of mass, momentum and energy, we can get *ρ*^{n+1}, **v**^{n+1} and *p*^{n+1} through
5.8
where the total energy *E*=*ρv*^{2}/2+3*p*/2 (the integration over the velocity spaces is carried out by direct sum of the integral kernel over the discrete velocity ordinates). Secondly, to get the value of , one needs the value of . The latter can be obtained by introducing the tensor . Multiplying equation (5.7) with *ξ*_{α}*ξ*_{β} and integrating it with respect to ** ξ**, we have
5.9
Once the value of is determined, one can obtain the value of , via numerical root finding methods, from the following nonlinear equation
5.10
and the value of from the equation
5.11
The values of , and

**v**

^{n+1}completely determine . Hence from equation (5.7), one get the distribution function

*f*

^{n+1}.

Note that at the final time step, the macroscopic quantities *z* and *T* can be obtained by solving the first two equations in equation (5.2).

### (c) Numerical results

The asymptotic preserving scheme is used to numerically simulate the Fermi gases in the near classical and degenerate limits, from the hydrodynamic flow regime to transition flow regime. In the following examples, we use the Neumann boundary condition ∂*f*/∂*x*=0.

### Example 5.1

The two-component Fermi gases in the near-classical limit with small Knudsen number, *Kn*=10^{−4}, so that the system is in the continuum regime where the quantum Euler equation can be derived from the semi-classical Boltzmann equation (Yang *et al.* 2007; Hu & Jin 2011). The initial distribution function at *t*=0 is given by the equilibrium distribution function *f*_{eq} with the following parameters
5.12
and *z*_{0}=0.0649, so the initial mass density *ρ* is 1 when *x*≤0.5, and 0.4 otherwise. Figure 3*a*,*b* shows that our results (circles) are in good agreements with the exact solutions (solid lines) of the quantum Euler equation (Hu & Jin 2011). Since the fugacity is very small, the behaviour of the quantum Fermi gases is similar to the classical gases.

### Example 5.2

The two-component Fermi gases in the nearly degenerate limit with *Kn*=10^{−4}. Similarly, the initial distribution function is set according to the equilibrium distribution function with the following parameters
5.13
and *z*_{0}=901.284. The initial mass density is 1 when *x*≤0.5, and 0.4 otherwise. Figure 3*c*,*d* illustrates that our results (circles) agree with the exact solutions (solid lines) of the quantum Euler equation well (Hu & Jin 2011). One could also find that the quantum Fermi gases in the near-classical limit and nearly degenerate limit behave quite differently.

### Example 5.3

The two-component Fermi gases with *Kn*=0.02. For the classical gases, the flow is in the slip-flow regime, where the Navier–Stokes equations can describe the flow well with appropriate velocity-slip and temperature-jump boundary conditions at the wall surfaces (Reese *et al.* 2003). The initial distribution function is given by the equilibrium one with the following parameters
5.14
and *z*_{0}=3.1887. The circles in figure 4 depict our numerical results of the two-component Fermi gases under the *s*-wave scattering. We also plot the results (solid lines) of the Fermi gas using the Maxwellian interaction potential with *Kn*=0.01 (Filbet *et al.* 2011) for comparison. Note that we consider two-component Fermi gases where the collisions happen between different components, while Filbet *et al.* (2011) considered one-component Fermi gas. This indicates that the viscosity and thermal conductivity in our system are half of those considered in Filbet *et al.* (2011). Therefore, the Knudsen number in our model should be scaled by a factor of 2 in comparison with that in Filbet *et al.* (2011). Indeed, one can see in figure 4 that our numerical results of the kinetic model equation agree well with the solution of the semi-classical Boltzmann equation. Since the calculation of the collision integral in the semi-classical Boltzmann is cumbersome and time-consuming, the kinetic model equation (3.2) provides a computationally efficient approach to study the flow dynamics of quantum Fermi gases.

### Example 5.4

The two-component Fermi gases in the nearly degenerate limit with *Kn*=0.5. For the classical gases, the flow is in the transition regime where the Navier–Stokes equations fail. The initial distribution function is given by the equilibrium one with the following parameters
5.15
and *z*_{0}=900. In this case, we have no available data to compare with, so we compare the dynamics of Fermi gas in the nearly degenerate regime with that in the near-classical regime (*z*≤0.065). From figure 5, we see that the dynamics of shock waves for the Fermi gases in the near-classical limit and in the nearly degenerate limit are quite different. In particular, we have observed two counter-gradient heat flux regions for the Fermi gas in the nearly degenerate limit (solid lines), where the direction of heat flux is from the cold region to the hot region. This phenomenon, which is in contrast to the Fourier’s heat transfer law, clearly indicates that the failure of the Navier–Stokes–Fourier equations in describing the quantum gas dynamics with large Knudsen numbers.

## 6. Concluding remarks

We have proposed a kinetic model equation for the quantum gas dynamics, which simplifies the complicated collision term of the semi-classical Boltzmann equation, but retains its main physical characters. The proposed kinetic model contains the quantum BGK model as a special case. Also, in the near-classical limit, the kinetic model reduces to the ellipsoidal statistical BGK model for classical gases. We have also derived the expressions of the shear viscosity and thermal conductivity from the kinetic model equation by means of the Chapman–Enskog expansion, and compared them with those derived from the semi-classical Boltzmann equation. It was found that our kinetic model equation yields the correct Prandtl number for the Fermi gases in the normal phase. Meanwhile, for the Bose gases, the Prandtl number given by our kinetic model is also closer to the realistic value than that of the quantum BGK model. Furthermore, to validate the proposed model, we have conducted numerical simulations on the one-dimensional shock wave problems, using the asymptotic-preserving scheme. It was found that the numerical results of the kinetic model equation agree well with the solutions of the quantum Euler and semi-classical Boltzmann equations. This implies that one can use the kinetic model equation to study dynamics of quantum Fermi gases for practical purposes.

## Acknowledgements

The research leading to these results has received funding from the European Communitys Seventh Framework Programme FP7/2007-2013 under grant agreement ITN GASMEMS no 215504. This work is also financially supported by the Engineering and Physical Sciences Research Council U.K. under Grants No. EP/EP/I036117/1/1.

## Appendix A. Monotonicity of

For the Bose gases, if is a monotonically decreasing function of *z*, its first derivative with respect to *z* should always be negative. Therefore, it is equivalent to show
A1
where *D*=2 or *D*=3. Obviously, one could find that equation (A1) holds if
A2

Equation (A2) can be proved as follows. According to the serial expansion of the Bose–Einstein function, i.e. , the coefficient of the term *z*^{k} (*k*≥2) in the left side of equation (A2) is
while the coefficient of the term *z*^{k} in the right side of equation (A2) is
Since the following inequality always holds,
we have *C*_{2}>*C*_{1}. Hence equation (A2) is proved. And for the Bose gases, is a monotonically decreasing function of *z*.

For the Fermi gases, however, we can only prove that is a monotonically decreasing function of *z* when *z* is very large. In this case according to the Sommerfeld expansion, the Fermi–Dirac function can be approximated as (the term with the under brace is zero since *z* is very large)
A3
On substituting this expression into equation (3.12), one can easily show that is a monotonically decreasing function of *z*, with and approaching 0.4836 and 0.5 as *z* approaches infinity. When *z* is not very large, our numerical calculations show that is a monotonically decreasing function of *z*, see figure 6 where *z* goes from zero to 10^{5}.

- Received November 11, 2011.
- Accepted February 7, 2012.

- This journal is © 2012 The Royal Society