Density fluctuations and random walks in an overdamped and supercooled simple liquid

TL;DR

This paper investigates the short-time dynamics of a simple liquid, focusing on how density fluctuations influence particle motion, and provides analytical and numerical insights into diffusion behavior in supercooled and overdamped states.

Contribution

It offers a combined analytical and numerical analysis of density fluctuations and particle diffusion in a Lennard-Jones liquid, linking microscopic behavior to macroscopic diffusion coefficients.

Findings

Analytical expression matches experimental diffusion data for liquid argon.

Density fluctuations are key to understanding reduced particle mobility.

Behavior is similar in supercooled and overdamped regimes.

Abstract

In this work, the short-time dynamics of simple liquid is explored both analytically and numerically with the focus on the interplay between the density fluctuations in a volume surrounding a chosen particle and its random walk motion. The particles interact via the Lennard-Jones potential with parameters corresponding to liquid argon. For large times, analytical calculations based on the fluctuation theory provides an explicit expression reproducing isothermal change of the self-diffusion coefficient in liquid argon corresponding to the experimental data. These results lead to the conclusion that such behavior is based on the reduced mobility of particles reflected in their density fluctuations that can be equivalently achieved in the cases of either low temperatures and pressures (supercooling) or moderate temperatures and high pressures (overdamping).

Figures

Tables

Table 1. Table 1: Parameters of the fitted thermodynamic quantities

	${\tilde{F}}_3$	${\tilde{F}}_2$	${\tilde{F}}_1$	${\tilde{F}}_0$	$\tilde{T}$	${\sigma }_T$
$\rho$	-0.05600	-0.63333	-34.387	1362.8	92.5	5.40
$\ln \left({\kappa }_T\right)$	$1.613\times {10}^{-4}$	$5.721\times {10}^{-2}$	0.1611	-19.82	92.5	5.40

Equations

$${\nu }^{-1}=\frac{⟨{\left(\Delta \rho \right)}^2⟩}{\rho }/{\left[\frac{⟨{\left(\Delta \rho \right)}^2⟩}{\rho }\right]}_{i.g.}=\frac{{\mu }_0}{RT}\rho {\kappa }_T,$$

$$\rho ={\rho }_0+\frac{1}{k}\ln \left[\frac{k{\mu }_0}{\nu RT}(P-P_0)+1\right],$$

$$\dot{{\mathbf{r}}_j}$$

$$\dot{{\mathbf{v}}_j}$$

$$U(r)=4ϵ\left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right],$$

$$\overline{\delta f^2}(\tau =mdt)=$$

$${\left(f_y(t_i+mdt)-f_y(t_i)\right)}^2+$$

$${\left(f_z(t_i+mdt)-f_z(t_i)\right)}^2],$$

$${\dot{\mathbf{r}}}_j$$

$${\dot{\mathbf{v}}}_j$$

$${\mathbf{F}}_j(\omega )=\frac{1}{3t_{obs}}{\int }_0^{t_{obs}}{\dot{\mathbf{v}}}_j{e}^{i\omega t}dt.$$

$${\mathbf{v}}_j(t)=\frac{1}{m_0}{\int }_0^t\xi (t)dt=\frac{t_{obs}}{2\pi }{\int }_{-\infty }^{+\infty }{\mathbf{F}}_j(\omega )e^{-i\omega t}d\omega ,$$

$$\left({\rho }^{1/3}\sqrt{{\mu }_0/RT)}\right)D\propto e^{-A\frac{S_{ex}}{N{k}_B}},$$

$$\frac{S_{ex}}{N{k}_B}=-\underset{0}{\overset{\rho }{\int }}\left[T{\left(\frac{\partial Z}{\partial T}\right)}_V+Z(\rho )-1\right]\frac{d\rho }{\rho }.$$

$$P_i=T{\left(\frac{\partial P}{\partial T}\right)}_V-P,$$

$$\frac{\Delta S_{ex}^{\prime }}{k_B}=-\underset{{\rho }_0}{\overset{\rho }{\int }}\left[\frac{{\mu }_0{P}_i}{\rho RT}+\frac{{\mu }_{0}P(\rho )}{\rho RT}-1\right]\frac{d\rho }{\rho }$$

$$\frac{\Delta S_{ex}^{\prime }}{k_B}=$$

$$+{\nu e^{-k{\rho }_s}\left(\mathrm{Ei}(k\rho )-\frac{e^{k\rho }}{k\rho }\right)]|}_{{\rho }_0}^{\rho },$$

