Self-diffusion in single-component Yukawa fluids

It was suggested in the literature that the self-diffusion coefficient of simple fluids can be approximated as a ratio of the squared thermal velocity of the atoms to the"fluid Einstein frequency,"which can thus serve as a rough estimate of the friction (momentum transfer) rate in the dense fluid phase. In this article we test this suggestion using a single-component Yukawa fluid as a reference system. The available simulation data on self-diffusion in Yukawa fluids, complemented with new data for Yukawa melts (Yukawa fluids near the freezing phase transition), are carefully analyzed. It is shown that although not exact, this earlier suggestion nevertheless provides a very sensible way of normalization of the self-diffusion constant. Additionally, we demonstrate that certain quantitative properties of self-diffusion in Yukawa melts are also shared by systems like one-component plasma and liquid metals at freezing, providing support to an emerging dynamical freezing indicator for simple soft matter systems. The obtained results are also briefly discussed in the context of the theory of momentum transfer in complex (dusty) plasmas.


I. INTRODUCTION
De Gennes in his seminal paper on liquid dynamics and inelastic scattering of neutrons 1 related the self-diffusion coefficient in classical atomic liquids to liquid structure, measured in terms of the radial distribution function (RDF) g(r), and the potential of inter-atomic interactions φ(r). The relation he put forward is where v T = T /m is the thermal velocity, T is the temperature, m is the atomic mass, and the characteristic frequency Ω E is defined as where n is the liquid number density. This characteristic frequency is called the Einstein frequency, because it represents the oscillation frequency of a single particle in the fixed environment of surrounding particles, characterized by a given RDF. 2 de Gennes also pointed out that since the interaction potential was very poorly known for liquids with measured self-diffusion coefficients, a direct comparison was not possible at that time. 1 The very same characteristic frequency Ω E was also earlier identified as the friction rate ξ in liquids by Kirkwood, Buff, and Green 3 and by Collins and Raffel 4 and a) Electronic mail: Sergey.Khrapak@dlr.de thus related to the self-diffusion coefficient by means of the Einstein relation, ξD ≃ T /m. This obviously leads to Eq. (1) to within the accuracy of a prefactor of order unity.
These early results from the transport theory of liquids are not expected to be exact, but they nevertheless suggest a sensible way of normalization of the self-diffusion constant, which may be expected to provide a slowly varying quantity in a wide parameter regime. The usefulness of this normalization can be easily tested using present day computational capabilities. In this article we demonstrate this by a detailed investigation of the behavior of D E in one-component Yukawa fluids. Yukawa fluids are particularly suitable for such a test because the Einstein frequency Ω E is trivially related to an important thermodynamics property, excess internal energy, which is relatively well known in a wide parameter regime. In addition, results related to transport properties of single component Yukawa systems are available in the literature and can be used for extensive tests of earlier predictions. In view of existing previous works on transport in Yukawa systems, which also considered self-diffusion in Yukawa fluids, it is important to carefully formulate the motivation for presenting still another article on this topic. Our main purpose here can be divided into the following four objectives: (i) to test quantitatively de Gennes's prediction of Eq. (1); (ii) to present new detailed data for the self-diffusion of Yukawa melts; (iii) compare the results for Yukawa melts with available data on diffusion in liquid metals at the melting temperature; and (iv) provide simple explanation on the observed quantitative similarity, which is likely to work for other related systems. To the best of our knowledge, these issues have not been addressed in previous studies.
The article is organized as follows. In Sec. II, we provide basic information on the physics of single component Yukawa systems. In Sec. III, we outline the procedure of numerical simulations performed in this work. In Sec. IV we summarize our main results and discuss their implications. This is followed by our conclusion in Sec. V.

