Transport Properties of the Simple Lennard-Jones/Spline Fluid I: Binary Scattering and High-Accuracy Low-Density Transport Coefficients

Inspired by recent work on the thermodynamic properties of the Lennard-Jones/spline (LJ/s) fluid [1], we have considered the transport properties of the simple LJ/s fluid. The binary scattering problem for LJ/s particles was solved numerically, and results were compared to the untruncated LJ fluid. The scattering dynamics are affected both by the restricted range of the LJ/s potential, and the stronger attraction between LJ/s particles at distances between the inflection point of the potential and the cutoff range. At small relative kinetic energies, it was found that the scattering cross section of the LJ/s particles is much smaller than that of the LJ particles. The shear viscosity, thermal conductivity, and the self-diffusion coefficient were calculated from the scattering cross sections by the Chapman-Enskog method, and a six-parameter equation with a worst case accuracy of roughly 1% over the temperature interval [0.1, 1000] in LJ units is provided. The smaller scattering cross section at low kinetic energies leads to transport coefficients of the LJ/s fluid to be greater than those of the LJ fluid at low temperatures, and were all found to be roughly 50% greater at T = 0.1, which is the lowest temperature considered.


INTRODUCTION
Knowledge of the transport properties of fluids is key to understanding a wide range of nonequilibrium behavior, such as thermal conduction, viscous flow, and diffusion. While the theory of non-equilibrium thermodynamics provides closure relations for describing fluxes in terms of thermodynamic driving forces, the functional dependence of the corresponding transport coefficients on thermodynamic variables must be supplied either from experiments, or calculated by means of statistical mechanics from a model interaction potential. The Lennard-Jones 12-6 potential (LJ) is a particularly ubiquitous and simple model potential that captures both the attraction between electrically neutral particles at long range, and their mutual repulsion at short range. The infinite range of this potential is, however, the cause of some problems in the application of periodic boundary conditions in direct molecular dynamics simulations. For this reason, the potential is often replaced in practical calculations by a truncated version with a finite range. One such potential is the Lennard-Jones/spline model (LJ/s) [2], which truncates the potential smoothly by means of a cubic spline. Recently, a detailed assessment was made on the effect of this type of truncation on the equilibrium thermodynamic properties of the resulting fluid [1]. Inspired by this work, it is of interest to see how this truncation affects non-equilibrium properties.
Transport properties of the untruncated Lennard-Jones fluid have been studied for nearly a century, and transport coefficients are known to a high degree of precision over a wide range of temperatures and densities. While the properties of the LJ and LJ/s fluids are expected to be similar, their difference may be of great importance in precise calculations of quantities of theoretical interest, such as the entropy generated by irreversible processes. In particular, the Chapman-Enskog values, which are expected to be accurate at low densities and low Knudsen numbers, are useful for several reasons-they supplement data from molecular dynamics in the low-density regime where the small number of particles gives low-quality statistics, and are also used as reference values in the modified Enskog theory [3], and for prediction through entropy scaling, as in e.g., [4]. Highaccuracy calculations are therefore important as the basis for further work, such as the prediction of transport coefficients at higher densities, which will be the topic for a future paper.

Binary Scattering
The binary scattering problem of the Lennard-Jones fluid is already well-studied, and detailed recent analyses can be found in [5] and [6]. We consider instead the Lennard-Jones/ spline potential which behaves like a Lennard-Jones potential up to a separation equal to the inflection point r = r s . At longer separation, the potential drops smoothly to uniform at r = r c by means of a cubic spline. The coefficients a and b are chosen such that V (r) and the force f (r) = −∇V (r) are continuous at r = r s . Exact values of these parameters can be found in [1]. We consider here the case that the masses of the scattering particles are equal, and all quantities are dimensionless (Lennard-Jones units). The shape of the potential and the corresponding force is shown in Figure 1 along with the standard LJ potential. In a large part of the spline region, the force of attraction between LJ/s-particles is significantly larger than the corresponding force between LJparticles. As we will see, this feature can lead to appreciable differences in the particle dynamics when collision trajectories are off-center. We are interested in the shape of trajectories for different possible initial relative velocities and positions. Let θ be the angle between the particle positions with θ = 0 parallel to the initial relative velocity. A useful quantity is the relative angular momentum which is a conserved quantity by virtue of V(r) being a central-force (θ -independent) potential. Since V(r) is also timeindependent, the total energy is also conserved. Let g be the initial relative velocity as the particles approach each other from infinity, and let the impact parameter ℓ be such that the initial (and constant) angular momentum is gℓ/2. The initial (and constant) total energy is then simply g 2 /4. Varying g and ℓ from 0 to ∞ then allow us to explore all possible trajectories during such an encounter. The trajectory may be found by solving the system Conservation of energy allows us to integrate the radial equation from second to first order However, the sign is ambiguous, and care must be taken to avoid unphysical solutions. Solving for dt and multiplying by p θ gives an equation for the trajectory dθ = ℓdr By spacetime inversion symmetry, the trajectory is symmetric about the line connecting the two particles when the interparticle distance is at the trajectory minimum r m . Let the angular coordinate at this point be . Since the particles approach from infinity, we have The initial angle at infinity is π, and we are interested in the scattering angle = π − 2 . The scattering angle as function of g and ℓ fully determines the contribution of binary collisions to the macroscopic hydrodynamical modes of the fluid in the low-density limit.