$$w(N,\theta )=a(N)\exp \left[N\theta -\kappa (\theta )\right].$$

$$\kappa (\theta )=\ln \mathbf{E}\left(e^{\theta N}\right)=\sum _{j=1}^{\infty }{M}_j\frac{{\theta }^j}{j!},$$

$$\bar{N}=\frac{d\kappa }{d\theta }$$

$$\mathrm{var}(N)=\overline{(\Delta N{)}^2}=\frac{d^2\kappa }{d{\theta }^2}=\frac{d\bar{N}}{d\theta }.$$

$$\frac{⟨(\Delta N{)}^2⟩}{N}=\frac{N}{V}{k}_{B}T{\kappa }_T,$$

$$\frac{⟨(\Delta N^2)⟩}{N}=\frac{RT}{{\mu }_0}\rho {\kappa }_T$$

$$\overline{(\Delta n{)}^2}=({\nu }_{0}V{)}^{-1}\bar{n}{e}^{-km\bar{n}},$$

$$\frac{d\bar{n}}{d\theta }=({\nu }_{0}V{)}^{-1}\bar{n}{e}^{-km\bar{n}},$$

$$\theta ={\nu }_{0}V\mathrm{Ei}(km\bar{n})+c_1.$$

$$\frac{d\kappa }{d\theta }=\frac{d\kappa }{d\bar{n}}\frac{d\bar{n}}{d\theta }=\frac{d\kappa }{d\bar{n}}({\nu }_{0}V{)}^{-1}\bar{n}{e}^{-km\bar{n}}=\bar{n},$$

$$\frac{d\kappa }{d\bar{n}}={\nu }_0\bar{n}{e}^{k\bar{n}}$$

$$\kappa (\bar{n})={\nu }_{0}V(km{)}^{-1}{e}^{km\bar{n}}+c_2.$$

$$\kappa (\bar{n})={\nu }_{0}V(km{)}^{-1}{e}^{km\bar{n}}=(km{)}^{-1}{e}^{km\bar{n}+b}=\frac{d\nu }{d\bar{n}}.$$

$$w(n,\bar{n})=a(n)\exp \left[n{\nu }_{0}V\mathrm{Ei}(km\bar{n})-{\nu }_{0}V(km{)}^{-1}{e}^{km\bar{n}}\right].$$

$$w(N,\rho )=a(N)\exp \left[N\left({\nu }_0\mathrm{Ei}(k\rho )-{\nu }_0(k\rho {)}^{-1}{e}^{k\rho }\right)\right],$$

Full text

I Introduction

In principle, numerical simulations of simple liquids by the method of molecular dynamics (MD) is a quite old topic, for example, quite accurate estimations of details of molecular motion and the self-diffusion coefficient in comparison with the actual experimental data for liquid argon are dated back to the seminal work by A. Rahman Rahman (1964); a review of the further developing of such simulations and their results can be found in Ref. Meier et al. (2004). However, their majority is limited by the vicinity of the liquid-vapor coexistence curve and moderately high supercritical pressures at low temperatures as well as time intervals corresponding to the stable normal diffusive regime. This fact may be conditioned by a limited number of experimental data available for simulations testing and discussion, especially liquified noble gases Barton and Speedy (1970); Suárez-Iglesias et al. (2015) as well as techniques for their obtaining. Thus, recent model studies utilizing MD simulations are shifted either to the region of supercritical and supercooled fluid states, where specific structure conditions lead to a variety of anomalous effects reflected in the transport coefficients even for relatively large characteristic times Baidakov and Kozlova (2010); Ediger and Harrowell (2012); Costigliola et al. (2016); Ohtori et al. (2017) or to a case of mixtures Thirumalai et al. (1989); Kob and Andersen (1995), where different local spatial scales of interacting particles lead to similar effects.

At the same time, the Lennard-Jones simple liquid can play the role of a model system even for discussing anomalous diffusion in more complex media mimicking the problems, which arise in biophysical systems Jeon et al. (2016); Ghosh et al. (2016). In this case, extremely short-time range dynamics may be crucial since it strongly depends on the microscopic surrounding of a moving particle leading to different effects, such as particle population splitting, non-ergodicity, etc. Cherstvy and Metzler (2013); Schulz et al. (2013, 2014). Due to the existence of inter-particle interactions, such walking processes and structural features should be taken into account not only in the context of trapping but also reaction binding, see e.g. Grebenkov et al. (2018) including such a hot topic as the direct MD simulation of forming mesoscopic objects like viruses Tarasova and Nerukh (2018). It also should be pointed out that dynamical features of particles motion in liquids are directly connected with their local microscopic structure of liquids, which determines rather complex transition from short- to long-times scales, as is has been revealed by considering a model system of hard spheres in Hopkins et al. (2010).

