## Abstract

The ideal quantum gas dynamics as manifested by the semiclassical ellipsoidal-statistical (ES) equilibrium distribution derived in Wu *et al.* (Wu *et al*. 2012 *Proc. R. Soc. A* **468**, 1799–1823 (doi:10.1098/rspa.2011.0673)) is numerically studied for particles of three statistics. This anisotropic ES equilibrium distribution was derived using the maximum entropy principle and conserves the mass, momentum and energy, but differs from the standard Fermi–Dirac or Bose–Einstein distribution. The present numerical method combines the discrete velocity (or momentum) ordinate method in momentum space and the high-resolution shock-capturing method in physical space. A decoding procedure to obtain the necessary parameters for determining the ES distribution is also devised. Computations of two-dimensional Riemann problems are presented, and various contours of the quantities unique to this ES model are illustrated. The main flow features, such as shock waves, expansion waves and slip lines and their complex nonlinear interactions, are depicted and found to be consistent with existing calculations for a classical gas.

## 1. Introduction

The kinetic Boltzmann equation has been commonly used to describe various transport phenomena in classical dilute gases covering a wide range of flow parameters, such as Reynolds numbers, Mach numbers and Knudsen numbers. Analogous to the classical Boltzmann equation, a semiclassical Boltzmann equation, which generalizes the collisions term to treat collisions of particles of quantum statistics, has been developed by Uehling & Uhlenbeck [1] (see [1,2]). In recent years, owing to the rapid advancements of nanotechnology, the device or structure characteristic length scales become comparable to the mean free path and the wavelength of energy and information carriers (mainly electrons, photons, phonons and molecules), and so some of the classical continuum transport laws are no longer applicable. It is generally believed that the microscopic description of the Boltzmann equation (classical or semiclassical) is adequate to treat transport phenomena in the mesoscale range. Hydrodynamic behaviour of quantum gases has been the subject of some prominent research [3–5], and the application of quantum Boltzmann hydrodynamic equations has been implemented in the analysis of electron flows in quantum semiconductor devices [6–8]. Owing to the different types of carriers involved simultaneously in a single problem, it is desirable to have a method that can allow them to be treated in a unified and parallel manner. With the semiclassical Boltzmann equation, it is possible to describe adequately the mesoscale transport of particles of arbitrary statistics.

The main difficulty encountered in solving the semiclassical Boltzmann equation is similar to that encountered in the classical counterpart, and is mainly due to the complicated integral nature of its collision term. The relaxation time approximation proposed by Bhatnagar, Gross and Krook (BGK) [9] for the classical Boltzmann equation provides a simpler form of the collision term, retains the principal effects of particle collisions and enables more tractable solution methods. The relaxation time approximation is quite general and is applicable to the semiclassical Boltzmann equation. The only change is that the equilibrium Maxwell–Boltzmann distribution in the classical case is replaced by the Bose–Einstein or Fermi–Dirac distribution, depending on the types of carrier particles considered. Indeed, the semiclassical Boltzmann–BGK equation has been widely used in electron carrier transport in semiconductor [10–17] and phonon energy transport in thermoelectric materials [18]. The main shortcoming of the classical BGK model is that the Prandtl number is equal to 1, while it is 2/3 from the full Boltzmann equation. In order to improve the classical Boltzmann–BGK model equation that does not produce correct Prandtl numbers, a new kinetic model called the ES model was proposed by Lowell & Holway [19]. Analogously, a semiclassical ES kinetic model for quantum gases in the normal phase has been recently developed by Wu *et al*. [20]. This new model mimics the classical ES model of Lowell & Holway [19], and was derived by maximizing the corresponding entropy function under some constraints. It possesses all the desirable features of a kinetic model and yields a correct Prandtl number for Fermi gases. In the semiclassical Boltzmann–BGK equation, if the equilibrium Bose–Einstein or Fermi–Dirac distribution is replaced by the semiclassical ES distribution, then one has the semiclassical Boltzmann–ES–BGK equation [20].