The Distance of Closest Approach
The piecewise nature of the LJ/s potential means that the equation determining r m depends on whether or not r m < r s . In any case, the physical root must be the largest real root of dr/dt = 0 satisfying r m ≤ r c , where we restrict our attention to ℓ ≤ r c . For ℓ > r c , the particles will always remain out of range, so r m = ℓ and trivially = 0. Whether or not r m < r s depends on whether or not the relative kinetic energy of the particles can overcome any energy barrier along the way to r = r s . The radial energy barrier is due to the centrifugal potential offset by the attractive potential. The height of this barrier is the maximum with r max ≥ r s . Then, r m < r s if and only if If r m ≥ r s , we need to find the largest real r m such that or, equivalently, the smallest real u m : = 1/r m such that with monomial coefficients We solve this equation in two stages-first, we locate all the real roots of the equation by diagonalizing the companion matrix such that |u m I − C| = 0 is equivalent to (11). The true turning point is characterized by a positive radial acceleration. A secure way of computing the physical solution is then to diagonalize C and pick the largest real diagonal element that is smaller than r c and gives a positive radial acceleration. The latter constraint is If no such eigenvalue exists, we infer that r m < r s , in which case we need to find the largest real r m such that or the smallest real u m such that subject to the constraint that r m < r s . The computation of eigenvalues is restricted in precision by the conditioning of the companion matrices. Therefore, we introduce as a second calculation stage that the values obtained by this method be refined further by using them as initial guesses to a more precise root-finding algorithm. For this, we choose Halley's method where P is either of the polynomials (11) or (16), and the prime denotes the derivative with respect to the argument. The convergence criterion is that the absolute value of P(u i ) is smaller than some tolerance. This way of computing r m ensures that we always obtain the physical solution. This can also be doublechecked by direct integration of (4), which is a comparably inefficient two-particle molecular dynamics approach, but serves as a way to check the consistency of the obtained value of r m .

The Scattering Angle
Having obtained accurate values of r m , the next step is to find the scattering angle by numerical integration of (7), where r m appears as one of the integration limits. In this limit, the integrand has a singularity, which will need to be dealt with some care. Let u : = r m r −1 , such that Let then, we apply a tanh-sinh type transformation, suitable for dealing with endpoint singularities (21) Applying the trapezoidal rule with spacing x, the quadrature is then which is further approximated by truncating the sum. The singularity is smeared out over an unbounded interval, and the quadrature exhibits exponential convergence to the analytical solution [7]. The scattering angle is then = π − 2 .

Transport Coefficients
The low-density, low-Knudsen number limit of the transport coefficients are obtained by solving the Boltzmann equation for the one-particle probability density, and then calculating the expectation values corresponding to the heat flux, pressure tensor, and the diffusion flux in order to obtain the thermal conductivity λ, the shear viscosity η, and the self-diffusion coefficient D, respectively. The bulk viscosity of the simple fluid is zero in this limit, because the particles have no internal degrees of freedom. We choose the Chapman-Enskog approach to solving the equation, which gives the following expressions for η and λ [3] where T is the reduced kinetic temperature, and the coefficients a 1 and b 1 satisfy the linear equations (24) where the matrix elements A pq and B pq are where c = δv/ √ 2T, and δv is the reduced peculiar velocity, while cc = cc − Ic 2 /3. The Sonine polynomials are where the subscript means r q = r (r − 1) (r − 2) . . . r − q + 1 .
The brackets in (25) denote integrals that involve the scattering angles obtained by solving the binary scattering problem, and can be expressed as linear combinations of the well-known collision integrals (n,s) [3].
where g ′ = g/2 √ T, with g the relative velocity, ℓ is the impact parameter defined in the preceding section, and is the scattering angle. General expressions for the coefficients in these linear combinations were derived by Tompson et al. [8,9], who also gave explicit calculations up to order 5.
In practice, the linear systems in (24) are solved approximately by truncating the sums at some finite number k and solving the resulting k × k linear systems, which is to truncate the series expansion of the Chapman-Enskog solution to the corresponding order in c 2 . We consider also the selfdiffusion coefficient, which can be obtained as the special case of the interdiffusion coefficient of a binary mixture when both components are actually the same. The interdiffusion coefficient D 12 of a binary mixture of particles with equal masses is [3] where d 0 is found from the solution to with n the particle density, and x i = n i /n, where n i is the number density of particles of component i. We are interested in the case of a simple fluid, where x 1 = x 2 . Let c 1 and c 2 be the dimensionless peculiar velocities of components 1 Frontiers in Physics | www.frontiersin.org and 2. Manipulation of the expressions given in [3] gives in that case and where the brackets are analogous to those in (25), with subscripts 1 and 2 denoting collisions between pairs of particles of types 1 and 2, respectively, and the subscript 12 denoting collision between the two different types. Treating x 1 x 2 d 0 as the coefficient to be determined, we need not worry about the exact values of x 1 and x 2 . Only that they are equal. This gives a result that is independent of the values of x 1 and x 2 , as it must be for the selfdiffusion coefficient. We see in particular that since [c 1 , c 2 ] 12 = −4 (1,1) , we obtain the first approximation which is a familiar result in the literature. Calculating higherorder collision integrals allows us to find correction terms to this result up to any desired order in the polynomial expansion. Such higher-order corrections are more important at higher temperatures, as the contributions from greater powers of c 2 become more significant. Again, expressions for the bracket integrals may be obtained from [8]. For convenience, the matrix elements A ′ pq /x 1 x 2 have been calculated from the bracket integrals and are given in the Supplementary Material up to order 5. Since the collision integrals all depend only on the temperature, the quantity nD is a function of the temperature only.