But in contrast to the studies mentioned above, which were primarily addressed to simplified model systems, this work is intended to consider an interplay of structural, fluctuational and diffusional properties of simple liquid via a case study of liquid argon in the range of a condition resembling the experimentally accessible states as close as possible to the latter. It is focused on the relatively unexplored in details of the short-time range of processes in this pure liquid under high pressure and low temperatures, where an influence of liquid’s microscopic structure is sufficient.

Respectively, changes in microstructure should be unavoidably reflected in as thermodynamic quantity such as the excess entropy, which determines the behavior of the macroscopic coefficient of self-diffusion Dyre (2018). Thus, the second part of this work deals with the large-time counterpart of the problem: an analytic predictive calculation of the self-diffusion coefficient in liquid argon for the same pressure-density-temperature conditions basing on an interplay between the density fluctuations and measurable thermodynamic quantities, the density and the isothermal compressibility. To assure physical relevance, the actual experimental information on thermodynamic and transport properties known in reference literature and databases are used for the direct comparison of calculated and measured data.

II Reduced density fluctuations and self-diffusion

II.1 Reduced fluctuations and self-diffusion along an isotherm and the saturation curve

To analyze the pressure- and density-dependent behavior of the self-diffusion in the considered LJ-fluid simulating liquid argon, it is worthy to address the relative density fluctuations (a ratio of the actual density fluctuations to the density fluctuations in a hypothetical medium with the same values of thermodynamic parameters but in the state of the ideal gas):

[TABLE]

which allows for analysis of an interplay between microscopic liquid structure and thermodynamics Goncharov et al. (2013) since $\kappa_{T}$, $\mu_{0}$, $R$ are the isothermal compressibility, the molar mass, and the gas constant, respectively.

It is known that the logarithm of this parameter has a practically linear dependence on the density Goncharov et al. (2013); Chora̧żewski et al. (2017), $\mathrm{ln}(\nu)=k\rho+b$, where $k$ and $b$ are substance-dependent parameters, the same for both the liquid-vapor coexisting curve and in the single-phase region (in the latter case this universality is violated under pressures higher than several hundred MPa only). Figure 1 illustrates this fact showing the logarithm of the reduced density fluctuations (1) as a function of the density simultaneously calculated along the saturation curve up to the triple point and along the isotherm $T=90~{}\mathrm{K}$ up to the vicinity of the freezing pressure at this temperature. The experiment-based database NIS was used as a source of thermodynamic data, i.e. they are completely independent of simulations.

At first, one can see that all the dots are placed along one practically straight line. At second, the sequence of blue circles corresponding to the simultaneous change of the density and the temperature along the saturation curve and the black squares corresponding to the pressure change along the isotherm overlap up to the density corresponding to the triple point, i.e. to the normal freezing density. But further, the set of squares continues this sequence along the same straight line calculated using the saturated data. Thus, from the point of view of such continuation, the isothermal fluctuations at high pressures should correspond to the fluctuations in metastable liquid argon under the normal saturated vapor, when argon it is accurately kept in a fluid non-frozen state, i.e. it represents a supercooled liquid.

Such behavior makes it possible to derive the so-called Fluctuation Theory-based Equation of State (FT-EoS)

[TABLE]

which is based on the mentioned universality that makes available predicting the density of liquids, from simple to quite complex substances, see Chora̧żewski et al. (2017) and references therein, using the data measured at normal conditions only. This possibility originates from the more physically relevant picture of elastic properties of liquids in comparison, say, with the empiric Tait equation widely used for a pure fitting approach, as discussed in Chora̧żewski et al. (2017).

Figure 4(A) shows the predicted density of argon along the isotherm $T=90~{}\mathrm{K}$ from the liquid-vapor coexistence curve to the vicinity of the freezing point ($25~{}\mathrm{MPa}$) at this temperature; the average absolute deviation between them is equal to 0.0065 %. Experimental data for the density under pressure as well as the saturated density, speed of sound and the heat capacity ratio used for computing the isothermal compressibility and, subsequently, $\nu$ and $k$ (see the description of the algorithm in Ref. Chora̧żewski et al. (2017)) are taken from the database of the National Institute of Standards and Technology (NIST) NIS , which is based on the high-accurate equation of state fitting a wide compilation of different experimental data Tegeler et al. (1999). The numerical values of parameters used to obtain plots and the fitting procedure are given in the Appendix.

