Density-Dependent Finite System-Size Effects in Equilibrium Molecular Dynamics Estimation of Shear Viscosity: Hydrodynamic and Configurational Study

TL;DR

This study investigates how finite system size affects the estimation of shear viscosity in molecular dynamics simulations, revealing oscillatory and scaling behaviors linked to hydrodynamic and configurational factors.

Contribution

It provides a detailed analysis of size effects on shear viscosity, distinguishing between hydrodynamic and configurational contributions through extensive simulations.

Findings

Oscillatory size effect at high density

L^{-1} scaling correction at low density

Long-time tail behavior confirms hydrodynamic origin

Abstract

We study the intrinsic nature of the finite system-size effect in estimating shear viscosity of dilute and dense fluids within the framework of the Green-Kubo approach. From extensive molecular dynamics simulations, we observe that the size effect on shear viscosity is characterized by an oscillatory behavior with respect to system size $L$ at high density and by a scaling behavior with an $L^{-1}$ correction term at low density. Analysis of the potential contribution in the shear-stress autocorrelation function reveals that the former is configurational and is attributed to the inaccurate description of the long-range spatial correlations in finite systems. Observation of the long-time inverse-power decay in the kinetic contribution confirms its hydrodynamic nature. The $L^{-1}$ correction term of shear viscosity is explained by the sensitive change in the long-time tail obtained from a finite system.

Equations

$$\eta ={\int }_0^{\infty }C(t)dt,C(t)=\frac{V}{k_{\mathrm{B}}T}⟨p_{xy}(0)p_{xy}(t)⟩.$$

$$p_{xy}=p_{xy}^{\mathrm{K}}+p_{xy}^{\mathrm{P}},$$

$$C(t)=C^{\mathrm{KK}}(t)+C^{\mathrm{PP}}(t)+2C^{\mathrm{KP}}(t).$$

$$p_{xy}=\frac{1}{V}\left[m\sum _i{v}_{i,x}{v}_{i,y}-\frac{1}{2}\sum _i\sum _{j\ne i}\frac{r_{ij,x}{r}_{ij,y}{\varphi }^{\prime }(r_{ij})}{r_{ij}}\right],$$

$$C^{\circ \star }(t)=\frac{V}{k_{\mathrm{B}}T}⟨p_{xy}^{\circ }(0)p_{xy}^{\star }(t)⟩,$$

$$\eta ={\eta }^{\mathrm{KK}}+{\eta }^{\mathrm{PP}}+2{\eta }^{\mathrm{KP}}\text{mbox},where{\eta }^{\circ \star }={\int }_0^{\infty }{C}^{\circ \star }(t)dt.$$

