Influence of intermolecular potentials on rarefied gas flows: fast spectral solutions of the Boltzmann equation

The Boltzmann equation with an arbitrary intermolecular potential is solved by the fast spectral method. As examples, noble gases described by the Lennard-Jones potential are considered. The accuracy of the method is assessed by comparing both transport coe  cients with variational solutions and mass/heat ﬂow rates in Poiseuille/thermal transpiration ﬂows with results from the discrete velocity method. The fast spectral method is then applied to Fourier and Couette ﬂows between two parallel plates, and the inﬂuence of the intermolecular potential on various ﬂow properties is investigated. It is found that for gas ﬂows with the same rarefaction parameter, di ↵ erences in the heat ﬂux in Fourier ﬂow and the shear stress in Couette ﬂow are small. However, di ↵ erences in other quantities such as density, temperature, and velocity can be very


I. INTRODUCTION
When the ratio of the molecular mean free path to the characteristic flow length becomes significant, the Boltzmann equation (BE) is the best tool to investigate the rarefied gas dynamics 1 . The BE employs a one-particle velocity distribution function (VDF) to describe the state of a macroscopic volume of gas consisting of a large number of molecules, where the linear streaming operator models the molecular transport and the nonlinear Boltzmann collision operator (BCO) describes the binary molecular collisions.
The intermolecular potential is incorporated into the BCO through the di↵erential crosssection (DCS). As the DCSs for realistic potentials such as Lennard-Jones (LJ) or the potentials from ab initio calculations are very complicated, the simple hard-sphere (HS) model with a constant value of DCS is widely adopted 2 . However, the viscosity and heat conductivity of the HS model are proportional to the square root of the gas temperature, which does not agree with experimental data for common gases. To overcome this drawback, variable HS 2 , variable soft-sphere 3 , generalized HS 4 , and generalized soft-sphere 5 models have been proposed for the direct simulation Monte Carlo (DSMC) simulation of the BE. Also, the µ-DSMC method has been proposed in order to reproduce an arbitrary viscosity variation with temperature 6 . In our recent fast spectral approximation of the BCO, some special forms of the DCS were used to recover Sutherland's formula for viscosity, as well as the viscosity of the LJ potential 7,8 .
Note that all these DCSs were proposed in order to match the viscosity, and sometimes the mass di↵usion coe cient, with experimental data or theoretical values, but they ignore or simplify the detailed dependence of the DCS on the deflection angle and the relative collision energy that are characteristic of realistic potentials. For gas mixtures, the use of simplified DCSs becomes problematic, since it is di cult to recover the mass di↵usion and thermal di↵usion coe cients simultaneously for general intermolecular potentials. As the intermolecular potential can strongly influence certain phenomena in rarefied gases 9,10 , a numerical method to solve the BE with realistic potentials is urgently needed.
Implementation of the LJ potential in DSMC has been reported previously 11 , but was time-consuming as the deflection angle was calculated for every binary collision. Recently, the LJ potential and some ab initio potentials were successfully implemented into the DSMC method by pre-calculating the deflection and storing the results in a table [12][13][14] . Alternatively, an LJ polynomial approximation model was proposed to represent the deflection angle as a polynomial expansion in non-dimensional collision parameters 15,16 . Realistic intermolecular potentials have also been used in some deterministic numerical methods for solutions of the BE 10,17,18 . However, the discrete velocity method developed for the linearized BE 10,17 has a very high computational cost, which means it can only be applied to simple geometries, while the accuracy of the projection-interpolation method 18 is not clear when the VDF has steep variations or large discontinuities.
The aim of the present paper is to implement realistic intermolecular potentials in the fast spectral method (FSM), which is a promising numerical method for solving the BE deterministically 7,8,19,20 . By testing the proposed method, we also demonstrate the conditions in which the variable HS model can be adopted.