II.2 Molecular dynamics simulations

Since density fluctuations in continuous media should be connected with the process of self-diffusion, consider some simulations related to the points of thermodynamics state discussed above.

The free available Python code git was used for simulations. It realizes the molecular dynamics algorithm considering Newton’s law dynamics

[TABLE]

for $60912$ particles interacting by forces $\mathbf{f}=-\nabla U$ determined via the Lennard-Jones pair-wise potential

[TABLE]

with the parameters $\epsilon/k_{B}=120~{}\mathrm{K}$, $\sigma=3.4~{}\mathrm{\AA}$ ($k_{B}$ is Boltzmann’s constant) that corresponds to the classic Rahman’s system Rahman (1964) but extended in size ($864\times 8$). Periodic boundary conditions were applied to the box of the side $L_{box}$, which was fixed and chosen in such a way that the mean density of particles within the box corresponded to the actual density of liquid argon at 90 K, i.e. $L_{box}=10.229\sigma$ respectively to the density $\rho_{sat}=1378.6~{}\mathrm{kg\,m^{-3}}$ on the saturation curve, and $L_{box}=10.0607\sigma$ respectively to the density $\rho_{20}=1432.1~{}\mathrm{kg\,m^{-3}}$ at 20 MPa on the isotherm. Note also that due to a limited number of experimental data on the self-diffusion coefficient of liquid argon, especially, under pressure for temperatures from a “normal range”, i.e. sufficiently below the critical point Barton and Speedy (1970); Suárez-Iglesias et al. (2015), the isotherm 90 K provides a most data-rich choice in a vicinity of the boiling point. In addition, the simulations were evaluated for the same density $\rho_{20}$ but the temperature $T=81.27~{}\mathrm{K}$, which was obtained as corresponding to this density in the supercooled (liquid) state along the saturation curve derived by extrapolating $\rho(T)$ dependence along this curve.

The potential’s cut-off was chosen as $r_{cut}=4.5\sigma$ that is twice of the conventional choice $r_{cut}=2.25\sigma$ to assure possible long-ranged interaction in a high-density compressed liquid: the chosen value corresponds to 8 van der Waals radii $r_{W}=1.88~{}\mathrm{\AA}$ of the argon atom, i.e. to the fourth coordination sphere (although comparative studies showed that the results do not differ significantly, i.e. nearest-neighbor interactions between atoms prevails, as expected). The time step of simulations was equal to $dt=10^{-14}~{}\mathrm{s}$.

Before measurements, the system was equilibrating up to 470 ps; the equilibration was controlled by the plots of the pressure and the temperature of the system. The co-ordinates of all atoms were recorded for the subsequent 4096 time iterations (40.96 ps). To exclude effects of periodization, all trajectories were corrected via the shift of co-ordinates on $L_{box}$ when a particle reached the box’s side and the respective co-ordinate jumps occurred.

The time-averaged mean square displacement (tMSD) for each trajectory was calculated for $m$ elementary time steps $dt$ as

[TABLE]

where $f_{x,y,z}$ are components either of the displacement $\mathbf{r}$ (for the true MSD) or the velocity $\mathbf{v}=(\mathbf{r}(t_{i}+mdt)-\mathbf{r}(t_{i}))(mdt)^{-1}$ (in this case, the mean square velocity, MSV, is considered).

Note that all trajectories, which crossed the boundary and were continued under periodicity condition, were shifted and “glued” in these points of jumps in such a way that they went out on the initial box keeping the continuity of trajectories and velocities. The ensemble-averaged eMSD and eMSV were calculated for all cases indicating that they coincide with tMSD (the walks are ergodic), thus further the notations MSD and MSV will be used.

