Stokes-Einstein relation and excess entropy scaling law in liquid Copper

We report an ab initio study of structural and dynamic properties of liquid copper as a function of temperature. In particular, we have evaluated the temperature dependence of the self-diffusion coefficient from the velocity autocorrelation function as well the temperature dependence of the viscosity from the transverse current correlation function. We show that LDA based results are in close agreement with experimental data for both the self-diffusion coefficient and the viscosity over the temperature range investigated. Our findings are then used to test empirical approaches like the Stokes-Einstein relation and the excess entropy scaling law widely used in the literature. We show that the Stokes-Einstein relation is valid for the liquid phase and that the excess entropy scaling law proposed by Dzugutov is legitimate only if a self-consistent method for determining the packing fraction of the hard sphere reference liquid is used within the Carnahan-Starling approach to express the excess entropy.


Introduction
The study of the relationship between structural and dynamic properties of metallic liquids is a longstanding open problem in condensed matter physics. For instance, transport properties such as the shear viscosity (η) and the diffusion coefficients (D) of liquid metals are key parameters in the study of the crystal nucleation and growth in metallic melts [1,2].The temperature dependences of the viscosity and the relaxation time also play a very important role in studying the liquid-glass transition in a glass-forming system. Therefore, finding such a link should be also of practical interest as it relates a quantity like the self-diffusion coefficient that can be difficult to measure accurately to an experimentally more accessible dynamic, thermodynamic or structural quantity.
One alternative is to approximate diffusion coefficients from viscosities through the Stokes-Einstein (SE) relation [3], D SE = k B T /2πRη, where k B is Boltzmann's constant and R is an effective diameter of particles. The interest of this relation lies on the fact that it is much more difficult to measure the diffusion coefficients of liquids than their viscosities since the influences of convective flow on the diffusion profile during annealing can be important using the long-capillary (LC) technique and its variation. At the same time, viscosity data are often available from standard data collections with a weak scatter between measurements made using different techniques. For instance, it is only ±6% for copper [4]. First derived for the motion of a macroscopic particle in a homogeneous viscous medium, its use for atomic diffusion in liquid alloys requires defining R at the atomistic scale. It is well accepted in the common literature [3] to use the first position in the pair-correlation function obtained either experimentally or by means of atomistic simulations.
its validity even for complex alloys [8] while others evidence its failure [9,10], indicating the importance of the chemical composition of the melt. Moreover, it seems that the temperature also plays an important role in its use [10]. These experiments raise many questions and among them the choice of the effective diameter of particles as well as in a more general aspect, the relation of transport properties with the melt structure.
Another alternative is to relate the dynamic properties to the structure and thermodynamics of melts using the universal scaling relationship like that proposed by Dzugutov [11]. Within this scheme, the pair excess entropy of a liquid, S 2 , is related to a dimensionless form of the diffusion coefficient, D * . The quantity S 2 can be easily computed from pair correlation functions. However, if this approach is legitimate for model fluids like Lennard-Jones or Hard-Sphere (HS) fluids, a coupled Monte Carlo/molecular dynamics study has shown that this two-body correlation approximation does not hold for liquid metals and alloys modelled by EAM potentials [12].
More generally, molecular-dynamics (MD) simulations represent a powerful means for getting a microscopic picture of structural and dynamic properties in the liquid state and for studying their relationships but such studies are fruitful at the condition that atomic interactions are accurately modelled.
For liquid Cu, structural and dynamic properties have been extensively studied using the semi-empirical embedded-atom model (EAM) formalism [13][14][15][16][17][18], as well as by the orbital free ab initio method [19]. However, if structural properties are in very good agreement with experimental data whatever the EAM potential, diffusion coefficients reported in classical MD simulations show an important scatter, depending on the specific implementation of EAM potentials. Another difficulty holds with the fact that the calculated diffusion coefficients are mainly referred to experimental LC data [20], with the consequence that they exhibit a scatter of about ±50% around the most recent QNS values [7]. Only simulated results provided by Alemany et al. [14] disagree with LC data and these authors suggest that the LC values may be erroneous. Let us mention that very recent MD simulations [18] refer to experimental QNS values and give a better agreement with them.
All these results show that it is necessary to revisit the determination of the structural and dynamic properties of liquid Cu using ab initio molecular dynamics simulations.
In a first step, we show that the determination of the self-diffusion coefficient as well as the viscosity depend on the exchange-correlation functional used in DFT calculations, namely the local-density approximation (LDA) and the generalized gradient approximation (GGA). Our study shows that GGA based calculations of the diffusion coefficient underestimates both LDA based values and experimental data while we obtain the opposite trend for the viscosity. We evidence that these results can be explained by a more important cage effect using GGA as compared to LDA, due to the increase of the icosahedral motifs.
In a second step, our results are used to discuss the Stokes-Einstein relation and its validity as a function of temperature. Finally, we test the universal scaling law relating the diffusion coefficient and the excess entropy of a liquid as proposed by Dzugutov [11]. We show that the scaling law is obeyed provided a self-consistent method for determining the packing fraction of the hard sphere reference fluid is used within the Carnahan-Starling approach to express the excess entropy. In order to strengthen the present discussion, we shall present our results together with those recently published for liquid Al [23].