Binary Scattering
Equations (11) and (16) were solved on a patchwork of uniform grids on the g, ℓ-plane numerically by implementation of a looping algorithm in C++, using the Eigen package [10] to calculate the eigenvalues of the companion matrices. Equation (7) was then solved for each grid point, by implementing the quadrature given in Equation (23), using the Boost C++ library [11]. Both r m and were calculated within an error tolerance of 10 −14 in reduced units. The distance of closest approach and the scattering angle obtained for both LJ/s and LJ are shown as functions of g for a few values of ℓ in Figure 2.
We observe that when g is small, the dynamics are dominated by the potential energy of interaction, and r m → 1 because that is the point where the mutual potential energy is equal to that at infinite separation. When g is large, r m → ℓ as the energy barrier between the particles is dominated by the centrifugal force. The difference between the LJ and the LJ/s potential can be observed when ℓ approaches the cutoff radius, r c , from below. The change in r m is more sudden as a function of g for the LJ/s potential, because the particles only interact in a short window of time as the particles pass one another, and increasing g decreases the length of this time interval. This has a more pronounced effect on the scattering angle, where we see that when g is small and ℓ is large, the LJ/s-particles only deflect each other slightly from their initial paths, as their window of interaction becomes very short. We observe the singularity where → −∞, where particles go into mutual orbit for an arbitrarily large number of revolutions before separating.
The total cross sections φ (n) for n = 1, . . . , 6 were calculated from the scattering angles, and some of them are given as functions of g in Figure 3. For the Lennard-Jones particles, an error was introduced by truncating the calculations at an upper bound for ℓ. Using the estimate given by Equation (24) in [5], this error was kept below machine precision by choosing the upper bound for ℓ large enough. A discretization error was also introduced, and will be assessed by its effect on the transport coefficients.
At low velocities, the cross sections of the LJ particles are orders of magnitude greater than those of the LJ/s-particles, due to the infinite range of the LJ-potential. At high velocities, their cross sections converge because the repulsive part of the potential, which is the same for both types, dominates the dynamics in this regime. The cross sections fall off asymptotically as the inverse cube root of the relative speed, due to scattering events becoming more similar to hard-sphere collisions with an effective radius proportional to r m . The distance of closest approach for head-on collisions is obtained by solving g 2 /4 = V(r m ), assuming that r m < r s r m = 1 2 (33) so that the scattering cross section at high kinetic energy is approximately proportional to r 2 m ∝ g −1/3 .