Figure 2(A) shows the MSD for three considered conditions, where the length variable is rescaled to the characteristic length included into the L-J potential; note also that the value $2^{1/6}\sigma\approx 2r_{W}$, i.e. to the characteristic inter-particle distance ($r_{W}$ is the van der Waals radius of argon). The time variable is kept dimensional for a comparison with times, which are accessible in real physical experiments for studying the self-diffusion. One can see that for short times the particle’s motion is ballistic that is highlighted by the green dashed parabola. The respective displacements are quite short, of order $0.15~{}\sigma$ that correspond to motions just in a small vicinity of the potential’s minimum, where forces acting on the particle is practically negligible. During the next time interval of around 1-2 ps, the MSD exhibits more sophisticated behavior corresponding to the crossover from quadratic to linear time dependence. The latter function corresponds to the normal diffusion and is detected for $\tau>2~{}\mathrm{ps}$. Different slopes for three different thermodynamic states indicate different values of the self-diffusion coefficient. To check the consistency of simulations, note that the diffusive regime for the saturated argon at $T=90~{}\mathrm{K}$ is found as $D_{sat}=2.3\,10^{-9}~{}\mathrm{m^{2}\,t^{-1}}$ that corresponds to the value obtained by other authors for the L-J system mimicking liquid argon at the same temperature Wei-Zhong et al. (2008) and belongs to the middle of the interval of its experimental values determined in different sources: $D_{sat}^{(exp)}=(1.89\pm 0.08)\,10^{-9}~{}\mathrm{m^{2}\,t^{-1}}$, $D_{sat}^{(exp)}=(2.10\pm 0.10)\,10^{-9}~{}\mathrm{m^{2}\,t^{-1}}$ Cini-Castagnoli and Ricci (1960), $D_{sat}^{(exp)}=2.43\,10^{-9}~{}\mathrm{m^{2}\,t^{-1}}$ Naghizadeh and Rice (1962). More compressed liquid at the same temperature shows smaller self-diffusion, and the least one is the self-diffusion at the least temperature but the same density as the previous one. This is quite expectable, but one can draw more interesting conclusions about details of molecular motions using MSVD shown in Fig. 2(B),(C).

First of all, Figure 2(B) indicates that the crossover from the ballistic to the diffusive motion in Figure 2(A) corresponds to the sufficiently non-monotonic behavior of MSVD in Figure 2(B). The boundary regimes are a ballistic parabola and a constant line, and there is a maximum in between. Physically it may be interpreted as a fast collision-reflection-permeation process when walking particles interact with particles, which comprise the first co-ordination shell. The value of this maximum is largest for the overdamped liquid at $T=90~{}\mathrm{K}$, i.e. particles move fast but they are sufficiently squeezed. The same squeezing for $T=81.27~{}\mathrm{K}$ results in a smaller “hump”, and the smallest one corresponds to the saturated case of fast-moving particles, for which a larger free volume is available. All cases result in the constant MSVD, i.e. in the thermalization of motion, when the dispersion of velocities is fixed and the random walk is Brownian.

All three MSVD curves are distinct after leaving the regime of the practically force-less ballistic flights in the vicinity of the potential well minimum. However, it is possible to rescale them to reveal a kind of universality, which reflects the basic physical premises of random displacements. First of all, magnitudes of the velocity should be determined by the system temperature as follows from the basics of statistical thermodynamics. Secondly, a particle gets a different velocity while it moves within an available free volume before the thermalization resulting in the diffusional motion. The respective statistics can be characterized by the reduced density fluctuation parameter (1), which is also equal to the reduced volume fluctuations. Since MSDV relates to the one-dimensional length units squared, the resulting dimensionless scaling factor takes the form $\varepsilon/k_{B}T\nu^{1/3}$, where $\varepsilon/k_{B}$ defines the characteristic L-J temperature.

Figure 2(C) shows that such scaling already results in merging TMSD curves and this picture is different from the one depicted in Figure 2(B). Especially demonstrative is the quite accurate coincidence of blue dashed and red dash-dotted curves, which correspond to velocities in the supercooled and overdamped liquids that confirms fluctuation theory-based conclusions discussed above in relation to Figure 1. Both the ballistic (highlighted by the green dotted parabola) and the diffusional regimes are characterized by the same universal curve for all three cases.

Such behavior of MSD and MSVD visible in Figure 2 for an individual trajectory can be discussed in terms of random walks when the deterministic process defined by Newton’s equations of motion (3)–(4) is replaced with the system

[TABLE]

where $\xi(t)$ is an appropriate random process.

Note that such representation differs from a conventional discussion of random walks with the transition from ballistic to diffusional MSD introducing the underdamped Langevin equation with a Gaussian (white or colored) random noise and velocity-dependent damping Bodrova et al. (2016); Fa (2018). The mentioned approach postulates a kind of random force and introduces the viscous damping from some macroscopic manifestations adjusted to the chosen noise character, while Eqs. (6)–(7) play a role of a stochastic counterpart to dynamical Eqs. (3)–(4), where the stochasticity originates from non-uniform distribution of nearest-neighbor particles acting on the observed one. This consideration provides a more direct interpretation for the behavior of MSD and MSDV shown in Fig 2, where the noise term in Eqs. (7) can be directly characterized by the recorded ensemble of trajectories. The most informative characteristic of such random forces (random accelerations) is its power spectral density (PSD) $\langle|\mathbf{F}_{j}(\omega)|^{2}\rangle$ over the period of observations $t_{obs}$ averaged over the ensemble of trajectories, where