Theoretically, the Chapman–Enskog expansion method is usually applied to the full Boltzmann or model Boltzmann equation (classical or semiclassical) to yield sets of governing hydrodynamic transport equations, such as Euler, Navier–Stokes and Burnett equations, which are, respectively, the zeroth-order, first-order and second-order expansions in the parameter Knudsen number [21]. The zeroth-order solution yields the equilibrium distribution in phase space and the Euler equations in physical space, and governs the ideal gas dynamics. The semiclassical ES equilibrium distribution is the zeroth-order solution of the semiclassical Boltzmann–ES–BGK equation. The main objective of this work is to investigate the ideal gas dynamics governed by the semiclassical ES equilibrium distribution, as it has not been well explored previously. The main numerical method adopted is a direct solver in phase space for the semiclassical ES equilibrium distribution function. The direct solver consists of two methods: one is the discrete ordinate method and the other the high-resolution shock-capturing method. In classical rarefied gas flow computation, the implementation of the discrete ordinate method to nonlinear model Boltzmann equations has been developed by Yang & Huang [22]. Such a direct method will allow one to examine the same physical flow problems, but with different gas particles. Also, if the classical limit situations of the same flow problem are considered, then one expects to obtain similar or identical flow structures for the three statistics. It is noted that even with the classical Maxwell–Boltzmann statistics, the present formulation includes fugacity, which is not seen in the original Boltzmann–BGK equation [9] and most other existing works based on it. Computer simulations of several two-dimensional Riemann problems in gas flows of arbitrary statistics are shown to illustrate the complex nonlinear wave manifestation of ideal quantum gas dynamics as governed by the semiclassical ES equilibrium distribution.

This paper is organized as follows. Elements of the semiclassical Boltzmann–ES–BGK equation are briefly described in §2. Its relationship with hydrodynamic equations is outlined. In §3, governing equations in two space dimensions and their non-dimensional forms are detailed. In §4, the present direct solution method in phase space is described. The implementation of the discrete ordinate method and total variation diminishing (TVD) schemes and the numerical quadrature to evaluate the macroscopic properties are detailed. A decoding procedure to determine the key parameters for determining the equilibrium ES distribution is given. In §5, numerical experiments of two-dimensional semiclassical Riemann problems of ideal gas dynamical flows are presented to illustrate the present method as well as the flow physics. Finally, concluding remarks are given in §6.

## 2. Semiclassical Boltzmann–ellipsoidal-statistical-Bhatnagar–Gross–Krook equation

In this section, the elements of the semiclassical Boltzmann kinetic equation are briefly described. Following [1,2], we consider the extension of the Boltzmann equation to quantum systems due to Uehling & Uhlenbeck, in which they took the Pauli exclusion principle into account,
2.1where *m* is the particle mass, *U* is the mean field potential and *f*(** p**,

**,**

*x**t*) is the distribution function that represents the average density of particles with momentum

**at the space–time point**

*p***,**

*x**t*. The parameter (

*δf*/

*δt*)

_{coll}. denotes the collision term, and, according to Uehling & Uhlenbeck [1], it takes the form 2.2where

*K*(

**p**,

**q**,

*Ω*) denotes the collision kernel,

*Ω*is the solid angle and

*θ*is a parameter that specifies the type of particle statistics. Here, for

*θ*=+1, Bose–Einstein particles are considered, for

*θ*=−1, Fermi–Dirac particles, and for

*θ*=0, Maxwell–Boltzmann classical particles are considered. It is well known that the collision integral of the semiclassical Boltzmann equation will automatically vanish when the gas distribution function is in equilibrium, and the equilibrium distribution function for general statistics can be expressed as 2.3where

**u**(

**x**,

*t*) is the mean velocity,

*T*(

**x**,

*t*) is temperature,

*k*

_{B}is the Boltzmann constant and is the fugacity, where

*μ*is the chemical potential. In (2.3), again,

*θ*=−1, denotes the Fermi–Dirac statistics,

*θ*=+1, the Bose–Einstein statistics, and

*θ*=0 denotes the Maxwell–Boltzmann statistics. To avoid the mathematical difficulty caused by the nonlinear integral collision term, the relaxation time approximation of BGK is generally applied to replace the collision term of Uehling & Uhlenbeck, thus the semiclassical Boltzmann–BGK equation reads 2.4where

*τ*is the relaxation time, which needs to be specified for each carrier scattering, and

*f*

^{eq}is given by equation (2.3). Inspired by the classical ES kinetic model [19] and using the maximum entropy principle, a new kinetic model equation has been proposed for dilute quantum gases in the normal phase [20]. This semiclassical ES kinetic model equation preserves the main properties of the semiclassical Boltzmann equation, including conservation of mass, momentum and energy, the entropy dissipation property and the rotational invariance. It also can provide correct Prandtl numbers. This semiclassical Boltzmann–ES–BGK equation can be expressed as 2.5where