Transport Coefficients
Having obtained the cross sections from solving the scattering problem, the collision integrals (n,s) needed for calculating the transport coefficients may be obtained. The collision integrals suffer a truncation error from calculating to an upper bound g s and a lower bound g i for g, which is easily assessed by considering the asymptotic behavior of the cross sections. Since the cross  sections fall off approximately as g −1/3 at large g, and falls off even slower at small g, the truncation error is, approximately where A (n) is some n-dependent constant. This becomes very small for g s >> 1 and g i << 1, and was kept below machine precision by choosing g s large enough and g i small enough for the temperature range of interest. Having control over the truncation errors in the ℓ, g-plane, the transport coefficients still suffer from two additional sources of error-the error from discretizing the ℓ, g-plane, and the error from truncating the momentumspace polynomial expansion of the solution to the Boltzmann equation. To assess the former, the transport coefficients of the LJ/s fluid were calculated using different grid sizes, controlled by a size parameter N. For the latter, the coefficients were calculated with different truncation orders in momentum space. These different solutions were compared to the most accurate solution, to obtain upper bounds for the relative errors. This is summarized in Figure 4, where we see that the upper bound for the relative error due to discretization and truncation are both on the order 10 −5 for the most accurate calculations that were carried out. The estimates are similar for the unmodified Lennard-Jones fluid.   The calculated values of η, λ, and nD of the LJ/s fluid are shown in Figure 5 as functions of the temperature. The lowest temperature considered here is 0.1, where the density has to be extremely low in order to stay away from the gas-liquid phase transition, see phase diagram in [1]. The highest temperature was chosen to be 1,000, to demonstrate the convergence to ordinary Lennard-Jones behavior at high temperatures. All transport coefficients appear to have two distinct power law regimes, separated by a sigmoidal region in the neighborhood of the critical temperature. This motivates a parameterization on the form ln (ψ) = a 0 + a 1 ln (T) + a 2 + a 3 ln (T) tanh ln (T) − a 4 a 5 (34) where ψ = η, λ, nD. The optimal parameters were found by a uniform search in the parameter space, applying Newton's method to approach the local optima. The best optimum found is assumed to be the global optimum. The relative errors of the optimal fitting are also given in Figure 5 as a function of temperature, and the numerical values of the optimal parameters for each of the transport coefficients are given in Table 1. The two asymptotical power law exponents are then a 1 ± a 3 , while a 4 and a 5 are the position and width of the sigmoidal region. The remaining coefficients a 0 and a 2 are simply bias-correction coefficients adjusting the positions of the optimal curves along the ln (ψ)-axis. We give also a comparison between the transport coefficients of the LJ/s fluid and those of the untruncated potential. As could be expected from the differences in the scattering cross sections at low velocities, the differences between the two fluids manifest themselves at low temperatures. At very low temperatures, the transport coefficients of the LJ/s fluid are significantly larger than those of the LJ fluid, around 50% at T = 0.1. In the Chapman-Enskog limit, that is, low density and low Knudsen number, the only mode of transport taken into account is that due to translation of particles, as opposed to transfer by collisions. Due to the truncated interaction range of the LJ/s particles, those particles generally move more freely at long range than LJ particles at the same temperature. At very low temperatures, the low kinetic energy of the average particles is more easily overcome by the infinite-ranged, though weak, attraction between LJ particles, while LJ/s particles are completely unaffected at distances greater than r c , no matter how low the kinetic energy is. This leads to an overall more efficient transport of particles and their momenta in a given direction. At high temperatures, the differences between the two fluids are negligible.

CONCLUSIONS
We have provided an account of the properties of binary scattering events between LJ/s particles and their untruncated LJ counterparts. The main difference is that for LJ/s particles, the scattering angle is always zero when the impact parameter is greater than the cutoff radius of the potential, whereas the untruncated potential gives rise to a significant region of deflecting and oribiting collisions even at large impact parameters when the kinetic energy is small. The attractive force experienced by the LJ/s particles between the cutoff radius and the inflection point of the potential is stronger than that of the LJ particles, which leads to a slightly different scattering angle landscape in the region of orbital collisions. These differences lead to the LJ/s potential providing a much smaller scattering cross section than that provided by the untruncated potential, which again leads to significant differences in transport properties at low temperatures.
The shear viscosity η, the thermal conductivity λ, and the self-diffusion coefficient D have been calculated, by application of the Chapman-Enskog method, for temperatures between 0.1 and 1000 in LJ units, and a six-parameter equation (34) has been provided, with a worst-case accuracy of roughly 1% across the entire temperature range by using the coefficient values listed in Table 1. Such a model is useful for making quick calculations where an error of up to 1% is acceptable. It was observed that all transport coefficients appear to obey a power law as a function of the temperature at both high and low temperatures, with an intermediate transition region close to the critical temperature.
Comparisons were made between the LJ/s fluid and the untruncated LJ fluid, see Figure 6. At high temperatures, the transport properties are dominated by the repulsive part of the potential, which is identical for the two potentials. Therefore, the transport coefficients of the two fluids converge to the same values as the temperature is increased. At lower temperatures, the differences in the transport coefficients between the two fluids become significant. At temperatures below unity, the values diverge rapidly away from one another, with transport coefficients of the LJ/s being roughly 50% larger than those of the LJ fluid at T = 0.1, the lowest temperature considered in this work. When the LJ/s potential replaces the LJ potential for the sake of convenience in molecular dynamics simulations, it is important to be aware of these differences when considering subcritical temperatures. Further work will assess transport properties at increased densities, using the values obtained here in the low-density limit.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and has approved it for publication.

FUNDING
Article processing charges are covered by the publishing fund of the Norwegian University of Science and Technology.