ORIGINAL RESEARCH article
Transport Properties of the Simple Lennard-Jones/Spline Fluid I: Binary Scattering and High-Accuracy Low-Density Transport Coefficients
- PoreLab, Department of Chemistry, Norwegian University of Science and Technology, Trondheim, Norway
Inspired by recent work on the thermodynamic properties of the Lennard-Jones/spline (LJ/s) fluid , 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.
Knowledge of the transport properties of fluids is key to understanding a wide range of non-equilibrium 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) , 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 . 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 , and for prediction through entropy scaling, as in e.g., . High-accuracy 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.
2.1. Binary Scattering
which behaves like a Lennard-Jones potential up to a separation equal to the inflection point r = rs. At longer separation, the potential drops smoothly to uniform at r = rc 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 = rs. Exact values of these parameters can be found in . 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 LJ-particles. As we will see, this feature can lead to appreciable differences in the particle dynamics when collision trajectories are off-center.
Figure 1. The Lennard-Jones/spline potential and the untruncated Lennard-Jones potential, plotted along with their corresponding forces.
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 time-independent, 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 g2/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
By spacetime inversion symmetry, the trajectory is symmetric about the line connecting the two particles when the interparticle distance is at the trajectory minimum rm. 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.
2.1.1. The Distance of Closest Approach
The piecewise nature of the LJ/s potential means that the equation determining rm depends on whether or not rm < rs. In any case, the physical root must be the largest real root of dr/dt = 0 satisfying rm ≤ rc, where we restrict our attention to ℓ ≤ rc. For ℓ > rc, the particles will always remain out of range, so rm = ℓ and trivially Θ = 0. Whether or not rm < rs depends on whether or not the relative kinetic energy of the particles can overcome any energy barrier along the way to r = rs. The radial energy barrier is due to the centrifugal potential offset by the attractive potential. The height of this barrier is the maximum of
with rmax ≥ rs. Then, rm < rs if and only if
If rm ≥ rs, we need to find the largest real rm such that
or, equivalently, the smallest real um : = 1/rm 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 |umI − 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 rc and gives a positive radial acceleration. The latter constraint is
If no such eigenvalue exists, we infer that rm < rs, in which case we need to find the largest real rm such that
or the smallest real um such that
subject to the constraint that rm < rs. 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 (ui) is smaller than some tolerance. This way of computing rm ensures that we always obtain the physical solution. This can also be double-checked 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 rm.
2.1.2. The Scattering Angle
Having obtained accurate values of rm, the next step is to find the scattering angle Θ by numerical integration of (7), where rm 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 , such that
then, we apply a tanh-sinh type transformation, suitable for dealing with endpoint singularities
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 . The scattering angle is then Θ = π − 2Φ.
2.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 λ 
where T is the reduced kinetic temperature, and the coefficients a1 and b1 satisfy the linear equations
where the matrix elements and are
where , and δv is the reduced peculiar velocity, while . The Sonine polynomials are
where the subscript means rq = 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) .
where , 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 c2. We consider also the self-diffusion 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 D12 of a binary mixture of particles with equal masses is 
where d0 is found from the solution to
with n the particle density, and xi = ni/n, where ni is the number density of particles of component i. We are interested in the case of a simple fluid, where x1 = x2. Let c1 and c2 be the dimensionless peculiar velocities of components 1 and 2. Manipulation of the expressions given in  gives in that case
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 x1x2d0 as the coefficient to be determined, we need not worry about the exact values of x1 and x2. Only that they are equal. This gives a result that is independent of the values of x1 and x2, as it must be for the self-diffusion coefficient. We see in particular that since , we obtain the first approximation
which is a familiar result in the literature. Calculating higher-order 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 c2 become more significant. Again, expressions for the bracket integrals may be obtained from . For convenience, the matrix elements 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.
3. Results and Discussion
3.1. 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  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 . Both rm 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.
Figure 2. The distance of closest approach, rm, and the scattering angle, Θ, from binary scattering of Lennard-Jones (solid lines) and Lennard-Jones/spline-particles (dashed lines) given as functions of the relative speed g for sampled values of the impact parameter ℓ.
We observe that when g is small, the dynamics are dominated by the potential energy of interaction, and rm → 1 because that is the point where the mutual potential energy is equal to that at infinite separation. When g is large, rm → ℓ 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, rc, from below. The change in rm 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 , 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.
Figure 3. Samples of the total cross sections ϕ(n) calculated for both Lennard-Jones (solid lines) and Lennard-Jones/spline (dashed lines) particles, as a function of the relative speed g. The stars are tabulated values of ϕ(1) for the LJ potential from .
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 rm. The distance of closest approach for head-on collisions is obtained by solving , assuming that rm < rs
so that the scattering cross section at high kinetic energy is approximately proportional to .
3.2. 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 gs and a lower bound gi 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 gs >> 1 and gi << 1, and was kept below machine precision by choosing gs large enough and gi 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 momentum-space 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.
Figure 4. Discretization and truncation errors of the calculated transport coefficients of the Lennard-Jones/spline fluid. Top left: the estimated discretization error of a grid with half the number of points as the largest grid used. Top right: The variation of the inf-norm of the discretization error with the number of grid points, proportional to the size parameter N. Bottom left: the estimated truncation error by neglecting 5th order terms and higher (truncation order 4) in the Sonine polynomial expansions. Bottom right: The variation of the inf-norm of the truncation error with the truncation order.
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 . 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
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 a1 ± a3, while a4 and a5 are the position and width of the sigmoidal region. The remaining coefficients a0 and a2 are simply bias-correction coefficients adjusting the positions of the optimal curves along the ln (ψ)-axis.
Figure 5. Logscale plot of transport coefficients against kinetic temperature. Bottom: Logscale plot of relative error in the optimized parameterization given by Equation (34).
Table 1. Optimized numerical values of the parameters defined in Equation (34) for each of the transport coefficients.
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 rc, 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.
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.
Figure 6. Comparison of the transport coefficients ψ = D, η, λ of the Lennard-Jones and Lennard-Jones/spline fluids at different temperatures T, where D is the self-diffusion coefficient, η is the shear viscosity, and λ is the thermal conductivity.
Data Availability Statement
The datasets generated for this study are available on request to the corresponding author.
The author confirms being the sole contributor of this work and has approved it for publication.
Article processing charges are covered by the publishing fund of the Norwegian University of Science and Technology.
Conflict of Interest
The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The author is grateful to the Research Council of Norway through its Centers of Excellence funding scheme, project number 262644, PoreLab. Prof. Bjørn Hafskjold is thanked for fruitful discussions.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2020.00271/full#supplementary-material
5. Kim S, Monroe C. High-accuracy calculations of sixteen collision integrals for Lennard-Jones (12-6) gases and their interpolation to parameterize neon, argon, and krypton. J Comp Phys. (2014) 273:358–73. doi: 10.1016/j.jcp.2014.05.018
8. Tompson RV, Tipton EL, Loyalka SK. Chapman-Enskog solutions to arbitrary order in Sonine polynomials IV: summational expressions for the diffusion- and thermal conductivity-related bracket integrals. Eur J Mech B Fluids. (2009) 28:695–721. doi: 10.1016/j.euromechflu.2009.05.002
9. Tompson RV, Tipton EL, Loyalka SK. Chapman-Enskog solutions to arbitrary order in Sonine polynomials V: summational expressions for the viscosity-related bracket integrals. Eur J Mech B Fluids. (2010) 29:153–79. doi: 10.1016/j.euromechflu.2009.10.002
10. Guennebaud G, Jacob B. Eigen v3 (2010). Available online at: http://eigen.tuxfamily.org
11. Boost C++ Libraries. Available online at: http://www.boost.org/
Keywords: single phase fluid, kinetic theory, transport properties, scattering, molecular dynamics
Citation: Kristiansen KR (2020) Transport Properties of the Simple Lennard-Jones/Spline Fluid I: Binary Scattering and High-Accuracy Low-Density Transport Coefficients. Front. Phys. 8:271. doi: 10.3389/fphy.2020.00271
Received: 09 January 2020; Accepted: 17 June 2020;
Published: 28 July 2020.
Edited by:Ramesh L. Gardas, Indian Institute of Technology Madras, India
Reviewed by:Sergey Khrapak, Helmholtz Association of German Research Centers (HZ), Germany
Alejandro Gil-Villegas, University of Guanajuato, Mexico
Copyright © 2020 Kristiansen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Kim R. Kristiansen, firstname.lastname@example.org