First principles determination of some static and dynamic properties of the liquid 3d transition metals near melting

We report an ab initio molecular dynamics simulation study of several static and dynamic properties of the liquid 3 d transition metals. The calculated static structure factors show qualitative agreement with the available experimental data, and its second peak displays an asymmetric shape which suggests a signiﬁcant local icosa-hedral short-range order. The dynamical structure reveals propagating density ﬂuctuations whose dispersion relation has been evaluated; moreover, its long wavelength limit is compatible with their respective experimental sound velocity. Results are reported for the longitudinal and transverse current spectral functions as well as for the respective dispersion relations. We also analyze the possible appearance of transverse-like low-energy excitations in the calculated dynamic structure factors. Several transport coeﬃcients have been evaluated and compared with the available experimental data.


Introduction
This paper presents recent ab initio molecular dynamics simulation studies for a range of bulk static, dynamic and electronic properties of the liquid 3d transition metals at thermodynamic conditions near melting. We review some of our previous calculations [1][2][3][4][5] along with several new results.
Despite the practical relevance of the 3d transition metals, it is remarkable that the experimental information concerning their static and dynamic properties in the liquid state is rather scarce.
The static structure factor, S(q), encompasses the information about the average short-range order of a liquid system, and was first obtained experimentally by Waseda and coworkers [6] in the decade of the 1970s for the 3d transition metals, by means of X-ray diffraction (XD) techniques. More recently, the S(q) of l-Ti, l-Fe, l-Ni and l-Zn have been determined again by either XD or neutron diffraction (ND) experiments [7][8][9][10][11]. In l-Zn, the new measurements showed a good agreement with the older ones, but for l-Ni and l-Fe, the comparison with Waseda and coworkers' XD (WXD) data for the S(q), revealed some discrepancies concerning the height of the main peak as well as the shape of the second peak; nevertheless, the positions of the maxima and minima of S(q) were coincident in all the experiments. On the other hand, the recent XD and ND data for the S(q) of l-Ti [9,10] yielded results well compatible between each other, but when compared with WXD data, they were displaced towards greater q-values by a significant amount of ≈ 0.20 Å −1 . Moreover, the recent XD data showed a second peak with a marked shoulder on the high-q side; actually, shoulders on the high-q side of the second peak of S(q) have also been found experimentally in other transition metals (Fe, Ni, Zr) [7,8,12]. We think that these disparities call for new precise structural data for the other 3d transition metals, both extracted from new additional experiments and/or from accurate simulations as the ab initio ones presented here.
To our knowlege, only l-Cu, l-Ni and l-Fe had been previously studied by other groups using ab initio simulation methods. For l-Cu, three ab initio studies evaluated its static properties at melting and some undercooled states [34,36,37]. As for l-Ni, Jakse et al. [38,39] studied a range of static and dynamic properties near melting and calculated the self-diffusion and shear viscosity coefficients. In the case of l-Fe, Ganesh and Widom [40] analyzed several static properties, getting also into the undercooled region.
In our recent studies of the liquid 3d transition metals [1][2][3][4][5], we were interested in the calculation and analysis of their dynamic properties. Besides its intrinsic interest, we were also motivated by the recent finding of some remarkable features such as low-energy transverse-like excitations in the dynamic structure factor and/or the appearance of a high-frequency branch in the transverse current dispersion relation, in addition to the usual one that appears at low frequencies. Indeed, the study of the microscopic processes behind the dynamic properties of liquid systems has been a long standing research field of Prof. Mryglod, and in this respect we spotlight his contributions to the development and application of the generalized collective modes (GCM) approach. As for the above mentioned features, there is still no clear explanation of their physical origin, although we have recently advanced some connection between the existence of a second branch in the transverse dispersion relation and the coupling of the transverse current with density fluctuations at all wavevectors [4,41].
Our studies of the bulk liquid 3d transition metals were performed with an ab initio molecular dynamics (AIMD) simulation method based on the density functional theory [42,43]. The liquid metal is modelled as an interacting system of ions and electrons where the ionic positions evolve according to classical mechanics while the electronic subsystem follows adiabatically. Table 1 provides information concerning the thermodynamic states as well as some technical details. The AIMD simulations were performed with the Quantum-ESPRESSO package [44,45], excepting l-Fe and l-Zn for which we used the VASP code [46][47][48][49]. The generalized gradient approximation of Perdew-Burke-Ernzerhof [50] was used to account for the electronic exchange-correlation energy, with the exceptions of l-Ti and l-Cu for which we used Perdew and Wang's approximation [51] and Perdew and Zunger's local density expression [52], respectively.
The ion-electron interaction was described by ultrasoft pseudopotentials [53], excepting l-Fe and l-Zn for which the projector augmented wave all-electron description [54,55] was used. For all the systems, the cutoff for the plane-wave representation of the orbitals was within the range 25.0-35.0 Ryd and the sampling of the Brillouin zone was performed by means of the single Γ point.
We recall that most previous studies of the liquid 3d transition metals were carried out by using semiempirical or more fundamental pair (or many-body) potentials coupled with liquid state theories or classical molecular dynamics simulations [56][57][58][59]. Moreover, in these approaches the electronic degrees Table 1. Input data of the liquid 3d transition metals considered in the present AIMD simulation study. ρ is the total ionic number density, (taken from [57]), T is the temperature, N part is the number of particles, ∆t is the ionic timestep, Z val is the number of valence electrons and N c is the total number of configurations. of freedom are hidden into the effective potential and, therefore, no electronic properties (i.e., density of states, conductivity, . . . ) can be consistently/reliably calculated. On the contrary, the AIMD simulation method yields parameter-free interatomic interactions derived from first principles, and the valence electrons are explicitly taken into account during the calculation.

