The effective hydrodynamic radius in the Stokes–Einstein relation is not a constant

Gan Ren (任淦)

School of Science,Civil Aviation Flight University of China,Guanghan 618307,China

Abstract Variants based on the assumption of effective hydrodynamic radius being a constant are usually adopted to test the Stokes–Einstein(SE)relation.The rationality of the assumption is examined by performing molecular dynamics simulations with the truncated Lennard-Jones-like (TLJ)model,Kob–Andersen model and ortho-terphenyl (OTP) model.The results indicate the assumption is generally not established except for special case.The effective hydrodynamic radius is observed to increase with decreasing temperature for TLJ model but is decreased for Kob–Andersen and OTP model;and which is almost a constant for TLJ particle with enough rigidity.The variant of SE relationD～T/ηis invalid for the three models except for the TLJ particle with enough rigidity.We propose similar inconsistency may be also existed in other liquids and the assumption should be critically evaluated when adopted to test the SE relation.

Keywords: Stokes–Einstein relation,effective hydrodynamic radius,Stokes’ formula

1.Introduction

The Stokes–Einstein(SE)relation[1]D=kBT/Cη rcombines the Einstein relationD=kBT/αand Stokes’ formulaα=Cη r,which correlates the diffusion coefficient D and the frictional coefficientαor the viscosityη,where kBis the Boltzmann constant,T is the temperature,ris 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 formulaD=kBT/Cη r,the three variants,D～T/η[5,8–11],D～τ-1[7,12–17],andD～T/τ[6,18–20],are usually adopted to test the SE relation,where τ is the structural relaxationtime.The variantD～T/ηis basedon the assumptionthatris a constant.Both the variantD～τ-1andD～T/τare based onD～T/η.TheD～τ-1is based on a further assumption that τ has a similar temperature dependence asη/T.TheD～T/τdepends on another approximate relationη=G∞τ,where G∞is the instantaneous shear modulus and is presumed to be a constant.Therefore the three variants are all based on the assumption thatris 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 thatris 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 leadsrto deviate from the crystallographic radius.Theris almost proportional to the volume fraction in a diluted organic aqueous solution[23].The gyration radius often adopted to estimaterin 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,rshould 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 thatris 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.Becauseris directly equal to the radius of hard sphere[2]and the variantD～T/η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

We chooseσ=0.34 nm,ε/kB=120 K,mass of molecule ism=39.95 g mol-1.The positive parameter n determines the degree of softness of the pair potential,and the limitn→∞corresponds to the hard sphere system.Asr=0.5σfor hard sphere,we have chosen eight ns within 1.0–6.0 to explore the variations ofrwith temperature and rigidity.The system contains 8192 particles in a cubic box with size 6.802 nm.

We adopt Lewis–Wahnstrom model of OTP [31],each phenyl ring is represented by a Lennard-Jones site(ε=5.276 kJ mol-1,σ=0.483 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 massm=78 g 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α,β∈{A,B}is given by

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.ax=A· cos(qz)with the maximum A is applied in the X

2.2.Analysis methods

The method proposed by Hess [37]is adopted to calculateηfor its reliability and fast convergence.An external force direction,whereq=2π/lwith l the box size.The shear viscosity is described by

where V is the maximum of the velocity in the X direction and ρ is the density.

Theαis determined by introducing a small forcefeto a part of particles in the linear response regime.The frictional force on an ionfr=αvis equal tofeafter reaching the nonequilibrium steady state.The frictional coefficient is thus determined by

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ηandα,theris calculated by

where we chooseC=4πin this work corresponding to the slip boundary conditions,for which the calculatedris equal to the radius of TLJ particle for low temperature or large n.

The diffusion coefficient is calculated via its asymptotic relation with the mean square displacement by

where Δr(t) is particle position displacement and〈〉denotes ensemble average.

3.Results and discussion

To examine the assumption and its possible influence on the SE relation,the viscosityη,frictional coefficientαand diffusion constant D for TLJ model at different temperature T are calculated and plotted in figure 1.If theris a constant,theηandαshould have the same changes with T.However,theηandαshow different variation with increasing T.Theαis increased with increasing T and n.Theη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η, α and D for TLJ model with different n as a function of T: (a)η versus T;(b)α versus T;(c) D versus T.

Figure 2.(a)Ther as a function of T for TLJ model with different n;the black dotted line is the reference line r = 0.17 nm.Verification of the validities of the SE relationD～T/αand its variantD～T/ηfor the TLJ model:(b)D～T/αand(c)D～T/η.The calculated data are represented by symbols and the dotted line is y=x.

Thercalculated 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.Theris 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 theris 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],theris equal to the radius of hard sphere.The minimum of T in our simulation is 7 K and the scale of the interaction isε/kB=120 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 smallerrespecially at a higher T and smaller n.

