Static structure, collective dynamics and transport coefficients in the liquid Li-Pb alloy. An ab initio molecular dynamics study

Static structure, collective dynamics and transport coefficients in the liquid Li-Pb alloy. An ab initio molecular dynamics study M.M.G. Alemany a , Jaime Souto-Casares a , Luis E. González b, ⇑ , David J. González b a Departamento de Física de Partículas, Área de Física de la Materia Condensada, Facultad de Física, Universidad de Santiago de Compostela, E-15706 Santiago de Compostela, Spain b Departamento de Física Teórica, Átomica y Óptica, Facultad de Ciencias, Universidad de Valladolid, E-47011 Valladolid, Spain article info Article history: Received 16 August 2021 Revised 1 October 2021 Accepted 3 October 2021 Available online 9 October 2021 Keywords: Liquid Li-Pb alloy Structure and dynamics Ab initio simulations abstract Several static and dynamic properties of the liquid Li-Pb alloy at diverse compositions, have been calculated by means of ab initio molecular dynamics simulation study. This alloy has attracted much attention because of the finding of fast sound at the Li 0:80 Pb 0:20 composition and also the technological interest of the the eutectic composition, Li 0:17 Pb 0:83 , as a component of the blanket in fusion reactors. Results are reported for total static structure factors, which are compared with the available experimental data. An additional analysis of the structure allows the quantification of the heterocoordinating tendencies in this alloy, which at the Li 0:80 Pb 0:20 composition are largest and lead to a Pb-centered polyhedral structure, where, however, Li 4 Pb units are not present. Regarding the collective dynamics, the calculated partial dynamic structure factors exhibit side peaks indicative of propagating density fluctuations, including density fluctuation modes with phase velocity greater than the hydrodynamic sound velocity. Also, the longitudinal and transverse dispersion relations have been calculated and its different branches analysed. We find all the high frequency branches to behave as optic-like modes, contrary to other interpretations in terms of an acoustic-like fast sound mode. Some transport coefficients such as self- and inter-diffusion coefficients, shear viscosities and adiabatic sound velocities, have also been calculated. Finally, the obtained results for the electronic density of states clearly indicate the metallic character of the liquid Li x Pb 1x alloy. Ó2021 The Authors. Published by Elsevier B.V. Thisisan openaccessarticleundertheCC BYlicense(http:// creativecommons.org/licenses/by/4.0/). 1. Introduction In the last four decades, liquid binary mixtures have been extensively investigated, both experimentally and by computer simulations, with the main aim focused towards understanding the microscopic mechanisms behind the collective excitations and the interdifussion processes. In this context, special attention has been devoted to those binary mixtures composed of particles with disparate masses because of the appearance of a high frequency mode with a phase velocity much greater than the value predicted by an extension of the hydrodynamic sound to large wavevectors. This has been the case in studies performed for Li 0:80 Pb 0:20 ,Na 0:50 Cs 0:50 and Li 0:30 Bi 0:70 using Molecular Dynamics (MD) simulations [1–4] and inelastic neutron scattering (INS) experiments [5–8]. In fact, kinetic theory calculations for binary mixtures with large atomic mass difference [2], have found that this high frequency mode was supported by the light atoms only and if the associated high frequency were interpreted in terms of a phase velocity of a propagating acoustic mode (also called ‘‘fast sound”) it would lead to phase velocities characteristic of the pure (light) component. Moreover, calculations using the revised Enskog theory (RET) [9] for binary hard sphere mixtures predicted, besides the hydrodynamic sound, two propagating collective modes. One was identified with the fast sound whereas the other, with an associated phase velocity smaller than the hydrodynamic one, was named ‘‘slow sound”. However, the very nature of this new high frequency excitation still remains a controversial point with differing opinions concerning it low qbehaviour. Thus, it has been suggested that as q!0 this new mode would merge into the hydrodynamic sound of the mixture whereas other authors have identified it with an opticlike mode that would take a non-zero frequency value as qgoes into the long-wavelength region. Molten Li 0:80 Pb 0:20 has been the system where this high frequency excitation has been more closely scrutinized. A first study was performed by Jacucci et al. [1] who carried out a classical https://doi.org/10.1016/j.molliq.2021.117775 0167-7322/Ó2021 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). ⇑ Corresponding author. E-mail address: [email protected] (L.E. González). Journal of Molecular Liquids 344 (2021) 117775 Contents lists available at ScienceDirect Journal of Molecular Liquids journal homepage: www.elsevier.com/locate/molliq molecular dynamics (CMD) simulation for the liquid Li 0:80 Pb 0:20 alloy at T = 1085 K. The liquid alloy was modelled as a mixture of positive Li and negative Pb ions and the interactions among the ions were represented by interatomic pair potentials consisting of a hard core repulsion plus a screened Coulomb interaction. They found that the Li-Li partial dynamic structure factor, S LiLi ðq; x Þ, displayed a high-frequency side peak, which dispersed linearly with a velocity (7500 m/s) which is substantially greater than the hydrodynamic sound velocity of the mixture (c s 2000 m/s). This feature was explained as a collective excitation supported by the light atoms and it was named the ‘‘fast sound” mode. Shortly after, additional INS experiments [6,7] on the liquid Li 0:80 Pb 0:20 alloy were performed in search of excitations with frequencies well above that of hydrodynamic sound. The obtained data confirmed the existence of a high-frequency mode supported by the light component (Li atoms) only, although its velocity ( 4500 m/s) was clearly smaller than the ‘‘fast sound” suggested by the previous CMD simulations [1]; in fact, this new estimate for the velocity was comparable to that of pure Li at similar thermodynamic conditions of number density and temperature. Moreover, the INS data of Alvarez et al. [7] suggested that this high frequency mode had features which clearly deviated from those expected for the propagation of a sound wave and pointed to the presence of rather localized, out-of-phase atomic motions, resembling those exhibited by Coulomb systems [10]. Actually, some theoretical models (viscoelastic model [4], generalized collective modes (GCM) theory [11]) envisage that a specific fast sound ‘‘mode” is not strictly necessary to account for the presence of a side peak in the partial dynamic structure factor of the light species. It may also appear as a consequence of a weighted sum of an ‘‘extended sound mode”, which goes into the hydrodynamic sound in the long-wavelength limit, and a kinetic (non-hydrodynamic) mode of much higher frequency which goes to a non-zero frequency for q!0 and whose weight in the sum decreases towards zero for decreasing wavevectors. Fernandez-Perea et al. [3] have also performed CMD simulations for the liquid Li 0:80 Pb 0:20 alloy by using the same interatomic pair potentials as Jacucci et al. [1]. However, their simulation box contained 6720 particles which allowed reaching q-values as small as 0.019 Å 1 and therefore to gather information well into the hydrodynamic region. Oddly enough, in constrast to the conclusions from the INS data of Alvarez et al. [7], their results showed the existence of two branches of collective excitations which merged, for q60:11 Å 1 into the usual hydrodynamic linear dispersion. Bryk and Mryglod [12] have also studied the static and dynamic properties of the liquid Li 0:80 Pb 0:20 alloy by combining CMD simulations (with the same pair potentials as in Ref.[1,3]) with the generalized collective modes (GCM) method. Basically, the time correlation functions yielded by the CMD simulations are analyzed within the framework of the GCM method which provides some additional insight into the hydrodynamic and non-hydrodynamic collective processes existing in the liquid with different spatial and time scales. Their results showed the existence of two branches of collective excitations which in the hydrodynamic limit were associated to hydrodynamic sound and optic-like excitations respectively. Interestingly, previous to the above mentioned developments associated with the liquid Li 0:80 Pb 0:20 alloy, the liquid Li-Pb alloy had already attracted much experimental and theoretical work. On the experimental side, we notice that most measured properties show a non-ideal dependence with composition and the maximum deviation from ideal behaviour is usually found around the composition, Li 0:80 Pb 0:20 . For instance, at this composition the stability function has a marked peak (suggesting strongly reduced concentration fluctuations) [13], the entropy of mixing shows a dip [14], the excess volume per atom has a maximum relative deviation of about 18% [15], the interdiffusion coefficient shows a maximum [16] and the electrical resistivity attains a maximum value of 500 lX cm, although the system still remains metallic [17]. It is believed that this behaviour is related to the existence of some electronic charge transfer from the Li atoms to the Pb atoms which is driven by the significant electronegativity difference between both type of atoms with a net effect of rendering the bond partially ionic. In fact, several studies showed that some of the thermodynamic [18,19], electronic [20], structural [21], and dynamic [16,22] properties of Li 0:80 Pb 0:20 could be at least qualitatively explained by assuming the existence, even if transient, of Li 4 Pb units in the melt. On the other hand, the static structure of the liquid Li-Pb alloy has been determined, for several temperatures and concentrations, by means of neutron scattering (NS) experiments [23,24]; moreover, the subsequent analysis pointed out to the existence of heterocoordinating tendencies which become more marked near the Li 0:80 Pb 0:20 composition. As for the dynamical properties, we mention that the first INS experiment was performed by Soltwisch et al. [25,22] for the liquid Li 0:80 Pb 0:20 alloy at T = 1023, 1098 and 1173 K. These experiments yielded the quasielastic part of the concentration-concentration Bathia-Thornton dynamic structure factor, S cc ðq; x Þ. This is so because if the 7 Li isotope is used in the alloy, then the coherent part of the scattered intensity (which depends on the momentum transfer, q, and energy transfer,  h x ) is proportional to S cc ðq; x Þ. Moreover, as the incoherent contribution is dominated by the 7 Li atoms, the analysis of the data was performed by fitting the measured total scattered intensity to a sum of a term proportional to the S cc ðq; x Þplus another one proportional to the incoherent contribution of the 7 Li atoms. By restricting the analysis to the quasi-elastic region and performing a fitting to the analytical hydrodynamic limits of both contributions, they obtained estimates of the interdiffusion and self-diffusion coefficients in this alloy. Although the liquid Li-Pb alloy has been intensively studied by means of MD simulations, specially after the detection of the high frequency mode, however it is startling that ab initio MD (AIMD) simulation studies of its dynamical properties are yet to be performed. Moreover, all previous CMD simulation studies have been carried out using exactly the same interatomic pair potentials, i.e. hard core plus screened Coulomb interaction. There have been, in fact, a few studies that used effective manybody embedded atom model potentials, that either did not [26], or did [27], include also screened Coulomb potentials between the components, but they did not delve into the analysis of dynamic properties. Therefore, it seems important to perform a study where the interactions among the atoms/ions are evaluated at a more fundamental level. To our knowledge, the only AIMD study of this alloy [28] was focused on two concentrations, namely Li 0:50 Pb 0:50 and Li 0:80 Pb 0:20 , and they reported results for static and electronic properties only. This paper describes an AIMD simulation study of the structural, dynamical and electronic properties of the liquid Li x Pb 1x alloy at several concentrations. Obviously, we have considered the Li 0:80 Pb 0:20 composition as we are interested in attaining a more detailled picture of its collective excitations by analysing magnitudes which are not yielded by experiment. Another interesting alloy is the eutectic composition, Li 0:17 Pb 0:83 , which due to its low melting point, high boiling point and appropriate absorption and activation cross section for neutrons has been proposed [29] as a blanket in nuclear fusion reactors. Moroever, it has the capability of acting as a tritium breeder and neutron multiplier and it could also be used as a coolant. Therefore, it is important to gather as much information as posible concerning the properties of the M.M.G. Alemany, J. Souto-Casares, L.E. González et al. Journal of Molecular Liquids 344 (2021) 117775 2 eutectic alloy. Finally, we have also studied two other concentrations for which some experimental data are available. In this way, it will be possible to analyze the way in which the various properties are influenced by the composition. The present AIMD simulations have been performed by using the PARSEC code [30] which has already been used to study some pure liquid metals [31–33] (including liquid Pb near its triple point) and alloys [34], as well as other types of liquids and glasses [35]. Nevertheless, the complex nature of the liquid Li x Pb 1x alloy poses a stringent test on the capability of this method to describe this type of systems. Computational details are briefly described in Section 2, and in Section 3we present and discuss our results, comparing them with available experimental data. Finally, in Section 4, we summarize our main conclusions. 2. Technical details Within the PARSEC method [30], the Density Functional Theorybased Kohn–Sham equations [36,37] were solved self-consistently on a rectangular three-dimensional real-space grid within a supercell geometry [38]. The core electrons of the Li and Pb atoms were represented by Troullier-Martins norm-conserving pseudopotentials [39], generated for the reference configurations [He]2s 1 2p 0 and [Xe](4f 14 5d 10 )6s 2 6p 2 6d 0 5f 0 , respectively, and with associated radial cutoffs of 2.8 a.u. and 3.2 a.u. The potentials were made separable by the procedure of Kleinman and Bylander [40], applied in real space, with the ppotential chosen to be the local component. A partial-core correction for nonlinear exchange correlation was included for Pb. The local density functional of Ceperley and Alder [41], as parametrized by Perdew and Zunger [42], was used and the single C-point was employed in sampling the Brillouin zone. A spacing of 0.64 a.u. was used for constructing the real-space grid. Calculations were performed for the thermodynamic states shown in Table 1. The total number of atoms considered in the simulation varied from 240 (for x Li = 0.50) to 330 atoms (for x Li = 0.80) and the simulation was started by placing the atoms, at random, in a cubic supercell with dimensions chosen so as to obtain the respective experimental number density [15]. The temperatures of the simulations were chosen as those of (or close to) the corresponding NS experiments. The AIMD simulation runs were performed with the ionic dynamics being generated using the Beeman algorithm [43] with Hellmann–Feynman forces [44]. After thermalization, the number of time steps that were used for calculating the static, dynamic and electronic properties reported below varied from 20000 (80 ps of simulation time for Li 0:80 Pb 0:20 )to 28000 (112 ps of simulation time for Li 0:17 Pb 0:83 ). 3. Results and discussion 3.1. Structural properties Fig. 1 shows the calculated partial pair distribution functions g ij ðrÞof the liquid Li x Pb 1x alloy at four concentrations. We observe that the main peak’s height of the g LiPb ðrÞis always similar (x Li = 0.17) or higher than those of the g LiLi ðrÞand g PbPb ðrÞ. This feature is indicative of heterocoordinating tendencies in the alloy and implies that each atom tends to be surrounded, in its first coordination shell, by atoms of the other species. Notice that when x Li = 0.80, the main peak in g PbPb ðrÞhas practically vanished and the second peak is unusually high. This confirms the results of the pioneering AIMD study of Senda et al. [28], although with increased statistical accuracy. In that study some similarity between the structure of the liquid Li 0:80 Pb 0:20 alloy and that of the crystalline compound Li 7 Pb 2 (with similar composition, x Li 0:78) was highlighted, including the lack of Pb atoms as near neighbors of a Pb atom, and the long Pb-Pb distance, which is somewhat smaller in the solid (4:75 Å 1 ) than in the liquid, where the maximun of g PbPb ðrÞis located around 5:2Å 1 . There is another solid compound in the phase diagram of the alloy, which is the Lirichest one, that after several revisions has been assigned as Li 17 Pb 4 [45] with a composition (x Li 0:81) even closer to the liquid one studied. Its atomic structure is rather complex, with 105 atoms in the primitive unit cell, and can be studied as composed of two Li 20 Pb 6 clusters, one Li 22 Pb 4 cluster and one Li 23 Pb 4 cluster. If full charge transfer is assumed from Li (becoming Li + ) to Pb (becoming Pb 4 ) the first two clusters would have negative charge (4) while the latter two would instead be positively charged (+6 and +7, respectively). Consequently some more electron charge must be localized outside the clusters (nominally 5 electrons) and in fact this is what has been found in a recent ab initio study of this solid compound [46], which therefore qualifies as an electride. In fact the charge transfer was found not to be complete but partial, so that the clusters retained the sign but not the magnitude of the charges, and the total number of electrons outside the clusters was indeed smaller and localized in 4 particular regions of the structure (see [46] for details). Alternatively, the atomic structure can be described as built from different types of Pb centered Li polyhedra [46]. Three of these building blocks contain 13 Li atoms while the fourth has 14 of them. Therefore, the average number of Li atoms around a Pb in the solid is 13:25, there are no Pb near neighbors of a Pb atom and the Pb-Pb distances are again quite long, in fact distributed within a range from approximately 4:75 to 5:50 Å [46]. Consequently, the liquid structure shares even more similarities with that of the Li 17 Pb 4 compound. Such polyhedral arrangements are common in other types of solid compounds, like oxides, and there is some consensus that the molten structures can also be described in terms of polyhedra, although distorted or even with changes in the number of faces of the polyhedra [47]. We will later analyze this possibility for the particular case of liquid Li 0:80 Pb 0:20 below. From the g ij ðrÞ, we can obtain the coordination numbers n ij describing the number of j-type particles around an i-type particle within a sphere of radius R ij , i.e. n ij ¼4 pq x j Z R ij 0 r 2 g ij ðrÞdr;ð1Þ where R ij is usually identified [48] with the position of the first minimum of the related partial radial distribution function, 4 p r 2 g ij ðrÞ.A quantitative estimate of the short range order (SRO) in the alloy and Table 1 Thermodynamic input data of the liquid LixPb1xalloy used in the present AIMD simulation study. The total ionic number density, q, was taken from [15]. x Li q(Å 3 )T(K) 0.17 0.0325 775 0.50 0.0390 775 0.62 0.0431 775 0.80 0.0436 1075 M.M.G. Alemany, J. Souto-Casares, L.E. González et al. Journal of Molecular Liquids 344 (2021) 117775 3 its dependence with concentration is provided by the Warren [49] SRO parameter for the first neighbour shell, a 1 . It is defined as a 1 ¼1n ij x j ðx i n j þx j n i Þðj–i¼1;2Þð2Þ where x j is the concentration of the j-type particles and n i ¼n ii þn ij (i;j= 1, 2). The computed values of n ij and a 1 are given in Table 2 where its is observed that a 1 takes negative values for all concentrations; this feature points to heterocoordinating tendencies in the alloy, which for the concentrations studied in this work are most marked in the case of x Li ¼0:80. Note also that the average number of Li atoms around a Pb atom in the alloy with x Li ¼0:80 is around 11, which is substantially smaller than what was found in solid Li 17 Pb 4 , pointing to an important change in the characteristics of the Li polyhedra around a Pb atom, if this picture is at all correct. In order to check if the polyhedral description of the liquid structure in Li 0:80 Pb 0:20 is a valid one (results are shown in figure 2) we have analyzed all of the configurations generated in the following way. For each Pb atom we have identified those Li atoms that are its near neighbors. Then we have checked if any of the other Pb atoms lies inside the convex hull of the set of these Li atoms, i.e. inside the polyhedron whose vertices are the Li atoms. In order to perform this check we have used a simplified version of the GJK (Gilbert-Johnson-Keerthi) algorithm, which is widely used to detect collisions of 2- or 3-dimensional objects in visualization programs. If the polyhedral description of the liquid were correct then no such cases should occur, as each Pb atom should be inside its own polyhedron. Note that this question is not just a simple comparison about some Pb-Pb distances being shorter than Pb-Li distances, which certainly happens as evidenced by the small, but not null, first peak of g PbPb ðrÞthat lies within the first peak of g PbLi ðrÞ. The problem is strictly a 3-dimensional geometric one: for example, the distance between the centers of two face-sharing regular tetrahedra is smaller than the distance between the center and the vertex of the tetrahedra. As a result of the check we have found that out of the 20000 configurations generated only in two of them there was just a single Pb atom that was inside the polyhedron associated to a different Pb atom. Therefore we can confidently say that the polyhedral picture is not perfect, but almost so. The average number of Li atoms that are near neighbors of a Pb atom, as already indicated above, is around 11. However, the number of Li atoms in the polyhedra are of course distributed, and Fig. 2 shows this distribution. We see the appearance of PbLi 7 to PbLi 15 polyhedra, with PbLi 11 being the most abundant one. Since the average concentration in the alloy implies 4 Li atoms per Pb atom, obviously there must be many Li atoms shared by more than one Pb atoms, i.e. Li atoms that belong to the polyhedra of several Pb atoms. But what can be definitely discarded is the presence of Li 4 Pb units that have been suggested by other authors to exist in order to justify the properties of the alloy at this particular composition [16,18–22]. In order to compare our AIMD results with the NS data of Ruppersberg and Reiter [23,24], we have first calculated the Ashcroft- Langreth (AL) partial static structure factors S ij ðqÞwhich are depicted in Fig. 3. The S LiLi ðqÞshow a prepeak whose position moves from q1:5Å 1 (for x Li ¼0:17) to q1:65 Å 1 (for x Li ¼0:80); moreover it also shows another peak at q2:45 Å 1 (for x Li ¼0:17) whose amplitude grows with increasing Li concentration so that it becomes its main peak while its position slightly moves to q2:55 Å 1 (for x Li ¼0:80). Remarkably, the prepeak exhibited by S LiLi ðqÞat x Li ¼0:80 had not been predicted by previous CMD simulations [1,4,3]. The S PbPb ðqÞhas for x Li ¼0:17 a main maximum at q2:35 Å 1 , which for increasing Li concentration diminishes its amplitude while, at the same time, another peak develops at a smaller q(q1:6Å 1 ) so that for x Li ¼0:80 has become the main peak. Note also that the position of the first peak or prepeak in the S LiLi ðqÞand S PbPb ðqÞare almost coincident with the position of the first minimum of S LiPb ðqÞ, which is another indication of a certain amount of chemical ordering. This type of ordering is clearly revealed in the Bhatia-Thornton (BT) partial structure factors [50–53] that we have also evaluated and are shown in Fig. 4, namely the number-number, S NN ðqÞ, the concentration-concentration, S CC ðqÞ, and the numberconcentration, S NC ðqÞ, partial structure factors. Their longwavelength limits provide microscopic information on the ordering tendencies of the liquid alloy, in particular S CC ðq!0Þ, for which we have obtained the following results: S CC ðq!0Þ=ðx Li x Pb Þ Fig. 1. Partial pair distribution functions g ij ðrÞof the liquid Li x Pb 1x alloy at x= 0.17, 0.50, 0.62 and 0.80. Full blue, dashed red and dotted lines correspond to g LiLi ðrÞ;g PbPb ðrÞ and g LiPb ðrÞ, respectively. Table 2 Calculated coordination numbers nij and the Warren-Cowley short-range order parameters aðiÞ 1for the liquid LixPb1xalloy at the thermodynamic states given in Table 1. x Li n LiLi n LiPb n PbPb n PbLi a 1 0.17 0.8 9.2 9.2 1.9 0.09 0.50 4.6 6.6 5.0 6.6 0.16 0.62 6.5 5.3 3.2 8.7 0.18 0.80 8.0 2.8 0.3 10.9 0.26 M.M.G. Alemany, J. Souto-Casares, L.E. González et al. Journal of Molecular Liquids 344 (2021) 117775 4 0.40, 0.35, 0.32 and 0.10 (20%) at x Li = 0.17, 0.50, 0.62 and 0.80, respectively. These estimates are smaller than unity, which clearly suggests heterocoordination tendencies, in agreement with the information provided by the Warren-Cowley SRO parameter. For all concentrations, S CC ðqÞshows a distinct main peak at around q1.6 Å 1 , which coincides with the positions of the first peak/ prepeak in S LiLi ðqÞand S PbPb ðqÞand minimum in S LiPb ðqÞ, pointing thus to their physical origin as due to chemical order. The total neutron weighted structure factor S T ðqÞis readily evaluated either from the AL or the BT partial structure factor. The expression for the latter case is S T ðqÞ¼ ðx Li b Li þx Pb b Pb Þ 2 S NN ðqÞþðx Li b Li þx Pb b Pb Þðb Pb b Li ÞS NC ðqÞ h þðb Pb b Li Þ 2 S CC ðqÞix Li b 2 Li þx Pb b 2 Pb hi 1 ; ð3Þ where b Li and b Pb denote the neutron scattering lengths of Li and Pb, respectively [54]. We recall that the liquid alloy samples for the NS measurements were prepared [23,24] from the 7 Li isotope and natural Pb; therefore, our evaluation of the S T ðqÞhas been made using the associated values b Li =2.22 fm and b Pb = 9.40 fm. In the case of x Li ¼0:80, due to the particular values of the scattering lengths, we have that x Li b Li þx Pb b Pb 0, and then S CC ðqÞx Li x Pb S T ðqÞ, so that it can be probed by a single NS experiment. Fig. 4 shows the calculated S T ðqÞof the liquid Li x Pb 1x alloy along with the available NS data at x Li ¼0:17;0:50;0:62 and 0:80 [23,24]. Notice that the shape of the experimental S T ðqÞundergoes significant changes as the Li concentration is increased, with its main peak moving to smaller q-values and its position being practically determined by that of the S PbPb ðqÞ. Nevertheless, despite those changes undergone by the S T ðqÞ, we highlight the excellent agreement with experiment achieved by the present AIMD calculations. 3.2. Dynamic properties 3.2.1. Self- and inter-diffusion The long-time behaviour of the motion of one particle in a liquid system is characterized by the self-diffusion coefficient, and in a binary alloy two different self-diffusion coefficients exist, one for each component, D 1 and D 2 . Additionally, in the case of mixtures one can consider the motion of the center of mass of the particles of each component, and the long time behaviour of the relative motion of the center of mass of particles of type 2 with respect to that of particles of type 1 comes characterized by the relative difussion coefficient, D 12 . In an ideal mixture in which all components are identical, the relative diffusion coefficient can be expressed directly in terms of the self-diffusion coefficients of each component, through the so called ideal formula D 0 12 ¼x 2 D s 1 þx 1 D s 2 . Any deviation from this rule is usually reported in terms of the quantity c 12 , defined so that D 12 ¼D 0 12 ð1þ c 12 Þ¼D 0 12 þx 1 x 2 D d 12 ; where we have also defined the ‘‘distinct” part of the relative diffusion coefficient, D d 12 . As its name suggests, it is a meassure of the influence of distinct particles, irrespective of their being of the same type or of different type, on the relative diffusion coefficient. In computer simulations all these diffusion coefficients can be obtained as time integrals of the corresponding velocity autocorrelation functions, namely, D s i from Z s i ðtÞ, which is the usual velocity autocorrelation function of a tagged i-type particle of the fluid, D 0 12 from Z 0 12 ðtÞ, the ideal (also called self-) relative velocity autocorrelation function, given as Z 0 12 ðtÞ¼x 2 Z s 1 ðtÞþx 1 Z s 2 ðtÞ, and D 12 from Z 12 ðtÞ, the relative velocity autocorrelation function, defined as the time autocorrelation function of the relative velocity of the center of mass of particles of type 2 with respect to that of particles of type 1. Finally, D d 12 is obtained as the time integral of the distinct relative velocity autocorrelation function, Z d 12 ðtÞ, defined so that Z 12 ðtÞ¼Z 0 12 ðtÞþx 1 x 2 Z d 12 ðtÞ:ð5Þ The sign of c 12 hints towards a dynamic tendency to homocoordination or heterocordination, namely, if c 12 >0 then particles of the same species have a greater tendency to diffuse together than particles of different species, and conversely in the case of c 12 <0. More details can be found in [55]. Finally, the interdiffusion coeffi- Fig. 3. Ashcroft-Langreth partial static structure factors S ij ðqÞof the liquid Li x Pb 1x alloy at x Li = 0.17, 0.50, 0.62 and 0.80. Full blue, dashed red and dotted lines correspond to S LiLi ðqÞ;S PbPb ðqÞand S LiPb ðqÞ, respectively. Fig. 2. Distribution of PbLi n polyhedra present in the Li 0:80 Pb 0:20 alloy. M.M.G. Alemany, J. Souto-Casares, L.E. González et al. Journal of Molecular Liquids 344 (2021) 117775 5 cient is given as D int ¼hD 12 , where h¼x 1 x 2 =S CC ðq!0Þ. For a nearly ideal mixture, h1; c 12 0 and, therefore, D int D 0 12 , which is also called Darken’s approximation. Fig. 5 shows the results obtained for the normalized self and relative VCFs. The self VCF of the heavier particles, Z s Pb ðtÞ, has the slower decay and smaller backscattering because of the velocity persistence of the heavy atoms when colliding with the lighter ones. However, when the amount of Li decreases (greater proportion of heavy atoms), Z s Pb ðtÞbecomes narrower and shallower. On the other hand, Z s Li ðtÞhas a steep decay and a marked backscattering which is more pronounced with increasing concentration of Pb ions. The calculated diffusion coefficients are given in Table 3. Notice that the D s Li is always greater than D s Pb because of the smaller mass of the Li atoms. Experimental estimates for the self-diffusion coefficients are only available for the liquid Li 0:80 Pb 0:20 alloy. Soltwisch et al. [25] fitted their measured INS data to a model that allowed to discriminate the incoherent quasielastic scattering contribution (which is dominated by 7 Li) and, by assuming an hydrodynamic aproximation, they obtained (for T = 1098 K) the following ‘‘experimental” estimates: D s Li = 2.13 0.18 and D s Pb = 0.33 0.03 (in 10 4 cm 2 /s units), to be compared with the present AIMD values of 1:64 and 0:45 respectively (see Table 3). As for the other concentrations, we are not aware of any other experimental data; nevertheless we stress that the application of this same AIMD method to other liquid metals, including liquid Pb near melting [31,32], has yielded estimates for the self-diffusion coefficients in very good agreement with the experimental data. As for the c LiPb , we obtain negative values (suggesting mild heterocoordinating tendencies) except for the x Li ¼0:80 concentration when a positive value is obtained. Such behavior at this concentration hints towards a transient persistance of the aforementioned Li polyhedra that would surround a Pb atom, since this would imply a collective displacement of at least several of the Li atoms of the polyhedron so that like particles tend to move together. The values of c LiPb , combined with those of h, which range from 2:5 to 10 for the concentrations considered, lead to interdiffusion coefficients that are also shown in Table 3. Here we observe a distinct deviation from the ideal mixture values, D 0 LiPb . The ratio D int =D 0 LiPb takes a value around two for the x Li ¼0:17;0:50 and 0:62 concentrations, where the relatively high values of hovercome the negative c LiPb , but it grows to a factor of fifteen when x Li ¼0:80, where both c LiPb and hpush up the value of D int . The interdiffusion coefficient in the liquid Li-Pb alloys has been measured for a range of concentrations and temperatures [16] and those values have also been included in Table 3. It is observed that the present AIMD calculations show a good agreement with experiment, including the substantial increase for x Li ¼0:80. 3.2.2. Collective dynamics The AL and the BT partial intermediate scattering functions (ISF), F ij ðq;tÞ, and their Fourier transforms, the AL and BT partial dynamic structure factors, S ij ðq; x Þ, describe the collective dynamics of the fluctuations in the partial densities, component-wise in the case of the AL functions, topological and chemical in the case of the BT ones (see Ref. [34] for more details). As a consequence of the continuity equation, that relates time derivatives of the densities with longitudinal currents, the dynamic structure factors are directly related to the Fourier transforms of Fig. 4. Bhatia-Thornton partial static structure factors and total structure factor of the liquid Li x Pb 1x alloy at x Li = 0.17, 0.50, 0.62 and 0.80. Continuous, dotted and green dashed lines correspond to S NN ðqÞ;S NC ðqÞand S CC ðqÞ=ðx Li x Pb Þrespectively. The red continuous lines stand for the calculated S T ðqÞ, whereas the open circles are the corresponding NS data of Ruppersberg et al. [23,24]. Note that for x Li ¼0:80 S CC ðqÞ=ðx Li x Pb Þpractically coincides with S T ðqÞ. Fig. 5. Normalized self, relative and ideal VACFs for the liquid Li x Pb 1x alloy at several x Li values. Full, dashed, thin (blue) and dotted (red) lines correspond to Z s Li ðtÞ;Z s Pb ðtÞ, Z LiPb ðtÞ, and Z 0 LiPb ðtÞ, respectively. M.M.G. Alemany, J. Souto-Casares, L.E. González et al. Journal of Molecular Liquids 344 (2021) 117775 6 the time correlation functions of the longitudinal components of the partial currents (AL and BT), C L ij ðq; x Þ, so that C L ij ðq; x Þ¼S ij ðq; x Þ x 2 =q 2 . We will consider below these longitudinal current correlation functions, which can be helpful in order to discern longitudinal modes not directly visible in S ij ðq; x Þ, as well as the transverse ones, C T ij ðq; x Þ, which provide additional information about collective dynamics in the system associated to the possible propagation of shear waves and enable the calculation of the shear viscosity of the liquid. Starting from the atomic positions and velocities obtained from the simulations we have evaluated the partial densities and partial currents (longitudinal and transverse components) according to their microscopic definition [34], and from them we have calculated the corresponding time correlation functions and finally their Fourier transforms numerically after applying a window function to alleviate the numerical inaccuracies that typically occur at long times (see appendix B of [56]). Fig. 6 shows, for some q-values, the calculated partial intermediate scattering functions, F ij ðq;tÞ, obtained for x Li ¼0:17;0:50 and 0:80, at similar q-values. At small q’s the partials F LiLi ðq;tÞ;F PbPb ðq;tÞ and F LiPb ðq;tÞare dominated by diffusive contributions which impose a slow decay and conceal the oscillations associated with the propagating density fluctuations. Nevertheless, for x Li ¼0:80 the decay is distinctly faster and the oscillations are more visible. The BT partial F NN ðq;tÞshows a behavior qualitatively similar to that of the AL partials. The existence of propagating density fluctuations can also be revealed by the presence, within some q-range, of side peaks/ shoulders in the partial dynamic structure factors S ij ðq; x Þ.Fig. 7 depicts our calculated S LiLi ðq; x Þ;S PbPb ðq; x Þand S NN ðq; x Þat different concentrations and their respective smallest attainable q-value, namely q min . Notice that S NN ðq¼q min ; x Þexhibits a shape similar to the hydrodynamic Rayleigh-Brillouin triplet, which for a binary system includes sound propagation contributions in the side peaks and contributions from thermal diffusion and also from interdiffusion in the central line. Adiabatic sound velocities have been estimated from the position of the Brillouin peak in S NN ðq min ; x Þ, denoted as x B ðq min Þ,asc s  x B ðq min Þ/q min . The values we have obtained for x Li ¼0:17;0:50;0:62 and 0:80 are c s 1500, 1350, 1600 and 1700 150 m/s, respectively. Comparison with experiment can only be performed for x Li ¼0:17 and 0:80, where measurements of the hydrodynamic sound velocity yielded  1720 m/s [57] and 2000 m/s [15], respectively. Side peaks/shoulders in S NN ðq; x Þappear also at higher wavevectors, and the qrange of their appearance has a weak concentration dependence. So, whereas for x Li ¼0:17 the calculated S NN ðq; x Þdisplay side peaks/shoulders up to q1.1 Å 1 , this range slightly increases with increasing x Li so that when x Li ¼0:80 side peaks/shoulders are found up to q1:3Å 1 . The S PbPb ðq; x Þfollow a similar trend. As for the S LiLi ðq; x Þ,we observe the appearance of small amplitude, but clearly discernible, side peaks/shoulders at frequencies much higher than those shown in Fig. 7, for q-values up to q1:3Å 1 in the case of x Li ¼0:17, and Table 3 Diffusion coefficients (in 10 4 cm 2 /s) and parameter cLiPb of the liquid LixPb1xalloy at the thermodynamicstates given in Table 1. The numbers in parenthesis are the experimental values for Dint obtained by Khairulin et al. [16]. x Li 0.17 0.50 0.62 0.80 D s Li 0.52 0.53 0.49 1.64 D s Pb 0.33 0.33 0.26 0.45 D LiPb 0.44 0.37 0.25 0.97 D 0 LiPb 0.49 0.43 0.35 0.69 D d LiPb 0.35 0.21 0.43 1.77 c LiPb 0.10 0.12 0.29 0.41 S CC ð0ÞD int 0.06 0.09 0.06 0.16 D int 1.0 1.0 0.8 10.0 (0.85) (1.2) (9.0) S CC ð0ÞD Darken int 0.07 0.11 0.08 0.11 Fig. 6. Partial intermediate scattering functions, F ij ðq;tÞ, for the liquid Li x Pb 1x alloy at three concentrations. The solid blue line corresponds to F LiLi ðq;tÞ, the dashed red line to F PbPb ðq;tÞ, and the dotted line is F LiPb ðq;tÞ. Circles denote the BT F NN ðq;tÞ. Fig. 7. Partial dynamic structure factors S ij ðq;xÞof the liquid Li x Pb 1x alloy at x Li ¼0:17;0:50;0:62 and 0:80 for q¼q min ¼0:32;0:34;0:35 and 0:32 Å 1 respectively. Full (blue) line, dashed (red) line and circles correspond to S LiLi ðq;xÞ;S PbPb ðq;xÞand S NN ðq;xÞ, respectively. The insets show 10 2 S ij ðq;xÞ. M.M.G. Alemany, J. Souto-Casares, L.E. González et al. Journal of Molecular Liquids 344 (2021) 117775 7 up to q1:4Å 1 for x Li ¼0:80. These results obtained for x Li ¼0:80 are qualitatively similar to those provided by the CMD simulations of Bosse et al. [1,2] which yielded S LiLi ðq; x Þwith side peaks up to q1:2Å 1 , although their amplitudes were more marked than the present AIMD ones. From the positions of the side peaks/shoulders in the S LiLi ðq; x Þ and S PbPb ðq; x Þ, we have evaluated the associated dispersion curves, x LiLi ðqÞand x PbPb ðqÞwhich are depicted in Fig. 8. For all concentrations we obtain two dispersion branches which are indicative of two processes clearly separated in frequency; moreover, their separation increases with increasing Li concentration. The small-qbehaviour of the lower frequency branch, x PbPb ðqÞ,is qualitatively similar to what has been obtained by other authors through CMD simulations or derived from experimental data. Typically, when approaching the hydrodynamic regime, the associated phase velocity is smaller than the hydrodynamic value but as q!0 this branch smoothly merges into hydrodynamic sound. The high-frequency branch shows up only in the Li-Li partials. In the case of x Li ¼0:80, the calculated high frequency branch is compatible with a phase velocity of 4500 m/s, which is very similar to the results derived from the INS data of Alvarez et al. [7] whose associated phase velocity was also around 4500 m/s. Interestingly, both values are much smaller than that of 7500 m/s obtained by means of CMD simulations. [2,3] The AIMD phase velocity obtained is not far from the hydrodynamic sound velocity that would correspond to pure liquid Li at the same conditions of density and temperature as those of the alloy, as the RET suggests: we have performed an additional AIMD simulation for such a system and obtained a sound velocity of 4150 m/s. Notwithstanding the importance of the agreement between experiment and our simulations, it has been the physical nature of these high frequency excitations that has sparked strong controversies. It is not clear whether or not when qdiminishes towards the hydrodynamic limit, this high frequency branch undergoes a continuous transition into the hydrodynamic sound mode and merges with it, as the RET proposes, or otherwise, the high frequency mode is really a non-hydrodynamic optic-like mode whose weight on the Li-Li partials decreases with decreasing qand of course vanishes at q¼0. This debate has mainly revolved around the liquid Li 0:80 Pb 0:20 alloy as this is the composition that has attracted the most attention and where this high frequency branch was first detected. This latter interpretation has been suggested by the authors of the INS measurements [7] and also by several theoretical approaches applied to the data obtained from CMD simulations, such as the viscoelastic theory [4], or the GCM approach [12]. In the present AIMD simulations the smallest attainable q-value is 0:35 Å 1 , which is probably outside the hydrodynamic region, so we anticipate that a definite conclusion concerning the low-q behavior of the high frequency branch will not be possible. However, a viable interpretation can be made by observing the trends in Fig. 8 and the additional data provided by the longitudinal current correlation functions to be discussed below. In other AIMD simulation studies of Li-based alloys, namely Li- Ba and Li-Bi [34,58], the evaluation of the respective dispersion relations from the peaks/shoulders of the partial dynamic structure factors yielded also two dispersion branches: a high frequency one associated to Li-Li partials, and a low frequency one. In both alloys when approaching the hydrodynamic region, the corresponding high frequency branch underwent a continuous decrease and merged with the hydrodynamic sound mode from above, while the low frequency brach joined it from below. Nevertheless, the results were also compatible with a scenario where there is not a specific fast-sound mode, but an optic-like non-hydrodynamic mode of higher frequency whose weight on the Li-Li partials decreases as qgoes into the hydrodynamic region. The trend obtained in the present study of the Li-Pb alloys, as observed in Fig. 8, points towards a clearer picture than the other Li alloys mentioned before. It does not suggest that the high frequency branch were to (continuously) change its slope and merge into hydrodynamic sound, but rather to tend to a finite non-zero value as qgoes to zero, which is the behavior exhibited by optic-like excitations in binary liquids. The case of x Li ¼0:80 is certainly the most suggestive of this interpretation, but for the other concentrations the picture also looks tenable. In fact, looking closely at the S LiLi ðq; x Þfor small q, and especially in the Li-rich alloys, a shoulder can be observed at a frequency similar to, or slightly larger than, that of the side peak/ shoulder corresponding to the Pb-Pb partials (see Fig. 7). Therefore, at least in this region of the phase diagram, it is possible that the true picture is in fact a combination of both theoretical scenarios considered, namely, the Li-Li partials show two modes, the opticlike non hydrodynamic high frequency one, and the extended sound mode with some positive dispersion, whereas the Pb-Pb partials basically exhibit only the extended sound mode, now with a negative dispersion. As indicated above, the spectra of the current correlation functions, C L ij ðq; x Þ, due to the apperance of the factor x 2 , show a depleted contribution coming from the smallx modes (in particular those of diffusive type), and therefore help to uncover longitudinal modes that are not visible in the partial dynamic structure factors either because of its small weight and/or because they are shielded by the diffusive modes. Fig. 9 shows the AIMD calculated Fig. 8. Dispersion relations, x LiLi ðqÞ(full circles) and x PbPb ðqÞ(full squares) of the maxima in the partials, S LiLi ðq;xÞand S PbPb ðq;xÞfor the liquid Li-Pb alloy at x Li ¼0:17;0:50 and 0:80. The slope of the dashed line is the calculated adiabatic sound velocity in the liquid alloy. For x Li ¼0:80 the slope of full line amounts to a sound velocity of 4150 m/s which corresponds to that of pure Li at the same density and temperature as the alloy, as obtained in an AIMD calculation (see text). Fig. 9. Partial longitudinal current correlation functions, C L ij ðq;xÞ, for the liquid Li x Pb 1x alloy at three concentrations. Full blue line: C L LiLi ðq;xÞ, dashed red line: C L PbPb ðq;xÞ, dot-dashed line: C L LiPb ðq;xÞ(multiplied by a factor of five for the lowest q-value), and open circles: C L NN ðq;xÞ. M.M.G. Alemany, J. Souto-Casares, L.E. González et al. Journal of Molecular Liquids 344 (2021) 117775 8 C L ij ðq; x Þfor three concentrations and q-values. The C L LiLi ðq; x Þand C L PbPb ðq; x Þalways exhibit only one peak each, located at a high, x L LiLi ðqÞ, frequency, and low, x L PbPb ðqÞ, one, respectively. The C L NN ðq; x Þmay however display two peaks for some q-ranges and concentrations. In Fig. 9 the cross term has been enhanced multiplying it by a factor of five for q¼q min , in order to better observe its behavior. It is seen that the C L LiPb ðqÞalso show extrema, that can be of positive or negative sign. It is interesting to note that for these small qvalues the positive peak in C L LiPb ðqÞis located near x L PbPb ðqÞ, and that also near this frequency a shoulder can be observed in the Li-Li partial. This is the behavior expected for an acoustic mode where all the atoms vibrate in phase. On the contrary, the negative peak in C L LiPb ðqÞis located near x L LiLi ðqÞ, and this suggests that this vibration mode is optic-like, with particles of different type oscillating in opposite phase. This gives further support to the interpreation of the high frequency branch that shows up in the Li-Li partials as an optic-like non-hydrodynamic mode, while the low frequency branch, which appears clearly (as peaks) in the Pb-Pb partials, and weakly (as shoulders) in the Li-Li partials, has an acoustic-like nature, i.e., sound propagation. The longitudinal dispersion relations, x L ij ðqÞ, are depicted for three concentrations in Fig. 10, where we have only included the maxima of the functions, and not the shoulders, in order not to overload the figure. First, notice that the x L PbPb ðqÞtakes smaller values than those of x L LiLi ðqÞ, because of the greater atomic mass of the Pb atoms/ions. For all concentrations, the x L LiLi ðqÞand x L PbPb ðqÞhave just one branch whereas x L NN ðqÞexhibits two. The q-variation of the dispersion relations (apart from their low-qbehavior, which obviously is linear for the acoustic mode and going to a non-zero value for the optic-like one) is very much influenced by the corresponding partial structure factors. Minima in x L PbPb ðqÞroughly coincide with maxima in S PbPb ðqÞ. And a similar argument applies to the x L LiLi ðqÞwhich, for x Li ¼0:50 and 0:80, have a clear minimum, whereas for x Li ¼0:17 its almost monotonously increasing shape can be related to the weak double-peak structure of the corresponding S LiLi ðqÞ. The x L NN ðqÞdispersion curve has two branches with low and high frequencies. The high frequency branch exists for the whole q-range and practically coincides with the x L LiLi ðqÞbranch whereas the low frequency one closely follows the x L PbPb ðqÞbranch. For the concentration x Li ¼0:17 the low frequency x L NN ðqÞbranch exists for the whole q-range, closely follows the x L PbPb ðqÞbranch and goes to zero as q!0. However, when x Li ¼0:50 and 0:80, the low frequency x L NN ðqÞbranch exists only in a very narrow q-range. We have also included in Fig. 10 for x Li ¼0:80, the results deduced from the INS experiments by Alvarez et al. [7], which are very close to our simulation data, and the two propagating modes obtained by Anento et al. [4] through the application of the viscoelastic theory to CMD simulations. We recall that the low frequency mode behaved similarly to sound propagation (without the contribution of thermal fluctuations ignored in the viscoelastic approach) while the high frequency mode was identified as non-hydrodynamic optic-like. Despite the fact that the CMD used the approximate pair potential of Jacucci et al. [1] and the present simulations use the much more accurate DFT method to describe the forces, we find a qualitatively similar behavior of the x L ij ðqÞand those viscoelastic modes. We end up the study of the collective dynamic properties of the system considering the spectra of the partial transverse current correlation functons, C T ij ðq; x Þ, in their AL and BT formulations, and the frequencies where they show peaks, x T ij ðqÞ, which are shown in Fig. 11. In a similar way as in the case of the longitudinal counterparts, we find that for all concentrations, x T LiLi ðqÞshows as highfrequency branch, and x T PbPb ðqÞdisplays a low-frequency branch, but now only for qlarger than the long-wavelength propagation gap for acoustic transverse excitations, that takes approximate values of 0:9;0:5 and 0:3Å 1 for x Li ¼0:17;0:50 and 0:80, respectively. C T NN ðq; x Þ, however, always shows a high frequency peak, which is near x T LiLi ðqÞ, and outside the acoustic propagation gap also a low frequency peak, close to x T PbPb ðqÞ, but only in a limited q-range in the case of x Li P0:50. Even though a theoretically correct description of the acoustic transverse dispersion relation near the propagation gap involves a non-analytic expression involving a square root, it has also been shown that for qjust above the propagation gap the dispersion is quasilinear [59], so that it can be reasonably well fitted by a linear expression, x T NN ðqÞc T ðqq c Þ, where q c is effectively somewhat smaller then the real propagation gap, whereas the slope, c T , yields an estimate of the velocity of propagation of the shear modes in the alloy. Such a fit has produced values of c T 2300 150 m/s for x Li ¼0:17, c T 3270 150 m/s for x Li ¼0:50, and c T 1000  150 m/s for x Li ¼0:80. For comparison, we mention that in the limiting case of pure liquid Pb at T¼650 K we have obtained c T 800 100 m/s [31]. Notice that in this linear region x T NN ðqÞis close to x T PbPb ðqÞ, which suggests that the acoustic shear modes are propagating through the heavier Pb atoms only. Concerning the optic-like transverse branch we just mention that its long-wavelength limit is comparable to that of the corresponding longitudinal one, as expected for systems where there are no unscreened Coulomb interactions between particles. We have also evaluated the alloy shear viscosity, g . This magnitude plays an important role in processes such as fluid transport, Fig. 10. Longitudinal dispersion relations x L LiLi ðqÞ(blue squares) x L PbPb ðqÞ(red circles) and x L NN ðqÞ(open and dashed triangles) for the Li x Pb 1x liquid alloy at several concentrations. The diamonds with error bars are the experimental INS data. The dashed curves stand for the viscoelastic low and high frequency modes obtained in previous CMD calculations [4]. M.M.G. Alemany, J. Souto-Casares, L.E. González et al. Journal of Molecular Liquids 344 (2021) 117775 9