Bridging the gap between molecular dynamics and hydrodynamics in nanoscale Brownian motions

TL;DR

This study uses molecular dynamics simulations to connect nanoscale Brownian motion with hydrodynamic behavior, revealing how solvation and molecular interactions influence particle dynamics and hydrodynamic properties.

Contribution

It bridges molecular dynamics and hydrodynamics in nanoscale Brownian motion, highlighting the effects of solvation structure and molecular interactions on hydrodynamic behavior.

Findings

Hydrodynamic radius varies with solvation free energy and potential type.

VACF matches Navier-Stokes analytical expressions at long times.

Long-time tail of VACF follows a t^{-3/2} decay.

Abstract

Through molecular dynamics simulations, we examined hydrodynamic behavior of the Brownian motion of fullerene particles based on molecular interactions. The solvation free energy and the velocity autocorrelation function (VACF) were calculated by using the Lennard-Jones (LJ) and Weeks-Chandler-Andersen (WCA) potentials for the solute-solvent and solvent-solvent interactions and by changing the size of the fullerene particles. We also measured the diffusion constant of the fullerene particles and the shear viscosity of the host fluid, and then the hydrodynamic radius $a_\mathrm{HD}$ was quantified from the Stokes-Einstein relation. The $a_\mathrm{HD}$ value exceeds that of the gyration radius of the fullerene when the solvation free energy exhibits largely negative values using the LJ potential. In contrast, $a_\mathrm{HD}$ becomes comparable to the size of bare fullerene, when the solvation free energy is positive using the WCA potential. Furthermore, the VACF of the fullerene particles is directly comparable with the analytical expressions utilizing the Navier-Stokes equations both in incompressible and compressible forms. Hydrodynamic long-time tail $t^{-3/2}$ is demonstrated for timescales longer than the kinematic time of the momentum diffusion over the particles' size. However, the VACF in shorter timescales deviates from the hydrodynamic description, particularly for smaller fullerene particles and for the LJ potential. This occurs even though the compressible effect is considered when characterizing the decay of VACF around the sound propagation time scale over the particles' size. These results indicate that the nanoscale Brownian motion is influenced by the solvation structure around the solute particles originating from the molecular interaction.

Tables

Table 1. Table 1: Radius of gyration of fullerene particle a 𝑎 a (in nm), mass density ratio ρ f / ρ p subscript 𝜌 f subscript 𝜌 p \rho_{\mathrm{f}}/\rho_{\mathrm{p}} between solvent fluid and fullerene particle, and size ratio L / a 𝐿 𝑎 L/a .

	C20	C60	C180	C240	C320	C540
$a$	0.20	0.34	0.60	0.69	0.80	1.04
${\rho }_f/{\rho }_p$	0.104	0.191	0.336	0.387	0.453	0.582
$L/a$	104.79	58.56	33.18	28.86	24.89	19.14

Table 2. Table 2: Physical parameters concerning host fluids for LJ and WCA potentials. η 𝜂 \eta denotes shear viscosity (in 10 − 4 superscript 10 4 10^{-4} kg m -1 s -1 ). η v subscript 𝜂 v \eta_{\mathrm{v}} denotes the bulk viscosity (in 10 − 4 superscript 10 4 10^{-4} kg m -1 s -1 ), and c 𝑐 c denotes the speed of sound (in nm/ps).

	$\eta$	${\eta }_v$	$c$
LJ	1.95	1.11	0.547
WCA	1.54	0.516	0.812

Table 3. Table 3: Diffusion constant D 𝐷 D (in 10 − 4 superscript 10 4 10^{-4} nm 2 /ps) and diffusion time τ D = a 2 / D subscript 𝜏 D superscript 𝑎 2 𝐷 \tau_{\mathrm{D}}=a^{2}/D (in ps) concerning fullerene particles for LJ and WCA potentials.

	C20	C60	C180	C240	C320	C540
$D$ (LJ)	8.75	5.05	2.87	3.13	2.52	1.88
$D$ (WCA)	15.8	9.24	6.52	5.64	5.62	4.04
${\tau }_D$ (LJ)	41	229	1246	1522	2549	5745
${\tau }_D$ (WCA)	23	125	548	845	1143	2680

Equations

$$M\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}=-{\int }_0^t\zeta (t-s)\mathbf{v}(s)\mathrm{d}s+\mathbf{R}(t),$$

$$M\frac{\mathrm{d}C(t)}{\mathrm{d}t}=-{\int }_0^t\zeta (t-s)C(s)\mathrm{d}s.$$

$$\tilde{C}(\omega )=\frac{MC(0)}{-iM\omega +\tilde{\zeta }(\omega )},$$

$$\mathbf{F}(t)=-6\pi \eta a\mathbf{v}-\frac{M_{\mathrm{f}}}{2}\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}-6a^2\sqrt{\pi \eta {\rho }_{\mathrm{f}}}{\int }_0^t\frac{\mathrm{d}\mathbf{v}/\mathrm{d}t}{\sqrt{t-s}}\mathrm{d}s,$$

$$\tilde{\zeta }(\omega )=6\pi \eta a-i\omega \frac{M_{\mathrm{f}}}{2}+6\pi a^2\sqrt{-i\omega \eta {\rho }_{\mathrm{f}}}.$$

$$C^{\nu }(t)=\frac{k_{\mathrm{B}}T}{3M}\frac{2{\rho }_{\mathrm{p}}}{3{\rho }_{\mathrm{f}}}\frac{1}{3\pi }{\int }_0^{\infty }\frac{e^{-st/{\tau }_{\nu }}{s}^{1/2}}{1+{\sigma }_{1}s+{\sigma }_2{s}^2}\mathrm{d}s,$$