[TABLE]

Figure 3 shows the averaged power spectral densities for three considered thermodynamic states (8) normalized to the unit area under each curve. All of them have specific qualitative shape features, which explains the behavior of MSVD and MSD shown in Figure 2.

The PSD curve has a clear maximum that indicates the existence of a leading oscillatory term; its period is equal to $2\pi/\omega_{max}=1.38~{}\mathrm{ps}$ for $\rho_{90}$ and $1.27~{}\mathrm{ps}$ for $\rho_{90}$, respectively. Note that these period values are coordinated approximately with the end of the transient process Figs. 2(B)–(C), i.e. larger time (and, respectively, spatial) scales are thermodynamic. In the parts of the PSD for the frequencies $\omega/2\pi<0.5~{}\mathrm{ps}^{-1}$, i.e. to times larger that $2\,\mathrm{ps}$, blue dashed and red dash-dotted curves in Figure 3 practically coincide (and they are quite close everywhere being distinct from the black solid line), that has correspondence to the equivalence of the slow density fluctuation processes under thermodynamic spatial and temporal scales as it is reflected in Figure 1. Finally, the PSD clearly tends to zero as $\omega\to 0$. Thus, the most probable particle’s motion has an isotropic oscillatory character within cages of radii corresponding to average interparticle distance but existing density fluctuations leads to a wide distribution with respect to frequencies instead of a localized peak.

This distribution, however, has a fast (exponential) decay (see the green guiding line in Figure 8). Practically, there are no frequency components for $\omega/2\pi>3.5~{}\mathrm{ps}^{-1}$, i.e. oscillations with periods less than $\tau_{bound}=0.28~{}\mathrm{ps}$. Looking at Figs. 2(B)–(C), one can see that it is accurately the right boundary of ballistic behavior for the mean square velocity displacements, i.e. a practical absence of oscillatory components with higher frequencies results in the principal absence of disturbances leading to deviations from the linear dependence of MSDV on time for the walking process with zero mean

[TABLE]

If one replaces the limits of the integral above by the values defined by $\tau_{bound}$, such velocity displacements will resemble so-called “smooth random process” extensively discussed recently in Filip et al. (2019) and in references therein. Such interpretation is well-fitted with the ideology of molecular dynamics since Eqs. (3)–(4) are dynamic, i.e. resulting trajectories are differentiable and their irregularity is based on the irregularity of surrounding particles acting on the observed particle. This is actually the case of Eqs. (7), (9). In turn, the integration of Eq. (9) as follows from Eq. (6) results in an irregular path. Taking into account times larger than $\tau_{bond}$, the velocities $\mathbf{v}_{j}(t)$ can be considered as effective random variables and their integration (a sequential summation) leads to the Brownian motion due to the Central Limit Theorem (there are no power-law tails here) with the MSD seen in Figure 2(A).

II.3 Self-diffusion under high pressures via the thermodynamic fluctuations route

Now it is possible to analyze the possibility of extending this fluctuation-based method from volumetric to transport properties, in particular, to self-diffusion in the thermodynamic limit of scales whose boundaries are revealed in the results of MD simulations above. The dimensionless reduced coefficient of self-diffusion in liquids in a majority of cases exhibit a universal exponential behavior Rosenfeld (1977); Dzugutov (1996); Rosenfeld (1999); Dyre (2018)

[TABLE]

where $S_{ex}/Nk_{B}$ is the excess entropy per particle, which defines the difference between the entropy of the system under study and the entropy of an equivalent ideal gas at the same temperature and density, i.e. refers to the same pair of states as Eq. (1); $A$ is a positive substance-dependent constant.

It should be pointed out that the excess entropy is not a quantity, which can be determined straightforwardly in an experiment and its value depends on fitting to a chosen model Dyre (2018). In particular, an approach relatively easily applicable to results of numerical simulations replaces the full excess entropy by the pair entropy Baranyai and Evans (1989); Saija and Giaquinta (1996), which is, in fact, the leading term in the expansion of the complete thermodynamic function into the Taylor series with respect to contributions of particles pairs, triplets, etc.

But the simplest possibility, valid not only in the low-density limit, is a usage of purely thermodynamic quantities, namely the compressibility factor $Z=\mu_{0}P/\rho RT$ and its isochoric derivative Vaz et al. (2012)