II. BOLTZMANN EQUATION
The state of a dilute monatomic gas is described by the VDF f (t, x, v) of the molecular The evolution of f is governed by the BE: where v@/@x is the streaming operator, while Q is the BCO defined by In the above equations, v, v ⇤ are the molecular velocities before the binary collision, while v 0 , v 0 ⇤ are the corresponding post-collision velocities. Conservation of momentum and energy yields v 0 = v + (|u|⌦ u)/2 and v 0 ⇤ = v ⇤ (|u|⌦ u)/2, where u = v v ⇤ is the relative pre-collision velocity and ⌦ is a vector in the unit sphere S 2 along the relative postcollision velocity v 0 v 0 ⇤ . The deflection angle ✓ between the pre-and post-collision relative velocities satisfies cos ✓ = ⌦ · u/|u|, with 0  ✓  ⇡. Finally, (✓, |v v ⇤ |) is the DCS.
For HS molecules with a molecular diameter d, it is d 2 /4, while for a general intermolecular potential the dependence of on |u| and ✓ is complicated and the numerical calculation of the DCS is necessary. Detailed information can be found in a recent publication 17 .
In this paper, we consider the following (6-12) LJ potential as an example: where ⇢ 0 is the intermolecular distance, ✏ is a potential depth, and d is the distance at which the potential is zero. As the interaction range of the LJ potential is ostensibly infinity, the total cross-section, i.e. the integral of the DCS with respect to the deflection angle, is infinity too. In practice, however, a finite cuto↵ either in the deflection angle 16,18,21,22 or in the radial potential 10,17,23 is introduced.

A. Normalizations
For practical calculations, it is convenient to introduce dimensionless variables. Here, the spatial location is normalized by the characteristic length`, temperature is normalized by T 0 , velocity is normalized by the most probable molecular speed v m = p 2k B T 0 /m, time is normalized by`/v m , molecular number density is normalized by n 0 , and the VDF is normalized by n 0 /v 3 m , where k B is the Boltzmann constant. Also, in the numerical evaluation of the DCS for the (6-12) LJ potential, the intermolecular distance ⇢ 0 is normalized by d.
Therefore, the BE becomes where d⌦dv ⇤ is the collision frequency and (✓, |u|v m ) is exactly the same as the DCS (✓, E) in Ref. 17 with the dimensionless relative collision The normalized density, flow velocity, and temperature are given by while the pressure tensor and heat flux, which are normalized by n 0 k B T 0 and n 0 k B T 0 v m , respectively, are given by where i, j = 1, 2, 3.