II. FORMULATION
Yukawa systems are characterized by the pairwise repulsive interaction potential between charged point-like particles of the form where Q is the particle charge, λ is the screening length, and r is the distance between two particles. The screening normally comes from the redistribution of electrical charges (electron and ions) in a plasma surrounding the particle. Fixed neutralizing background corresponds to the absence of screening and the pure Coulomb interaction potential, the limit known as the one-component plasma (OCP) model. Yukawa potential is considered as a reasonable starting point to model interactions in threedimensional isotropic complex (dusty) plasmas and colloidal dispersions, 5-8 although it is also well recognized that considerable deviations can occur, especially at long interparticle separations. [9][10][11][12][13][14][15][16] The phase state of Yukawa systems is conventionally described by the two dimensionless parameters, 17 which are the coupling parameter Γ = Q 2 /aT , and the screening parameter κ = a/λ, where a = (4πn/3) −1/3 is the Wigner-Seitz radius. The screening parameter κ determines the softness of the interparticle interaction. It can be varied from the very soft long-ranged Coulomb interaction at κ = 0 to the hard-sphere interaction limit at κ → ∞. In the context of complex plasmas and colloidal suspensions low-and moderate-κ regime is of particular interest and this regime will be addressed below.
Yukawa systems are conventionally called "strongly coupled" when the coupling parameter Γ is large, that is when the potential interaction energy stored in the system dominates over the kinetic (thermal) energy. In this case the so-called Yukawa fluid regime is realized. When Γ increases further (i.e. temperature decreases), Yukawa fluids can crystallize into the body-centered-cubic (bcc) or face-centered-cubic (fcc) lattice (bcc lattice is thermodynamically favorable at week screening). The critical coupling values at which the fluid-solid phase transition occurs are usually denoted Γ m , where the subscript "m" refers to melting. This critical coupling strongly depends on the strength of screening. Accurate values Γ m (κ) have been tabulated 17,18 and fitted by simple expressions. [19][20][21] If Yukawa fluid is quickly cooled down to even lower temperatures, a Yukawa glass may form. The glass-transition line appears to be almost parallel to the melting line in the extended region of the phase diagram. 22,23 For the Yukawa potential the Einstein frequency is trivially related to the internal excess energy, which is relatively well known analytically. This relation is the consequence of the identity ∆φ(r) = φ(r)/λ 2 and can be expressed as 24,25 where ω p = 4πQ 2 n/m is the plasma frequency and u ex is the ecxess energy per particle in units of temperature. Very accurate analytical expressions for the quantity u ex (κ, Γ) are available in the literature. [26][27][28][29] Particularly simple practical expressions are based on the Rosenfeld-Tarazona scaling, 30,31 which states that for a wide class of soft repulsive potentials the thermal component of the internal energy exhibits a quasi-universal scaling on approaching the fluid-solid phase transition. The expression for u ex used in this work is that from Ref. 26.
Self-diffusion in one-component Yukawa fluids has been extensively studied using molecular dynamics (MD) simulations by several authors. In particular, Ohta and Hamaguchi tabulated the values of the reduced quantity D/a 2 ω E in a wide range of κ and Γ corresponding to the fluid regime. 32 The characteristic Einstein frequency ω E corresponds to an ideal crystalline lattice and is thus different from that used in this paper [Eq.
(2)]. Daligault investigated the diffusion properties of one-component plasmas and binary ionic mixtures from the weakly to the strongly coupled regimes. 33 He then presented a practical interpolation formula for the self-diffusion coefficient in Yukawa one-component plasmas valid for a wide range of inverse screening lengths and over the entire fluid region. 34 This practical formula operates with the reduced quantity D/a 2 ω p . Rosenfeld used one more normalization D R = Dn 1/3 /v T to demonstrate that D R (and also similarly reduced shear viscosity) exhibits quasi-universal behavior as a function of reduced excess entropy for various simple model systems, 35 including Yukawa fluids. 31 The usefulness of the "fluid Einstein frequency" Ω E to normalize self-diffusion data was not explored previously.

III. METHODS
In this paper we use the data for the self-diffusion coefficient in Yukawa fluids tabulated by Ohta and Hamaguchi in a wide range of κ and Γ values corresponding to the fluid state. 32 In addition, we performed complementary molecular dynamics (MD) simulations for Yukawa fluids at the melting temperature (Yukawa melts). The melting temperatures or, equivalently, coupling parameters Γ m are those tabulated for different κ in Ref. 17. Our MD simulations were performed on graphics processing unit (NVIDIA GeForce GTX 960) using the HOOMDblue software. [36][37][38] We used N = 50653 particles interact-ing according to the pairwise potential (4) and placed in a cubic box with periodic boundary conditions. There was no need to introduce hard cores since the potential is repulsive and diverges at short distances. At long distances the potential was cut off at L cut ≃ 14.5λ. The numerical time step was set to ∆t = 5 × 10 −3 ma 3 /Q 2 . Simulations were performed in the canonical ensemble (N V T ) with the Langevin thermostat at a temperature corresponding to the desired target coupling parameter Γ. The dissipation rate used was small enough to ensure that the system dynamics is not affected by dissipation (i.e. we are dealing with the Newtonian dynamics limit, opposite to the over-damped Brownian limit). 20,39,40 Equilibrium of the system was first achieved by running the simulation for at least two million time steps (depending on the value of κ we also used longer equilibration runs). Then, the particle positions and velocities were saved for 300Ω −1 E with a resolution of (1/20)Ω −1 E to get sufficient statistics.
From MD simulation with Yukawa melts the following information was collected: We calculated the RDF g(r) and the normalized velocity autocorrelation function (VAF) as well as the mean squared displacement (MSD), MSD(t) = ∆r 2 (t) .
Here v(t) is the velocity of a given particle at a time t, ∆r = r(t) − r(0) is its displacement from the initial position, and · · · denotes an average over all particles. The time-dependent diffusion coefficients can be evaluated employing either MSD: D(t) = MSD(t)/6t, or VAF: D(t) = v 2 T t 0 Z(t ′ )dt ′ (Green-Kubo relation). These definitions are not equivalent, but in the limit t → ∞ they should approach one single value, the conventional self-diffusion constant D.

