1. Introduction
The Stokes–Einstein (SE) relation [1] $D={k}_{{\rm{B}}}T/C\eta r$ combines the Einstein relation $D={k}_{{\rm{B}}}T/\alpha $ and Stokes' formula $\alpha =C\eta r,$ which correlates the diffusion coefficient D and the frictional coefficient $\alpha $ or the viscosity $\eta ,$ where kB is the Boltzmann constant, T is the temperature, $r$ is the effective hydrodynamic radius and is equal to the radius for a hard sphere particle, C is a constant depending on the boundary condition, and which is 6π for no-slip boundary condition and 4π for the slip [2].
The SE relation has been successfully applied to many situations, such as the protein diffusion [3] and the oxygen transport [4] in solutions. However, it is proposed to be invalid for liquid undergoes supercooling [5–7]. Instead of the original formula $D={k}_{{\rm{B}}}T/C\eta r,$ the three variants, $D\sim T/\eta $ [5, 8–11], $D\sim {\tau }^{-1}$ [7, 12–17], and $D\sim T/\tau $ [6, 18–20], are usually adopted to test the SE relation, where τ is the structural relaxation time. The variant $D\sim T/\eta $ is based on the assumption that $r$ is a constant. Both the variant $D\sim {\tau }^{-1}$ and $D\sim T/\tau $ are based on $D\sim T/\eta .$ The $D\sim {\tau }^{-1}$ is based on a further assumption that τ has a similar temperature dependence as $\eta /T.$ The $D\sim T/\tau $ depends on another approximate relation $\eta ={G}_{\infty }\tau ,$ where G∞ is the instantaneous shear modulus and is presumed to be a constant. Therefore the three variants are all based on the assumption that $r$ is constant under different conditions. Its rationality is important to the conclusion drawn on the validity or not of the SE relation.
However, there exists some evidences show that $r$ is varied with conditions. A classical case is the deviations from SE relation for alkali ions in aqueous solutions [21, 22]. It is shown that the dielectric polarization leads $r$ to deviate from the crystallographic radius. The $r$ is almost proportional to the volume fraction in a diluted organic aqueous solution [23]. The gyration radius often adopted to estimate $r$ in aqueous macromolecule solutions is also dependent on conditions [24]. Another example is the smaller and lighter ion has a smaller D than the larger and heavier ion in ionic systems at a certain temperature, such as in NaCl aqueous solution [25], [EMI+][NO3−] [26] and [EMI+][BF4−] [27] ionic liquids at room temperature. The larger and heavier ion at a lower temperature may even have a greater D than the smaller and lighter ion at a higher temperature. The phenomena is different from the estimation according to the SE relation and the gas kinetic theory [28]. If the SE relation is valid, $r$ should be changed with conditions, and otherwise the SE relation is breakdown.
The SE relation is a result of the combination of Einstein relation and Stokes's law. Einstein relation as a special case of the fluctuation-dissipation theorem and which is established in equilibrium and near equilibrium state. The supercooled liquids are still in equilibrium and the Einstein relation should be valid [29]. Moreover, the viscosity increases with decreasing temperature and the Reynolds number is decreased, so the Stokes's law applicable for liquids at the normal temperature should be also valid for supercooled liquids. In addition, the Stokes's law is derived for hard sphere particle in fluids [2], there is no reasons that $r$ is always or directly equal to the radius of soft sphere under different conditions.
Based on the above facts, we have enough reasons to suspect the assumption especially when conditions have large changes, such as liquid undergoes supercooling. Therefore it is necessary to examine the rationality of the assumption. Because $r$ is directly equal to the radius of hard sphere [2] and the variant $D\sim T/\eta $ is observed to be breakdown in the two classical supercooled liquid models, Kob–Andersen model and ortho-terphenyl (OTP) [6], in this work, we performed molecular dynamics (MD) simulations with the truncated Lennard-Jones-like (TLJ) model with different rigidity, OTP and Kob–Andersen model to test the rationality of the assumption and to examine its influence on the validity of SE relation.
2. Simulation details and analysis methods
2.1. Simulation details
The TLJ model potential [30] is given by$\sigma =0.34\,{\rm{nm}},$ $\varepsilon /{k}_{{\rm{B}}}=120\,{\rm{K}},$ mass of molecule is $m=39.95\,{\rm{g}}\,{{\rm{mol}}}^{-1}.$ The positive parameter n determines the degree of softness of the pair potential, and the limit $n\to \infty $ corresponds to the hard sphere system. As $r=0.5\sigma $ for hard sphere, we have chosen eight ns within 1.0–6.0 to explore the variations of $r$ with temperature and rigidity. The system contains 8192 particles in a cubic box with size 6.802 nm.
$\begin{eqnarray}u\left(r\right)=\left\{\begin{array}{l}\varepsilon \left[{\left(\displaystyle \frac{\sigma }{r}\right)}^{2n}-2{\left(\displaystyle \frac{\sigma }{r}\right)}^{n}+1\right],\,0\lt r\lt \sigma ,\\ 0,\,r\geqslant \sigma .\end{array}\right.\end{eqnarray}$
We choose We adopt Lewis–Wahnstrom model of OTP [31], each phenyl ring is represented by a Lennard-Jones site ($\varepsilon =5.276\,{\rm{kJ}}\,{{\rm{mol}}}^{-1},\,\sigma =0.483\,{\rm{nm}}$ ), and the three sites form a rigid isosceles triangle with an angle of 75° and the bond length is 0.483 nm. Each site has a mass $m=78\,{\rm{g}}\,{{\rm{mol}}}^{-1}.$ The system contains 3072 molecules with a constant density 1.0746 g cm−3.
The Kob–Andersen system [32] contains a binary (80:20) mixture of 8000 Lennard-Jones particles consisting of two species of particles, A and B, in a cubic box with size 6.392 nm. The interaction between two particles of type $\alpha ,\,\beta \in \left\{\left.A,\,B\right\}\right.$ is given by${\varepsilon }_{AB}/{\varepsilon }_{AA}=1.5,$ ${\varepsilon }_{BB}/{\varepsilon }_{AA}=0.5,$ ${\sigma }_{{AB}}/{\sigma }_{{AA}}=0.8,$ and ${\sigma }_{{BB}}/{\sigma }_{{AA}}=0.88;$ we choose ${\varepsilon }_{AA}/{k}_{{\rm{B}}}=120\,{\rm{K}},$ ${\sigma }_{AA}=0.34\,{\rm{nm}}$ and ${m}_{A,B}=39.95\,{\rm{g}}\,{{\rm{mol}}}^{-1}.$
$\begin{eqnarray}{V}_{\alpha \beta }\left(r\right)=4{\varepsilon }_{\alpha \beta }\left[{\left(\displaystyle \frac{{\sigma }_{\alpha \beta }}{r}\right)}^{12}-{\left(\displaystyle \frac{{\sigma }_{\alpha \beta }}{r}\right)}^{6}\right],\end{eqnarray}$
where All simulations were performed in NVT ensemble with the GROMACS package [33, 34]. The simulated temperature range is within 7–138 K with 12 temperatures for TLJ model, 260–400 K with 13 temperatures for OTP and 66–500 K with 24 temperatures for Kob–Andersen model. The variant of SE relation is observed to be invalid in the chosen simulated temperature range for OTP and Kob–Andersen model [6]. The system temperature was kept a constant by the Nosé–Hoover thermostat [35, 36]. The periodic boundary conditions were applied to all three dimensions. The interactions were calculated directly with the cutoff of 0.35 nm, 2.0 nm and 0.85 nm for TLJ, OTP and Kob–Andersen model, respectively.
2.2. Analysis methods
The method proposed by Hess [37] is adopted to calculate $\eta $ for its reliability and fast convergence. An external force ${a}_{x}=A\cdot \,\cos \left(qz\right)$ with the maximum A is applied in the X direction, where $q=2\pi /l$ with l the box size. The shear viscosity is described by
$\begin{eqnarray}\eta =A\rho /V{q}^{2},\end{eqnarray}$
where V is the maximum of the velocity in the X direction and ρ is the density.The $\alpha $ is determined by introducing a small force ${f}_{e}$ to a part of particles in the linear response regime. The frictional force on an ion ${f}_{r}=\alpha v$ is equal to ${f}_{e}$ after reaching the non-equilibrium steady state. The frictional coefficient is thus determined by$\eta $ and $\alpha ,$ the $r$ is calculated by$C=4\pi $ in this work corresponding to the slip boundary conditions, for which the calculated $r$ is equal to the radius of TLJ particle for low temperature or large n.
$\begin{eqnarray}\alpha =\displaystyle \frac{{f}_{e}}{v}.\end{eqnarray}$
In this work, 64 TLJ particles, 160 OTP molecules and 400 B particles are separately chosen for each model to keep enough statistical accuracy and avoid too much disturbance on the system. After getting the $\begin{eqnarray}r=\displaystyle \frac{\alpha }{C\eta },\end{eqnarray}$
where we choose The diffusion coefficient is calculated via its asymptotic relation with the mean square displacement by${\rm{\Delta }}{\boldsymbol{r}}\left(t\right)$ is particle position displacement and $\langle \rangle $ denotes ensemble average.
$\begin{eqnarray}D=\mathop{\mathrm{lim}}\limits_{t\to \infty }\displaystyle \frac{\left\langle {\rm{\Delta }}{{\boldsymbol{r}}}^{{\rm{2}}}\left(t\right)\right\rangle }{6t},\end{eqnarray}$
where 3. Results and discussion
To examine the assumption and its possible influence on the SE relation, the viscosity $\eta ,$ frictional coefficient $\alpha $ and diffusion constant D for TLJ model at different temperature T are calculated and plotted in figure 1. If the $r$ is a constant, the $\eta $ and $\alpha $ should have the same changes with T. However, the $\eta $ and $\alpha $ show different variation with increasing T. The $\alpha $ is increased with increasing T and n. The $\eta $ is also increased with increasing T for each n, but it is decreased with increasing n when n is within 1.0–2.0 at a certain temperature. The D is increased with increasing T for all n, and which is almost increased with increasing n at a certain temperature except for n = 6.
Figure 1. The |
The $r$ calculated by equation (5 ) is plotted in figure 2(a). It is not a constant but is increased as T decreases when n is within 1.0–2.0. Meanwhile, it is also increased with increasing n at a certain temperature. The $r$ is approaching 0.17 nm with decreasing T as well as the increasing of n, and is around 0.17 nm for n = 3.0 and 6.0 at all temperatures. The results indicate the $r$ is varied with conditions except for the large rigidity, and the assumption is not valid for the TLJ model. As point out in the Stokes' formula [2], the $r$ is equal to the radius of hard sphere. The minimum of T in our simulation is 7 K and the scale of the interaction is $\varepsilon /{k}_{{\rm{B}}}=120\,{\rm{K}}.$ The particle has no enough energy to penetrate the surrounding particles when n is large or T is small, and which looks like a hard sphere. On the contrary, it looks like a soft sphere with a smaller $r$ especially at a higher T and smaller n.
Figure 2. (a) The |
The SE relation $D\sim T/\alpha $ and its variant $D\sim T/\eta $ are tested with the data shown in figure 1 and the related results are plotted in figures 2(b) and (c). The data are rescaled by the value at T = 7 K. The data for $D\sim T/\alpha $ with different n are almost all fallen onto the same reference line $y=x,$ which indicate that the SE relation is definitely established for all n. However, the variant $D\sim T/\eta $ deviates from $y=x$ for all ns except for n = 3.0 and 6.0. The deviations are decreasing with increasing n as shown by figure 2(c), which is consistent with the $r$ shown in figure 2(c). The results indicate $D\sim T/\eta $ is only valid when the $r$ can be considered as a constant and otherwise is not equivalent to the SE relation $D\sim T/\alpha .$
The $\eta ,$ $\alpha $ and D for OTP and Kob–Andersen model are plotted in figure 3. Both $\eta $ and $\alpha $ are increased with decreasing T and show awfully different changes compared with the TLJ model. The D is also increased with increasing T similar as the TLJ model. To examine the assumption, we rescaled the $r$ by the value r = 0.402 nm at T = 400 K for the OTP and r =0.133 nm at T = 500 K for Kob–Andersen model, respectively. As the figure 4(a) shown, the rescaled effective hydrodynamic radius ${r}_{s}$ is not a constant for both OTP and Kob–Andersen model, and both are almost decreased with decreasing temperature. Similar results as shown by $r\sim T/D\eta $ or $D\sim {\left(\eta /T\right)}^{-\xi }$ with $\xi \lt 1$ are also observed in other supercooled liquids, including supercooled water [12], supercooled binary Lennard-Jones liquids [11, 38], water/methanol solutions [10], and tris-Naphthylbenzene [39]. However, the changes of $r$ show opposite trend compared with the TLJ model.
Figure 3. The |
Figure 4. (a) The rescaled effective hydrodynamic radius |
As discussed above, the $r$ of the TLJ model can be explained by the rigidity of the particle under different conditions. However, particles are interacting through Lennard-Jones potential in both OTP and Kob–Andersen model. Both repulsive and attractive interactions are present at the same time; and the negative attractive interaction usually plays a more important role. To give a unified picture, the trend of $r$ changing with temperature can be understood as follows. A molecule is usually not freely moving in liquid but partially drags the effective shells composed of surrounding molecules to move together. The first coordination shell usually plays a more important role [25, 40] and forms a composite particle along with the central molecule. The interaction between the central molecule and its first solvation shell is different under different conditions, which can be described by the coordination number n and a factor p. A larger p corresponds to a larger probability for molecules in the first shell to move together with the central molecule. Therefore the average number of molecules in the composite particle is $1+np,$ where 1 denotes the central molecule. By assuming both free molecule and the composite particle are spheres and adopting a mean field approach, the effective hydrodynamic radius of the composite particle is roughly ${\left[1+np\right]}^{1/3}{r}_{0},$ where ${r}_{0}$ is the effective hydrodynamic radius of a free molecule. Because the frictional force applied to the composite particle $\sim C\eta v{\left[1+np\right]}^{1/3}{r}_{0}$ should be equal to the sum of the frictional force applied to each molecule in the composited particle $\sim C\eta v\left[1+np\right]r,$ where v is the mean velocity of the composite particle. Therefore the average effective hydrodynamic radius $r$ of a molecule can be described by $r={r}_{0}/{\left(1+np\right)}^{2/3}.$ The p can be estimated by $p\,=1/\left({{\rm{e}}}^{{\rm{\Delta }}E/{k}_{{\rm{B}}}T}+1\right)$ if we roughly separate the molecules into two parts, in the shell and out of the shell, and ${\rm{\Delta }}E$ is the energy difference of molecule in the shell and out of the shell. Then the $r$ is $r={r}_{0}/{\left[1+n/\left({{\rm{e}}}^{{\rm{\Delta }}E/{k}_{{\rm{B}}}T}+1\right)\right]}^{2/3}.$
If ${r}_{0}$ and ${\rm{\Delta }}E$ are known, we can calculate the $r.$ Or if we know ${\rm{\Delta }}E,$ we can calculate the ratio of $r$ at different temperature. However, ${\rm{\Delta }}E$ is not a constant as the conditions vary [41]. Its sign can be found from the correlation. Due to the pure repulsive interaction, TLJ molecules are more likely to stay out of the shell than in the shell; however, Kob–Andersen and OTP molecules are more likely to stay in the shell because of the attractive interaction. Therefore ${\rm{\Delta }}E\gt 0$ for TLJ model but ${\rm{\Delta }}E\lt 0$ for Kob–Andersen and OTP model. With formula $r={r}_{0}/{\left[1+n/\left({{\rm{e}}}^{{\rm{\Delta }}E/{k}_{{\rm{B}}}T}+1\right)\right]}^{2/3}$ and the sign of ${\rm{\Delta }}E,$ the trend of $r$ changing can be qualitatively explained. At lower temperatures or lager n, ${{\rm{e}}}^{{\rm{\Delta }}E/{k}_{{\rm{B}}}T}$ is large; TLJ molecules have a weaker correlation and behave like a free molecule. On the contrary, ${{\rm{e}}}^{{\rm{\Delta }}E/{k}_{{\rm{B}}}T}$ is small for Kob–Andersen and OTP model at lower temperatures; molecules have a stronger correlation and show a collective motion. The correlation is increased for TLJ model and is decreased for Kob–Andersen and OTP model as temperature increases. Therefore the trend of $r$ changing with temperature for TLJ model is reverse with Kob–Andersen and OTP model. In summary, the changes of $r$ is a collective effect, and which is caused by correlation between the central molecule and its surrounding shells.
The SE relation $D\sim T/\alpha $ and the variant $D\sim T/\eta $ for OTP and Kob–Andersen model are tested by $D\sim {\left(\alpha /T\right)}^{-\xi }$ and $D\sim {\left(\eta /T\right)}^{-\chi }.$ If $\xi $ or $\chi =1.0,$ the SE relation or the variant is established and otherwise invalid. The logarithm of D and $T/\alpha ,$ $T/\eta $ for OTP and Kob–Andersen are plotted in figure 4(b). The variant $D\sim T/\eta $ is breakdown and behaves as a fractional form $D\sim {\left(\eta /T\right)}^{-\chi }$ with $\chi \approx 0.9$ for the two liquids. The results are similar as previous studies [6] as well as the observed in other liquids such as supercooled water [8], aqueous NaCl solutions [42] and ionic liquids [43]. However, the SE relation $D\sim T/\alpha $ is definitely valid for the two liquids, because the exponent $\xi $ s in $D\sim {\left(\alpha /T\right)}^{-\xi }$ are so close to 1.0.
Combined the results given by the TLJ model, OTP and Kob–Andersen model, it indicates that the assumption of $r$ is constant is usually invalid except for some special case such as the particle having enough rigidity. The SE relation given by variant $D\sim T/\eta $ is only established when $r$ can be considered as a constant. The SE relation $D\sim T/\alpha $ is valid for OTP and Kob–Andersen model even in the supercooled region.
4. Conclusions
In summary, we have examined the assumption of $r$ being a constant by performing MD simulations with TLJ model, OTP and Kob–Andersen model as well as explored its influence on the SE relation. Our results indicate the assumption is usually invalid. The trend of $r$ changing is increased with decreasing T for TLJ model for small rigidity and is almost a constant for large rigidity. It is decreased with decreasing T for OTP and Kob–Andersen model. The changes of $r$ can be qualitatively explained by the collective effect caused by the correlation of the molecule with its surrounding shells. Molecules have negative correlations with surrounding shells to move together due to the pure repulsive interaction in the TLJ model, especially at low temperatures; however, positive correlations are present in OTP and Kob–Andersen model because of the attractive interaction. The different correlations lead reverse changes of $r$ with temperature comparing TLJ model with OTP and Kob–Andersen model. The SE relation given by variant $D\sim T/\eta $ is invalid for TLJ model, OTP and Kob–Andersen model, and which shows a fractional form for OTP and Kob–Andersen model. The SE relation is definitely established for the three models. Our simulations indicate that $r$ is an important parameter to the conclusion drawn on the validity of SE relation and we propose the assumption should be carefully evaluated when used to test the SE relation.
Our results only simulate three models, similar inconsistency may also appear in other supercooled liquids, especially for strongly correlated systems, such as ionic liquids and liquids at much lower temperatures. The temperature ranges of our MD simulations are limited to 260–400 K for OTP model and 66–500 K for Kob–Andersen model, respectively, it is still far from the glass transition point, the changes of $r$ needs to be further explored. Moreover, the same assumption is also adopted when testing the Stokes–Einstein–Debye relation for the molecular rotation; its rationality is still need to be justified. Our future work could be, by employing more computer resources, extending our simulations to low temperatures to see the trend of $r$ changing and explore the assumption in the Stokes–Einstein–Debye relation.