$$C^{\mathrm{c}}(t)=\frac{k_{\mathrm{B}}T}{M}\frac{e^{-{\alpha }_{1}t/{\tau }_{\mathrm{c}}}}{1+2{\rho }_{\mathrm{p}}/{\rho }_{\mathrm{f}}}\left[\cos \left(\frac{{\alpha }_{2}t}{{\tau }_{\mathrm{c}}}\right)-\frac{{\alpha }_1}{{\alpha }_2}\sin \left(\frac{{\alpha }_{2}t}{{\tau }_{\mathrm{c}}}\right)\right],$$

$$\tilde{\zeta }(\omega )=\frac{4\pi }{3}\eta a{x}^2\frac{(1+x)(9-9iy-2y^2)+x^2(1-iy)}{2x^2(1-iy)-(1+x)y^2-x^2{y}^2},$$

$$\tilde{c}={\left[c^2-\frac{i\omega }{{\rho }_{\mathrm{f}}}\left(\frac{4}{3}\eta +{\eta }_{\mathrm{v}}\right)\right]}^{1/2},$$

$${\xi }_{\mathrm{E}}=\frac{8}{3}{\left(\frac{2\pi k_{\mathrm{B}}T{m}_{\mathrm{Ar}}M}{m_{\mathrm{Ar}}+M}\right)}^{1/2}\frac{{\rho }_{\mathrm{N}}g(r_{\max })r_{\max }^2}{m_{\mathrm{Ar}}}\frac{1+2\chi }{1+\chi },$$

$$\frac{k_{\mathrm{B}}T}{6\pi \eta a_{\mathrm{HD}}}=\frac{1}{3}{\int }_0^{\infty }C(t)\mathrm{d}t.$$

$${\xi }_{\mathrm{E}}=\frac{8}{3}{\left(\frac{\pi k_{\mathrm{B}}T}{m_{\mathrm{Ar}}}\right)}^{1/2}{\rho }_{\mathrm{N}}{r}_{\max }^{2}g(r_{\max }),$$

Full text

I Introduction

In colloids, the macroscopic solute particles are dispersed in a liquid solvent. To predict flow behaviors of colloidal dispersions, not only the motions of the solute particles but also their coupling with the solvent dynamics must be considered. Note that the spatial and temporal scales of solute particles are in orders of magnitudes larger than those of solvent molecules. Thus, the dynamics of colloidal dispersions are mostly governed by the coupling effects of Brownian motion of the colloidal particles and the hydrodynamics of the solvent Russel (1981); Russel et al. (1992); Dhont (1996); Bian et al. (2016).

This indicates that the colloidal system is a typical example of multi-scale physics, including hierarchical scales, and offers a good target to be solved in computational science. In recent years, many coarse-grained methods have been developed utilizing the scale separation between solute and solvent particles. These include stochastic rotation dynamics/multiple-particle collision dynamics methods Malevanets and Kapral (2000); Padding and Louis (2004); Padding et al. (2005); Padding and Louis (2006); Gompper et al. (2009); Huang et al. (2012); Theers et al. (2016), the lattice Boltzmann method Ladd (1993); Lobaskin and Dunweg (2004); Cates et al. (2004); Chatterji and Horbach (2005); Poblete et al. (2014), Stokesian dynamics Ermak and McCammon (1978); Brady and Bossis (1988), direct numerical simulations using the immersed boundary method Peskin (2002); Atzberger et al. (2007); Sharma and Patankar (2004), Fluid Particle Dynamics Tanaka and Araki (2000); Kodama et al. (2004); Tanaka and Araki (2006); Furukawa and Tanaka (2010); Furukawa et al. (2018), and the Smoothed Profile Method Nakayama and Yamamoto (2005); Kim et al. (2006); Yamamoto et al. (2007); Nakayama et al. (2008); Iwashita et al. (2008); Yamamoto et al. (2008, 2009); Nakayama et al. (2010); Tatsumi and Yamamoto (2012); Luo et al. (2009).

The aforementioned methods enable consistent simulation of fluctuating hydrodynamic descriptions for colloidal particles. Their reliability has been conventionally tested by calculating the velocity autocorrelation function (VACF) of a single solute particle and comparing it with the analytical solution obtained by solving the generalized Langevin equation, which considers the hydrodynamic memory effect. It utilizes the general expression for the frequency-dependent hydrodynamic friction coefficient $\tilde{\zeta}(\omega)$ of a rigid spherical particle suspended in a viscous fluid, Here, $\tilde{\zeta}(\omega)$ is derived from the Navier–Stokes (NS) equation (see more details in Section II). In particular, the hydrodynamic memory effect induces a long-time tail in the VACF owing to algebraic decay $t^{-3/2}$, which was discovered by Alder and Wainwright using molecular dynamics (MD) simulations Alder and Wainwright (1970). Recently, the hydrodynamic memory effect was measured through experiments using a particle tracking technique Franosch et al. (2011); Jannasch et al. (2011); Huang et al. (2011); Kheifets et al. (2014); Mo and Raizen (2019).

However, when the solute particle is of nanoscale and is comparable in size to the solvent particles, the separation of spatial and temporal scales becomes unclear and the validity of the continuum description becomes questionable. Moreover, the most generally used assumption is the incompressible condition for the host fluid; however, its validity also becomes unreasonable because of the sound effect propagating over the molecular length scale. That is, the particle momentum is transferred by sound waves at short time intervals and by vortex formation around the particles at long time intervals. Therefore, the whole aspect of the multi-scale hierarchy must be clarified by using all-atom MD simulations. The question that arises is how the molecular interactions are relevant to the hydrodynamics interactions occurring by sound propagation and momentum diffusion via the kinematic viscosity.