B. Linearized Boltzmann equation
When the system state deviates only slightly from equilibrium, the BE (4) can be linearized. We express the VDF around the global equilibrium state as where h is the deviation function satisfying |h/f eq | ⌧ 1. The evolution of h is governed by the linearized BE: @h @t with the linearized BCO where ⌫ eq = n 0 d 2`R S 2 |u| f eq (v 0 ⇤ )d⌦dv ⇤ is the equilibrium collision frequency.

OPERATOR
In this section, we focus on the numerical approximation of the BCO; the approximation of the linearized collision operator (9) can be performed according to the relation L(h) = . For simplicity, the coe cient n 0 d 2`i s regarded as 1.
We rewrite the BCO in the Carleman representation as 7 where is Dirac's delta function, and the DCS becomes The VDF is periodized on the truncated velocity domain D L = [ L, L] 3 . For simplicity, we adopt uniform discretization in velocity space: v k (j k ) = 2j k L/N k with k = 1, 2, 3, where j k 2 [ N k /2, N k /2 + 1, · · · , N k /2 1] and N k is the number of velocity grid points in the k-th velocity direction, although in the simulation of highly rarefied gas flows the velocity space would be better discretized non-uniformly 8,20 . The VDF is approximated by a truncated Fourier series: is the Fourier spectrum of the VDF, i is the imaginary unit, and ⇠ j = j⇡/L are the frequency components.
, and B S (a sphere of radius S centered on the origin) is the support of the VDF 19 . Our numerical experience suggests that is a good choice 7,8 . The truncated BCO is also expanded by the Fourier series, where its j-th Fourier mode is related to the Fourier coe cientf of the VDF as follows: with e, e 0 being the vectors in the unit sphere S 2 .
The integration with respect to ⇢ in Eq. (13) can be approximated by a numerical quadrature. Suppose ⇢ r and ! r (r = 1, 2, · · · , M r ) are the abscissas and weights of a quadrature for ⇢ where ✓ p (' q ) and ! p (! q ) are the p (q)-th point and weight in the Gauss-Legendre quadrature with J 0 being the zeroth-order Bessel function of first kind.
Thus, combining Eqs. (12) and (14), b Q can be calculated through FFT-based convolution, with a computational cost of O(M 2 M r N 3 log N ). Since M and M r can be far smaller than N , the FSM proposed here is faster than conventional spectral methods that have a cost of O(N 6 ). Note that in our previous works 7, 8 , a special form of DCS was proposed to approximate the DCS for the LJ potential, and since that special DCS can be decomposed into the form of 0 (|y|, |z|) = 0 1 (|y|) 0 2 (|z|), the integration with respect to ⇢ in Eq. (13) can be expressed analytically, resulting in a computational cost of O(M 2 N 3 log N ) for the BCO. Here, for general DCSs, one must approximate the integration with respect to ⇢ or ⇢ 0 by a numerical quadrature to get a computational gain; and this approximation extends the applicability of the FSM.
As with the FSM that was developed for specific forms of DCS, this new FSM conserves mass, while momentum and energy are conserved at spectral accuracy.
To obtain the kernel mode (l, m), ⇢ is first discretized and then 0 (⇢ r , s) is calculated.
For the (6-12) LJ potential, for each relative collision energy E, the DCS is a continuous When ⇢ r is determined, the integral given by Eq. (15) is calculated numerically, where is uniformly discretized into 8000 sections. The key part is to calculate the DCS 0 (⇢ r , ⇢ 0 ). We first check the continuity of the DCS as ⇢ 0 goes from 0 to R. If 0 (⇢ r , ⇢ 0 ) is continuous, then Eq. (15) is approximated by the Gauss-Legendre quadrature of order is discretized non-uniformly by 60 points, with most of the points located near ⇢ 0 d , while the remaining region ⇢ 0 2 [⇢ 0 d , R] is approximated by the Gauss-Legendre quadrature of order 60. In the numerical integration of 0 , a DCS with deflection angle less than 0.05 radians

IV. NUMERICAL ACCURACY
To assess the accuracy of the proposed FSM, we run two test cases. The first is the calculation of the transport coe cients of five noble gases, and the second is the calculation of mass/heat flow rates in Poiseuille/thermal transpiration flows. We compare our results with those from the variational method 24 and the discrete velocity method 10,17 .

A. Transport coe cients
The shear viscosity µ 0 and thermal conductivity  0 are calculated as where µ and  are the reduced shear viscosity and thermal conductivity, respectively. The two functions h µ (v) and h  (v) satisfy the following integral equations: To find h µ and h  , Eq. (17) is solved by the following iterative scheme (with k the iteration step): The molecular velocity space [ 6,6] 3 is discretized by 64 ⇥ 24 ⇥ 24 uniform grid points, while M = 8 is chosen in the discretization of the solid angle, see Eq. (14). Potential depths for the five noble gases are adopted from Ref. 17 Table I, where we see that the di↵erence between the FSM results and those from the variational and discrete velocity methods 17 is small: the maximum relative error is less than 0.5%.
It is interesting to see how the inverse Schmidt number, defined as the ratio of mass di↵usivity to momentum di↵usivity (viscosity), changes between the various noble gases.
Here, the mass-di↵usion coe cient is calculated as where D is the reduced mass-di↵usion coe cient and h D (v) satisfies the following equation Similar to Eq. (17), Eq. (20) is solved in the following iterative scheme: Numerical results from the FSM for noble gases and the HS gas at T = 300 K are shown in Table II, together with those from the variational method 1,25 . We find that the relative error between the two methods is about 2%.

B. Poiseuille and thermal transpiration flows
We now consider a monatomic gas confined between two parallel infinite plates located at x 2 = ±`/2. In Poiseuille flow, the wall temperature is fixed at T 0 , and a uniform pressure gradient is imposed on the gas in the x 3 direction: the pressure is given by n 0 k B T 0 (1 + ⇠ P x 3 /`) with |⇠ P | ⌧ 1. In thermal transpiration flow, the pressure is fixed at n 0 k B T 0 , but a temperature gradient is imposed on both walls: the wall temperature is T = T 0 (1 where subscripts P and T stand for the Poiseuille and thermal transpiration flows, respectively.
We assume a di↵use gas-wall interaction, so h ↵ is zero for gas molecules entering the computational domain. Due to symmetry, only half of the spatial domain is considered: the normalized x 2 varies from 1/2 to 0. The dimensionless mass and heat flow rates are These flow rates are a function of the rarefaction parameter, defined as In the numerical simulations, the spatial domain 0.5  x 2  0 is divided into 100 non-uniform sections, with most of the discrete points placed near the wall: We use the following iterative scheme to solve Eq. (22): where k is the iteration step and the spatial derivative is approximated by a second-order upwind finite di↵erence. Iterations are terminated when the relative di↵erence in mass and heat flow rates between two consecutive steps is less than 10 6 .
Tables III, IV, and V compare our numerical results for G P , G T , and Q T with those by Sharipov and Bertoldo 10 . The mass flow rate G T is not shown, as G T = Q P according to the Onsager-Casimir relation, and our numerical results show that the relative di↵erence between G T and Q P is less than 0.2%. For 0.025, the di↵erence between our results and those of Sharipov and Bertoldo is . 1%, which increases to about 2% at = 0.01.
These di↵erences are small, as the numerical accuracy of the discrete velocity method itself is about 0.8%. 10

V. APPLICATIONS
We now apply the FSM for the BE with LJ potentials to solve Couette and Fourier flows between two parallel plates. The five noble gases He, Ne, Ar, Kr, and Xe, as well as the variable HS gas, are considered and the e↵ect of the intermolecular potential on the flow properties is investigated. Note that for the variable HS gas, the DCS is proportional to |u| 1 2! , where ! is the viscosity index (i.e. the gas viscosity is proportional to T ! ). For the HS gas, ! = 0.5, while for He and Xe at T = 300 K, !=0.66 and 0.85, respectively.

A. Planar Fourier flow
The geometry is the same as that of the Poiseuille flow in Section IV B, except that the plate at x 2 = 1/2 has a temperature T 0 T /2, while the plate at x 2 = 1/2 has a temperature T 0 + T /2. Also, there is no pressure gradient along the x 1 and x 3 directions.
We first assume that the temperature di↵erence T is negligible compared to T 0 , so that the BE (4) can be linearized to Eq. (8) by expressing the VDF as f = f eq + h T /T 0 .
The spatial region 1/2  x 2  0 is discretized by 100 non-uniform grid points, with most of the grid points located near the wall. The three-dimensional molecular velocity domain [ 6,6] 3 is discretized by 32 ⇥ 128 ⇥ 32 grid points, and the number of frequency components is 32 ⇥ 48 ⇥ 32. Assuming di↵use gas-wall interaction, the boundary condition while at ) is used, and the iterations are terminated when the maximum relative di↵erence in the density n = R hdv, temperature T = 2 R hv 2 dv/3 n, and heat flux q 2 = R hv 2 v 2 dv between two consecutive steps is less than 10 5 .
Typical density and temperature profiles are shown in Fig. 1 for a rarefaction parameter of = 0.1 and T 0 = 300 K. Although they have the same rarefaction parameter, the macroscopic properties of the six gases are quite di↵erent. The di↵erences are summarized   Table VI for = 0.1, 1, and 10. At the wall, when = 0.1, the relative di↵erence in density between He and the HS gas is 22%. This di↵erence between LJ and HS potentials increases as k B T 0 /✏ decreases: the density of Xe at the wall is 40% larger than that of the HS gas. For the temperature, the largest di↵erence between the noble gases and the HS gas reaches 25%. As increases, relative di↵erences in the densities and the temperature decrease: when = 1, relative di↵erences in the density and temperature of the HS gas and Xe at the wall are reduced to 11.4% and 8.7%, respectively, while they are 1.8% and 1.5% by = 10. As further increases, the hydrodynamic flow regime is reached and there is no di↵erence between the various gases. Interestingly, the di↵erences in the heat flux between the various gases are small and first increase and then decrease with . At = 0.1, the relative heat flux di↵erence between the HS gas and Xe is only 1.1%; this increases to 2.5% at = 1, and then decreases to 2% by = 10. Therefore, if only the heat flux is of interest, the HS gas model can be safely used, with a numerical error of less than two percent.
We also consider the variable HS gas with viscosity index ! = 0.66 and 0.85, at = 0.1: at x 2 = 0.5, the gas densities are 0.051 and 0.058, respectively. When compared to that of He and Xe, we find that the variable HS model does not produce significant improvement on the HS model, i.e. there are still about 20% and 14% relative di↵erences in gas density and temperature between the variable HS model and the LJ potential, respectively.
We then consider the nonlinear heat transfer between the two parallel plates by reducing the temperature of the plate at x 2 = 0.5 to T 0 /2, while that at x 2 = 0.5 remains at T 0 = 300 K. The BE (4) is solved in an iterative manner: with the following di↵use boundary conditions where n w1 = 2 p 2⇡/T 0 R v 2 <0 fv 2 dv and n w2 = 2 p ⇡/T 0 R v 2 >0 fv 2 dv. We compare He and Kr with the HS gas, as the results for Xe are very close to Kr, while the results for Ne and Ar lie between those for He and Kr. Figure 2 shows the density and temperature profiles when = 0.1 and 1. As in the linear case, the variations in the density

