Role of three-body interactions in formation of bulk viscosity in liquid argon

With the aim of locating the origin of discrepancy between experimental and computer simulation results on bulk viscosity of liquid argon, a molecular dynamic simulation of argon interacting via ab initio pair potential and triple-dipole three-body potential has been undertaken. Bulk viscosity, obtained using Green-Kubo formula, is different from the values obtained from modeling argon using Lennard-Jones potential, the former being closer to the experimental data. The conclusion is made that many-body inter-atomic interaction plays a significant role in formation of bulk viscosity.


I. INTRODUCTION
Argon above its melting temperature is a typical simple fluid.Consisting of spherical atoms that interact via short-range repulsion and long-range attraction, and are heavy enough for the quantum effects to be small, fluid argon and heavier noble gases are an excellent choice of a real system to be used for testing various approaches in classical theory of fluids.
An inter-particle interaction in argon is commonly represented by a well known 12-6 Lennard-Jones pair potential 1 , The two parameters, σ LJ and ǫ LJ , are usually determined by fitting thermodynamic properties, derived from the potential (1) by theoretical or computational methods, to corresponding experimental data.
It is known that Lennard-Jones potential is only an approximation to real interaction in argon.Several experimental results obtained for argon at large pressures are better explained if a larger steepness, compared to Lennard-Jones, of argon-argon interaction potential at small inter-atomic separation distances is taken into account 2,3 .Accurate argon-argon interatomic potentials have been calculated by direct ab initio quantum chemical calculations [4][5][6] or obtained by inversion of experimental data 7 .Moreover, many-body dispersion, exchange and induced polarization contributions to inter-atomic interactions are not small and noticeably influence thermodynamic properties of argon 8,9 .The most widely used of these contributions is triple-dipole dispersion interaction, derived by Axilrod and Teller 10,11 and Muto 12 , and account of this contribution in addition to ab initio pair potential is sufficient to describe thermodynamic properties of argon with good accuracy [13][14][15][16][17] .
By virtue of Henderson theorem 18,19 , which states that, for fluids with only pairwise interactions, and under given conditions of temperature and density, the pair potential which gives rise to a given radial distribution func-tion g(r) is unique up to a constant, the thermodynamic properties of the system with many-body interactions can be described by a model system with an appropriate effective pair potential.Generally, the effective potential depends on the thermodynamic state of the system and thermodynamic property to be described [20][21][22] .Van der Hoef and Madden 21 have demonstrated that the account of triple-dipole and dipole-dipole-quadrupole dispersion interactions moves the effective potential of argon towards Lennard-Jones form (1).Moreover, the possibility of consistent description of many thermodynamic properties of argon, using Lennard-Jones potential in a wide domain of thermodynamic states [23][24][25] , suggests that the state dependence of the effective potential is weak.
There is no analogous reason for kinetic properties of a system with many-body interactions to be equivalent to those of a system with a corresponding effective pair potential.Nevertheless, experimental data on self-diffusion, shear viscosity and thermal conductivity coefficients of argon have been shown to be accurately described by Lennard-Jones model with the parameters obtained by fitting thermodynamic data 26,27 .
Bulk viscosity is a noticeable exception.Bulk viscosity of argon has been measured experimentally [28][29][30][31][32][33][34][35] , and its behavior can be qualitatively described by the results of a molecular dynamics simulation of a Lennard-Jones system 36 .However, when results of simulations with Lennard-Jones potential are rescaled in an attempt to describe experimental data liquid argon, bulk viscosity, contrary to other kinetic properties, appears strongly underestimated (e.g. up to 50% in Ref. [27)].
In view of the above, I propose that the source of this discrepancy may lie in neglect of many-body interactions.Previous molecular dynamics simulations of systems consisting of 108 particles interacting via ab initio pair potential and Axilrod-Teller-Muto (ATM) interaction indicated that a triple-dipole interaction does not affect the bulk viscosity of liquid xenon near its triple point 37 and dense gaseous krypton 38 .However, the error in the values of bulk viscosity obtained from molecular dynamics simulation of the systems with such a small number of particles can be quite large.For example, the values of the reduced bulk viscosity of the Lennard-Jones systems consisting of 128 and 256 particles at the reduced temperature T * = 0.722 and the reduced density ρ * = 0.8442, reported in Refs 36, 39-41, range from 0.89 to 1.47, with the ratio of the latter to the former of 1.65.
This paper presents the results of more accurate molecular dynamics simulations of a liquid consisting of 1372 argon atoms with ab initio+ATM interaction, which demonstrate that bulk viscosity, determined from Green-Kubo formulae, significantly changes with the account of three-body interaction, moving results towards experimental data.

II. INTERACTION
Nasrabad et al 16 undertook a Monte Carlo simulation of argon using combination of ab initio pair interaction 4 and ATM triple-dipole dispersion interaction 10 to test their ability to predict vapor-liquid equilibrium.Although more accurate ab initio pair potentials for argon have become available recently 5,6 , and other manybody contributions to inter-atom interaction can be calculated 8 , we use the same interaction as Nasrabad et al because, being able to predict accurately the phase diagram of argon 16 , it is computationally more efficient.
Specifically, the ab initio pair interaction potential used in the present work is described by a function 16 where and numerical values of the parameters A, α, β, b, and C 2n are given in Ref. [16].The ATM triple-dipole interaction has form 10 where the r ik are the lengths of the sides, α, β, and γ are the angles of the triangle formed by three argon atoms, and ν = 7.32 • 10 −108 J•m 9 for argon 13,14 .For simulations of argon using Lennard-Jones potential (1) the values σ LJ = 3.3952 Å and ǫ LJ = 116.79K are used 25 .