Results and discussion
The evaluation of static and dynamic properties of the bulk liquid metal was performed by using the total number of equilibrium configurations listed in table 1.

Static properties
The calculated static structure factors, S(q), are plotted in figure 1 along with the corresponding WXD data [6]. The figure also includes the more recent XD and ND data [7][8][9][10][11] for l-Ti, l-Fe, l-Ni and l-Zn.
Aside from l-Sc and l-Ti, the calculated S(q) shows a good qualitative agreement with the WXD data. There are some slight disparities, namely: (i) a very small phase shift in l-V, l-Mn and l-Co, (ii) the calculations predict an asymmetric shape for the second peak with a more or less marked shoulder, whereas the XD data display, for all systems, a symmetric second peak. Notice that shoulders on the high-q side of the second peak of S(q) have also been experimentally observed in l-Ti, l-Fe, l-Ni and l-Zr [7][8][9][10]12]; moreover, they become more marked upon undercooling [60][61][62].
As for l-Ti, comparison with the recent XD and ND data [9,10] for the S(q), reveals an excellent agreement with the present resuts, including the shape of the second peak as well as the shoulder on its high-q side. Consequently, we believe that the WXD data for l-Ti are somewhat unreliable and that something similar might also happen with l-Sc, where an analogous (although opposite) dephasage is visible when compared with the present AIMD result.
From the calculated S(q), within the range q 1.2 Å −1 , we used a least squares fit, S(q) = s 0 + s 2 q 2 , to obtain an estimate for S(q → 0) and the results are given in table 2. Then, we evaluated the isothermal compressibility, κ T , by resorting to the relation S(q → 0) = ρk B T κ T , where k B is Boltzmann's constant. Table 2 lists the obtained results along with the available experimental data.  present AIMD calculations. Open circles: WXD data from [6]. The blue circles are XD data [7][8][9][10][11] and red ones are ND data [7][8][9][10]. The inset shows the corresponding pair distribution function.  [13] In figure 1, we also depicted the associated pair distribution functions, g(r), which are compared with their respective XD data [6]. The associated coordination number (CN) was evaluated by integrating the radial distribution function, 4πρr 2 g(r), up to the position of its first minimum, r min . Table 2 shows the obtained values for the CN's which are typical of simple liquid metals near their respective triple points [35]. A more detailed three-dimensional description of the short range order in these liquid metals is achieved by applying the common neighbour analysis [63] (CNA) method. This method allows one to characterize the local environment surrounding each atomic pair that contributes to the peaks of the g(r), in terms of the number and properties of the common nearest neighbours of the pair under consideration. Each bonded pair of atoms is characterized by four indices, with the first index being 1 if the pair belongs to the first peak of g(r), the second index stands for the number of common first neighbours, the third index is the number of bonds that connect those shared first neighbours and the fourth index discriminates among those configurations having the same first three indices but a different topology. The CNA method allows one to discern different local structures such as FCC, HCP, BCC and ICOS. For example, the close-packed (CP) structures, FCC and HCP, are composed of 142x-type pairs, the BCC is typified by 144x and 166x pairs whereas the perfect ICOS is characterized by 155x pairs and the distorted ICOS is described by the 154x and 143x pairs.
The CNA study is performed on inherent structures in which the ions are brought to local minima of the potential energy surface. For each metal, this calculation was made for four/five ionic configurations, well separated in time (≈ 15.0−20.0 ps), and the results were averaged. These are plotted in figure 2 where it is noticed that the five-fold symmetry dominates, as the sum of perfect and distorted ICOS structures varies between ≈ 59% (l-Sc) and ≈ 72% (l-Cr, l-Co) of the pairs. The amount of local BCC-type pairs is also remarkable as it ranges from ≈ 32% (l-Ti) to ≈ 6% (l-Zn) of the pairs. Finally, the CP-type pairs are almost absent in the early transition metals but become significant for Mn and for the late elements (Co to Zn). This behaviour of the CP vs BCC local structures correlates with the corresponding phase diagrams of the elements, since Sc to Fe melt from a BCC phase, while Co to Zn do so from CP structures, either FCC (Co, Ni, Cu) or HCP (Zn). The phase diagram of Mn is rather complex, and the high temperature BCC phase only exists down to 100 K below the melting point, where it transforms into FCC, and this latter phase appears to somehow partially survive locally also after melting.