B. Planar Couette flow
The geometry is the same as that for the Poiseuille flow in Section IV B, except that the plate at x 2 = `/2 moves in the x 3 direction with a speed V wall , while the other plate moves in the opposite direction with the same speed. Also, there is no pressure gradient along the We are interested in the gas velocity and the shear stress. The velocity, which is normalized by the wall speed, is V 3 = R hv 3 dv; the shear stress, which is normalized by     where The profiles of macroscopic quantities and reduced VDFs when = 0.1 are shown in Figs. 4 and 5, respectively, where we see that the relative di↵erence in gas velocity is close to that in the linearized Couette flow, and the use of the variable HS model only slightly improves the accuracy. The reduced VDF also has a relatively large di↵erence between the noble gas and the HS gas at v 2 ⇠ 0. The shear stresses in the various gases are, however, very close to each other in nonlinear Couette flow. This is because the gas temperature is around T 0 so the rarefaction parameters are nearly the same. Therefore, if only the shear stress is of interest, the HS model can be used. is not shown, since it nearly overlaps that of the HS model.

VI. CONCLUSIONS
We have presented a general fast spectral method to solve the Boltzmann equation with arbitrary intermolecular potentials. Specifically, through comparison with results from the variational and discrete velocity methods, we have demonstrated the accuracy of the FSM for the realistic (6-12) LJ potential. As an application, the FSM has been applied to planar in Couette flow are required, the hard-sphere model can be safely adopted. Otherwise, the di↵erential cross-section of a realistic intermolecular potential should be adopted when the molecular mean free path is comparable to, or smaller than, the characteristic flow length.
Although we have only considered one-dimensional flows here, the computational time required for the Boltzmann collision operator remains unchanged for two-and threedimensional flows, as the molecular velocity space is always three-dimensional. Our proposed numerical method can also be applied to mixtures of monatomic gases using ab initio potentials 13 .