MD simulations have been intensively performed for observing the long-time tail in the VACF for pure Lennard–Jones (LJ) and Weeks–Chandler–Andersen (WCA) fluids Levesque and Ashurst (1974); Erpenbeck and Wood (1985); McDonough et al. (2001); Dib et al. (2006); Lesnicki et al. (2016); Han et al. (2018); Ignatyuk et al. (2018). Very recently, the velocity field generated through MD simulations was directly compared with that described by the linearized NS equations Han et al. (2018). The validity of the combined Langevin and hydrodynamic model for the Brownian motion has also been examined using MD simulations Bocquet and Barrat (1994); Bocquet et al. (1997); Nuevo et al. (1998); Ould-Kaddour and Levesque (2000); Schmidt and Skinner (2003); Sokolovskii et al. (2006); McPhie et al. (2006); Ould-Kaddour and Levesque (2007); Jung et al. (2017). Many efforts have been devoted to discussions of the microscopic origin of the hydrodynamic radius $a_{\mathrm{HD}}$ and the boundary condition at the solute-solvent interface with respect to the Stokes–Einstein (SE) relation, $D=k_{\mathrm{B}}T/(c\pi\eta a_{\mathrm{HD}})$, between the diffusion constant $D$ of the solute particle and the shear viscosity of the fluid $\eta$ Hansen and McDonald (2013). Here, $k_{\mathrm{B}}$ and $T$ are the Boltzmann constant and the temperature, respectively. The constant $c$ is determined by stick ($c=6$) or slip ($c=4$) boundary conditions imposed at the particle-fluid interface. Note that the concepts of the hydrodynamic radius $a_{\mathrm{HD}}$ and its association with the hydrodynamic boundary condition become ambiguous at molecular scales. In fact, $a_{\mathrm{HD}}$ is actually “defined” by the SE relation and its value is influenced by the choice of $c$. Given that the (macroscopic) hydrodynamics is implemented with the stick boundary condition, $c=6$ is a natural choice when the solute size is to be varied continuously over a wide range by bearing in mind the macroscopic limit.

In Ref. Li (2009), Li Eqs. (6) + (7)demonstrated that the van der Waals interaction between nanoscale LJ clusters and solvent particles plays a crucial role in determining the hydrodynamic radius $a_{\mathrm{HD}}$. In addition, the effect of solvation free energy on the hydrodynamic transport of the fullerene particles suspended in a water solvent has been investigated by Morrone et al. Morrone et al. (2012) Other MD simulation study has been reported for a system composed of a LJ cluster suspended in LJ fluids Chakraborty (2011). Remarkably, Chakraborty used MD simulations to demonstrate the crossover of hydrodynamic effects from compressible to incompressible fluids Chakraborty (2011). However, the interplay between the hydrodynamic behavior and solvation free energy has not been thoroughly elucidated yet.

In this study, we used MD simulations to comprehensively examine both hydrodynamic and thermodynamic properties of nanoscale fullerene particles dispersed in two types of solvents by using the LJ and WCA potentials. The contributions of the present study are threefold. First, we analyzed the solvation free energy of a fullerene particle to investigate how its solvation structure depends on the molecular interaction. Second, we quantified the hydrodynamic radius $a_{\mathrm{HD}}$ from the diffusion constant and the SE relation by assuming the stick boundary condition. We examined the impact of intermolecular interactions on the hydrodynamic radius $a_{\mathrm{HD}}$ and discussed the results in terms of the solvation free energy. Third, we investigated the VACF of the fullerene particle to characterize the hydrodynamic long-time tail. The sound propagation effect on the VACF is then discussed. The VACF in MD simulations was compared with the analytic expressions utilizing the frequency-dependent friction $\tilde{\zeta}(\omega)$, which was obtained by solving the NS equation.

The remainder of this paper is organized as follows. Section II introduces the hydrodynamic model for the VACF in the Brownian motion. We explain the MD simulation details in Section III, and present the numerical results regarding the solvation free energy, hydrodynamic radius, and VACF in Section IV. Our conclusions are drawn in Section V, before presenting an Appendix that provides the numerical results for the VACF in pure LJ and WCA fluids.

II Overview of hydrodynamic descriptions of VACF

Here, we briefly review the theoretical descriptions of the hydrodynamics for the VACF of a colloidal particle. The generalized Langevin equation for a spherical particle with mass $M$ suspended in a fluid exhibiting fluctuating hydrodynamics has been analyzed in various studies Zwanzig and Bixon (1970); Chow and Hermans (1972, 1973); Hauge and Martin-Löf (1973); Bedeaux and Mazur (1974); Hinch (1975); Metiu et al. (1977); Español (1995); Felderhof (2005). Moreover, Bian et al. reviewed the recent progress on the Brownian motion Bian et al. (2016).

The equation of motion is written as

[TABLE]

where $\zeta(t)$ and $\bm{R}(t)$ represent the memory kernel of the friction coefficient and the random force acting on the particle, respectively. $\bm{R}(t)$ satisfies the fluctuation-dissipation theorem $\langle\bm{R}(t)\cdot\bm{R}(0)\rangle=3k_{\mathrm{B}}T\zeta(t)$ with the zero-mean value $\langle\bm{R}(t)\rangle=\bm{0}$. Here, $\langle\cdots\rangle$ represents the ensemble average. The VACF of the particle is defined as $C(t)=\langle\bm{v}(t)\cdot\bm{v}(0)\rangle$, the time evolution of which is given by

[TABLE]

The Laplace transform into the frequency ($\omega$) domain reduces to

[TABLE]

where $C(0)=k_{\mathrm{B}}T/M$ according to the equipartition theorem. Note that the zero-frequency limit $\tilde{C}(0)$ corresponds to the diffusion constant $D=k_{\mathrm{B}}T/\zeta_{0}$ with $\zeta_{0}=\tilde{\zeta}(0)$. This is equivalent to the Einstein relation, where $D$ is determined via the mean square displacement at long times.

For an incompressible fluid, the linearized NS equation results in the Basset–Boussinesq–Oseen equation,

[TABLE]