[TABLE]

Using the standard thermodynamic definition of the internal pressure

[TABLE]

the derivative mentioned above can be easily calculated and Eq. (11) rewritten as

[TABLE]

if one considers a difference between excess entropies of two states with the densities $\rho_{0}$ and $\rho$ placed on one isotherm. After substituting into the integrand the pressure expressed from FT-EoS (2) as an explicit function of the density, the integral (12) is taken analytically that results in

[TABLE]

where $\mathrm{Ei}(k\rho)$ is the exponential integral, and $P_{i}=\mathrm{const}$ as it follows from the basic construction of FT-EoS (2); although it fulfils only approximately in average, the actual changes of $P_{i}$ with respect to its value within the considered range of pressures can be neglected. Accuracy of the density prediction shown in Figure 4(A) confirms this.

Figure 4(B) shows the curve of the relative coefficient of self-diffusion change calculated using Eqs. (13), (10) with $A=3.2$ in comparison with the raw experimental ratio of the coefficient of self-diffusion under pressure to its saturated value taken from Naghizadeh and Rice (1962). One can see that the proposed density fluctuation-based model quite accurately reproduces the non-linear character of the self-diffusion coefficient’s diminishing with the growing pressure.

II.4 Particle number distribution as an exponential dispersion model

Now let us consider this model from the point of view of the probability density functions (p.d.f.) for particles located within a small (but not too small, see below) volume surrounding a chosen particle. This analysis can be done following the approach called an exponential dispersion model (EDM) Jørgensen (1997) describing statistical distributions for which the variance of a random variable is a function of its mean value.

Such property is fulfilled, when the p.d.f. has a form

[TABLE]

Here it is written with respect to the number of particles $N$; $\kappa(\theta)$ is the cumulant function defined in a standard way as the natural logarithm of the moment-generating function (moments are denoted as $M_{j}$ here)

[TABLE]

i.e. coefficients of its expansion into the Taylor series are the moments of distribution; the parameter $\theta$, with respect to which this expansion is written, is called the the canonical parameter; and $a(N)$ is some suitable function, which assures the norm of this distribution.

Following this definition, the mean number of particles within the chosen volume and its variance are determined as

[TABLE]

and

[TABLE]

Consider the volume as small but sufficiently macroscopic, i.e. one can express the number fluctuations in a fixed volume via the standard statistical thermodynamics expression Landau and M. (2013)

[TABLE]

where $k_{B}$ is Boltzmann’s constant. Since we operate with the macroscopic volume, for which one can determine the isothermal compressibility, the medium can be considered as uniform at such scales, i.e. $N/V\cong\langle N\rangle/V\equiv\bar{n}=\rho/m$, where $m$ is the mass of one particle. From the microscopical point of view, this condition is fulfilled when the volume’s radius exceeds 3-4 mean inter-particle distances, when the radial distribution function (r.d.f.) approaches a horizontal line and, respectively, the number of particles inside a selected volume will grow with the growth of the latter linearly as proportional to the thermodynamic density.

Under this assumption, Eq. (17) takes the form

[TABLE]

coinciding with Eq. (1). For ideal gas the right-hand side is equal to 1, otherwise it has a form of the exponential function $\exp\left[-(k\rho+b)\right]=\exp\left[-(km\bar{n}+b)\right]$. Respectively, multiplying the nominator and denominator in the left-hand side by the fixed volume squared $V^{2}$ and applying the same assumption that actual number density is equal to the mean number density, it is possible to conclude that the system satisfies the conditions of the EDM:

[TABLE]

where the notation $\nu_{0}=\exp(kb)$ is introduced.

Therefore, the mean value satisfies the ordinary differential equation

[TABLE]

which can be easily solved by the variable separation method:

[TABLE]

Within the same way, it is possible to find the cumulant function using Eq. (15):

[TABLE]

i.e.

[TABLE]

with the solution

[TABLE]

Note that $c_{1}$ and $c_{2}$ are additive constants, and a substitution of the expressions (19), (20) into Eq. (14) will change the latter by a constant factor $\exp(c_{2})$ and the still indefinite function $a(n)$ is the $n$-only dependent multiplier. Therefore, one can put both $c_{1}=0$ and $c_{2}=0$ without loss of generality.

Thus, the conclusion is that the cumulant function in this case is the derivative of the reduced bulk modulus with respect to the density: with the solution

[TABLE]

The resulting p.d.f. has a form

[TABLE]