IV. RESULTS AND DISCUSSION
The reduced diffusion coefficient, D E , as a function of the reduced coupling parameter Γ/Γ m = T m /T is shown in Fig. 1 First, we observe that the data points for the dependence of D E on Γ/Γ m are grouping around a single universal curve. No systematic dependence on κ (that is on the interaction potential softness) is observed. This should be expected in view of the isomorph theory of systems with hidden scale invariance. 42 isomorphs, too. Application of the isomorph theory to Yukawa systems has been recently discussed in detail. 44 Second, we see that the quantity D E decreases from ≃ 3 to ≃ 0.6 upon increase in coupling from 0.01Γ m to Γ m . This should be compared with almost three orders of magnitude variation of the quantities D/ω E a 2 and D/ω p a 2 , 32,34 and thirty times variation in D R in a similar range of Γ. 31 Thus, we observe a rather weak dependence of D E on the relative coupling strength Γ/Γ m . This suggests a very rough estimate for the strongly coupled fluid regime, which is accurate to within a factor of two in an extended region of the phase diagram. Third, more careful examination of the data plotted in Fig. 1 identifies three regimes of self-diffusion. In the weakly coupled regime, Γ 0.03Γ m the diffusion coefficient slowly decreases with coupling. This is the regime where the binary collision approximation for a gaseouslike phase should be appropriate (we did not test this quantitatively, however). Then we observe a plateaulike behavior in a wide regime of coupling strength 0.03Γ m Γ 0.3Γ m , corresponding to the fluid phase, where the reduced diffusion constant is almost constant, D E ≃ 2. Finally, on approaching the fluid-solid phase transition the reduced diffusion constant roughly scales as D E ∝ (Γ m /Γ) and reaches the value of approximately D E ≃ 0.6 at freezing (Γ = Γ m ).
In Figures 2 and 3 we show the results related to Yukawa melts obtained in this work. In figures 2 dt ′ , the inset shows the corresponding velocity autocorrelation functions, Z(t). In both cases the diffusion coefficient is expressed using the Rosenfeld normalization, DR = Dn 1/3 /vT . time-dependent diffusion coefficients calculated either via MSD (a) or Green-Kubo relation (b) are plotted. The diffusion coefficients are reduced according to the Rosenfeld prescription, D R = Dn 1/3 /v T . It is observed that they tend to D R ≃ 0.03 as t → ∞, similar to several other model systems near the fluid-solid phase transition studied recently (Hertzian, Gaussian-core, and inversepower-law interactions). 39 The convergence is somewhat better when using Green-Kubo formula [see figure 2 (b)]. Interestingly, the time-dependent self-diffusion coefficient evaluated at a time point at which the VAF crosses zero [at this point D(t) reaches its first maximum] is roughly twice the actual diffusion constant. Figure 3 shows essentially the same information as Fig. 2 (b), but with a different normalization, D E from Eq. (3). In this case the self-diffusion coefficients approach D E ≃ 0.6 as t → ∞. Somewhat more scattering in D E values is observed compared to D R [compare Fig. 2 (b) and Fig. 3]. The oscillations in time-dependent diffusion coefficients are related to the oscillations in the velocity autocorrelation function, which are characteristic for the caged short-time dynamics of particles in the fluid state. The values of differently normalized self-diffusion constants at freezing can be related using the following simple empirical consideration. Near the melting point of a crystalline lattice, an average deviation of particles from their lattice sites is roughly related to the temperature via The Einstein frequencies of a dense fluid near freezing and of a crystalline solid near melting (at the melting temperature) are close. According to the Lindemann melting rule we expect at melting δr 2 ∼ L 2 ∆ 2 , where L ∼ 0.1 is the Lindemann parameter and ∆ = n −1/3 is the mean interparticle separation. Using this to estimate the Einstein frequency at the melting temperature we obtain This is exactly what we have observed for Yukawa melts.
Since Yukawa interaction has not been assumed in this consideration, neither explicitly, nor implicitly it would make sense to test whether the derived relations above can operate in real systems as well. We have estimated the values of the reduced self-diffusion coefficient of some liquid metals at the freezing point using the data summarized by March and Tosi 45 (see Table 5.2 of their book). These results are shown in Table I. Except the two extreme cases of Cu (D R ≃ 0.040) and Zn (D R ≃ 0.027) all the other data are scattered in a rather narrow range (D R ≃ 0.032±0.03). Moreover, using a more recent value for the self-diffusion in Cu melt (3.27 × 10 −5 cm 2 s −1 according to Ref. 46) we obtain D R ≃ 0.033, which also belongs to the specified range. Similar relationships for the reduced coefficient of self-diffusion in liquid metals at the melting temperature were proposed in Refs. 47 and 48 on the basis of quite different arguments. In case of extremely soft Coulomb interaction, corresponding to the OCP limit, the self-diffusion coefficient at the fluidsolid phase transition is D R ≃ 0.031. This value has been obtained using the fitting formula by Daligault,34 adopting Γ m ≃ 172 at the fluid-solid phase transition. 17 The hard-sphere result for the self-diffusion near freezing is somewhat smaller D R ≃ 0.02, 49 possibly signaling the crossover between soft sphere and hard sphere regimes of collective motion and transport. 50 All this implies that certain tendencies identified in this work are unlikely specific to the Yukawa interaction potential, but can be applicable to other model and real systems, provided the interactions are soft enough.
To conclude the discussion, we point out one more potential application of the obtained results. The discovered dependence of D E on relative coupling strength, combined with the elementary expression for the diffusion coefficient D = T /mξ, allows us to estimate the characteristic momentum transfer frequency in Yukawa fluids, We get in the simplest rough approximation ξ ∼ Ω E , improvements can be made by making use of the data plotted in Fig. 1. This can be relevant for the momentum transfer frequency in particle-particle collisions in complex (dusty) plasmas. Previously, a unified theory of momentum transfer in complex plasmas has been developed using the binary collision approach, assuming the interaction potential between charged species of the Yukawa type. [51][52][53][54] This is normally adequate for electron-particle and ion-particle collisions, but the binary collision approximation for particle-particle collisions can only be relevant for extremely rarefied particle clouds with weak electrostatic coupling between them. The results presented in this article suggest a useful method to evaluate momentum transfer rates beyond the weak coupling regime.