which describes a force acting on a spherical particle with instantaneous velocity $\bm{v}$ and acceleration $\mathrm{d}\bm{v}/\mathrm{d}t$ in low-Reynolds-number regimes Landau and Lifshitz (1987). $\eta$ and $\rho_{\mathrm{f}}$ denote the shear viscosity and the mass density of the solvent fluid, respectively. In addition, $a$ and $M_{\mathrm{f}}=4\pi a^{3}\rho_{\mathrm{f}}/3$ are the particle radius and the added mass due to the replacement of the fluid by the particle, respectively. That is, the particle is considered to move with the mass $M+M_{\mathrm{f}}/2$ in the incompressible fluid, where the sound is assumed to propagate with the infinite speed.

According to Eq. (4), the frequency-dependent friction coefficient $\tilde{\zeta}(\omega)=-\tilde{F}(\omega)/\tilde{v}(\omega)$ is expressed as

[TABLE]

Thus, the zero-frequency limit $\zeta_{0}=6\pi\eta a$ corresponds to the Stokes drag force, resulting in the SE formula, $D=k_{\mathrm{B}}T/(6\pi\eta a)$, under the stick boundary condition. The third term, which is proportional to $\sqrt{\omega}$, causes the decay of $\zeta(t)$ to $t^{-3/2}$, which is the source of the long-time tail in VACF. The expression of VACF can be written through Eqs. (3) and (5) to

[TABLE]

with the fullerene particle mass density $\rho_{\mathrm{p}}=M/(4\pi a^{3}/3)$ and the kinematic time of the momentum diffusion over the particle size $\tau_{\nu}=a^{2}/\nu$ Paul and Pusey (1999); Chakraborty (2011). Here, the kinematic viscosity is defined as $\nu=\eta/\rho_{\mathrm{f}}$. Factors $\sigma_{1}$ and $\sigma_{2}$ are defined as $\sigma_{1}=(1/9)(7-4\rho_{\mathrm{p}}/\rho_{\mathrm{f}})$ and $\sigma_{2}=(1/9)^{2}(1+2\rho_{\mathrm{p}}/\rho_{\mathrm{f}})^{2}$, respectively. The asymptotic behavior of $C^{\mathrm{\nu}}(t)$ is expressed as $C^{\mathrm{\nu}}(t)\simeq(2k_{\mathrm{B}}T/3\rho_{\mathrm{f}})(4\pi\nu t)^{-3/2}$ for long times. Note that the zero-time value of the VACF becomes $C^{\mathrm{\nu}}(0)=k_{\mathrm{B}}T/(M+M_{\mathrm{f}}/2)$ owing to the effect of added mass $M/2$, which deviates from the result of the equipartition theorem.

To describe the short-time relaxation of the VACF appropriately, a correction term is introduced as

[TABLE]

where $\tau_{\mathrm{c}}=a/c$ denotes the sound propagation time over the particle size with the speed of sound in a compressible fluid Chow and Hermans (1973); Zwanzig and Bixon (1975); Chakraborty (2011). In addition, $\alpha_{1}=(1+\rho_{\mathrm{f}}/2\rho_{\mathrm{p}})$ and $\alpha_{2}=(1-\rho_{\mathrm{f}}^{2}/4\rho_{\mathrm{p}}^{2})^{1/2}$. The initial value of the VACF given by $C(t)=C^{\mathrm{\nu}}(t)+C^{\mathrm{c}}(t)$ eventually recovers the result of the equipartition theorem, $C(0)=k_{\mathrm{B}}T/M$. In deriving Eq. (7), the time separation as $\epsilon=\tau_{\mathrm{c}}/\tau_{\nu}\ll 1$ was assumed, in which the contribution of sound propagation to the decay of the VACF is much faster than that of momentum diffusion owing to the fluid viscosity. Here, $\epsilon$ represents the non-dimensional factor required to characterize the fluid incompressibility Tatsumi and Yamamoto (2012). This linear combination formula has been examined via stochastic rotation dynamics Padding and Louis (2006) and MD simulations Chakraborty (2011).

Previous studies have also analyzed the frequency-dependent hydrodynamic friction coefficient, $\tilde{\zeta}(\omega)$, in the compressible fluid by using the linearized NS equations and the relationship between the pressure and density fields, $\nabla p=c^{2}\nabla\rho_{\mathrm{f}}$ Zwanzig and Bixon (1970); Bedeaux and Mazur (1974); Metiu et al. (1977); Felderhof (2005). $\tilde{\zeta}(\omega)$ is expressed as

[TABLE]

with $x=a(-i\omega\rho_{\mathrm{f}}/\eta)^{1/2}$ and $y=a\omega/\tilde{c}$ Felderhof (2005). Here, the frequency-dependent speed of sound $\tilde{c}$ is given by

[TABLE]

with the bulk viscosity $\eta_{\mathrm{v}}$. Equation (8) can be applied to a high-compressibility fluid exhibiting $\epsilon>1$, where the sound propagation is slower than the momentum diffusion. High-compressibility factor $\epsilon>1$ causes a peculiar “backtracking,” which corresponds to the negative contribution in the VACF. This is due to the inversion of the particle at short time scales; this in turn is induced by the non-uniform fluid density field around the moving solute particle Felderhof (2005). Note that a sufficiently small ratio of $\eta_{\mathrm{v}}/\eta$ is another important factor for a slower sound propagation because larger bulk viscosity of a fluid attenuates sound propagation. Comparisons with the simulated VACF have been made through multi-particle collision dynamics Belushkin et al. (2011); Poblete et al. (2014) and direct numerical simulation of fluctuating hydrodynamics Tatsumi and Yamamoto (2012). However, to the best of our knowledge, a thorough examination of the VACF obtained from Eq. (8) has not yet been performed via MD simulations .

III Model and simulation methods

The Gromacs package was used to conduct MD simulations for one fullerene particle suspended in a solvent consisting of $N=160,000$ Ar molecules Hess et al. (2008); Abraham et al. (2015). This simulation setup was similar to that in a previous MD simulation study Ishii and Ohtori (2016). For the fullerene particles, Cn ($n=20$, $60$, $120$, $240$, $320$, and $540$) were used. These fullerenes are good models of spherical particles. We utilized the geometrical coordinates provided by Tomanek Tomanek (2014). All the C-C distances in the fullerene were constrained with the LINCS algorithm. The radii of gyration of the fullerenes $a$ are listed in Table 1.