Note that it can be rewritten back to the number of particles $N=nV$ and the thermodynamic mass density $\rho=m\bar{n}=mn$, in the form independent on a particular value of the volume $V$ (keeping it larger then the region of r.d.f.’s oscillations)

[TABLE]

i.e. the exponential term of the EDM’s p.d.f. can be associated with a part in the thermodynamic entropic contribution into the self-diffusion (13), which contains the same combination of the exponential integral and the exponential function.

However, it can be shown that Eq. (21) has a more wide applicability and can be used not only in the case of high (liquid) densities that can be demonstrated via a transition to the limiting case of extremely small number densities, i.e., $|k|\rho<<1$. This may be achieved considering the vapor branch at almost ideal gas conditions, when $\nu_{0}\approx 1$, and $k<0$, see Goncharov et al. (2013), where the inverse reduced density fluctuations under such conditions were studied based on actual experimental data for gaseous argon.

In the considered case $\exp(k\rho)\approx 1+k\rho=1+km\bar{n}$, and the cumulant function

[TABLE]

since $\bar{n}=\bar{N}/V$.

The series representation of the exponential integral function in this case is

[TABLE]

where $\gamma$ is Euler’s constant, and $\left[\mathrm{ln}(km/V)-i\pi\right]=\mathrm{Re}\mathrm{ln}(km/V)$ is a real number.

Where,

[TABLE]

The exponential of the first term is equal to

[TABLE]

the last term can be neglected as small, the second term does not depend on $\bar{N}$ and, therefore, can be combined with $a(N)$. Thus, the resulting p.d.f. is

[TABLE]

Taking the first factor as $1/N!$ [it is possible since $a(N)$ is a still indefinite function and it is required to norm the p.d.f. to unity], this expression reduces to the Poisson distribution

[TABLE]

that is known for the ideal gas (or in the case of small volumes containing an extremely small number of particles inside) v Smoluchowski (1904); Landau and M. (2013).

Thus, the revealed non-trivial dependence on the liquid density, which can be connected with average inter-particle distances, may serve for a future microscopic probabilistic background for the studying transport processes in liquids under pressure, which exhibit a non-trivial behavior.

III Conclusion

The main results of this work can be summarized as follows. Although the self-diffusion in simple liquids in the thermodynamic limit behaves like a normal diffusion process, its detailed picture for short times is more complex and exhibits properties of a random mixture of oscillatory modes whose frequencies are defined by free volume available for particles caged within their first coordination shells. The properties of such oscillations can be highlighted in the most direct way by the exploration of the mean square velocity displacement and connected with the thermodynamic density fluctuations. Note also that such behavior can be associated with the recent topic of non-stationary transient modes in model Ornstein-Uhlenbeck and related stochastic processes Thiel et al. (2016); Cherstvy et al. (2018); Kuśmierz et al. (2018). In addition, the approach considered in this work is not limited by the particular L-J potential only, and may be applied to systems with other potentials, e.g. in colloidal systems Pellicane (2012). This is because it addresses the density (or local volume) fluctuations, and, respectively, the liquid’s structure factor, which can be calculated as a statistical quantity for liquid with various kinds of interactions.

Finally, it should be pointed out that some kind of statistical equivalence of liquids behavior under elevated pressures or lowered temperatures can be described in a uniform way within a frame of the generalized linear model. This follows from a response of the density fluctuations on the growing particles packing that lead to the same value of the reduced density fluctuations. As a result, an analytic predictive expression for calculating the density and the self-diffusion coefficient with an accuracy corresponding to the experimental one is available. Thus, this open perspective for analytical prediction of this parameter not only for argon but for more complex substances, e.g. organic molecular liquids, since the basic reference data can be extracted from thermodynamic databases.

Acknowledgment

This work is partially supported by the National Science Centre (Poland) Grant No. 2016/23/B/ST8/02968.

Appendix

The density and the natural logarithm of the isothermal compressibility are fitted by polynomials using the standard MATLAB’s function polyfit pol as

[TABLE]

where $\tilde{f}$ is either $\rho$ measured in $\mathrm{kg\,m^{-3}}$ or $\mathrm{ln}(\kappa_{T})$ measured in $\mathrm{Pa}^{-1}$, the temperature is in Kelvins. The respective coefficients are given in Table. 1.

For the key dimensionless fluctuation parameter, Eq. (1), its natural logarithm was linearly fitted as $\mathrm{ln}(\nu)=k\rho+b$ that gave values of the coefficients as $k=0.00564~{}\mathrm{m^{3}\,kg^{-1}}$ and $b=-4.94$.