A Numerical Scheme for the Quantum Boltzmann Equation Efficient in the Fluid Regime
aa r X i v : . [ m a t h . NA ] S e p A NUMERICAL SCHEME FOR THE QUANTUM BOLTZMANN EQUATIONEFFICIENT IN THE FLUID REGIME
FRANCIS FILBET , JINGWEI HU AND SHI JIN
Abstract.
Numerically solving the Boltzmann kinetic equations with the small Knudsen numberis challenging due to the stiff nonlinear collision term. A class of asymptotic preserving schemeswas introduced in [6] to handle this kind of problems. The idea is to penalize the stiff collisionterm by a BGK type operator. This method, however, encounters its own difficulty when appliedto the quantum Boltzmann equation. To define the quantum Maxwellian (Bose-Einstein or Fermi-Dirac distribution) at each time step and every mesh point, one has to invert a nonlinear equationthat connects the macroscopic quantity fugacity with density and internal energy. Setting a goodinitial guess for the iterative method is troublesome in most cases because of the complexity of thequantum functions (Bose-Einstein or Fermi-Dirac function). In this paper, we propose to penalizethe quantum collision term by a ‘classical’ BGK operator instead of the quantum one. This is basedon the observation that the classical Maxwellian, with the temperature replaced by the internalenergy, has the same first five moments as the quantum Maxwellian. The scheme so designed avoidsthe aforementioned difficulty, and one can show that the density distribution is still driven towardthe quantum equilibrium. Numerical results are present to illustrate the efficiency of the new schemein both the hydrodynamic and kinetic regimes. We also develop a spectral method for the quantumcollision operator. Introduction
The quantum Boltzmann equation (QBE), also known as the Uehling-Uhlenbeck equation, describesthe behaviors of a dilute quantum gas. It was first formulated by Nordheim [13] and Uehling andUhlenbeck [16] from the classical Boltzmann equation by heuristic arguments. Here we mainly considertwo kinds of quantum gases: the Bose gas and the Fermi gas. The Bose gas is composed of Bosons,which have an integer value of spin, and obey the Bose-Einstein statistics. The Fermi gas is composedof Fermions, which have half-integer spins and obey the Fermi-Dirac statistics.Let f ( t, x, v ) ≥ t , position x andparticle velocity v , then the quantum Boltzmann equation reads:(1.1) ∂f∂t + v · ∇ x f = 1 ǫ Q q ( f ) , x ∈ Ω ⊂ R d x , v ∈ R d v . Here ǫ is the Knudsen number which measures the degree of rarefaction of a gas. It is the ratiobetween the mean free path and the typical length scale. The quantum collision operator Q q is(1.2) Q q ( f )( v ) = Z R dv Z S dv − B ( v − v ∗ , ω ) [ f ′ f ′∗ (1 ± θ f )(1 ± θ f ∗ ) − f f ∗ (1 ± θ f ′ )(1 ± θ f ′∗ )] dωdv ∗ where θ = ~ d v , ~ is the rescaled Planck constant. In this paper, the upper sign will always correspondto the Bose gas while the lower sign to the Fermi gas. For the Fermi gas, we also need f ≤ θ bythe Pauli exclusion principle. f , f ∗ , f ′ and f ′∗ are the shorthand notations for f ( t, x, v ), f ( t, x, v ∗ ), f ( t, x, v ′ ) and f ( t, x, v ′∗ ) respectively. ( v, v ∗ ) and ( v ′ , v ′∗ ) are the velocities before and after collision.They are related by the following parametrization: v ′ = v + v ∗ | v − v ∗ | ω,v ′∗ = v + v ∗ − | v − v ∗ | ω, (1.3) This work was partially supported by NSF grant DMS-0608720 and NSF FRG grant DMS-0757285. FF was supportedby the ERC Starting Grant Project NuSiKiMo. SJ was also supported by a Van Vleck Distinguished Research Prizeand a Vilas Associate Award from the University of Wisconsin-Madison. where ω is the unit vector along v ′ − v ′∗ . The collision kernel B is a nonnegative function that onlydepends on | v − v ∗ | and cos θ ( θ is the angle between ω and v − v ∗ ). In the Variable Hard Sphere(VHS) model, it is given by(1.4) B ( v − v ∗ , ω ) = C γ | v − v ∗ | γ where C γ is a positive constant. γ = 0 corresponds to the Maxwellian molecules, γ = 1 is the hardsphere model.When the Knudsen number ǫ is small, the right hand side of equation (1.1) becomes stiff andexplicit schemes are subject to severe stability constraints. Implicit schemes allow larger time step,but new difficulty arises in seeking the numerical solution of a fully nonlinear problem at each timestep. Ideally, one wants an implicit scheme allowing large time steps and can be inverted easily. In[6], for the classical Boltzmann equation, Filbet and Jin proposed to penalize the nonlinear collisionoperator Q c by a BGK operator:(1.5) Q c = [ Q c − λ ( M c − f )] + λ [ M c − f ]where λ is a constant that depends on the spectral radius of the linearized collision operator of Q c around the local (classical) Maxwellian M c . Now the term in the first bracket of the right hand sideof (1.5) is less stiff than the second one and can be treated explicitly. The term in the second bracketwill be discretized implicitly. Using the conservation property of the BGK operator, this implicitterm can actually be solved explicitly. Thus they arrive at a scheme which is uniformly stable in ǫ ,with an implicit source term that can be inverted explicitly. Furthermore, under certain conditions,one could show that this type of schemes has the following property: the distance between f andthe Maxwellian will be O ( ǫ ) after several time steps, no matter what the initial condition is. Thisguarantees the capturing of the fluid dynamic limit even if the time step is larger than the mean freetime.Back to the quantum Boltzmann equation (1.1), a natural way to generalize the above idea is topenalize Q q with the quantum BGK operator M q − f . This means we have to invert a nonlinearalgebraic system that contains the unknown quantum Maxwellian M q (Bose-Einstein or Fermi-Diracdistribution) for every time step. As mentioned in [7], this is not a trivial task compared to theclassical case. Specifically, one has to invert a nonlinear 2 by 2 system (can be reduced to one nonlinearequation) to obtain the macroscopic quantities, temperature and fugacity. Due to the complexity ofthe quantum distribution functions (Bose-Einstein or Fermi-Dirac function), it is really a delicate issueto set a good initial guess for an iterative method such as the Newton method to converge.In this work we propose a new scheme for the quantum Boltzmann equation. Our idea is basedon the observation that the classical Maxwellian, with the temperature replaced by the (quantum)internal energy, has the same first five moments as the quantum Maxwellian. This observation wasused in [7] to derive a ‘classical’ kinetic scheme for the quantum hydrodynamical equations. Therefore,we just penalize the quantum collision operator Q q by a ‘classical’ BGK operator, thus avoid theaforementioned difficulty. At the same time, we have to sacrifice a little bit on the asymptotic property.Later we will prove that for the quantum BGK equation, the so obtained f satisfies:(1.6) f n − M nq = O (∆ t ) for some n > N, any initial data f , i.e. f will converge to the quantum Maxwellian beyond the initial layer with an error of O (∆ t ).Another numerical issue is how to evaluate the quantum collision operator Q q . In fact (1.2) can besimplified as(1.7) Q q ( f )( v ) = Z R dv Z S dv − B ( v − v ∗ , ω ) [ f ′ f ′∗ (1 ± θ f ± θ f ∗ ) − f f ∗ (1 ± θ f ′ ± θ f ′∗ )] dωdv ∗ so Q q is indeed a cubic operator. Almost all the existing fast algorithms are designed for the classicalBoltzmann operator based on its quadratic structure. Here we will give a spectral method for theapproximation of Q q . As far as we know, this is the first time to compute the full quantum Boltzmanncollision operator with the spectral accuracy.The rest of the paper is organized as follows. In the next section, we give a brief introduction to thequantum Boltzmann equation: the basic properties, the quantum Maxwellians and the hydrodynamiclimits. In section 3, we present the details of computing the quantum collision operator by the spectralmethod as well as the numerical accuracy. Our new scheme to capture the hydrodynamic regime isgiven in section 4. In section 5, the proposed schemes are tested on the 1-D shock tube problem of NUMERICAL SCHEME FOR THE QUANTUM BOLTZMANN EQUATION EFFICIENT IN THE FLUID REGIME 3 the quantum gas for different Knudsen number ǫ ranging from fluid regime to kinetic regime. Thebehaviors of the Bose gas and the Fermi gas in both the classical regime and quantum regime areincluded. Finally some concluding remarks are given in section 6.2. The Quantum Boltzmann Equation and its Hydrodynamic Limits
In this section we review some basic facts about the quantum Boltzmann equation (1.1). • At the formal level, Q q conserves mass, momentum and energy.(2.1) Z R dv Q q ( f ) dv = Z R dv Q q ( f ) vdv = Z R dv Q q ( f ) | v | dv = 0 . • If f is a solution of QBE (1.1), the following local conservation laws hold: ∂∂t Z R dv f dv + ∇ x · Z R dv vf dv = 0 ,∂∂t Z R dv vf dv + ∇ x · Z R dv v ⊗ vf dv = 0 ,∂∂t Z R dv | v | f dv + ∇ x · Z R dv v | v | f dv = 0 . (2.2) Define the macroscopic quantities: density ρ , macroscopic velocity u , specific internal energy e as ρ = Z R dv f dv, ρ u = Z R dv vf dv, ρe = Z R dv | v − u | f dv (2.3) and stress tensor P and heat flux q P = Z R dv ( v − u ) ⊗ ( v − u ) f dv, q = Z R dv
12 ( v − u ) | v − u | f dv, (2.4) the above system can then be recast as ∂ρ∂t + ∇ x · ( ρu ) = 0 ,∂ ( ρu ) ∂t + ∇ x · ( ρu ⊗ u + P ) = 0 ,∂∂t (cid:18) ρe + 12 ρu (cid:19) + ∇ x · (cid:18)(cid:18) ρe + 12 ρu (cid:19) u + P u + q (cid:19) = 0 . (2.5) • Q q satisfies Boltzmann’s H-Theorem,(2.6) Z R dv ln (cid:18) f ± θ f (cid:19) Q q ( f ) dv ≤ , moreover,(2.7) Z R dv ln (cid:18) f ± θ f (cid:19) Q q ( f ) dv = 0 ⇐⇒ Q q ( f ) = 0 ⇐⇒ f = M q , where M q is the quantum Maxwellian given by(2.8) M q = 1 θ z − e ( v − u )22 T ∓ , where z is the fugacity, T is the temperature (see [7] for more details about the derivation of M q ). This is the well-known Bose-Einstein (‘-’) and Fermi-Dirac (‘+’) distributions. FRANCIS FILBET , JINGWEI HU AND SHI JIN
Hydrodynamic Limits.
Substituting M q into (2.3) (2.4), the system (2.5) can be closed,yielding the quantum Euler equations: ∂ρ∂t + ∇ x · ( ρu ) = 0 ,∂ ( ρu ) ∂t + ∇ x · (cid:18) ρu ⊗ u + 2 d v ρeI (cid:19) = 0 ,∂∂t (cid:18) ρe + 12 ρu (cid:19) + ∇ x · (cid:18)(cid:18) d v + 2 d v ρe + 12 ρu (cid:19) u (cid:19) = 0 . (2.9)With the macroscopic variables ρ , u and e , they are exactly the same as the classical Euler equations.However, the intrinsic constitutive relation is quite different. ρ and e are connected with T and z (used in the definition of M q (2.8)) by a nonlinear 2 by 2 system: ρ = (2 πT ) dv θ Q dv ( z ) ,e = d v T Q dv +22 ( z ) Q dv ( z ) , (2.10)where Q ν ( z ) denotes the Bose-Einstein function G ν ( z ) and the Fermi-Dirac function F ν ( z ) respec-tively, G ν ( z ) = 1Γ( ν ) Z ∞ x ν − z − e x − dx, < z < , ν > z = 1 , ν > , (2.11) F ν ( z ) = 1Γ( ν ) Z ∞ x ν − z − e x + 1 dx, < z < ∞ , ν > , (2.12)and Γ( ν ) = R ∞ x ν − e − x dx is the Gamma function.The physical range of interest for a Bose gas is 0 < z ≤
1, where z = 1 corresponds to the degeneratecase (the onset of Bose-Einstein condensation). For the Fermi gas we don’t have such a restrictionand the degenerate case is reached when z is very large. For small z (0 < z < z , G ν ( z ) = ∞ X n =1 z n n ν = z + z ν + z ν + . . . , (2.13) F ν ( z ) = ∞ X n =1 ( − n +1 z n n ν = z − z ν + z ν − . . . . (2.14)Thus, for z ≪
1, both functions behave like z itself and one recovers the classical limit.On the other hand, the first equation of (2.10) can be written as(2.15) Q dv ( z ) = ρ (2 πT ) dv θ where ρ (2 πT ) dv is just the coefficient of the classical Maxwellian, which should be an O (1) quantity.Now if θ →
0, then Q dv ( z ) →
0, which means z ≪ Q ν . Thisis consistent with the fact that one gets the classical Boltzmann equation in QBE (1.1) by letting θ → ρ and e [1]. NUMERICAL SCHEME FOR THE QUANTUM BOLTZMANN EQUATION EFFICIENT IN THE FLUID REGIME 5 Computing the Quantum Collision Operator Q q In this section, we discuss the approximation of the quantum collision operator Q q . The methodwe use is an extension of the spectral method introduced in [12, 5] for the classical collision operator.We first write (1.2) as(3.1) Q q = Q c ± θ ( Q + Q − Q − Q ) , where(3.2) Q c ( f )( v ) = Z R dv Z S dv − B ( v − v ∗ , ω )[ f ′ f ′∗ − f f ∗ ] dωdv is the classical collision operator. The cubic terms Q – Q are Q ( f )( v ) = Z R dv Z S dv − B ( v − v ∗ , ω ) f ′ f ′∗ f ∗ dωdv, Q ( f )( v ) = Z R dv Z S dv − B ( v − v ∗ , ω ) f ′ f ′∗ f dωdv, Q ( f )( v ) = Z R dv Z S dv − B ( v − v ∗ , ω ) f f ∗ f ′ dωdv, Q ( f )( v ) = Z R dv Z S dv − B ( v − v ∗ , ω ) f f ∗ f ′∗ dωdv. (3.3)In order to perform the Fourier transform, we periodize the function f on the domain D L =[ − L, L ] d v ( L is chosen such that L ≥ √ R , R is the truncation of the collision integral whichsatisfies R = 2 S , where B (0 , S ) is an approximation of the support of f [14]). Using the Carlemanrepresentation [2], one can rewrite the operators as (for simplicity we only consider the 2-D Maxwellianmolecules), Q c ( f )( v ) = Z B R Z B R δ ( x · y )[ f ( v + x ) f ( v + y ) − f ( v + x + y ) f ( v )] dxdy (3.4)and Q ( f )( v ) = Z B R Z B R δ ( x · y ) f ( v + x ) f ( v + y ) f ( v + x + y ) dxdy, Q ( f )( v ) = Z B R Z B R δ ( x · y ) f ( v + x ) f ( v + y ) f ( v ) dxdy, Q ( f )( v ) = Z B R Z B R δ ( x · y ) f ( v + x ) f ( v + x + y ) f ( v ) dxdy, Q ( f )( v ) = Z B R Z B R δ ( x · y ) f ( v + y ) f ( v + x + y ) f ( v ) dxdy. (3.5)Now we approximate f by a truncated Fourier series,(3.6) f ( v ) ≈ N − X k = − N ˆ f k e i πL k · v , ˆ f k = 1(2 L ) d v Z D L f ( v ) e − i πL k · v dv. Plugging it into (3.4) (3.5), one can get the k -th mode of ˆ Q q . The classical part is the same as thosein the previous method [12]. We will mainly focus on the cubic terms.Define the kernel modes(3.7) β ( l, m ) = Z B R Z B R δ ( x · y ) e i πL l · x e i πL m · y dxdy. FRANCIS FILBET , JINGWEI HU AND SHI JIN
Following [12], β ( l, m ) can be decomposed as(3.8) β ( l, m ) = πM M − X p =0 α p ( l ) α ′ p ( m )with(3.9) α p ( l ) = φ ( l · (cos θ p , sin θ p )) , α ′ p ( m ) = φ ( m · ( − sin θ p , cos θ p )) , where φ ( s ) = Lπs sin( πL Rs ), M is the number of equally spaced points in [0 , π ] and θ p = π pM . Then • The k -th coefficient of ˆ Q is N − X l,m,n = − N l + m + n = k β ( l + n, m + n ) ˆ f l ˆ f m ˆ f n = πM M − X p =0 N − X n = − N N − X l,m = − N l + m = k − n α p ( l + n ) α ′ p ( m + n ) ˆ f l ˆ f m ˆ f n = πM M − X p =0 N − X n = − N ˆ g k − n ( n ) ˆ f n . (3.10) Terms inside the bracket is a convolution (defined as ˆ g k − n ( n )), which can be computed bythe Fast Fourier Transform (FFT). However, the outside structure is not a convolution, sinceˆ g k − n ( n ) itself depends on n . So we compute this part directly. • The k -th coefficient of ˆ Q is N − X l,m,n = − N l + m + n = k β ( l, m ) ˆ f l ˆ f m ˆ f n = πM M − X p =0 N − X n = − N N − X l,m = − N l + m = k − n α p ( l ) α ′ p ( m ) ˆ f l ˆ f m ˆ f n . (3.11) In this case, both inside and outside are convolutions. The FFT can be implemented easily. • The k -th coefficient of ˆ Q is N − X l,m,n = − N l + m + n = k β ( l + m, m ) ˆ f l ˆ f m ˆ f n = πM M − X p =0 N − X n = − N α p ( l + m ) N − X l,m = − N l + m = k − n α ′ p ( m ) ˆ f l ˆ f m ˆ f n . (3.12) Factoring out α p ( l + m ), both inside and outside are convolutions again. • The k -th coefficient of ˆ Q is N − X l,m,n = − N l + m + n = k β ( m, l + m ) ˆ f l ˆ f m ˆ f n = πM M − X p =0 N − X n = − N α ′ p ( l + m ) N − X l,m = − N l + m = k − n α p ( m ) ˆ f l ˆ f m ˆ f n . (3.13) This term can be evaluated similarly as ˆ Q . Remark 3.1.
The computational cost of this quantum solver is O ( M N log N ) , which mainly comesfrom computing Q . This cost is higher than O ( M N ) of the discrete velocity model. But takinginto account the high accuracy and small value of log N ( N is not very big in the real simulation),our method is still more attractive than the quadrature method. The fast algorithm for the quantumcollision operator remains an open problem. Numerical Accuracy.
To illustrate the accuracy of the above method, we test it on a steadystate, namely, we compute Q q ( M q ) and check its max norm. In all the numerical simulations, theparticles are assumed to be the 2-D Maxwellian molecules.Let ρ = 1, T = 1, from (2.10) one can adjust θ to get z that lies in different physical regimes.When θ = 0 .
01 ( ~ = 0 . z Bose = 0 . z Fermi = 0 . θ , say θ = 9 ( ~ = 3), z Bose = 0 . z Fermi = 3 . NUMERICAL SCHEME FOR THE QUANTUM BOLTZMANN EQUATION EFFICIENT IN THE FLUID REGIME 7
Figure 1.
The Maxwellians at ρ = 1, T = 1, θ = 0 .
01. Left: Bose gas; Center:classical gas; Right: Fermi gas.
Figure 2.
The Maxwellians at ρ = 1, T = 1, θ = 9. Left: Bose gas; Center:classical gas (same as in Fig.1); Right: Fermi gas.In Table 1, we list the values of k Q c ( M c ) k L ∞ and k Q q ( M q ) k L ∞ computed on different meshesN=16, 32, 64 (number of points in v direction), M=4 (number of points in angular direction θ p ; itis not necessary to put too many points since M won’t effect the spectral accuracy, see [12]). Thecomputational domain is [ − , × [ − ,
8] ( L = 8).These results confirm the spectral accuracy of the method, although the accuracy in the quantumregime is not as good as that in the classical regime. This is because the regularity of the quantumMaxwellians becomes worse when θ is increasing, or strictly speaking, the mesh size ∆ v is not smallenough to capture the shape of the Maxwellians. To remedy this problem, one can add more gridpoints or more effectively, shorten the computational domain. For the Bose-Einstein distribution, wealso include the results computed on [ − , × [ − ,
6] in Table 1. One can clearly see the improvements.
FRANCIS FILBET , JINGWEI HU AND SHI JIN ×
16 32 ×
32 64 ×
64 convergence rateclassical gas 2.1746e-04 3.8063e-12 1.9095e-16 20.0253Bose gas θ = 0 .
01 2.1084e-04 2.5512e-10 1.9080e-16 20.0036 θ = 9 0.4891 0.0310 1.3496e-04 5.9117 θ = 9, L = 6 0.1815 0.0052 4.0278e-06 7.7298Fermi gas θ = 0 .
01 2.2397e-04 1.6485e-10 1.9152e-16 20.0445 θ = 9 8.9338e-04 2.0192e-06 1.5962e-10 11.2081 Table 1.
Comparison of the quantum collision solver on different Maxwellians ( L = 8unless specified).3.2. Relaxation to Equilibrium.
Let us consider the space homogeneous quantum Boltzmann equa-tion for the 2-D Maxwellian molecules. As already mentioned, this equation satisfies the entropy con-dition, and the equilibrium states are the entropy minimizers. Hence, we first consider the quantumBoltzmann equation for a Fermi gas with an initial datum 0 ≤ f ≤ θ and observe the relaxationto equilibrium of the distribution function. Then, we take a Bose gas for which the entropy is nowsublinear and fails to prevent concentration, which is consistent with the fact that condensation mayoccur in the long-time limit.Fermi gas. The initial data is chosen as the sum of two Maxwellian functions(3.14) f ( v ) = exp (cid:18) − | v − v | (cid:19) + exp (cid:18) − | v + v | (cid:19) ; v ∈ R , with v = (2 , T end = 0 .
5, which is very close to the stationarystate.In the spatially homogeneous setting, Pauli’s exclusion principle facilitates things because of theadditional L ∞ bound 0 ≤ f ( t ) ≤ θ . In this case, the convergence to equilibrium in a weak sensehas been shown by Lu [10]. Later Lu and Wennberg proved the strong L stability [9]. However,no constructive result in this direction has ever been obtained, neither has any entropy-dissipationinequality been established.In Fig.3 we report the time evolution of the entropy and the fourth and sixth order moments ofthe distribution with respect to the velocity variable. We indeed observe the convergence to a steadystate of the entropy and also of high order moments when t → ∞ . Figure 3.
Fermi gas. Time evolution of the entropy, fourth and sixth order moments.In Fig.4 we also report the time evolution of the level set of the distribution function f ( t, v x , v y )obtained with N = 64 modes at different times. Initially the level set of the initial data correspondsto two spheres in the velocity space. Then, the two distributions start to mix together until the NUMERICAL SCHEME FOR THE QUANTUM BOLTZMANN EQUATION EFFICIENT IN THE FLUID REGIME 9 stationary state is reached, represented by a single centered sphere. It is clear that the sphericalshapes of the level sets are described with great accuracy by the spectral method.
Figure 4.
Fermi gas. Time evolution of the distribution function f ( t, v x , v y ) with N = 64 modes at times t = 0 , . , .
04 and 0 . f ( v ) = 14 π T exp (cid:18) − | v − v | T (cid:19) + exp (cid:18) − | v + v | T (cid:19) ; v ∈ R , with v = (1 , /
2) and T = 1 / t → ∞ in Fig.5.In Fig.6 we report the time evolution of the level set of the distribution function f ( t, v x , v y ) obtainedwith N = 64 modes at different times and observe the trend to equilibrium.4. A Scheme Efficient in the Fluid Regime
So far we have only considered spatially homogeneous quantum Boltzmann equations, now whathappens for spatially inhomogeneous data? Due to the natural bound 0 ≤ f ( t ) ≤ θ , the Boltzmann-Fermi model seems to be well understood mathematically [17]. The situation is completely differentfor the Boltzmann-Bose model, since singular measures may occur [17]. Figure 5.
Bose gas. Time evolution of the entropy, fourth and sixth order moments.
Figure 6.
Bose gas. Time evolution of the distribution function f ( t, v x , v y ) with N = 64 modes at times t = 0 , . , .
04 and 0 . ∂f∂t + v · ∇ x f = 1 ǫ Q c ( f ) . The first-order scheme reads:(4.2) f n +1 − f n ∆ t + v · ∇ x f n = Q c ( f n ) − λ ( M nc − f n ) ǫ + λ ( M n +1 c − f n +1 ) ǫ , where λ is some appropriate approximation of |∇Q c | (can be made time dependent). To solve f n +1 explicitly, we need to compute M n +1 c first. Since the right hand side of (4.2) is conservative, it vanishes NUMERICAL SCHEME FOR THE QUANTUM BOLTZMANN EQUATION EFFICIENT IN THE FLUID REGIME11 when we take the moments (multiply by φ ( v ) = (1 , v, v ) T and integrate with respect to v ). Then(4.2) becomes(4.3) U n +1 − U n ∆ t + Z φ ( v ) v · ∇ x f n dv = 0 , where U = ( ρ, ρu, ρe + ρu ) T is the conserved quantities. Once we get U n +1 , M n +1 c is known. Now f n +1 in (4.2) is easy to obtain.When generalizing the above idea to the quantum Boltzmann equation (1.1), the natural idea is toreplace Q c and M c in (4.2) by Q q and M q respectively. However, as mentioned in section 2, one hasto invert the nonlinear system (2.10) to get z and T . Experiments show that the iterative methodsdo converge when the initial guess is close to the solution (analytically, this system has a solution [1]).But how to set a good initial guess for every spatial point and every time step is not an easy task,especially when ρ and e are not continuous.Here we propose to use a ‘classical’ BGK operator to penalize Q q . Specifically, we replace thetemperature T with the internal energy e in the classical Maxwellian using relation e = d v T (true forclassical monatomic gases) and get(4.4) M c = ρ (2 πT ) dv e − ( v − u )22 T = ρ (cid:18) d v πe (cid:19) dv e − dv e ( v − u ) . An important property of M c is that it has the same first five moments as M q .Now our new scheme for QBE (1.1) can be written as(4.5) f n +1 − f n ∆ t + v · ∇ x f n = Q q ( f n ) − λ ( M nc − f n ) ǫ + λ ( M n +1 c − f n +1 ) ǫ . Since the right hand side is still conservative, one computes M n +1 c the same as for (4.2).It is important to notice that z and T are not present at all in this new scheme, thus one does notneed to invert the 2 by 2 system (2.10) during the time evolution. If they are desired variables foroutput, one only needs to convert between ρ , e and z , T at the final output time.4.1. Asymptotic Property of the New Scheme.
In this subsection we show that the new scheme,when applied to the quantum BGK equation, has the property (1.6). Consider the following timediscretization:(4.6) f n +1 − f n ∆ t + v · ∇ x f n = ( M nq − f n ) − λ ( M nc − f n ) ǫ + λ ( M n +1 c − f n +1 ) ǫ . Some simple mathematical manipulation on (4.6) gives(4.7) f n +1 −M n +1 q = 1 + ( λ − ∆ tǫ λ ∆ tǫ ( f n −M nq ) − ∆ t λ ∆ tǫ v ·∇ x f n +( M nq −M n +1 q )+ λ ∆ tǫ λ ∆ tǫ ( M n +1 c −M nc ) . Assume all the functions are smooth. When λ > ,(4.8) | f n +1 − M n +1 q | ≤ α | f n − M nq | + O ( ǫ + ∆ t ) , where 0 < α = | λ − ∆ tǫ | / | λ ∆ tǫ | < ǫ and ∆ t . The O ( ǫ ) term comes from thesecond term of the right hand side of (4.7). The O (∆ t ) term is from the third and fourth terms. Then(4.9) | f n − M nq | ≤ α n | f − M q | + O ( ǫ + ∆ t ) . Since ∆ t is taken bigger than ǫ , this implies the property (1.6). It is interesting to point out that f approaches M q , not M c , with (4.6). Remark 4.1.
The first order (in-time) method can be extended to a second order by an Implicit-Explicit (IMEX) method (see also [6] ): (4.10) f ∗ − f n ∆ t/ v · ∇ x f n = Q q ( f n ) − λ ( M nc − f n ) ǫ + λ ( M ∗ c − f ∗ ) ǫ ,f n +1 − f n ∆ t + v · ∇ x f ∗ = Q q ( f ∗ ) − λ ( M ∗ c − f ∗ ) ǫ + λ ( M nc − f n ) + λ ( M n +1 c − f n +1 )2 ǫ . This scheme can be shown to have the same property (1.6) on the quantum BGK equation. Numerical Examples
In this section, we present some numerical results of our new scheme (4.5) (a second order finitevolume method with slope limiters [8] is applied to the transport part) on the 1-D shock tube problem.The initial condition is(5.1) (cid:26) ( ρ l , u l , T l ) = (1 , ,
1) if 0 ≤ x ≤ . , ( ρ r , u r , T r ) = (0 . , , .
25) if 0 . < x ≤ . The particles are again assumed to be the 2-D Maxwellian molecules and we adjust θ to getdifferent initial data for both the Bose gas and the Fermi gas.In all the regimes, besides the directly computed macroscopic quantities, we will show the fugacity z and temperature T as well. They are computed as follows. First, (2.10) ( d v = 2) leads to(5.2) Q ( z ) Q ( z ) = θ π ρe . We treat the left hand side of (5.2) as one function of z , and invert it by the secant method. Once z is obtained, T can be computed easily using for example the first equation of (2.10). To evaluate thequantum function Q ν ( z ), the expansion (2.13) is used for the Bose-Einstein function. The Fermi-Diracfunction is computed by a direct numerical integration. The approach adopted here is taken from [15](Chapter 6.10).When approximating the collision operator Q q , we always take M = 4, N = 32 and L = 8, except L = 6 for the Bose gas in the quantum regime.5.1. Hydrodynamic Regime.
We compare the results of our new scheme (4.5) with the kineticscheme (KFVS scheme in [7]) for the quantum Euler equations (2.9). The time step ∆ t is chosen bythe CFL condition, independent of ǫ . Fig.7 shows the behaviors of a Bose gas when θ = 0 .
01. Fig.8shows the behaviors of a Bose gas when θ = 9. The solutions of a Fermi gas at θ = 0 .
01 are verysimilar to Fig.7, so we omit them here. Fig.9 shows the behaviors of a Fermi gas when θ = 9. Allthe results agree well in this regime, which exactly implies the scheme (4.5) is asymptotic preserving(when the Knudsen number ǫ goes to zero, the scheme becomes a fluid solver).5.2. Kinetic Regime.
We compare the results of our new scheme (4.5) with the explicit forwardEuler scheme. The time step ∆ t for the new scheme is still chosen by the CFL condition. When theKnudsen number ǫ is not very small, 10 − or 10 − , the above ∆ t is also enough for the explicit scheme.Fig.10 shows the behaviors of a Bose gas when θ = 0 .
01. Fig.11 shows the behaviors of a Bose gaswhen θ = 9. The solutions of a Fermi gas at θ = 0 .
01 are very similar to Fig.10, so we omit themhere. Fig.12 shows the behaviors of a Fermi gas when θ = 9. Again all the results agree well whichmeans the scheme (4.5) is also reliable in the kinetic regime. To avoid the boundary effect, all thesimulations in this subsection were carried out on a slightly larger spatial domain x ∈ [ − . , . Conclusion
A novel scheme was introduced for the quantum Boltzmann equation, starting from the scheme in[6]. The new idea here is to penalize the quantum collision operator by a ‘classical’ BGK operator soas to avoid the difficulty of inverting the nonlinear system ρ = ρ ( z, T ), e = e ( z, T ). The new scheme isuniformly stable in terms of the Knudsen number, and can capture the fluid (Euler) limit even if thesmall scale is not numerically resolved. We have also developed a spectral method for the quantumcollision operator, following its classical counterpart [12, 5].So far we have not considered the quantum gas in the extreme case. For example, the Bose gasbecomes degenerate when the fugacity z = 1. Many interesting phenomena happen in this regime.Our future work will focus on this aspect. Acknowledgments.
The second author would like to thank Mr. Bokai Yan for helpful discussionson the spectral method of the collision operator.
NUMERICAL SCHEME FOR THE QUANTUM BOLTZMANN EQUATION EFFICIENT IN THE FLUID REGIME13
Figure 7.
Bose gas. ǫ = 1 e − θ = 0 . z l = 0 . z r = 7 . e −
04. Density ρ , velocity u , fugacity z and temperature T at t = 0 .
2. ∆ t = 0 . x = 0 . ◦ : New scheme (4.5)for QBE (1.1). References [1] L. Arlotti and M. Lachowicz,
Euler and Navier-Stokes limits of the Uehling-Uhlenbeck quantum kinetic equations ,J. Math. Phys., (1997), pp. 3571–3588,[2] , T. Carleman, Sur la th´eorie de l’´equation int´egrodiff´erentielle de Boltzmann , Acta Math., , (1933) pp. 91–146,[3] , C. Cercignani, The Boltzmann Equation and Its Applications , Springer-Verlag, (1988)[4] , M. Escobedo and S. Mischler,
On a quantum Boltzmann equation for a gas of photons , J. Math. Pures Appl., , (2001), pp. 471-515,[5] F. Filbet and C. Mouhot and L. Pareschi, Solving the Boltzmann equation in NlogN , SIAM J. Sci. Comput., ,(2006), pp. 1029–1053[6] F. Filbet and S. Jin, A class of asymptotic-preserving schemes for kinetic equations and related problems with stiffsources , J. Comput. Phys., , (2010), 7625-7648[7] J. Hu and S. Jin,
On kinetic flux vector splitting schemes for quantum Euler equations preprint.[8] , R. J. LeVeque,
Numerical Methods for Conservation Laws , Birkh¨auser Verlag, (1992)[9] X. Lu and B. Wennberg,
On stability and strong convergence for the spatially homogeneous Boltzmann equationfor Fermi-Dirac Particles , Arch. Ration. Mech. Anal., , (2003), pp. 1-34[10] X. Lu,
On spatially homogeneous solutions of a modified Boltzmann equation for Fermi-Dirac particles , J. Stat.Phys., , (2001), pp. 353-388[11] X. Lu,
A modified Boltzmann equation for Bose-Einstein particles: isotropic solutions and long-time behavior , J.Stat. Phys., , (2000), pp. 1335-1394,[12] C. Mouhot and L. Pareschi, Fast algorithms for computing the Boltzmann collision operator , Math. Comput., ,(2006), pp. 1833–1852[13] L. W. Nordheim, On the kinetic method in the new statistics and its application in the electron theory of conduc-tivity , Proc. R. Soc. London, Ser. A, , (1928), pp. 689–698,[14] , L. Pareschi and G. Russo,
Numerical solution of the Boltzmann equation I. Spectrally accurate approximation ofthe collision operator , SIAM J. Numer. Anal., , (2000), pp. 1217–1245. Figure 8.
Bose gas. ǫ = 1 e − θ = 9, z l = 0 . z r = 0 . ρ , velocity u , fugacity z and temperature T at t = 0 .
2. ∆ t = 0 . x = 0 .
01. Solid line:KFVS scheme [7] for quantum Euler equations (2.9); ◦ : New scheme (4.5) for QBE(1.1). [15] W. H. Press and S. A. Teukolsky and W. T. Vetterling and B. P. Flannery, Numerical Recipes: The Art of ScientificComputing , Cambridge University Press, (2007)[16] E. A. Uehling and G. E. Uhlenbeck,
Transport phenomena in Einstein-Bose and Fermi-Dirac gases. I , Phys. Rev., , (1933) pp. 552–561,[17] C. Villani, A review of mathematical topics in collisional kinetic theory , North-Holland, S. Friedlander and D.Serre, (2002) pp. 71-305.
Francis FilbetUniversit´e de Lyon,Universit´e Lyon I, CNRSUMR 5208, Institut Camille Jordan43, Boulevard du 11 Novembre 191869622 Villeurbanne cedex, FRANCEe-mail: fi[email protected]
Jingwei HuDepartment of Mathematics,University of Wisconsin-Madison,480 Lincoln Drive, Madison,WI 53706, USAemail: [email protected]
NUMERICAL SCHEME FOR THE QUANTUM BOLTZMANN EQUATION EFFICIENT IN THE FLUID REGIME15
Figure 9.
Fermi gas. ǫ = 1 e − θ = 9, z l = 3 . z r = 1 . ρ ,velocity u , fugacity z and temperature T at t = 0 .
2. ∆ t = 0 . x = 0 .
01. Solidline: KFVS scheme [7] for quantum Euler equations (2.9); ◦ : New scheme (4.5) forQBE (1.1). Shi JinDepartment of Mathematics,University of Wisconsin-Madison,480 Lincoln Drive, Madison,WI 53706, USAemail: [email protected]
Figure 10.
Bose gas. ǫ = 1 e − θ = 0 . z l = 0 . z r = 7 . e −
04. Density ρ , velocity u , fugacity z and temperature T at t = 0 .
2. ∆ t = 0 . x = 0 . ◦ : New scheme (4.5) for QBE (1.1). NUMERICAL SCHEME FOR THE QUANTUM BOLTZMANN EQUATION EFFICIENT IN THE FLUID REGIME17
Figure 11.
Bose gas. ǫ = 1 e − θ = 9, z l = 0 . z r = 0 . ρ ,velocity u , fugacity z and temperature T at t = 0 .
2. ∆ t = 0 . x = 0 .
01. Solidline: Forward Euler scheme for QBE (1.1); ◦ : New scheme (4.5) for QBE (1.1). Figure 12.
Fermi gas. ǫ = 1 e − θ = 9, z l = 3 . z r = 1 . ρ ,velocity u , fugacity z and temperature T at t = 0 .
2. ∆ t = 0 . x = 0 .
01. Solidline: Forward Euler scheme for QBE (1.1); ◦◦