The interaction is described by the LJ potential, $U_{\mathrm{LJ}}(r)=4\epsilon_{\alpha\beta}[(\sigma_{\alpha\beta}/r)^{12}-(\sigma_{\alpha\beta}/r)^{6}]$, where $r$ is the distance between two atoms and $\alpha,\beta=$ Ar, C. The Lorentz–-Berthelot combination rule of $\sigma_{\alpha\beta}=(\sigma_{\alpha}+\sigma_{\beta})/2$ and $\epsilon_{\alpha\beta}=\sqrt{\epsilon_{\alpha}\epsilon_{\beta}}$ was utilized for the interactions between Ar and C atoms. Furthermore, the parameters, $\sigma_{\mathrm{Ar}}$ = $\sigma_{\mathrm{C}}$ = 0.34 nm, $\epsilon_{\mathrm{Ar}}/k_{\mathrm{B}}$ = 117.8 K and $\epsilon_{\mathrm{C}}/k_{\mathrm{B}}$ = 43.3 K were used, and the cutoff distance $r_{\mathrm{c}}$ was chosen as 1.2 nm or 0.382 nm. The value 1.2 nm corresponds to the conventional value in the LJ potential, whereas 0.382 nm corresponds to $r_{\mathrm{c}}=2^{1/6}\sigma_{\mathrm{Ar}}$. This generates a purely repulsive potential, which is the so-called WCA potential, $U_{\mathrm{WCA}}(r)=U_{\mathrm{LJ}}(r)+\epsilon_{\alpha\beta}$ $(r<r_{\mathrm{c}})$. In this study, potentials with these two cutoff lengths are referred to as LJ and WCA potentials, respectively. Note that the Ar-Cn and Ar-Ar interactions are of the same type; both of them are chosen from either LJ or WCA potential.

The linear dimension of the simulation box was $L=19.91$ nm, and the mass and number densities of the solvent were $\rho_{\mathrm{N}}=20.27$ nm*-3* and $\rho_{\mathrm{f}}=\rho_{\mathrm{N}}m_{\mathrm{Ar}}=1,345$ kg/m3, respectively. Here, $m_{\mathrm{Ar}}$ represents the atomic mass of Ar. This number density corresponds to $\rho_{\mathrm{N}}{\sigma_{\mathrm{Ar}}}^{3}=0.797$ in the LJ units. The mass density of the fullerene is denoted by $\rho_{\mathrm{p}}=M/(4\pi a^{3}/3)$, considering the mass of the fullerene as $M=nm_{\mathrm{C}}$, where $m_{\mathrm{C}}$ is the atomic mass of carbon. Then, the mass density ratios between the Ar solvent fluid and the fullerene particles, i.e., $\rho_{\mathrm{f}}/\rho_{\mathrm{p}}$, are presented in Table 1. Each system was first equilibrated with the $NVT$ ensemble at the temperature $T=95$ K, corresponding to $k_{\mathrm{B}}T/\epsilon_{\mathrm{Ar}}=0.806$ in the LJ units. Then, the $NVE$ ensemble simulations were performed for 10 ns to generate 20 independent trajectories at each system. In all simulations, periodic boundary conditions were utilized with a time step of 1 fs.

The parameters of the host fluid were determined beforehand through MD simulations for a pure solvent particle system as follows: Shear viscosity $\eta$ and bulk viscosity $\eta_{\mathrm{v}}$ were determined using the Green–Kubo formula for the off-diagonal and diagonal stress tensor, respectively. In addition, the speed of sound $c$ was quantified from the numerical calculations of $\sqrt{(\partial p/\partial\rho_{\mathrm{f}})_{T}}$. The obtained parameters are presented in Table 2.

Periodic images can influence the hydrodynamic behavior in MD simulations due to its long-range interaction Yeh and Hummer (2004). The larger size ratio $L/a$ is thus required to characterize the spatial extent of the momentum transfer with respect to viscosity in MD simulations. The size ratios $L/a$ of our system are presented in Table 2, ranging from 19.14 (C540) to 104.79 (C20). Note that MD simulations of LJ cluster dispersions have been performed to demonstrate long-time tails with the size ratios $L/a=22.58$ and $13.54$ in Ref. Chakraborty (2011). To check the finite-size effects on the VACF, we simulated another system with $N=20,000$ Ar particles. The corresponding linear dimension was $L=9.96$ nm. In this smaller system, the desired hydrodynamic long-time tail was masked by the finite-size artifact and was hardly observed, particularly in a larger fullerene particle system (for example, $L/a=9.57$ for C540). Thus, in the following sections, we show the simulation results for $N=160,000$ Ar particle systems.

IV Results and discussion

IV.1 Solvation free energy

We first examined the solvation free energy $\Delta G$, which is the transfer free energy of a solute from a vapor to a solvent, for the fullerene particles both in LJ and WCA potential systems. Note that $\Delta G$ was calculated using the Bennett acceptance ratio method Benett (1976). Figure 1 shows the results for $\Delta G$ as a function of the bare fullerene size, $a$. $\Delta G$ largely decreases with increasing fullerene size for the solvent using the LJ potential; this serves as a good solvent for a larger sized fullerene. By contrast, the $\Delta G$ of the WCA potential becomes positive, resulting in a solvation structure that differs from that obtained in the LJ potential case. These results lead to the conclusion that the change in solvation free energy of the fullerene is largely negative owing to the van der Waals attraction between the carbons and the solvent Ar molecules.