The SE relationD～T/αand its variantD～T/η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 forD～T/αwith different n are almost all fallen onto the same reference liney=x,which indicate that the SE relation is definitely established for all n.However,the variantD～T/ηdeviates fromy=xfor 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 thershown in figure 2(c).The results indicateD～T/ηis only valid when thercan be considered as a constant and otherwise is not equivalent to the SE relationD～T/α.

Theη,αand D for OTP and Kob–Andersen model are plotted in figure 3.Bothηandα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 therby 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 radiusrsis not a constant for both OTP and Kob–Andersen model,and both are almost decreased with decreasing temperature.Similar results as shown byr～T/Dη or D～(η/T)-ξwithξ＜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 ofrshow opposite trend compared with the TLJ model.

As discussed above,therof 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 ofrchanging 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 [1+np]1/3r0,wherer0is the effective hydrodynamic radius of a free molecule.Because the frictional force applied to the composite particle～Cη v[1+np]1/3r0should be equal to the sum of the frictional force applied to each molecule in the composited particle～Cη v[1+np]r,where v is the mean velocity of the composite particle.Therefore the average effective hydrodynamic radiusrof a molecule can be described byThe p can be estimated byp=if we roughly separate the molecules into two parts,in the shell and out of the shell,andΔEis the energy difference of molecule in the shell and out of the shell.Then theris

Figure 3.Theη,α and D as a function of T for OTP and Kob–Andersen model:(a)η versus T;(b)α versus T;(c)D versus T.The data for OTP is colored in black and the red is for Kob–Andersen model.

Figure 4.(a) The rescaled effective hydrodynamic radiusrs as a function of T for OTP and Kob–Andersen model.(b)Verification of the validities of the SE relationD～T/αand its variantD～T/η for OTP and Kob–Andersen model;the calculated data are represented by symbols and the solid lines are fitted by D～(α/T)-ξand D～(η/T)-χ,respectively.The data for OTP is colored in black and the red is for Kob–Andersen model.

Ifr0andΔEare known,we can calculate ther.Or if we knowΔE,we can calculate the ratio ofrat different temperature.However,ΔEis 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ΔE＞ 0 for TLJ model butΔE＜ 0 for Kob–Andersen and OTP model.With formulaand the sign ofΔE,the trend ofrchanging can be qualitatively explained.At lower temperatures or lager n,is large;TLJ molecules have a weaker correlation and behave like a free molecule.On the contrary,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 ofrchanging with temperature for TLJ model is reverse with Kob–Andersen and OTP model.In summary,the changes ofris a collective effect,and which is caused by correlation between the central molecule and its surrounding shells.

The SE relationD～T/αand the variantD～T/ηfor OTP and Kob–Andersen model are tested byD～(α/T)-ξandD～(η/T)-χ.Ifξorχ=1.0,the SE relation or the variant is established and otherwise invalid.The logarithm of D andTα,Tηfor OTP and Kob–Andersen are plotted in figure 4(b).The variantD～T/ηis breakdown and behaves as a fractional formD～(η/T)-χwithχ≈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 relationD～T/αis definitely valid for the two liquids,because the exponentξs inD～(α/T)-ξare so close to 1.0.

Combined the results given by the TLJ model,OTP and Kob–Andersen model,it indicates that the assumption ofris constant is usually invalid except for some special case such as the particle having enough rigidity.The SE relation given by variantD～T/ηis only established whenrcan be considered as a constant.The SE relationD～T/αis valid for OTP and Kob–Andersen model even in the supercooled region.

4.Conclusions

In summary,we have examined the assumption ofrbeing 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 ofrchanging 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 ofrcan 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 ofrwith temperature comparing TLJ model with OTP and Kob–Andersen model.The SE relation given by variantD～T/η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 thatris 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 ofrneeds 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 ofrchanging and explore the assumption in the Stokes–Einstein–Debye relation.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (No.12104502) and the Science Foundation of Civil Aviation Flight University of China(No.J2021-054).The computations of this work were conducted on the Tian-2 supercomputer,and the author thanks Yanting Wang (ITP CAS) for supporting.