Dynamic properties
We evaluated several dynamic magnitudes and their associated time correlation functions. Due to the macroscopically isotropic behaviour of the fluid, those correlation functions with a ì q-dependence are transformed into | ì q|-dependent magnitudes.

Single-particle dynamics
We report results for the self-diffusion coefficient, D, which, along with the shear viscosity, plays an important role in the study of crystal nucleation and growth in metallic melts. We evaluated the velocity autocorrelation function, Z(t), and the mean square displacement, δR 2 (t), of a tagged ion in the fluid, for all the 3d liquid metals. The diffusion coefficient can be determined through both of these functions, and both approaches yielded practically the same result for D, which is given in table 3 along with the available experimental/semiempirical data.

Collective dynamics
In a liquid, the dynamics of density fluctuations is usually described by the intermediate scattering function, F(q, t), which is defined as autocorrelation function of the microscopic wavevector-dependent density [35]. The Fourier transform (FT) of F(q, t) into the frequency domain leads to the dynamic structure factor, S(q, ω), which can be directly measured by either INS and/or IXS experiments.
Another important collective magnitude is the current due to the overall motion of the particles. It is a vectorial magnitude which is usually split into its longitudinal and transverse components with respect to ì q. Then, the longitudinal, J L (q, t), and transverse J T (q, t), current correlation functions are obtained as autocorrelation functions of the corresponding quantities [35]. Their time FT yields the associated spectra, J L (q, ω) and J T (q, ω), respectively.
In figures 3-4, we plotted, for some metals and for some range of q-values, the obtained AIMD results for the F(q, t). For small q-values, the F(q, t) show an oscillatory behaviour, indicative of wave propagation, which is superposed on another decaying component, indicative of relaxation modes. The oscillations become weaker with increasing q-values and disappear at around q ≈ 4/5q p , where q p denotes the position of the main peak of S(q). For larger q values, the decaying component becomes dominant leading at q ≈ q p to a very slow monotonous decay of the F(q, t). The relative strength of both 23606-5  components depends on the system and we observe that l-Mn has the weaker oscillations superposed on a very strong diffusive component. Figures 5-6 show the associated S(q, ω). The behaviour of the calculated S(q, ω) is qualitatively similar for all the 3d metals studied in this paper. Specifically, the S(q, ω) show side-peaks up to q-values ≈ (3/5)q p , which thereafter evolve into shoulders lasting up to q ≈ (4/5)q p . Therefrom, the S(q, ω) show a monotonous decreasing behaviour. Notice that for small q-values, the S(q, ω) of l-Mn takes relatively large values when ω → 0; this is due to the important diffusive component of its F(q, t) at those small q-values. From the positions of the side-peaks, ω m (q), in the S(q, ω), a dispersion relation for the density excitations was obtained and is plotted in figure 7. In the hydrodynamic region (small q), the slope of the dispersion relation curve is the q-dependent adiabatic sound velocity c s (q) = v th γ/S(q), with v th = k B T/m being the thermal velocity and γ being the ratio of specific heats. In the q → 0 limit, the c s (q) reduces to the bulk adiabatic sound velocity c s . However, the small size of the simulation box  implies that the smallest attainable q value, namely q min , is not small enough so as to permit an accurate determination of c s ; nevertheless, an approximate estimation can be obtained from value of ω m (q min )/q min . We also evaluated the longitudinal and transverse currents, J L (q, t) and J T (q, t), as well as their respective spectra, J L (q, ω) and J T (q, ω). The obtained J L (q, ω) show one peak for each q-value and the frequencies associated to those peaks, ω L (q), constitute the longitudinal dispersion relation for the associated collective modes. Figure 7 represents the obtained dispersion relations, which show the typical behaviour found in other liquid systems [35].
As for the J T (q, ω), it may exhibit peaks within some q-range, which are related to propagating shear waves. For all the metals studied in this paper, the associated J T (q = q min , ω) already showed a peak and with increasing q-values the associated frequency of the peak, ω T (q) also increased and reached a maximum value for q ≈ (2/3)q p ; therefrom it decreased and vanished at ≈ 3.0q p . An example of this behaviour is given in figure 8 which shows the calculated J T (q, ω) for l-V and l-Cr, as a function of ω, and for a range 0 < q 4.0 Å −1 . Moreover, for some metals (l-V, l-Cr, l-Co, l-Cu, l-Ni and l-Zn) we found that the J T (q, ω) exhibited, within the q-region 0.8 q/q p 1.2, another peak with a higher frequency. Figure 8 shows, for l-V and l-Cr, a closer view of that specific region where the two maxima are visible. For the other metals, (i.e., l-Sc, l-Ti, l-Mn and l-Fe) the results are not conclusive because, besides the lower frequency peak, we also found higher frequency shoulders instead of peaks.
The appearance of another, high frequency branch in the transvese dispersion relation, had been reported first in some liquid metals at high pressure, i.e., Li, Na, Fe, Al and Pb, [2,[64][65][66]. Subsequently, it was also found in l-Tl and l-Pb [64,66,67] at ambient pressure.
From the calculated J T (q, t), we also estimated [35] the shear viscosity coefficient η, as follows. The memory function representation of J T (q, t), namelỹ where the tilde denotes the Laplace transform, introduces a generalized shear viscosity coefficientη(q, z).
The area under the normalized J T (q, t) gives mJ T (q, z = 0)/(k B T), from whichη(q, z = 0) ≡η(q) can be obtained. Extrapolation to q → 0 gives the standard shear viscosity coefficient η. Different extrapolation functions lead to different values of η, which are portrayed in the uncertainty assigned to the results. Table 3 shows these results along with the available experimental/semiempirical data.
Recently, it was suggested [68] that transverse-like low-energy excitations might be observed as weak shoulders located in the region between the quasielastic and the longitudinal inelastic peaks of the S(q, ω). These excitations are usually visible within a small q-range around q p /2, because for smaller/greater q-values they are overcome by the quasielastic/inelastic peaks. These transverse-like excitations were first detected in the IXS spectra of l-Ga [68] at 313 K and subsequently in l-Cu, l-Sn, l-Na, l-Fe, l-Ni and l-Zn [3,29,30,[69][70][71]. We analyzed the results for S(q, ω) and found that within the range 1.00 Å −1 q 1.45 Å −1 , some weak shoulders appear in the ω-region located between the quasielastic and the inelastic peaks. This is shown in figure 9; moreover, the energies associated to these shoulders are close to those corresponding to the peaks in the transverse current spectra. The real physical origin of these excitations is still a matter of debate, and some alternative interpretations, for instance heat waves, were proposed within the framework of the GCM approach [72].

Electronic properties: density of states
The partial and total electronic density of states, n(E), were obtained from the self-consistently determined eigenvalues; they were averaged over four/six ionic configurations well separated in time (≈ 15.0 ps), and the sampling of the Brillouin zone was performed using an 8 × 8 × 8 Monkhorst-Pack set. Figure 10 shows the obtained results for the electronic total n(E) associated to the outer valence electrons, which is dominated by the 3d states, except in the early transition elements, where the d band    occupation is low. Notice that for the metals Sc to Mn we are not showing the contribution from the 3s and 3p atomic states, which were included in the self-consistent calculation but lie way below the Fermi energy. In the figure we observe the progresive filling of the d band, which is complete at Cu, together with a decrease of its width and an increase of its height as the occupation rises.

Conclusions
An ab initio molecular dynamics simulation method was used to calculate several static and dynamic properties of the bulk liquid 3d transition metals near their respective melting points.
The results for the static structure, as characterized by the static structure factor S(q), show a good agreement with the available experimental data. The only exception is l-Sc where a clear difference is observed in the phase of the oscillations, but this may be due to inaccuracies in the old experimental data, similar to the case of the old data for l-Ti. The S(q) display, in most cases, an asymmetric shape in the second peak which qualitatively agrees with the recent ND and/or XD data. This feature is related to the appearance of icosahedral short-range order in the liquid metal. This is confirmed by a more detailed study of the liquid static structure by using the CNA method, which revealed a marked abundance of five-fold structures. On the other hand, secondary local structures also present in the liquid correlate rather well with the solid phase from where the metals melt.
The calculated dynamic structure factors, S(q, ω), show side-peaks which are indicative of collective density excitations. A detailed analysis of the S(q, ω) revealed some type of excitations which have similar features as the transverse-like excitation modes found in IXS and INS data for several liquid metals near their respective melting points.
From the calculated longitudinal and transverse current correlation functions, we evaluated the associated spectral functions and the corresponding dispersion relations. For some metals, we obtained two branches of transverse modes, with the high frequency branch appearing over a limited q-range which is always located in the second pseudo-Brillouin zone.