The solvation structure around the fullerene particle was investigated with respect to the radial distribution function (RDF), $g(r)$, between the center of mass of the fullerene and solvent particles. The results are shown in Fig. 2. As demonstrated in Fig. 2(a), for the solvent using the LJ potential, the maximum peak of $g(r)$ increases with increasing the fullerene size, and correspondingly the $g(r)$ exhibits the intense oscillation. This strong solvation structure, in which the fullerene particle is presumably bounded by the solvent particles, is consistent with the negative $\Delta G$ value resulting from the van der Waals attraction between the fullerene and the solvent molecules. In contrast, the peak of $g(r)$ decreases in the case of the solvent of the WCA potential, as shown in Fig. 2(b). This solvophobic property of the fullerene particle also agrees with the positive value of $\Delta G$. Furthermore, we examined $g(r-(a+\sigma_{\mathrm{Ar}}/2))$ to take into account the peak position shift with increasing the fullerene size $a$. The profiles of $g(r-(a+\sigma_{\mathrm{Ar}}/2))$ are shown in Fig. 2(c) and (d) for the LJ and WCA potentials, respectively. It is demonstrated that the peak positions are scaled using the distance $r-(a+\sigma_{\mathrm{Ar}}/2)$ both in the LJ and WCA systems. The first maximum positions were observed to be located around $0.17$ nm in all RDFs. This size corresponds to $\sigma_{\mathrm{C}}/2$ because the C atom has a collision diameter $\sigma_{\mathrm{C}}$ at the spherical surface with the gyration radius $a$. Then, the effective size of the fullerene particle can be defined from the difference between the first maximum position $r_{\mathrm{max}}$ of $g(r)$ and the radius of the solvent particle, which is expressed as $a_{\mathrm{RDF}}=r_{\mathrm{max}}-\sigma_{\mathrm{Ar}}/2$.

IV.2 Hydrodynamic radius

We calculated the mean square displacement of the fullerene particle, $\langle\delta r(t)^{2}\rangle=\langle|\bm{r}(t)-\bm{r}(0)|^{2}\rangle$, where $\bm{r}(t)-\bm{r}(0)$ represents the displacement vector of the center of mass of the fullerene particle during the time interval $t$. Diffusion constant $D$ was determined from the Einstein relation, $D=\lim_{t\to\infty}\langle\delta r(t)^{2}\rangle/(6t)$, as shown in Fig. 3. Table 3 presents diffusion constant $D$ and diffusion time $\tau_{\mathrm{D}}=a^{2}/D$, during which the fullerene particle diffuses over the radius. Note that the diffusion time $\tau_{\mathrm{D}}$ is more than an order of magnitude larger than $\tau_{\mathrm{c}}$ and $\tau_{\nu}$. These results indicate that the diffusion constant $D$ decreases with increasing $a$, and correspondingly, diffusion time $\tau_{\mathrm{D}}$ increases. Furthermore, diffusion constant $D$ seems to reduce owing to the attraction of the LJ potential compared with the value obtained using the WCA potential.

Hydrodynamic radius $a_{\mathrm{HD}}$ was determined through the SE relation assuming the stick boundary condition of $a_{\mathrm{HD}}=k_{\mathrm{B}}T/(6\pi D\eta)$. Figure 4 shows the comparison between either $a_{\mathrm{HD}}$ or $a_{\mathrm{RDF}}$ and the bare fullerene size $a$. We observed that the increasing manner of $a_{\mathrm{RDF}}$ was akin to that of $a$ both for the LJ and WCA potentials, by exhibiting the constant difference $a_{\mathrm{RDF}}-a\approx\sigma_{\mathrm{C}}/2$. These behaviors are consistent with the scaled RDF profiles, $g(r-(a+\sigma_{\mathrm{Ar}}/2))$ (see Fig. 2(c) and (d)). However, $a_{\mathrm{HD}}$ of the LJ potential increases more than bare radius $a$, whereas that of the WCA potential is comparable to $a_{\mathrm{RDF}}$. In particular, for large fullerene particles in the LJ solvent, $a_{\mathrm{HD}}$ was larger than $a_{\mathrm{RDF}}$ by a value corresponding to several solvation shells. Note that the hydrodynamic radius $a_{\mathrm{HD}}$ will become a larger value, if we assume the slip boundary condition for the SE relation, $a_{\mathrm{HD}}=k_{\mathrm{B}}T/(4\pi D\eta)$.

This apparent deviation of $a_{\mathrm{HD}}$ from $a$ and $a_{\mathrm{RDF}}$ in the LJ potential, which increases more than $a$, is explained by the negative value of the solvation free energy $\Delta G$ and the associated strong solvation structure around the fullerene particle, as demonstrated in Figs. 1 and 2. However, it is reasonable to assume that the hydrodynamic radius $a_{\mathrm{HD}}$ will merge into the bare size $a$ at the macroscopic regime ($a_{\mathrm{HD}}/a\to 1$), where the size of the solute particle becomes many orders of magnitude larger than the solvent particle. This is due to the fact that the spatial resolution for the molecular size is completely lacked and the hydrodynamic description becomes justified with the stick boundary condition at the macroscopic regime.

IV.3 VACF

Numerical results pertaining to the center of mass VACF of fullerene particles with regard to LJ and WCA potentials are depicted in Fig. 5 and Fig. 6, respectively. Additionally, numerical results obtained for VACF in pure LJ and WCA fluids are reported in Appendix A.

In each plot depicted in Fig. 5 and Fig. 6, results obtained from MD simulations have been compared against hydrodynamic descriptions previously explained in Section II. In addition, the short time decay of the VACF has also been compared against the Enskog theory, yielding the exponential decay relation, $C^{\mathrm{E}}(t)=(k_{\mathrm{B}}T/M)\exp(-\xi_{\mathrm{E}}t)$, using the Enskog friction coefficient given by

[TABLE]

with $\chi=I/Ma^{2}$ obtained using the moment of inertia $I$ of the fullerene particles Subramanian and Davis (1975). Moreover, $g(r_{\mathrm{max}})$ denotes the peak height of the solute-solvent RDF, $g(r)$, at $r_{\mathrm{max}}$ (refer Fig. 2). It must be noted that the Enskog type exponential decay has a physical origin different from that of the NS equation.