III. SIMULATION
Meier et al 36 undertook a systematic study of the influence of the number of particles and the cutoff radius for pair interaction on the bulk viscosity of Lennard-Jones system.In view of their results, simulations were performed in a cubic box containing N = 1372 particles, and the cutoff radius for pair interactions was set to 5σ LJ .Three-body interactions were cut off when the distance between any pair of the atoms in the triplet exceeded one quarter of the simulation box length (around 3σ LJ for the densities studied in this work).Usual periodic boundary conditions and minimum image convention were applied.The simulations were started with the particles in a face-centered-cubic lattice, with randomly assigned velocities.Forces arising from three-body interactions were calculated using formulas given by Allen and Tildesley 42 , and an expression for forces due to ab initio pair interaction was obtained by applying gradient operator to Eq. ( 2).Newton's equations of motion were solved using velocity-Verlet algorithm with the time step ∆t • ǫ LJ /m/σ LJ = 0.003.
The runs were made at the experimental densities at various temperatures along the 40 atm isochore, taken from Ref. [33].Every simulation was initiated in the NVT ensemble and run for at least 2•10 5 time steps to attain thermodynamic equilibrium.After equilibration the thermostat was turned off and the NVE ensemble was invoked to calculate bulk and shear viscosities.The length of the production period was 4•10 6 time steps for the system interacting via Lennard-Jones potential, and between 10 6 and 3•10 6 time steps for the system with ab initio + ATM interaction, depending on the state point.
Bulk viscosity, ζ, and shear viscosity, η, were calculated using Green-Kubo formulas 43 : where V is volume, k B is Boltzmann constant, T is temperature, t is time, δp = p − p is the deviation of the instantaneous pressure p from its average value p , σ αβ is an off-diagonal element of the stress tensor, the angular brackets denote equilibrium ensemble averages over short trajectory sections of the phase-space trajectory of the system with multiple (every time step) time origins t 0 .The stress tensor was calculated using formulae given by Lee and Cummings 44 .The integration in Eqs ( 5) and ( 6) was carried out up to τ L = L/c, where L is simulation box length and c is sound velocity taken from Ref. [33].Depending on the state point, the value of τ L was between 4.80 and 11.25 ps.The statistical error in time correlation functions was estimated using formula given by Frenkel and Smit 45 , where t run is the length of the simulation, and the correlation time τ X was approximated as the time during which time correlation function decays e ≈ 2.718 times.
Time-correlation functions C(t) used for calculation of bulk viscosity at density 1.258 g/cm 3 .Solid and dashed lines correspond to the simulation results with ab initio + ATM and Lennard-Jones interaction, respectively.

IV. RESULTS
Fig. 1 and Table I present simulation results for the bulk viscosity obtained using ab initio + ATM (Eqs (2) and ( 4)) and Lennard-Jones (Eq.( 1)) interaction, respectively.Bulk viscosity, determined from Green-Kubo formulas, changes with the account of three-body interaction, moving towards experimental data.However, this change is not sufficient to obtain numerical agreement with experiment, especially at lower densities.Typical behavior of time correlation functions C(t) = δp(t)δp(0) is shown in Fig. 2. Fernandez et al 27 demonstrated that, contrary to bulk viscosity, the values of shear viscosity of argon obtained from molecular dynamics simulation of a Lennard-Jones system agree with experimental data.Lee and Cummings 44 and Marcelli et al 46 found that the influence of triple-dipole interaction on shear viscosity of argon is small.The results of the present simulation, shown in Fig. 3 and Table I, agree with these findings.

V. CONCLUSION
The message of this paper is that many-body interactions play a more substantial role in determining the value of the bulk viscosity than other transport coefficients.The present results from the molecular dynamic simulation of liquid argon demonstrate that even account of a single many-body contribution, ATM triple-dipole interaction, shifts the values of the bulk viscosity of argon towards experimental data.Larger sensitivity of the bulk viscosity to many-body interaction, compared to other transport coefficients, can be intuitively explained in the case of gaseous state.Bulk viscosity of a non-relativistic monoatomic gas calculated from the Boltzmann equation, which takes into account only pair collisions of atoms, appears to be zero, in contrast to heat conductivity and shear viscosity which have non-zero values in the same approximation 47 .A non-zero value of bulk viscosity appears in the approximations corresponding to higherorder terms in the virial expansion 48,49 , which correspond to the explicit account of at least three-atom collisions which, in turn, are sensitive to three-body inter-atomic interaction.
T , K ρ, g/cm 3 Bulk  I. Bulk and shear viscosities of argon obrained from molecular dynamics simulations using Lennard-Jones (LJ) and ab initio pair + Axilrod-Teller-Muto three-body (AI+ATM) interaction, and corresponding experimental data 33,50 .Error in the simulation data is calculated using Eq. ( 7).

3 Figure 1 .
Figure 1.Bulk viscosity of liquid argon at T = (90−140) K.Error bars connected with solid and dashed lines correspond to the simulation results with ab initio + ATM and Lennard-Jones interaction, respectively.Experimental points are taken from Refs[33] (circles, pressure 40 atm) and[29] (square with error bar, pressure 40 kg/cm 2 ).

3 Figure 3 .
Figure 3. Shear viscosity of liquid argon at T = (90−140) K. Error bars connected with solid and dashed lines correspond to the simulation results with ab initio + ATM and Lennard-Jones interaction, respectively.Dotted line corresponds to the interpolation data for pressure 40 atm taken from Ref. [50].