*f*

^{ES}is an anisotropic reference equilibrium distribution function given by 2.6In equation (2.6),

*λ*

_{αβ}is a matrix and is its inverse, and

*C*

_{α}=(

*p*

_{α}/

*m*)−

*u*

_{α}are the components of the thermal velocity. The reference equilibrium distribution function

*f*

^{ES}is derived by the maximum entropy principle under some constraints that the mass, momentum and energy are conserved and the following relations are satisfied: 2.7a 2.7b and 2.7cwhere

*D*is the dimension of momentum space, and

*W*

_{αα}=

*P*

_{αα}is required for the conservation of energy.

In principle, after solving the semiclassical Boltzmann–ES–BGK equation for the (non-equilibrium) distribution function *f*, all the macroscopic dynamic variables of interest such as number density, momentum density and energy density are low-order moments of the distribution function, and are defined by
2.8a
2.8b
and
2.8cOther higher order moments such as the pressure tensor *P*_{αβ} and the heat flux vector *Q*_{α}(**x**,*t*) can also be defined accordingly as
2.9and
2.10Lastly, to determine the semiclassical ES distribution *f*^{ES}, one usually needs to calculate *z* and the matrix *λ*_{αβ} simultaneously through
2.11and
2.12where ∥⋅∥ denotes the determinant of a matrix. To obtain *W*_{αβ} to completely determine *f*^{ES} and require that *W*_{αα}=*P*_{αα} to ensure conservation of energy, a natural choice leading to rotational invariance of the kinetic model equation is [19]
2.13where *p* is the hydrodynamic pressure defined as the average of the diagonal components of the pressure tensor *P*_{αβ}, i.e. *p*=*P*_{αα}/*D*, and *b* is a parameter. Requiring *W*_{αβ} to be positive definite, the parameter *b* satisfies
2.14for *D*=2 and *D*=3.

The conservation laws of macroscopic properties can be obtained by multiplying equation (2.5), respectively, by 1, **p** and **p**^{2}/2*m* and integrating the resulting equations over all **p**. Consequently, the integrals of the collision terms in all three cases vanish automatically, resulting in the conservation laws in differential equation form for the conserved macroscopic quantities, i.e. number density *n*(**x**,*t*), momentum density *mn***u**(**x**,*t*) and energy density *ϵ*(**x**,*t*) as follows:
2.15a
2.15b
and
2.15cIn this work, we will focus our effort on the case of the *semiclassical ES Euler* limit, in which the particle distribution function is always in reference equilibrium, i.e. *f*=*f*^{ES}, the collision term of equation (2.5) vanishes automatically and the set of governing equations, equations (2.15), will be termed semiclassical ES Euler equations in order to distinguish it from the semiclassical Euler equations governed by *f*=*f*^{eq}. It is noted that *f*^{ES} is not isotropic, it has quantities like *W*_{αβ} and is quite different from *f*^{eq}, the standard Bose–Einstein and Fermi–Dirac distributions, thus it is worth some detailed study to explore the gas dynamical flow features inherent in this peculiar equilibrium distribution before we take on the non-equilibrium gas flows as governed by the semiclassical Boltzmann–ES–BGK equation, i.e. equation (2.5). In §5, we will present the comparison of results with these two semiclassical Euler solutions. In recent years, notable numerical methods to describe the ideal quantum gas flows have been developed [23–25].

## 3. Governing equations in two space dimensions

Here, we give a detailed description of the theory and governing equation in two space dimensions. Without losing generality, we also neglect the influence of the externally applied field *U*(**x**,*t*). It is noted that quantum statistics of ideal gases in two dimensions has striking features such as the specific heat for an ideal gas of Fermi particles is identical with that for an ideal Bose gas [26]. The general method in three space dimensions can also be similarly formulated. When *f*=*f*^{ES}, the semiclassical Boltzmann–ES–BGK equation in two space dimensions can be expressed as
3.1with *υ*_{x} and *υ*_{y} as particle velocities. The macroscopic properties are the various moments of *f*^{ES}(*υ*_{x},*υ*_{y},*x*,*y*,*t*),
3.2a
3.2b
3.2c
and
3.2dAdditional information from *f*^{ES} can also be defined accordingly to
3.3a
3.3b
and
3.3cWe also have the gas pressure *p*(*x*,*y*,*t*)=(*W*_{xx}+*W*_{yy})/2. The pressure tensor *P*_{αβ}(*x*,*y*,*t*) can be obtained through
3.4For the conservation of energy, *W*_{xx}+*W*_{yy}=*P*_{xx}+*P*_{yy} is required. To obtain the new *f*^{ES}, one needs *z* and *λ*_{xx},*λ*_{xy}*λ*_{yy}, which can be determined by solving the following equations simultaneously:
3.5aand
3.5bComputationally, this requires a root finding procedure, and either the Newton–Ralphson method or the bisector method can be employed.

It is emphasized here that in semiclassical ES equilibrium flows, the main concern here is the *f*^{ES} and the macroscopic properties generated from the moments of *f*^{ES}, thus we first have *W*_{αβ} instead of *P*_{αβ}, and the pressure tensor quantities are computed from *W*_{αβ} and specified by parameter *b*. In the semiclassical Boltzmann–ES–BGK non-equilibrium flows, our concern will be *f*, and we will have the pressure tensor *P*_{αβ} instead of *W*_{αβ}.

### (a) Governing equations in non-dimensional form

Before proceeding to discretize the equation in this section, we introduce the characteristic properties of *V* _{0}, *t*_{0} and *n*_{0} (or *z*_{0}) for the purpose of normalization. Let us choose some reference parameters as
3.6with *L* defined as the characteristic length of the problem. Hence, the definitions of non-dimensional variables are introduced as
3.7After non-dimensionalizing the equation, we have the normalized semiclassical Boltzmann–ES equation in two space dimensions as follows:
3.8with and as particle velocities. The normalized two-dimensional semiclassical reference distribution function becomes
3.9where . From here on, we shall omit all the hat signs, and all the equations are in non-dimensional form.

### (b) Fermi and Bose functions in relation to hydrodynamic properties

For the sake of completeness and for later comparison, we also describe the case for *f*=*f*^{eq} in two space dimensions (*D*=2) in the non-dimensional form
3.10The macroscopic moments, i.e. number density *n*(*x*,*y*,*t*), momentum *j*(*x*,*y*,*t*) and energy density *ϵ*(*x*,*y*,*t*), can now be expressed in terms of quantum functions and are given by
3.11a
3.11b
and
3.11cand
3.11dwhere is the thermal wavelength and *β*=1/*k*_{B}*T*(*x*,*t*). The functions of order *ν* are, respectively, defined for either the Fermi–Dirac or the Bose–Einstein statistics as
3.12and
3.13Here, applies for Fermi–Dirac integral and for the Bose–Einstein integral, where *Γ*(*ν*) is the gamma function. The definition of macroscopic quantities in terms of Fermi or Bose functions applies for both cases of quantum distributions. One only needs to replace the Fermi function with the Bose function or vice versa while maintaining the same procedure.

## 4. A direct solver in phase space

### (a) Implementation of the discrete ordinate method

In the two-dimensional case, we may apply the Gauss–Hermite quadrature rule over the interval for each direction. The Gauss–Hermite quadrature rule reads
4.1or
4.2The discrete points *υ*_{α} and the corresponding weights *W*_{α}, with *α*=*σ*,*δ*, are tabulated in a table of Gauss–Hermite quadrature [27]. Applying the discrete ordinate method to equation (3.8) with (*υ*_{α}), then we have a set of partial differential equations in physical space (*x*,*y*,*t*) instead of a single equation in phase space,
4.3where is the value of *f*^{ES}(*υ*_{x},*υ*_{y},*x*,*y*,*t*) at the discrete velocity point (*υ*_{x}=*υ*_{σ},*υ*_{y}=*υ*_{δ}). Once the discrete distribution functions *f*_{σ,δ}(*x*,*y*,*t*) are solved, for every time level, we can acquire and update the macroscopic moments in the physical space using a selected quadrature method as described below,
4.4a
4.4b
4.4c
4.4d
4.4e
4.4f
and
4.4g

Once the values of *W*_{αβ} have been obtained, the gas pressure can be computed by using *p*(*x*,*y*,*t*)=(*W*_{xx}+*W*_{yy})/2, and the pressure tensor components by following relations (3.4).

Here, note that the selection criteria of suitable quadratures are to guarantee that the macroscopic moments can be accurately computed by requiring the accurate representation of the distribution function with suitable velocity ranges being covered. In the present work, Gauss–Hermite quadrature was employed; however, the equally spaced Newton–Cotes formula can also be used if one needs to cover the high-energy tail of the distribution function.

### (b) Implementation of high-resolution shock-capturing schemes

In this section, we describe the numerical algorithm to solve the set of equation (4.3). A class of TVD schemes is implemented to solve the set of conservation laws. We firstly discretize the space *x*,*y* and time *t* into a number of cells centred at *i*,*j* at time *n*, hence we approximate *f*^{ES}_{σ,δ}(*x*,*y*,*t*) by . In terms of numerical fluxes, the evolution from the *nth* time level to the (*n*+1)th time level is expressed by
4.5where the numerical fluxes are defined by
4.6The split flux vectors in equation (4.6) are defined by
4.7Here, and . The numerical flux for a TVD scheme can be given by
4.8with *ϕ*(*θ*_{x}_{i,j}) chosen to be the van Leer limiter [28] given by
4.9

The expression of can be derived similarly.

### (c) Determination of *z* and *λ*_{αβ}

After acquiring the macroscopic moments, we need to update the values of *z* and *λ*_{αβ} in order to obtain the *f*^{ES} distribution function for the next time step. To solve *z*(*x*,*y*,*t*), which is the root of equation
4.10we can use numerical root finding methods. After the new value of *z* is determined, we can immediately find *λ*_{αβ} before advancing the solution to the next time step,
4.11We can then substitute the above quantities to the *f*^{ES} distribution function.

## 5. Computational results

In this section, we present the results of selected computational tests. Here, we test our algorithm to solve two-dimensional Riemann problems of semiclassical gas dynamics.

### (a) Two-dimensional Riemann problems

We computationally study the two-dimensional Riemann problems for ideal quantum gas dynamics using the present direct solver for the semiclassical Boltzmann–ES–BGK equation in the *Kn*=0 limit, i.e. *f*=*f*^{ES}. Following earlier work [29–31], we selected several configurations to be tested among the 19 configurations classified. We also follow the notations used there, namely, *J* denotes the slip line, *R* the rarefaction wave, *S* the shocks, the forward arrow indicates forward waves and the backward arrow indicates backward waves. In the Riemann problems of two-dimensional gas dynamics, the initial state is in equilibrium and kept constant in each of the four quadrants. See the illustration in figure 1. The set-up is made so that only one elementary wave, a one-dimensional shock, ** S**, a one-dimensional rarefaction wave,

**, or a two-dimensional slip line (or contact discontinuity),**

*R**J*, appears at each interface. The computational domain is (

*x*,

*y*)∈[0,1]×[0,1], and is divided into equal intervals with

*Δx*=

*Δy*. After some mesh refinement tests, we found that 200×200 grids with

*Δx*=

*Δy*=0.005 and the Gauss–Hermite quadrature rule with 20 abscissas in both

*υ*

_{x}and

*υ*

_{y}directions will yield good resolution of the major flow features, and all the results below are obtained with these meshes and discrete velocity system. Specifically, we consider two configurations, namely

*Configuration 5*and

*Configuration 17*among the 19 genuinely different configurations classified in [29–31]. To satisfy these constraints, the set initial data must satisfy the Rankine–Hugoniot relations, i.e. for given a pair of left and right states (denoted by the indices

*l*and

*r*), we define the two configurations as follows.

*Configuration 5*. This is the four slip lines case with motion in opposite directions. We have

Here, *E*_{ij} denotes an elementary wave *E* between the *i*th and *j*th quadrant, and *E* can be either *J*, ** S** or

**. The initial constant states at the four quadrants are specified as After some time**

*R**t*>0, the slip lines

*J*

_{21}and

*J*

_{32}will meet the sonic circle of the constant state in the second quadrant and continue as almost straight lines so that a quarter of the sonic circle lies in between. Identical symmetric flow patterns for the slip lines

*J*

_{34}and

*J*

_{41}occur in the fourth quadrant. Inside the subsonic area, the slip lines bend and end in spirals. Centred about the point (1/2(

*u*

_{x1}+

*u*

_{x3}),1/2(

*u*

_{y1}+

*u*

_{y2})), there is an oval part of the subsonic area bounded by shocks. These shocks form simple Mach reflections near the slip lines.

*Configuration 17*. In this configuration, we have two slips, one shock and one rarefaction,
Here, the functions *Φ*_{41} and *Ψ*_{32} are the same as those defined in [30]. The initial conditions at the four quadrants are specified, respectively, by
It is noted that with the macroscopic quantities (*z*,** u**,

*T*) specified for each quadrant at time

*t*=0, one can then define the equilibrium distribution

*f*

^{ES}=

*f*

^{eq}for each quadrant at

*t*=0 as the initial condition. Thus, the two-dimensional Riemann problem is completely specified on the distribution function and one starts marching in time the solution of the distribution function. This is rather different from the traditional two-dimensional Riemann problems [29,30].

In this configuration, the shock, slip line and rarefaction fan are conjectured to be propagating from interfaces of the four quadrants. A shock shall be seen propagating through quadrants 3 and 2, a rarefaction fan through quadrants 4 and 1 and slip discontinuities through quadrants 2 and 1 and through quadrants 4 and 1.

We first report the results for *Configuration 5* for three statistics with parameter *b*=0.5. In figure 2, the number density, pressure and fugacity contours are shown for the three statistics at about the time instant similar to that reported in [29,30]. From the density and pressure contours, one can identify the interesting and complicated wave patterns. At the right top corner, the two slip lines *J*_{21} and *J*_{32} meet the sonic circle of the constant state in the second quadrant and continue as almost straight lines so that a quarter of the sonic circle lies in between. Our results are consistent with calculations in [29,30] and with the theoretical conjecture in [31] for the same configuration and with the same initial data. There are detectable differences among the three statistics, although the overall wave patterns are similar. Quantitatively, comparing the same contours for the three statistics, the numerical values of number density and pressure for the Bose–Einstein statistics are the largest among the three statistics, and the Maxwell–Boltzmann statistics always lie between the other two as dictated by the *θ* values.

Next, we report the results for *Configuration 17*. The same contours of flow properties as in figure 2 are shown in figure 3. After some time, as shown in the figure, the two slip lines *J*_{21} and *J*_{34} part the wave patterns into a left and right part and connect with a vortex inside the subsonic region. For *Configuration 17*, again, the wave patterns shown in figure 3 are consistent with calculations in [29,30] and with the theoretical conjecture in [31] for the same configuration and with the same initial data. Notable differences among the three statistics can be observed, although the overall wave patterns are similar. Quantitatively, the numerical values of number density and pressure contours for the Bose–Einstein statistics are the largest among the three statistics, and the flow properties obtained by using Maxwell–Boltzmann statistics are always greater than those obtained with Fermi–Dirac statistics and lower than those with Bose–Einstein statistics, which is as it should be.

### (b) Mesh refinement tests

Here, we investigate the consistency and convergence of our numerical implementation by increasing the number of grids in the computational domain and therefore reducing their size. We expect a decrease in the truncation error, and therefore observe a convergence pattern to the solution of the chosen configuration under test. We report in figure 4 our findings in *Configuration 5* and *Configuration 17* using a grid of 400×400, i.e. four times the total of cells used in our initial test grid. We observe in figure 4*a*,*b* that slip line structures exhibit much sharper profiles with a refined mesh. This is also true in the vicinity of the rarefaction fans and contact discontinuities by observing the structures in figure 4*c*,*d*. Both results in figure 4 are computed with a Courant–Friedrichs–Lewy condition equal to 0.6 and corresponds to an output time of *t*=0.25. Most of our results here presented are obtained using a grid system of 200×200, and they exhibit slightly smeared profiles compared to those using the grid system of 400×400. The trend of grid convergence of the solutions can be clearly observed.

### (c) Effect of parameter *b* in the ellipsoidal-statistical equilibrium distribution on flow properties

Here, we investigate quantitatively and qualitatively the effects of the dimensionless parameter *b* in the ES distribution model, *f*^{ES}, on the wave flow patterns and properties. From qualitative observations on different configurations, it is shown that parameter *b* influences the results by increasing or decreasing the sharpness of the profiles of every configuration. Quantitatively, we note that it is not deterministic how this parameter may affect the overall magnitude of flow properties. We report this observation by using the number density contours for *Configuration 5* evaluated with both quantum statistics in figure 5. We calculate three cases for the discrete values of *b*={−0.5,0,0.75}; and the number density wave patterns for Bose–Einstein and Fermi–Dirac statistics with *b*=0.5 are shown in figure 2*a*,*g*, respectively. Also it is noted that by setting *b*=0, we recover the traditional isotropic semiclassical equilibrium distribution, i.e. *f*^{ES}=*f*^{eq}. By observing figure 2*a*, one can find differences in the density number contours and their numerical values; however, they are still consistent with calculations in [29,30].

### (d) Observations on *W*_{αβ} and pressure tensor components

In figures 6 and 7, the contours of *W*_{αβ} and *P*_{αβ} tensor components of the ES model are depicted for *Configuration 5* with *b*=0.5 using ideal gases following Bose–Einstein and Fermi–Dirac statistics, respectively. It is readily noted that contours of diagonal tensor components, *W*_{xx} and *W*_{yy}, are very similar, but a closer look reveals that both differ slightly numerically. By contrast, the *W*_{xy} components, as depicted in figures 6*b* and 7*b* are numerically small compared with the former ones. However, this component contribution cannot be neglected in the present anisotropic model.

By looking at pressure tensor components in figures 6 and 7, it is evident that *W*_{αβ} and *P*_{αβ} tensor quantities are closely related. This was expected by recalling relations (3.4); moreover, it is also easy to note that the contribution of the main diagonal tensor components will generally be greater than the rest of the tensor components.

Comparing the pressure components, as before, we observe that the main diagonal components, *P*_{xx} and *P*_{yy}, exhibit similar flow patterns, but again, slight numerical differences can be identified. Note that the main diagonal components will be equal when the dimensionless parameter *b*=0, i.e. *f*^{eq} is used. Based on the quantitative observations on several tested configurations, it is safe to conjecture that due to the slight differences in *P*_{xx} and *P*_{yy} and the contribution of *P*_{xy} components, as triggered by the dimensionless parameter *b* of the ES distribution function, the anisotropic ES model contributes to relax the fixed Prandtl number condition of the traditional equilibrium distribution model, *f*^{eq}. Note the similar patterns between pressure contours in figures 6 and 7, and the results also show that numerical pressure values are greater for Bose–Einstein statistics, as expected. Finally, although the kinds of gases used in these examples are Bose and Fermi gases, the results are comparable to the earlier works [29,30] (classical gas) in terms of the expected flow patterns.

## 6. Concluding remarks

The semiclassical ES equilibrium distribution of Wu *et al.* [20] has been derived using the maximum entropy principle and conserves the mass, momentum and energy, but differs from the standard Bose–Einstein or Fermi–Dirac distribution. This ES distribution is anisotropic, thus it can possess additional high-order moments, and therefore, its gas dynamical features are not well known and are worth exploring. In this work, the unsteady quantum gas dynamical flow features as dictated by the semiclassical ES equilibrium distribution are numerically studied. The computational method treats the governing equation in phase space and employs the discrete ordinate method and high-resolution shock-capturing schemes. Specifically, we describe the solution method in detail for the equation in two space dimensions. A decoding procedure is devised for the semiclassical ES distribution which is different from that for standard Bose–Einstein or Fermi–Dirac distributions. Computations of two-dimensional Riemann problems for gas flows of arbitrary particle statistics are presented, and all our results are in good agreement with the calculations in [29,30] and consistent with the theoretical conjecture in [31]. The effect of the parameter *b* of the ES model is studied. Results showed how this parameter interacts with the pressure tensor, which can influence the Prandtl number in the non-equilibrium flows, and therefore the transport properties of classical and quantum gases when the full semiclassical Boltzmann–ES–BGK equation is considered. These computational examples serve the purposes of exploring the nonlinear manifestation of shock waves, contact lines and rarefaction waves, and testing the robustness of the present method. All the expected flow profiles comprising shock and rarefaction waves and contact discontinuities of semiclassical ideal gases and their nonlinear interactions can be observed with considerably good detail and are consistent with available results. The present work emphasizes building the unified and parallel framework for treating the semiclassical gas dynamics of three statistics. Extensions to three space dimensions and more complex geometries in general coordinates are straightforward, and can be carried out following the same framework presented here.

## Funding statement

This work was done under the auspices of the National Science Council, Taiwan, through grant nos. NSC-99-2221-E002-084MY3, NSC-100-2221-E002-106MY3 and CASTS 10R80909-4. The work for Z.L. was supported by the National Nature Science Foundation of China under grant no. 91016027.

- Received June 24, 2013.
- Accepted October 2, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.