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 LJ-particles. 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 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 g²/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 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 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₁ and b₁ satisfy the linear equations

where the matrix elements Apq and Bpq are

where c=δv/2T, and δv is the reduced peculiar velocity, while cc=cc-Ic2/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/2T, 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². 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 D₁₂ of a binary mixture of particles with equal masses is [3]

where d₀ 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₁ = x₂. Let c₁ and c₂ be the dimensionless peculiar velocities of components 1 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₁x₂d₀ as the coefficient to be determined, we need not worry about the exact values of x₁ and x₂. Only that they are equal. This gives a result that is independent of the values of x₁ and x₂, as it must be for the self-diffusion coefficient. We see in particular that since [c1,c2]12=-4Ω(1,1), 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 c² become more significant. Again, expressions for the bracket integrals may be obtained from [8]. For convenience, the matrix elements Apq′/x1x2 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.

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⁻¹⁴ 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 ϕ⁽ⁿ⁾ 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 g2/4=V (rm), assuming that r_m < r_s

so that the scattering cross section at high kinetic energy is approximately proportional to rm2∝g-1/3.

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⁽ⁿ⁾ 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 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⁻⁵ 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

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₁ ± a₃, while a₄ and a₅ are the position and width of the sigmoidal region. The remaining coefficients a₀ and a₂ 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.