V. CONCLUSION
To summarize, we have tested de Gennes approximation for the self-diffusion coefficient in simple fluids [Eq. (1)] using the data for self-diffusion in onecomponent Yukawa systems. Equation (1) is shown to be inexact, but nevertheless useful. It turns out to be accurate to within a factor of two in the entire strongly coupled fluid regime.
From the simulations of Yukawa melts we have observed that the properly normalized self-diffusion coefficient is approximately constant. The value D R ≃ 0.03 is consistent with the values reported for other simple model systems at freezing, such as Hertzian, Gaussiancore, and inverse-power-law fluids. 39 It is also consistent with the self-diffusion coefficient measured in liquid metals at the melting temperature, giving strong support to the emerging dynamic freezing criterion for simple atomistic systems. This can potentially complement the well-known empirical dynamic freezing criterion for over-damped colloidal fluids. 55 Simple explanation for the observed universality of the self-diffusion coefficient has been suggested.
Finally, it has been pointed out that the obtained results can be of some use in estimating the momentum transfer rates in particle-particle collisions in complex (dusty) plasmas. The existing approaches have been often limited to the binary collision approximation, only appropriate for the weakly coupled gaseous regime. The observed trends in self-diffusion can be useful in estimating momentum transfer rates in the strongly coupled fluid regime. The reported results can also be of some use in the context of self-diffusion in strongly coupled plasmas, experimental measurements of which have been recently reported. 56