Figure 5 demonstrates that the VACF of all fullerene particles using the LJ potential system exhibits a long-time tail beyond $\tau_{\nu}$. MD simulation results obtained for $t\gg\tau_{\nu}$ were observed to be consistent with those obtained using hydrodynamic descriptions of Eqs. (6) + (7) or Eq. (8), the long time asymptote of which can be expressed as $(2k_{\mathrm{B}}T/3\rho_{\mathrm{f}})(4\pi\nu t)^{-3/2}$.

In the initial time region, VACF results obtained from MD simulations demonstrated good agreement with analytical hydrodynamics predictions. This might be puzzling since the VACF is expected to be governed by the Enskog kinetic theory, yielding the exponential decay, $C^{\mathrm{E}}(t)=(k_{\mathrm{B}}T/M)\exp(-\xi_{\mathrm{E}}t)$, owing to limitation pertaining to the continuum description of solvent fluids. Note that the decay time of $C^{\mathrm{c}}(t)$ can be expressed as $\tau_{\mathrm{c}}/(1+\rho_{\mathrm{f}}/2\rho_{\mathrm{p}})$ in Eq. (7). From $\rho_{\mathrm{f}}/\rho_{\mathrm{p}}$ values in Table 1, this time scale was observed to be relatively close to the Enskog time $\xi_{\mathrm{E}}^{-1}$ determined from the MD simulations. Due to its construction, the hydrodynamic description by Eqs. (6) + (7) should agree with the MD results at $t$ close to [math], and the notable point is that the short time decay of $C^{\mathrm{c}}(t)$, i.e., $\tau_{\mathrm{c}}/(1+\rho_{\mathrm{f}}/2\rho_{\mathrm{p}})$, for the fullerene is close to $\tau_{\mathrm{c}}$ according to the $\rho_{\mathrm{f}}/\rho_{\mathrm{p}}$ values in Table 1. As shown in Appendix A, this kind of agreement does not hold in pure LJ and WCA fluids, where the mass density ratio is estimated as $\rho_{\mathrm{f}}/\rho_{\mathrm{p}}=\rho_{\mathrm{N}}/(\pi/6)\approx 1.522$. It has been demonstrated in Fig. 5 that the Enskog theory provides a reasonable explanation for short time VACF decays observed over small time instants, $t\lesssim\xi_{\mathrm{E}}^{-1}$. It must be noted that VACF of pure LJ and WCA solvents, for which the tagged solvent particle could be considered as a consolidated solute, could be well described using the Enskog theory, as demonstrated in Appendix A.

Deviations from the theoretical expressions described in Eqs. (6) + (7) and Eq. (8) become noticeable during the intermediate time period prior to commencement of the kinematic time $\tau_{\nu}$ over which velocity diffuses the radius of the fullerene particles. VACF obtained from MD simulations were observed to be less compared to those obtained from hydrodynamic descriptions involving linearized NS equation. This decrease in VACF is directly related to the hydrodynamic radius $a_{\mathrm{HD}}$ which was observed to be larger compared to the bare fullerene radius $a$ in accordance with the following VACF integral,

[TABLE]

It is also of interest to observe in Fig. 5 that the replacement of $a$ by $a_{\mathrm{HD}}$ in Eq. (6) results in better characterization of VACF obtained from MD simulations, particularly with regard to larger fullerene particles, e.g., Cn ($n\geq 120$). In this expression, the mass density of the solute particle correspondingly $\rho_{\mathrm{p}}$ changes in accordance with the relation $\rho_{\mathrm{p}}+4\pi\rho_{\mathrm{f}}(a_{\mathrm{HD}}^{3}-a^{3})/3$, incorporating mass of the solvent particle within the hydrodynamic radius, $a_{\mathrm{HD}}$. This observation implies that the fullerene particle transport occurs in conjunction with that of the surrounding solvation structure, the size of which is characterized by $a_{\mathrm{HD}}$.

When the size of the fullerene particles becomes comparable with that of the solvent particles, the observed value of the compressibility factor given by $\epsilon=\tau_{\mathrm{c}}/\tau_{\mathrm{\nu}}$ increases and finally exceeds unity in the C20 case, as described in Fig. 5(a). As already mentioned in Section. II, high fluid compressibility may result in VACF backtracking owing to sound propagation. In fact, a negative contribution of VACF obtained from MD simulations can be observed in Fig. 5(a), whereas results of hydrodynamic descriptions obtained using Eqs. (6) + (7) and Eq. (8) never yield negative VACF values. In general, the backtracking effect requires a sufficiently small value of bulk viscosity compared to shear viscosity Felderhof (2005). However, MD simulations provide finite values of bulk viscosity $\eta_{\mathrm{v}}$, thereby giving rise to an attenuation of the sound propagation in accordance with Eq. (9).

Finally, VACF results obtained for the WCA potential system have been plotted in Fig. 6, which demonstrates the overall VACF agreement between MD simulations and analytical expressions of hydrodynamics. That is, both the sound propagation effect due to fluid compressibility and long-time tail caused by kinematic viscosity can be thoroughly emulated at the molecular interaction level. As depicted in Fig. 4(b), the hydrodynamic radius $a_{\mathrm{HD}}$ approximately equals that of bare fullerene $a$, and this coincidence is in line with the agreement between MD and hydrodynamics observed in Fig. 6. However, small deviations were observed for time scales around $\tau_{\nu}$, especially with regard to smaller fullerene particles, C20 and C60, the hydrodynamic radius of which $a_{\mathrm{HD}}$ slightly exceeds $a$. Furthermore, over shorter time scales, VACF values were observed to be well described by the exponential decay predicted using the Enskog theory; this observation agrees well with that corresponding to the LJ potential system.

V Conclusions and Final Remarks

By performing MD simulations, thermodynamic and hydrodynamic properties of a single fullerene particle suspended in Ar fluids have been investigated in this study. The solvation free energy and the VACF were calculated to reveal the hydrodynamic behavior of said particles from the viewpoint of the molecular interaction by using LJ and WCA potentials.