$$\phi(r)=\left\{\begin{array}[]{cc}4\varepsilon\left[{\left(\frac{\sigma}{r}\right)}^{12}-{\left(\frac{\sigma}{r}\right)}^{6}+\frac{1}{4}\right],&r\leq 2^{1/6}\sigma,\\ 0,&r>2^{1/6}\sigma.\end{array}\right.$$

$$t_0={\left[4n{\sigma }^2{\left(\frac{\pi k_{\mathrm{B}}T}{m}\right)}^{1/2}g(\sigma )\right]}^{-1},$$

$$C^{\mathrm{KK}}(t)=\frac{2k_{\mathrm{B}}T}{5n{L}^3}\sum _{\mathbf{k}\ne 0}\exp \left(-\frac{2k^2\eta }{n}t\right),$$

$$C^{\mathrm{KK}}(t)=\frac{2k_{\mathrm{B}}T}{5n}{\left(\frac{8\pi \eta }{n}t\right)}^{-3/2}.$$

$$\begin{split}C^{\mathrm{KK}}_{\mathrm{I}}(t)&=\frac{m^{2}}{Vk_{\mathrm{B}}T}\bigg{\langle}\sum_{i}v_{i,x}(0)v_{i,y}(0)v_{i,x}(t)v_{i,y}(t)\bigg{\rangle},\\ C^{\mathrm{KK}}_{\mathrm{IJ}}(t)&=\frac{m^{2}}{Vk_{\mathrm{B}}T}\bigg{\langle}\sum_{i}\sum_{j\neq i}v_{i,x}(0)v_{i,y}(0)v_{j,x}(t)v_{j,y}(t)\bigg{\rangle}.\end{split}$$

$$C^{\mathrm{PP}}(0)=\frac{2\pi n^2}{15}{\int }_0^{\infty }{r}^4\left[{\varphi }^{\prime \prime }(r)+\frac{4}{r}{\varphi }^{\prime }(r)\right]g(r)dr.$$

$$C_{\mathrm{I},L}^{\mathrm{KK}}(t)$$

$$C_{\mathrm{IJ},L}^{\mathrm{KK}}(t)$$

$${\eta }_L^{\mathrm{KK}}={\eta }_{\infty }^{\mathrm{KK}}+\frac{F}{L^3}+\frac{G}{L},$$

Full text

I Introduction

The shear viscosity coefficient in dense fluids, along with diffusion coefficient and thermal conductivity, constitute important ingredients of hydrodynamic theory. These transport coefficients are connected to the corresponding time-correlation functions of microscopic fluctuating variables via the Green–Kubo relations. At equilibrium, or in small deviations from equilibrium, a systematic connection between the correlation functions and the hydrodynamic equations has been established in the long-wavelength and small-frequency hydrodynamic limit HansenMcDonald2013 ; BoonYip1980 ; BalucaniZoppi1995 ; Forster1975 .

The molecular dynamics (MD) simulation technique provides several methods to evaluate the shear viscosity coefficient of a molecular fluid. First, the shear viscosity coefficient can be estimated by comparing flow patterns generated from nonequilibrium MD simulations with those predicted by hydrodynamic theory Erpenbeck1984 ; EvansMorriss1989 ; BackerLowe2005 . While these direct methods can, in principle, be used to calculate nonlinear as well as linear transport coefficients, they have an inherent arbitrariness in producing the nonequilibrium flow fields with respect to the thermostat, barostat, and other factors that influence the motion of individual particles by an external applied force. The second type of method examines the long-time behavior of the long-wavelength correlations in the equilibrium fluctuations of the transverse momentum field Palmer1994 . It is based on the observation that in the hydrodynamic limit, these correlations decay exponentially with the exponent proportional to the shear viscosity. The third approach is to compute the shear viscosity $\eta$ using the Green–Kubo integral LevesqueVerletKurkijarvi1973 ; Erpenbeck1988 ; MaginnMesserlyCarlsonRoeElliott2018 ,

[TABLE]

Here, $C(t)$ is the normalized shear-stress autocorrelation function (SACF), which is calculated from an equilibrium MD simulation of a system with volume $V$ at temperature $T$. Contrary to the other two approaches, detailed information about the dynamics of the microscopic fluctuating variable (i.e., the off-diagonal pressure tensor component $p_{xy}$) is available.

I.1 Density Dependence

Structural relaxation in dense fluids is a complex process and coupled with momentum relaxation. Since the microscopic expression of $p_{xy}$ [See Section II.1] is written as the sum of kinetic (K) and potential (P) contributions,

[TABLE]

these relaxation processes can be investigated from the resulting decomposition of the SACF into kinetic-kinetic (KK), potential-potential (PP), and kinetic-potential (KP) components,

[TABLE]

While kinetic theory has been successful in explaining the behavior of kinetic term ErnstHauge1971 ; ErpenbeckWood1982 , it has not given a complete answer for theoretical substantiation of the asymptotic long-time behavior of the SACF and prediction of density dependence of the shear viscosity. This is because the behaviors of the potential and cross terms exhibit distinctively different features from the kinetic term. In particular, in the high density region near the fluid-solid transition point, high-frequency viscoelastic modes play an important role but their influence on the shear viscosity is still not clearly understood.

As illustrated in Figure 1, the SACF at high density is dominated by the potential term, whereas the long-time behavior of the SACF at low density is governed by the kinetic term. A significant contribution by the potential term explains why the prediction of the shear viscosity using kinetic theory fails in the high density region. At low densities, a substantial increase in the kinetic contribution implies that the dynamics is no longer dominated by structural relaxation. Competition between the two contributions, in addition to that of the cross term, makes the physical nature of the density dependence of the shear viscosity difficult to understand.

Previous MD simulation results IsobeAlder2009 ; IsobeAlder2010 ; IsobeAlder2012 suggest an intimate connection between shear viscosity and its configurational nature at high density. The configurational nature includes aspects of the dependence of the shear stress on the atomic configuration. For example, the long-time tail of the SACF observed in high density region, often referred to as the “molasses tail,” has a fundamentally different physical origin from the hydrodynamic vortex flow of the velocity autocorrelation function ChoiHan2017 ; HanKimTalknerKarniadakisLee2018 and also from the kinetic fluxes related to shear viscosity ErnstHauge1971 . This molasses tail effect becomes more pronounced as fluid density approaches the fluid-solid transition point and results in a markedly enhanced shear viscosity near solidification. MD simulation results on the molasses tail reveals that it is due to the orientational correlations of bonds connecting colliding pairs of particles at long times and high densities near solidification. The existence of algebraic tails, followed by the stretched exponential, is understood as having a hydrodynamic origin which includes much shorter wavelength fluctuations LaddAlder1989 ; LeegwaterBeijeren1989 . The long-range many-body ‘excluded volume’ effects showing up in the potential component $C^{\mathrm{PP}}(t)$ at high density have been noted for the Lennard-Jones fluid in the past StassenSteele1995a ; StassenSteele1995b and an empirical equation of state of shear viscosity that is accurate in a broad range of density and temperature has been proposed Woodcock2005 .

I.2 Finite System-Size Effects

In the process of analyzing the physical meaning of the simulation results, the effects of finite simulation system size should be effectively removed. These effects are unavoidable and cause serious difficulty in the interpretation of the simulation results. The periodic boundary conditions conventionally employed in MD simulations as an effort to correct physical properties observed in finite systems, and mimic the infinite system, still differ from material properties of the bulk system due to neglect of long-wavelength fluctuations of corresponding dynamical variables.

Compared to the self-diffusion coefficient, for which there have been a significant number of MD studies on the finite system-size effects since the celebrated Yeh–Hummer system-size correction YehHummer2004 , there have been only a few studies that investigate the particular issue for the shear viscosity. Most of the recent MD studies reported a weak system-size dependence DaivisEvans1995 ; YehHummer2004 ; MoultosZhangTsimpanogiannisEconomouMaginn2016 . However, it has not been clearly confirmed whether this system-size effect exhibits a scaling law. Moreover, very little has been known about its density dependence, not to mention its physical interpretation.

Our previous work KimHan2018 showed how configurational degrees of freedom influence values of the shear viscosity as the simulation system size changes at high density. The MD simulation results for simple and complex fluids have identified the presence of system size-dependent behaviors of the shear viscosity at high densities, especially for the case of small system sizes owing to limited configurational rearrangements in finite systems Petravic2004a ; Petravic2004b ; KimHan2018 . This restriction leads to complex oscillatory behavior of the shear viscosity (in the plot of $\eta_{L}$ versus the system size $L$), which however quickly disappears for system sizes larger than a critical size.

It is, however, expected that the effect of finite system size in estimating the shear viscosity will have a density dependence if the potential and kinetic components of the SACF exhibit different system-size effects and compete with each other. The behavior of the kinetic part has been theoretically studied using the mode-coupling theory with the linearized Navier–Stokes equations, and finite system-size effects on the algebraic tail of $C^{\mathrm{KK}}(t)$ have been investigated ErnstHauge1971 ; ErpenbeckWood1982 . The simulation system-size dependence of the SACF and the shear viscosity makes it even harder to understand the nature of their density-dependent behavior. Hence, it is crucial to understand the density-dependent behavior and relative weight of each contribution to the SACF to discern the actual behavior of the SACF and shear viscosity in the thermodynamic limit.

In this paper, we investigate the behavior of the three contributors to the SACFs. Our study focuses on the shear viscosity of simple fluids interacting via the Weeks–Chandler–Andersen (WCA) potential at three typical densities as a function of system size. The SACFs obtained from equilibrium MD simulation exhibit a variety of different patterns depending on system density. We relate system-size effects with the hydrodynamic behavior (induced by the kinetic component) as well as the configurational restriction (induced by the potential component). It is noted that a comprehensive understanding of the behavior of the SACF has, until now, been a demanding task because the short- and long-time behavior of the SACF exhibit qualitative differences, and moreover, the long-time behavior sensitively depends on simulation system size. The knowledge acquired from this study is theoretically important in understanding the hydrodynamic and configurational features in collective dynamics. Our study also provides practical suggestions on how to evaluate the shear viscosity in thermodynamic limit.

The rest of the present paper is organized as follows. In Section II, we present a brief description of our MD simulation as well as the decomposition of the SACF and the shear viscosity. In Section III, we discuss the density-dependent behavior of the SACF by examining the physical origins of characteristic long-time behaviors of the kinetic and potential components. In Section IV, we investigate two types of finite system-size effect on the shear viscosity coefficient, which have hydrodynamic and configurational origins. Our conclusions are presented in Section V.

II Method

II.1 Decomposition of SACF

For a simple fluid system under periodic boundary conditions, the molecular expression of the off-diagonal pressure tensor component $p_{xy}$ is written as ThompsonPlimpton2009

[TABLE]

where $\phi(r)$ is the interaction pair potential, $m$ is the mass of a fluid particle, $v_{i,\alpha}$ is the $\alpha$-component of the velocity vector of the $i$th particle, and $r_{ij}$ and $r_{ij,\alpha}$ are the interatomic distance and the $\alpha$-component of the displacement vector between particles $i$ and $j$, respectively. The first and second terms in the square brackets correspond to the kinetic and potential components (i.e., $p_{xy}^{\mathrm{K}}$ and $p_{xy}^{\mathrm{P}}$), respectively, see Eq. (2).

The three contributors to the SACF $C(t)$, denoted by $C^{\mathrm{KK}}(t)$, $C^{\mathrm{PP}}(t)$, and $C^{\mathrm{KP}}(t)$ (see Eqs. (1) and (3)), are correspondingly defined as

[TABLE]

where $\circ$ and $\star$ are either K or P. Note that the time-correlation functions are multiplied by a prefactor of $V/{k_{\textnormal{B}}T}$ to normalize with respect to the inherent system-size dependence of the microscopic stress tensor. Then we represent the shear viscosity $\eta$ by the following decomposition:

[TABLE]

II.2 MD Simulation

As a microscopic model of a simple fluid, we consider fluid particles moving in a three-dimensional cubic domain of side length $L$ with periodic boundary conditions. The interaction between particles is given by the WCA potential:

[TABLE]

We use reduced units of mass, length, and energy, i.e., $m=\sigma=\varepsilon=1$, along with the Boltzmann constant $k_{\textnormal{B}}=1$. The cutoff radius of the potential is set as $r_{\mathrm{c}}=2^{1/6}$. MD simulations of $N=128,\ldots{},8192$ particles were performed at number densities $n=0.45$, 0.65, and 0.85 and temperature $T=1$. Under this ambient temperature condition, the three number densities describe typical fluid states; the high density $n=0.85$ corresponds to a dense fluid near the fluid-solid transition point IsobeAlder2010 , whereas the low density $n=0.45$ corresponds to a fluid with strong hydrodynamic character. The simulation box size $L$ is accordingly determined as $L=(N/n)^{1/3}$. Simulations were conducted using the velocity Verlet algorithm implemented in LAMMPS Plimpton1995 with timestep size $\Delta t=0.002$. Each equilibrium sample was obtained through equilibration period $\mathcal{T}_{\mathrm{equil}}=10^{5}\Delta t=200$ with subsequent production run for period $\mathcal{T}=10^{6}\Delta t=2000$. For each set of simulation parameters, a total of $\mathcal{N}=524288$ sample trajectories were calculated.

We computed the SACF $C(t)$ and the shear viscosity $\eta$ as well as their components $C^{\circ\star}(t)$ and $\eta^{\circ\star}$ as follows. From each sample trajectory, instantaneous values of $p_{xy}^{\mathrm{K}}$ and $p_{xy}^{\mathrm{P}}$ (see Eq. (4)) were collected at every five timesteps to calculate $C^{\circ\star}(t)$ using a standard time-averaging procedure KimBorodin2015 . The results at different densities are shown in Figures 1 and 2, which will be discussed in Section III. We use a scaled time $\tau=t/t_{0}$ with $t_{0}$ being the average collision time

[TABLE]

where $\sigma$ is defined as the position of the first maximum of the radial distribution function $g(r)$. Numerical time-integration of the correlation functions, $C(t)$ and $C^{\circ\star}(t)$, was conducted using the trapezoidal rule to obtain the time interals up to time $t$, e.g., $\eta(t)=\int_{0}^{t}C(t^{\prime})dt^{\prime}$. For each set of simulation parameters, the $t^{*}$ value that satisfies $\eta(t^{*})\approx\eta=\lim_{t\rightarrow\infty}\eta(t)$ was chosen a posteriori from the time profile of $\eta(t)$ and was used to determine the values of $\eta$ and $\eta^{\star\circ}$. For each physical quantity, we computed its sample mean and standard deviation over $\mathcal{N}=524288$ samples.

III Shear-stress Autocorrelation Function

We begin with an overall description of the time profiles of the SACF and its components at low, intermediate, and high densities displayed in Figure 2. The short-time behavior of the SACF is dominated by the potential term regardless of number density, whereas its long-time behavior is influenced by both kinetic and potential terms resulting in distinctive density-dependent decay patterns. The slow algebraic decays observed at low and intermediate densities are mainly due to the kinetic correlation $C^{\mathrm{KK}}(\tau)$, whereas the potential correlation $C^{\mathrm{PP}}(\tau)$ is responsible for the faster non-algebraic decay at high density. The cross correlation $C^{\mathrm{KP}}(\tau)$ is as influential as the other two components at intermediate timescales, but displays negligible impact on the asymptotic long-time behavior of the SACF.

The decay patterns of the kinetic and potential correlations exhibit distinct density dependencies. At all three densities, the slope in the log-log plot of $C^{\mathrm{KK}}(\tau)$ at large $\tau\gtrsim 10$ is insensitive to change in number density, suggesting that the exponent $\alpha$ of the inverse-power decay $C^{\mathrm{KK}}(\tau)\sim\tau^{-\alpha}$ remains the same as $\alpha\approx 1.5$. The existence of such a universal exponent is strong evidence that this relaxation process has a hydrodynamic origin ErnstHauge1971 . On the other hand, the potential correlation decays overall faster than algebraically, implying that it has a different physical origin. With increasing number density, $C^{\mathrm{PP}}(\tau)$ increases, leading to an exponential tail at the high density limit. Our previous MD study on the shear viscosity at high density KimHan2018 suggests that the relaxation process involved in the potential correlation has a configurational origin.

To understand the intrinsic nature of kinetic and potential correlations, in the following sections we systematically examine their density dependence along with the system-size dependence. In Section III.1, we introduce an analytic expression of $C^{\mathrm{KK}}(t)$ based on molecular hydrodynamics to explain the long-time algebraic decay of the kinetic correlation and examine its validity by comparing with MD simulation results. In Section III.2, we examine the configurational nature of the potential correlation.

III.1 Hydrodynamic Nature of Kinetic Contribution

From the mode-coupling theory with the linearized Navier–Stokes (LNS) equations ErnstHauge1971 ; ErpenbeckWood1982 , an analytic expression for the long-time behavior of $C^{\mathrm{KK}}(t)$ in finite systems is given as

[TABLE]

where $k$ denotes the magnitude of wavevector $\bm{k}=\frac{2\pi}{L}(n_{x},n_{y},n_{z})$ with integers $n_{x}$, $n_{y}$, $n_{z}$. In the thermodynamic limit $L\rightarrow\infty$, the summation of exponential decay functions over wavevectors becomes an inverse-power decay:

[TABLE]

Since this approach assumes that the motion of a molecular fluid is well described by the LNS equations, the validity of Eqs. (9) and (10) is not guaranteed for short time $t$ comparable to or shorter than the average collision time $t_{0}$. However, our previous MD study on molecular hydrodynamics HanKimTalknerKarniadakisLee2018 showed that the dynamics of a molecular fluid is fairly well described even for several collision times. Longitudinal velocity field contributions are neglected and thus the effect of sound wave propagation under periodic boundary conditions is not included in Eq. (9).

In Figure 3, we first compare the time profiles of $C^{\mathrm{KK}}(t)$ obtained from MD simulations with the theoretical prediction given in Eq. (9). For large $t$, LNS predicts the asymptotic behavior of $C^{\mathrm{KK}}(t)$ regardless of number density and system size. As expected, the oscillations observed in MD results, which are caused by sound wave propagation across periodic boundaries ChoiHan2017 ; HanKimTalknerKarniadakisLee2018 , are not reproduced by the theoretical prediction. In addition, the discrepancy between MD and LNS is evident at short times, where dynamics is governed by detailed molecular interactions rather than by hydrodynamic laws. As observed in Figure 2, while the increase in number density weakens the kinetic correlation, it does not alter the inverse-power decay form of the long-time decay. It is also observed in Figure 3b that the time region exhibiting this characteristic decay gradually extends as the system size increases. Figure 4 clearly shows the transition from the short-time exponential decay to the long-time inverse-power decay observed from the $N=8192$ system at low density $n=0.45$.

To investigate further the hydrodynamic origin of the kinetic correlation, we subdivide $C^{\mathrm{KK}}(t)$ into the atomic (or self) correlation $C^{\mathrm{KK}}_{\mathrm{I}}(t)$ and the interatomic correlation $C^{\mathrm{KK}}_{\mathrm{IJ}}(t)$:

[TABLE]

Figure 5a shows the temporal changes in the contributions of these correlations to the total kinetic correlation. While $C_{\mathrm{I}}^{\mathrm{KK}}(t)$ dominates at short times, $C_{\mathrm{IJ}}^{\mathrm{KK}}(t)$ persists at long times. Since any term in $C^{\mathrm{KK}}$ at $t=0$ that is odd in velocity component vanishes due to symmetry, $C_{\mathrm{IJ}}^{\mathrm{KK}}(0)=0$ and thus $C^{\mathrm{KK}}(0)=C_{\mathrm{I}}^{\mathrm{KK}}(0)$. However, $C_{\mathrm{IJ}}^{\mathrm{KK}}(t)$ soon outgrows $C_{\mathrm{I}}^{\mathrm{KK}}(t)$ as the interatomic correlation starts to have nonzero values. This transition corresponds to the physical picture that momentum carried by a fluid particle is dissipated to the surrounding particles over time. In fact, as shown in Figure 5b, the long-time algebraic tail of $C^{\mathrm{KK}}(\tau)$ is well described solely by the interatomic correlation $C_{\mathrm{IJ}}^{\mathrm{KK}}(\tau)$.

III.2 Configurational Nature of Potential Contribution

Here we study the configurational aspect of the potential-potential component of the SACF. We first observe the influence of number density on $C^{\mathrm{PP}}(\tau)$. Figures 6a shows that the increase in number density strengthens the potential correlation. At low and intermediate densities, $C^{\mathrm{PP}}(\tau)$ is observed to decay rapidly. However, owing to numerical imprecision, their functional forms cannot be determined. On the other hand, an exponential tail is clearly observed at high density $n=0.85$. The decay pattern appearing after $\tau\approx 5t_{0}$ is well-described by an exponential decay with exponent 0.218. We believe that the different decay patterns of dilute and dense fluids are caused by the difference in the timescales of structural relaxation. The significance of structural relaxation in $C^{\mathrm{PP}}(t)$ has been repeatedly noted by Isobe and Alder IsobeAlder2009 ; IsobeAlder2010 ; IsobeAlder2012 . For dense two- and three-dimensional simple fluids, they have demonstrated the existence of the molasses regime in SACF and related its development to the reorganization of atomic clusters. The higher the number density, the longer the time it takes for the atomic cluster to dissociate, thereby result in an extensive molasses regime near the fluid-solid transition point.

Figure 6c compares the radial distribution functions (RDFs) at low, intermediate, and high densities. As a static equilibrium quantity, the RDF does not provide complete information to compute the potential correlation $C^{\mathrm{PP}}(t)$ except for the value at $t=0$ ZwanzigMountain1965 ,

[TABLE]

However, from the intensity of peaks in the RDF, a much slower structural relaxation is expected for high density than for the case of low density.

We next discuss the finite system-size effect on $C^{\mathrm{PP}}(\tau)$. The influence of system size is not noticeable in Figure 6b (note, however, that the log scale is used for the values of $C^{\mathrm{PP}}(\tau)$ here). Contrary to the kinetic correlation (see Figure 3b), the pattern of the tail is not sensitive to system size and a change in system size does not alter the dynamic feature of $C^{\mathrm{PP}}(t)$. However, perceptible system-size dependency is observed in the values of $C^{\mathrm{PP}}$ at $t=0$ (see Figure 9a). In fact, similar size dependency is observed consistently at different times, leading to similar size dependency in $\eta^{\mathrm{PP}}$ to be discussed in Section IV.2. At high density, $C^{\mathrm{PP}}(0)$ values show an oscillatory behavior with respect to system size, where increase in system size dampens the magnitude of the oscillation. At low density, $C^{\mathrm{PP}}(0)$ monotonically increases with increasing system size but at a very slow rate.

In principle, $C^{\mathrm{PP}}(0)$ can be computed using information on the equilibrium distribution such as the RDF, see Eq. (12). Figure 6d compares the RDFs obtained from various system sizes. Essentially the same $g(r)$ is reproduced as for small systems for the range of $r\lesssim L$. However, peaks appearing beyond this range are missing in the $g(r)$ for small systems. A more pronounced system-size dependence is expected to be observed at higher densities until all peaks in $g(r)$ are reproduced for sufficiently large systems. This induces a system-size dependence on $C^{\mathrm{PP}}(0)$.

IV Shear Viscosity

Depending on the number density of the simulation system, the shear viscosity exhibits two types of system-size behavior that are different in nature. Figure 7 compares the system-size effects on shear viscosity and its components at low and high densities. While an $L^{-1}$ correction with respect to system size $L$ is observed at low density, shear viscosity at high density shows an oscillatory behavior that dampens with increasing system size. The latter observation is very similar to that of our previous MD study of the Lennard-Jones fluid near the triple point KimHan2018 . The shear viscosity at intermediate density overall exhibits an $L^{-1}$ scaling behavior but with minor deviations at small simulation systems (not shown).

The effect of the system size on the shear viscosity is determined by the competition between kinetic and potential contributions in the SACF. At low density, where the portion of the kinetic contribution is significant (as in Figure 2a), the shear viscosity follows the overall behavior of the kinetic-kinetic component, denoted as the hydrodynamic system-size effect, resulting in an $L^{-1}$ scaling behavior. On the other hand, the SACF of a dense fluid is governed by the potential-potential component (as in Figure 2c). Consequently, the shear viscosity is mainly influenced by the system-size dependence of the potential-potential component, denoted as the configurational system-size effect. The difference between the viscosity values of smaller systems and the largest system are smaller than 3 percent for the low density and 1 percent for the high density. This is consistent with the weak system-size dependence observed in previous MD studies DaivisEvans1995 ; YehHummer2004 ; MoultosZhangTsimpanogiannisEconomouMaginn2016 .

In the following subsections, we study the hydrodynamic and configurational system-size effects on shear viscosity based on the observations made in Section III. Empirical formulas that describe the scaling behavior of the hydrodynamic system-size effect are proposed in Section IV.1. The relation between the complex size-dependent behavior of shear viscosity and that of $C^{\mathrm{PP}}(t=0)$ value is discussed in Section IV.2.

IV.1 Hydrodynamic System-Size Effect

In Section III.1, we separated $C^{\mathrm{KK}}(t)$ into the atomic contribution $C_{\mathrm{I}}^{\mathrm{KK}}(t)$ and the interatomic contribution $C_{\mathrm{IJ}}^{\mathrm{KK}}(t)$ and verified that the latter is responsible for the hydrodynamic long-time tail in $C^{\mathrm{KK}}(t)$ while the other dominates at shorter times. Likewise, we separate $\eta^{\mathrm{KK}}$ into two components, $\eta_{\mathrm{I}}^{\mathrm{KK}}$ and $\eta_{\mathrm{IJ}}^{\mathrm{KK}}$, and study their respective system-size dependence. Based on scaling relations observed in the MD simulation data of $C_{\mathrm{I}}^{\mathrm{KK}}(t)$ and $C_{\mathrm{IJ}}^{\mathrm{KK}}(t)$, their empirical formulas are proposed. From the relation $\eta_{\star}^{\mathrm{KK}}=\int_{0}^{\infty}C^{\mathrm{KK}}_{\star}(t)dt$ with $\star=\mathrm{I}$ or $\mathrm{IJ}$, the corresponding finite size effects of $\eta_{\star}^{\mathrm{KK}}$ are deduced.

In Figures 8a and 8b, we confirm that $C_{\mathrm{I},L}^{\mathrm{KK}}(t)$ and $C_{\mathrm{IJ},L}^{\mathrm{KK}}(t)$ for finite system size $L$ satisfy the following asymptotic formulas:

[TABLE]

where $C_{\mathrm{I},\infty}^{\mathrm{KK}}(t)$ and $C_{\mathrm{IJ},\infty}^{\mathrm{KK}}(t)$ denote the corresponding correlations in the thermodynamic limit $L\rightarrow\infty$. By assuming MD data from the largest $N=8192$ system at $n=0.45$ as bulk data, the time profiles of $f(t)$ and $g(t/L)$ were computed for each small system size using Eqs. (13) and (14). The remarkable coincidence of the time profiles obtained from various small system sizes indicates that our empirical asymptotic expressions are valid.

By applying the Green–Kubo relation to Eqs. (13) and (14), we obtain the finite system-size correction of $\eta^{\mathrm{KK}}_{L}$:

[TABLE]

where $F=\int_{0}^{\infty}f(t)dt$ and $G=\int_{0}^{\infty}g(t)dt$. Figures 8c and 8d display the actual scaling behaviors of $\eta_{\mathrm{I}}^{\mathrm{KK}}$ and $\eta_{\mathrm{IJ}}^{\mathrm{KK}}$ observed in MD simulations, which are consistent with Eq. (15).

Different finite system-size effects on $C_{\mathrm{I}}^{\mathrm{KK}}(t)$ and $C_{\mathrm{IJ}}^{\mathrm{KK}}(t)$ lead to different system-size corrections of $\eta^{\mathrm{KK}}$. The leading $L^{-1}$ correction is obtained from the interactomic correlation $C_{\mathrm{IJ}}^{\mathrm{KK}}(t)$. As discussed in Section III.1 using Figure 5, $C_{\mathrm{IJ}}^{\mathrm{KK}}(t)$ is responsible for the long-time tail of $C^{\mathrm{KK}}(t)$. Hence, the main finite system-size effect in $\eta^{\mathrm{KK}}$ is attributed to the disturbance in the long-time tail.

IV.2 Configurational System-Size Effect

In our previous study on the shear viscosity of a dense fluid KimHan2018 , we observed the remarkable resemblance between the size-dependent behaviors of $C(0)$ and the shear viscosity $\eta$ and drew attention to the significance of the configurational contribution in the estimation of finite system-size effects on shear viscosity. However, we were not certain whether the relation holds at intermediate and low densities and therefore did not claim the applicability of the uncertainty quantification method to other liquid states. Here we confirm that the system-size dependencies of $C^{\mathrm{PP}}(0)$ and $\eta^{\mathrm{PP}}$ are in consistent agreement regardless of number density, thereby giving validity to our quantification method for configurational system-size effect.

Figure 9 compares the system-size dependencies of $C^{\mathrm{PP}}(0)$ and $\eta^{\mathrm{PP}}$ at $n=0.45$ and 0.85. While the patterns of each quantity at the two densities are significantly different, the patterns of the two quantities at each density are remarkably similar. The damped oscillatory behavior of $\eta^{\mathrm{PP}}$ at high density is also observed in $C^{\mathrm{PP}}(0)$. On the other hand, a weak monotonic increase of $\eta^{\mathrm{PP}}$ at low density with respected to increasing system size is consistent with that of $C^{\mathrm{PP}}(0)$. In other words, the finite system-size effect of $\eta^{\mathrm{PP}}$ is well reproduced by $C^{\mathrm{PP}}(0)$.

While $C^{\mathrm{PP}}(t)$ and $\eta^{\mathrm{PP}}$ are related via the Green-Kubo relation, it is important to note that this does not automatically guarantee the remarkable resemblance of the finite system-size effects on $C^{\mathrm{PP}}(0)$ and $\eta^{\mathrm{PP}}$. In fact, such a coincidence is not expected to happen for the pair of kinetic components, $C^{\mathrm{KK}}(0)$ and $\eta^{\mathrm{KK}}$. In the potential case, the system-size dependence of $C^{\mathrm{PP}}(t)$ at each time $t$ remains similar resulting in a similar behavior for $\eta^{\mathrm{PP}}$, whereas in the kinetic case the main finite system-size effect on $\eta^{\mathrm{KK}}$ results from the change in the long-time tail. Hence, while $\eta^{\mathrm{PP}}$ is a dynamic property, its finite system-size effect can be roughly estimated from a static quantity $C^{\mathrm{PP}}(0)$, the accurate value of which can be much more easily computed.

V Conclusions

The calculation of transport coefficients has long been of interest and importance, see Ref. RiceGray1965 . We have studied the density-dependent finite system-size effects on the shear-stress autocorrelation function $C(t)$ and the shear viscosity coefficient $\eta$ by investigating hydrodynamic and configurational nature in the relaxation process of the off-diagonal shear-stress tensor component $p_{xy}$. Systematic MD simulations for a three-dimensional simple fluid at low, intermediate, and high densities have revealed that the shear viscosity of a dense fluid exhibits an oscillatory behavior that dampens with increasing system size $L$, whereas that of a dilute fluid has an $L^{-1}$ finite system-size correction. The former finite system-size effect at high density was identified to originate from configurational nature of the potential-potential component of the SACF, $C^{\mathrm{PP}}(t)$, whereas the latter at low density was shown to have a hydrodynamic origin arising from the kinetic-kinetic component, $C^{\mathrm{KK}}(t)$. Competition between these two contributions explains not only the strong density dependence of the shear viscosity but also that of its finite system-size effects.

Using analytic results obtained from the mode-coupling approach with the linearized Navier–Stokes equations ErnstHauge1971 ; ErpenbeckWood1982 , we performed a quantitative analysis on the kinetic correlation $C^{\mathrm{KK}}(t)$ and the kinetic component $\eta^{\mathrm{KK}}$ of the shear viscosity. The inverse-power decay of $C^{\mathrm{KK}}(t)$ was clearly observed in MD simulations at low and intermediate densities, which confirms the applicability of molecular hydrodynamic theory ChoiHan2017 ; HanKimTalknerKarniadakisLee2018 . Using the scaling behavior of the finite system-size effects of $C^{\mathrm{KK}}(t)$, we showed that the $L^{-1}$ correction of $\eta^{\mathrm{KK}}$ is due to the change in the long-time tail of $C^{\mathrm{KK}}(t)$. As seen Figure 3b, the latter change includes an oscillatory disturbance caused by traveling of sound waves across periodic boundaries and a cross-over from the algebraic decay to an exponential regime ascribable to the cut-off at low wavenumbers introduced by periodic boundaries. Similar behaviors have been observed in the velocity autocorrelation function ChoiHan2017 ; AstaLevesqueVuilleumierRotenberg2017 and the same form of correction has been derived for the self-diffusion coefficient DunwegKremer1993 ; YehHummer2004 . Contrary to the kinetic correlation, an analytic approach to determine the potential correlation $C^{\mathrm{PP}}(t)$ has not been successful due to the complicated structure of a nonlinear four-particle correlation function which plays a role in the structural relaxation. Motivated by our previous MD simulation study KimHan2018 , we focused on the remarkable resemblance between the finite system-size effects on $C^{\mathrm{PP}}(0)$ and $\eta^{\mathrm{PP}}$ regardless of density. This observation provides a practical suggestion of using $C^{\mathrm{PP}}(0)$ to roughly estimate the finite system-size effect on $\eta^{\mathrm{PP}}$. As a static equilibrium quantity, the accurate value of $C^{\mathrm{PP}}(0)$ can be much more easily computed than $\eta^{\mathrm{PP}}$ and thus the critical system size for an accurate value of $\eta^{\mathrm{PP}}$ can be efficiently estimated. We also point out that $C^{\mathrm{PP}}(0)$ can be investigated analytically using the radial distribution function.

Acknowledgements.

This work was supported in part by the U.S. Army Research Laboratory and was accomplished under Cooperative Agreement No. W911NF-12-2-0023, “Alliance for the Computationally-guided Design of Energy Efficient Electronic Materials (CDE3M)”.