Simulation details
The AIMD simulations have been performed by using the Vienna ab initio simulation package [24] using the projected augmented-wave method to describe the electron-ion interaction [25]. A plane-wave basis set with an energy cutoff of 273 eV is used. The exchange-correlation energy is described using either the local-density approximation [26,27] or the generalized gradient approximation in the Perdew-Burke-Ernzerhof form [28]. Only the Γ-point sampling is considered to represent the supercell Brillouin zone.

43603-2
Equations of motion are solved using Verlet's algorithm in the velocity form with a time step of 1 fs. Atomic motions are carried out in the NVT ensemble by means of a Nose thermostat to control temperature. For each temperature, we use a cubic cell with periodic boundary conditions containing 256 atoms. The liquid sample is first equilibrated at 2000 K (well above the experimental melting temperature of 1358 K) and then cooled to 1800, 1600, and 1398 K, successively with a rate of 10 13 K/s. At each temperature, we adjust the volume V of the supercell to reproduce the experimental densities [4] and after an equilibration of 10 ps, the run is continued further during 80 ps for production of the structural and dynamic quantities.
A number of 2000 configurations are collected to calculate the static structure factor, the pair-correlation function, as well as the diffusion coefficient and the viscosity. Among these configurations, we select ten independent ones, to compute their inherent structures [29] for the purpose of studying the local ordering. This is done numerically by carrying out a conjugate gradient energy-minimization in order to suppress the kinetic energy.

Self-diffusion coefficient
In a first step, we consider the single-atom dynamics through the velocity auto-correlation function, and its time integration that gives access to the self-diffusion coefficient, D [4]. Let us mention that this method requires shorter simulation times than that based on the mean-square displacement to have statistically meaningful results [4,30].
We report the temperature dependence of D using the LDA and GGA approximations in figure 1.
We can see that LDA calculations are close to the experimental QNS values [7] while GGA calculations underestimate LDA calculations by about 20% above the melting point. We also compute an activation energy assuming an Arrhenius-type behaviour for the diffusion process in liquid copper. As shown in figure 1, the data are well fitted in the whole range of temperatures and the derived activation energy is E D = 370 ± 5 and 475 ± 5 meV respectively for the LDA and GGA approximations. Let us mention that the LDA value is in close agreement with the experimental value, namely 337 ± 5 meV [7].  [7] along with their Arrhenius fit (dash-dotted lines). AIMD results of references [21] (full circle) and [22] (triangle up) are also plotted for comparison.

Viscosity
To compute the viscosity, we use a direct method based on the transverse current correlation function (see reference [31] for a detailed description of this technique), which has the advantage of yielding a generalized q-dependent shear viscosity from which the hydrodynamic limit can be evaluated [32,33].
In figure 2 we compare LDA and GGA results with the assessed values [3]. As for the temperature dependence of the self-diffusion coefficient, we observe that the temperature dependence of the viscosity obtained by LDA is close to the temperature dependence of the assessed viscosity while GGA values overestimate LDA and the assessed values. More particularly, LDA values are consistent with experimental data of Kehr et al. [34] and Brillo et al. [35], both experimental sets being used in the assessment.
Both computed self-diffusion coefficients and viscosities are then used to test the validity of the Stokes-Einstein relation. In figure 2, we report viscosities computed via the Stokes-Einstein relation using R given by the position of the first peak in the computed pair-correlation function. At T = 1600 K, the LDA value afforded by the Stokes-Einstein relation is 2.78 mPa·s, close to the direct one, 2.71 mPa·s and the assessed value of 2.81 mPa·s. We obtain a similar agreement at higher temperature, i.e., T = 1800 K, while at lower temperature, namely T = 1400 K, we observe a small discrepancy, the value afforded by the Stokes-Einstein relation being 10% smaller than the direct one. GGA calculations follow the same trend. For instance at T = 1600 K, the value obtained from the SE relation is 3.45 mPa·s which can be compared to that obtained by the direct method, namely 3.55 mPa·s. Note that the present study does not support the conclusions of the EAM based study [18] claiming that the SE relation is not obeyed in liquid Cu.