As observed, the solvation free energy $\Delta G$ demonstrated the strong dependence on the intermolecular potential. As regards LJ potential, the attraction energy between fullerene and solvent particles was observed to overwhelms the entropy loss owing to the exclusion of solvent particles, contributing to more negative value of $\Delta G$ for larger fullerene particle. Correspondingly, the solvation was highly structured around the fullerene particle, as observed in RDF. In contrast, $\Delta G$ was observed to become positive with regard to the WCA potential, only utilizing the short-range, repulsive part of the LJ potential.

The hydrodynamic radius $a_{\mathrm{HD}}$ was quantified from the SE relation using the shear viscosity of the pure solvent and the diffusion constant of the fullerene particles. Remarkably, $a_{\mathrm{HD}}$ of LJ potential was observed to exceed the bare size of fullerene $a$, whereas the comparable relationship between $a_{\mathrm{HD}}$ and $a$ was observed with regard to the WCA potential. This difference of $a_{\mathrm{HD}}$ could be attributed to the strength of solvation quantified by $\Delta G$. There still exists a difference between $a_{\mathrm{HD}}$ and $a$ of an order of a molecular length scale corresponding to several solvation shells. When the difference $a_{\mathrm{HD}}-a$ remains at the molecular level, the ratio $a_{\mathrm{HD}}/a$ converges to unity for macroscopic values of $a$ (or $a_{\mathrm{HD}}$), ensuring that the continuum description remains valid. We also note the direct evidence of the stick boundary condition cannot be directly assessed from our MD simulations. It is still natural to assume the stick boundary condition at the macroscopic regime, where those hydrodynamics descriptions become valid. Furthermore, it is speculated that $\Delta G$ and $a_{\mathrm{HD}}$ will depend on the examined thermodynamic condition by changing density and temperature at molecular scales.

VACF results obtained from MD simulations were directly compared against those obtained using analytical expressions based on the generalized Langevin equation and hydrodynamics involving shear and bulk viscosities as well as the speed of sound of pure solvents. As observed, VACF decay demonstrates a long-time tail $t^{-3/2}$, which is purely governed by the kinematic viscosity $\nu$ for time scales larger compared to the kinematic time $\tau_{\nu}$. For time scales shorter than $\tau_{\nu}$, the sound propagation effect is expected to be observed in the VACF. However, VACF for the LJ potential could not be appropriately predicted using the hydrodynamic description, albeit the NS equation of compressible fluids was employed. In contrast, VACF for the WCA potential system was observed to be in more accord with the corresponding hydrodynamic description even at approximately the sound propagation time $\tau_{\mathrm{c}}$, particularly for larger fullerene particles. Note that the origin of the difference of the VACF results between MD simulations and hydrodynamics remains elusive. To address this concern, it is essential to include not only sound propagation but also frequency- and wave-number dependent formalism related to viscoelastic properties of host fluids Grimm et al. (2011); Puertas and Voigtmann (2014).

In summary, the proposed study demonstrated the impact of the intermolecular interaction on the hydrodynamic behavior in the Brownian motion in all-atom MD simulations. For a real colloidal particle measuring a radius 1 $\mu$m, it is still difficult to simulate macroscopic hydrodynamics with molecular descriptions via MD simulations. In contrast, the proposed simulation system involving nanoscale fullerene particles enabled to resolve time scales up to the microscopic level. In particular, results obtained from MD simulations performed in this study were observed to bridge hierarchical time scales, the Enskog time $\xi_{\mathrm{E}}^{-1}$, the sound propagation time $\tau_{\mathrm{c}}$, and the kinematic time $\tau_{\nu}$, and the diffusion time $\tau_{\mathrm{D}}$.

Acknowledgements.

The authors thank Rei Tatsumi, Takuya Iwashita, Hideyuki Mizuno, and Kazuo Yamada for helpful discussions. This work was supported by JSPS KAKENHI Grant Numbers, JP17J01006 (Y.I.), JP18H01188 (K.K.), and JP15K13550 (N.M.). This work was also supported in part by the Post-K Supercomputing Project and the Elements Strategy Initiative for Catalysts and Batteries from the Ministry of Education, Culture, Sports, Science, and Technology. Y. I. is supported by the JSPS fellowship. The numerical calculations were performed at Research Center of Computational Science, Okazaki Research Facilities, National Institutes of Natural Sciences, Japan.

Appendix A VACF of pure LJ and WCA fluids

Figure 7 demonstrates VACF of pure LJ and WCA fluids using MD simulations. Similarly to Figs. 5 and 6, VACF values were compared against those obtained using analytical expressions described in Eqs. (6) + (7) and Eq. (8). Furthermore, values of the exponential decay were plotted in accordance with the Enskog theory, $\exp(-\xi_{\mathrm{E}}t)$. Here, the Enskog friction coefficient was given by

[TABLE]

where $g(r)$ and $r_{\mathrm{max}}$ represent the RDF and its first peak position within the system, respectively Hansen and McDonald (2013). Hydrodynamic radii were quantified as $a_{\mathrm{HD}}\approx 0.348\sigma_{\mathrm{Ar}}$ (LJ) and $0.346\sigma_{\mathrm{Ar}}$ (WCA), respectively, values of which were obtained from the SE relation involving the diffusion constant and shear viscosity.

As demonstrated in Fig. 7, the long-time tail is perfectly characterized through use of the hydrodynamic description, $C(t)\sim(2k_{\mathrm{B}}T/3\rho_{\mathrm{f}})(4\pi\nu t)^{-3/2}$. This is true for cases involving both LJ and WCA potentials. However, analytical expressions described in Eqs. (6)

• (7) and Eq. (8) demonstrate little ability to reproduce the MD results over a short time regime. Alternatively, the Enskog theory demonstrates better agreement with results obtained using MD simulations for both LJ and WCA potentials, as illustrated in Fig. 7.