Discussion
In order to understand the significant discrepancies between GGA and LDA calculations of the selfdiffusion coefficient and the viscosity, we inspect structural properties of liquid copper using both functionals.
In figure 3, we report the pair-correlation function g (r ) for all the three temperatures in the liquid range using GGA and LDA. As already obtained for liquid aluminum [23], we do not see any appreciable 43603-4 The curves for 1600 K, and 1398 K are shifted upwards by an amount of 1 and 2, respectively. difference between the two sets of calculations. Let us mention that our g (r ) at T = 1600 K compares well with g (r ) interpolated from X-ray diffraction experiments [36], while our g (r ) at T = 1398 K compare also well with the results from neutron diffraction experiments at T = 1393 K [37]. The calculated g (r ) from previous ab initio [21,22] or EAM model [13][14][15][16][17][18] are also consistent with our curves.
Coordination numbers are obtained from the integration of the computed pair-correlation functions g (r ) up to its first minimum. We find an increase of coordination number from N c = 11.8 ± 0.1 at T = 1800 K to N c = 12.4 ± 0.1 at T = 1398 K while GGA calculations give N c = 12.0 ± 0.1 at T = 1800 K and N c = 12.5 ± 0.1 at 1398 K. We obtain that LDA calculations give coordination numbers smaller than those obtained using GGA ones, but the differences remain small and cannot explain the differences obtained in the dynamic properties.
A deeper insight into the presence and the nature of the local atomic environments can be obtained from the common-neighbor analysis introduced in reference [38] and described in some detail in reference [39]. In figure 4, we display the most abundant pairs, averaged over the ten inherent configurations

43603-5
for each temperature. For both LDA and GGA calculations, 142x (sum of 1422 and 1421), 1431 and 15xx (sum of 1551 and 1541) pairs are found to be the main ones. We can also put forward two important remarks. 15xx pairs are the most abundant ones and they are the only pairs that display a linear increase as the temperature decreases. Let us mention that our results are in agreement with the previous ab initio based analysis [40,41]. The 15xx bonded pairs characterize the icosahedral order [42][43][44] while the 142x bonded pairs are related to close packed structures like the face-centered cubic and hexagonal close-packed structures [39]. The 1431 pairs can be considered as distorted icosahedra or distorted close-packed structures [45]. However, they cannot be responsible for the differences between dynamic properties using either GGA or LDA since they are similar in both approximations. The same comment holds for the 142x pairs. On the contrary, a significant increase of the number of 15xx pairs using the GGA approximation is a strong indication that the icosahedral short range order is more pronounced when using the GGA approximation.
To study the influence of ISRO on the dynamic properties, we inspect the velocity auto-correlation function from which the self-diffusion coefficient is obtained. In figure 5 we report this function computed at T = 1398 K using the LDA and GGA approximations. The function decays quickly first and gets to a minimum value after nearly 0.1 ps, which is known as the cage effect due to the temporary trapping of atoms by their neighbours. We can note that this first minimum becomes deeper using GGA, which reveals an increase of the cage effect. This minimum is also less pronounced when the temperature increases, as shown in the inset for LDA (a similar behaviour is seen for the GGA). This behaviour can be related to the increase of the icosahedral motifs using GGA as compared to LDA for the same temperature or when the temperature decreases as seen in figure 5. Indeed, as discussed in reference [41], the backscattering regime, that is predominant for liquid Cu, is more pronounced in the phase which presents the highest ISRO and consequently the diffusivity is smaller. We now discuss the universal scaling relationship proposed by Dzugutov [11], that relates a dimensionless form of the diffusion coefficient, D * to the excess entropy of the liquid phase, S E . The relationship can be written as: where D * is reduced with respect to the diffusion coefficient D using uncorrelated binary collisions described by the Enskog theory as:  the number density ρ within the framework of the hard-sphere fluid as: where m and σ are the atomic mass and the hard-sphere diameter. Parameter g (σ) is the value of the HS pair-correlation function g (r ) at the contact distance. In a first step, we study the scaling law assuming Dzugutov's idea that the excess entropy can be represented by the two-body correlation approximation denoted by: Note that σ and g (σ) are approximated by the position and height of the first peak in the ab initio computed pair-correlation function. Results for liquid Cu and previous ones for Al [23] are displayed in figure 6 (a) and show that the S 2 approximation is also not sufficiently accurate when utilizing either LDA or GGA calculations. It is worth mentioning that the truncation errors in equation (3.4), due to the limited range of g (r ) calculated by AIMD, do not exceed 2%, and, therefore, do not affect the results for Dzugutov's law.
To go beyond the two-body correlation approximation, we use a more reliable determination of S E but still keeping the framework of the HS reference fluid. It is based on the Carnahan-Starling equation which is known to give a quasi-exact equation of states for HS fluids in a range of values of packing fraction characteristic of melts [46]. In this case, the excess entropy can be written as:

43603-7
where the packing density is given by ξ = πρσ 3 /6. All entropy calculations can be done using the packing fraction ξ and its evolution with temperature. To determine ξ, we use the experimental densities [4] and σ values obtained previously. Note that the temperature dependence of ξ is mainly related to that of the density, as σ does not vary in the investigated temperature range. From figure 6 (a), one can see that the agreement with the original fit of Dzugutov is better since the scatter in data becomes smaller.
To explain the remaining discrepancy, we turn back to the determination of σ and g (σ). It is well known from perturbation theories [4] that the effective HS diameter adjusted to represent the structure and thermodynamics of a real system is different from the simple choice for σ and g (σ) made above. In particular, this method leads to a lack of thermodynamic consistency between the HS reference fluid and the real one. To evidence this failure, we compare in table 1 the isothermal compressibility calculated by the Carnahan Starling expression, i.e., ρk 4 ] (values in parenthesis) with that obtained from the extrapolation of S(q) at q = 0, S(q) being the structure factor obtained from ab initio calculations. Table 1 indicates for both Al and Cu that the thermodynamic consistency, i.e.
ρk B T χ T = S(0), via calculations of the isothermal compressibility is not respected. We present only LDA calculations for comparison as we obtain the same trend for GGA based results.  . To obtain S HS 2 , we compute the pair-correlation function of the hard sphere model, g HS (r ), using accurate integral equation method [47]. This method is based on the exact Ornstein-Zernike convolution equation h(r ) = c(r ) + ρ h(r )c(|r − r |)dr = c(r ) + γ(r ), (3.6) where h(r ), c(r ) and γ(r ) are respectively the total, direct and indirect correlation functions between two atoms, separated by a distance r in a liquid composed of N atoms with number density ρ [4]. In order to determine the pair-correlation function g (r ) = h(r ) + 1, equation (3.6) should be solved together with a closure relation whose formal expression is as follows: In equation (3.7), B (r ) is the so-called bridge function composed of an infinite series of elementary diagrams [4], and u HS (r ) is the hard-sphere potential. For B (r ) we use the efficient approximate formulation for the hard-sphere model proposed by Rodgers and Young [48]: which ensures an accurate description of thermodynamic properties of the HS fluid through the optimization of parameter A [49]. The set of equations (3.6), (3.7), and (3.8) are solved numerically using the

43603-8
algorithm of Labik et al. [50], which combines the Newton-Raphson method and the traditional iterative technique. Following the authors of reference [51,52], we use the tangent linear differentiation technique that is an essential ingredient for improving the accuracy of the integral equation method and for computing thermodynamic functions of the HS fluid that involve the derivative of g (r ).
From  4 ], or the extrapolation of S(q) at q = 0, i.e. ρk B T χ T = S(0), S(q) being obtained from LDA. This close correspondence is a strong indication of the thermodynamic consistency of our approach to determine the excess entropy.
Taking advantage of this simple analytical expression of the excess entropy, we plot in figure 6 (c) the relationship between the excess entropy and the dimensionless form of the diffusion coefficient, D * for both Cu and Al data. It is found that the scaling law proposed by Dzugutov is legitimate with the excess entropy derived from the Carnahan-Starling approach. Let us mention that this approach does not require any adjustable parameter.

Conclusion
In summary, we have computed the temperature dependence of the self-diffusion coefficient as well the viscosity of liquid copper by ab initio molecular dynamics using two different exchange and correlation potentials, LDA and GGA. The comparison of the calculated self-diffusion coefficients with the most recent QNS experimental data and the comparison of the calculated viscosities with the assessed values favor the LDA approximation to compute the dynamic properties of liquid copper. We show that the origin of the discrepancy using GGA is due to an enhancement of the icosahedral short range order (ISRO) that directly impact the dynamic properties of liquid copper. Our results allow us to discuss the applicability of the SE relation to the description of the relation between diffusivity and viscosity and to the firm establishment of its validity for liquid copper over the temperature range investigated. Finally, we show that the proposed scaling law of Dzugutov which relates the scaled diffusivity to the excess entropy can be successfully applied if a self-consistent method for determining the pair-correlation function of the HS fluid is used as well as the Carnahan Straling approach for expressing the excess energy is introduced. Note that the excess entropy based on the Carnahan Starling equation is given by an analytical expression and does not contain any adjustable parameter. How this formulation can be applied to liquids such as silicon, which presents loosely packed structures, is an important open question to be addressed in a future work.