Exploring Challenging Properties of Liquid Metallic Systems through Machine Learning: Liquid La and Li4Pb Systems

In this machine learning (ML) study, we delved into the unique properties of liquid lanthanum and the Li4Pb alloy, revealing some unexpected features and also firmly establishing some of the debated characteristics. Leveraging interatomic potentials derived from ab initio calculations, our investigation achieved a level of precision comparable to first-principles methods while at the same time entering the hydrodynamic regime. We compared the structure factors and pair distribution functions to experimental data and unearthed distinctive collective excitations with intriguing features. Liquid lanthanum unveiled two transverse collective excitation branches, each closely tied to specific peaks in the velocity autocorrelation function spectrum. Furthermore, the analysis of the generalized specific heat ratio in the hydrodynamic regime investigated with the ML molecular dynamics simulations uncovered a peculiar behavior, impossible to discern with only ab initio simulations. Liquid Li4Pb, on the other hand, challenged existing claims by showcasing a rich array of branches in its longitudinal dispersion relation, including a high-frequency LiLi mode with a nonhydrodynamic optical character that maintains a finite value as q → 0. Additionally, we conducted an in-depth analysis of various transport coefficients, expanding our understanding of these liquid metallic systems. In summary, our ML approach yielded precise results, offering new and captivating insights into the structural and dynamic aspects of these materials.


INTRODUCTION
Computer simulation techniques play an important role in providing new insights into many problems concerning condensed matter physics and materials science as well as an additional explanation and confirmation of experimental results.
Starting with the groundbreaking work of Alder and Wainwright 1 on the fluid to solid phase transition of a model hard sphere system, the molecular dynamics (MD) simulation technique has undergone a huge development.The next step into this process was to account for the interactions beyond simple elastic collisions so as to provide a more realistic description of the dynamic behavior of atoms and/or molecules.This was achieved by Rahman 2 who implemented a classical molecular dynamics (CMD) simulation, where the classical equations of motion were solved for a realistic model of a physical system, made up of real (argon) atoms interacting via a Lennard-Jones-type pair potential.Subsequently, this CMD approach was extended to more complex systems by introducing more elaborated types of interatomic potentials intended to describe the physics behind the interactions among the atoms.Nevertheless, in certain cases, such as systems where the electronic and atomic structures are closely interwoven, relying solely on interatomic empirical potentials is insufficient.It becomes imperative to incorporate an explicit description of the electronic properties into the simulation.
−5 This breakthrough enabled the direct determination of the potential responsible for the interactions in a many-body system, starting directly from the electronic configuration of its components, thereby making possible the simulation of covalently bonded and metallic systems.This approach has proved to be so accurate that the ab initio molecular dynamics (AIMD) methods based on DFT have now become the usual technique for the study of a wide range of condensed matter systems.Most AIMD methods are based on the Kohn−Sham orbital representation of DFT; however, its application poses heavy computational demands and therefore AIMD simulations are severely limited by the finite system size (∼100 s of atoms) and short time scales (∼10 s of ps).However, there are some physical magnitudes whose accurate determination would involve the use of length and time scales which are not currently accessible by these AIMD simulation methods.For example, and within the realm of liquid systems, we mention the long wavelength limit of a range of static and time-dependent correlation functions as well as some transport coefficients.A proper account of these physical magnitudes requires the use of large number of atoms (tens of thousands) along with much longer simulation times (∼100− 1000 s of ps), impossible nowadays with AIMD methods, but achievable by resorting to CMD simulations based on interatomic potentials at the cost of a significantly reduced accuracy.
Over the past decade, machine learning (ML) interatomic potentials have been developed as a way of addressing the abovementioned problems.By employing a database of quantum-mechanical (commonly DFT-based) calculations, it is possible to construct interatomic potentials that allow system sizes and simulation times typical of CMD simulations but with quantum mechanical accuracy.
We have chosen two distinct systems, liquid lanthanum (l-La) and liquid Li 4 Pb, to show the versatility and capability of ML potentials to study various aspects of liquid metal behavior inaccessible to AIMD.
Two basic magnitudes related to the structural short-range order in a liquid metal are the static structure factor, S(q) and the pair distribution function, g(r).
In the case of l-La near melting, the available experimental structural data were obtained more than 40 years ago by Breuil and Tourand 6 by means of neutron diffraction (ND) and, a few years later, by Waseda and Tamaki 7 through X-ray diffraction (XD).Since then, no additional measurements of S(q) have been reported for this metal.As for the thermophysical magnitudes of l-La, there are experimental data for the adiabatic velocity of sound 8 and shear viscosity. 9On the theoretical side, few works have studied the structural properties of l-La, and among them, we mention the study by Waseda et al. 10 who used effective interionic pair potentials, derived within the pseudopotential perturbation theory, combined with the Percus−Yevik integral equation to calculate static structure, i.e., S(q) and g(r), of some liquid rare-earth metals.More recently, Patel et al. 11−13 used a pseudopotential combined with the hard sphere reference system to evaluate the S(q) and g(r) as well as some thermodynamic and electronic properties of the liquid lanthanide metals.
The Li 4 Pb liquid alloy has already attracted a great deal of both theoretical and experimental work.In common with other binary alloys with a large mass disparity between both components, it exhibits a high-frequency mode with an associated phase velocity much greater than the extended hydrodynamic sound as found in both MD studies 14−18 and inelastic neutron scattering (INS) experiments. 19,20In fact, this high-frequency mode was first found in a CMD study of Li 4 Pb 14,15 using interionic pair potentials modeled with a hard core repulsion plus a screened Coulombic interaction.These results showed that the high-frequency mode produced a highfrequency side peak appearing in the Li−Li partial dynamic structure factor.It was interpreted as an extended hydrodynamic mode with an associated phase velocity of ≈7500 m/s and was named the "fast sound" mode.
Additional CMD simulations using the same interionic interaction but with a larger number of particles 16 confirmed the previous results and were able to discern two branches of collective excitations, as obtained from the maxima in the longitudinal current correlation spectra, that were identified as the fast sound mode and a new slow sound mode.It was concluded that these two modes merged, at very low q-values, into the usual hydrodynamic sound.Both modes showed a sharp transition in their phase velocities around q c ≈ 0.08 Å −1 , with an abrupt increase for q > q c toward a "fast velocity" around 8000 m/s in the case of the high-frequency mode and a decrease toward a "slow velocity" around 1200 m/s in the case of the low frequency one.However, another study using CMD simulations in combination with the generalized collective modes (GCM) method 17 found two branches of collective excitations which in the q → 0 limit were identified as a hydrodynamic (lowfrequency) branch and an optic-like (high frequency) branch that tends to a nonzero frequency instead of merging with hydrodynamic sound.Similar conclusions were also obtained from CMD simulations analyzed within the viscoelastic model by Anento et al. 18 On the other hand, the INS experiments, 19,20 that were limited to q larger than 0.6 Å −1 due to kinematic conditions, confirmed the existence of this high-frequency mode but with an associated phase velocity around 4500 m/s.Moreover, they suggested that this branch had features that pointed to a kinetic (nonhydrodynamic) mode with a nonzero value in the q → 0 limit.
There are therefore two contradicting interpretations on the behavior of the high-frequency mode as the wavelength increases (either it merges with hydrodynamic sound or it behaves as a kinetic mode).The comparison with the experiment is not conclusive due to the wavevector region explored and moreover is somewhat hindered by the use of an effective potential in the CMD simulations that may not capture all the details of the true interactions among atoms in the alloy.Recent first-principles simulations with realistic interactions derived from electronic structure calculations 21,22 appear to support the optic-like view, but in this case, the problem lies in the value of the smallest wavevector allowed by the periodic boundary conditions inherent to the calculations, which is 0.32 Å −1 .This value is a direct consequence of the small number of particles affordable in such high cost simulations and is quite large so as to make conclusive statements on the behavior of the modes for q → 0.
This paper reports a simulation study of several static and dynamic properties of l-La at thermodynamic conditions just above melting at ambient pressure and the liquid Li 4 Pb alloy just above its liquidus temperature.The study has been performed first using AIMD simulations, then creating the corresponding ML interatomic potentials based on the aforementioned AIMD simulations, and finally using these ML potentials in CMD runs with much larger time and length scales.
An important aim of this study is to calculate, using ML simulation methods, the spectrum of collective excitations in the liquid Li 4 Pb alloy and analyze them, with special emphasis in its controversial high-frequency branch, and in particular its behavior in the long wavelength region.For l-La, we provide the first AIMD study and use the ML derived interatomic potential to better analyze the results in the hydrodynamic Journal of Chemical Theory and Computation regime, uncovering an unconventional behavior in the generalized specific heat ratio.
The review is structured as follows: the next section summarizes the basic ideas underlying the AIMD and ML simulation methods along with some technical details.In Section 3, we report the results of the calculations for both l-La and the Li 4 Pb alloy, which are compared with the available experimental data along with some discussion.Finally, a brief summary and conclusions are given in Section 4.

COMPUTATIONAL METHOD
The reference database of l-La was generated by AIMD simulations performed using the DFT-based QUANTUM ESPRESSO (QE) package. 23,24Within this framework, the electronic exchange-correlation energy was described by the generalized gradient approximation (GGA) of Perdew−Burke− Ernzerhof. 25The ion−electron interaction was accounted for by means of an ultrasoft pseudopotential generated from a scalarrelativistic calculation including nonlinear core corrections.In this pseudopotential, the 5s, 5p, 5d, 6s, and 6p orbitals are explicitly considered as valence states, and 11 valence electrons per atom were employed during simulations.The kinetic energy cutoff for the plane-wave expansion of the wave function was 530 eV, and the single Γ point was used for sampling the Brillouin zone.
−29 We used the projector augmented wave allelectron description of the electron−ion-core interaction provided by the VASP distribution, which considers one and four valence electrons for Li and Pb, respectively; moreover, the exchange correlation functional used was again GGA, while the kinetic energy cutoff for the plane-wave expansions was taken as 480 eV.
We have employed the SIMPLE-NN package 30 to develop deep neural network potentials (DNNPs) which, due to their high flexibility, are a suitable model to represent the challenging problem of the potential energy surface of a liquid metal or alloy.To transform the atomic coordinates into atomic descriptors, SIMPLE-NN employs Behler−Parrinello Gaussian functions. 31,32For l-La, we used the Gastegger method 33 to select the parameters of the 35 total Gaussian functions employed, whereas for Li 4 Pb, the Imbalzano method 34 provided better parameters for the 36 Gaussian functions of each element.
The choice of either method (Gastegger or Imbalzano) is based on the better performance obtained when applying the learning process to the training data sets.The main difference between both methods lies in the expressions employed to select the widths and positions of the Gaussians.In particular, the Gastegger method uses a fixed width for the shifted Gaussians that is equal to the distance between the centers of the two adjacent functions, whereas the Imbalzano method increases the width as the distance increases, effectively creating a finer grid closer to the central atom, where small variations of the position can have a greater impact on the potential.
The neural network employed for the DNNP of l-La consisted of three layers, with 50 neurons per layer.The loss function included the energy per atom, atomic forces, and stress.The DNNP was trained on more than 11,700 configurations, using a training and validation split of 80:20.
In the case of the DNNP of liquid Li 4 Pb, the neural network architecture was composed of four layers with 50 neurons each.The training configurations included data from four different concentrations of the Li−Pb alloy: Li 0.17 Pb, Li 0.50 Pb, Li 0.62 Pb, and Li 0.80 Pb.We selected 200 configurations from each concentration, taken from a previous AIMD study, 21 and the corresponding energies and forces were computed with VASP.However, due to extrapolation errors that appeared when testing the DNNP, additional AIMD runs were performed for liquid Li 0.17 Pb and liquid Li 0.80 Pb, maintaining the temperature but increasing the atomic density by a 20% and by a 40% so as to allow for configurations where shorter interatomic distances were present.From these new AIMD runs, a total of 190 additional configurations were included in the training of the final DNNP.Due to the small number of total configurations, a training and validation split of 90:10 was employed.The loss function included the energy per atom and the atomic forces.
While in principle, it would seem unnecessary to include additional compositions apart from the 80% Li in the training data set in order to develop the NN potential, in fact it turns out to be profitable to do so for two reasons: first because the instantaneous local atomic arrangements in the alloy can fluctuate departing from the nominal alloy concentration, especially when considering samples with a large number of atoms; and second, from a more practical viewpoint, because such strategy will allow the future study of the alloy at different concentrations using the same DNNP.
Note that the inclusion of the stress in the loss function for l-La, as well as the larger number of configurations included in the data set for that system, as compared to the case of Li 4 Pb, is justified by the higher degree of complexity of the interactions in this system, as evinced by the much higher number of valence electrons needed in the AIMD study in order to produce accurate forces and energies.
Both DNNPs employed fully connected layers, with the hyperbolic tangent as the activation function, were trained using a mini-batch of 30, used random sampling along with Adam optimizer with a learning rate of 0.001, and employed an L 2 regularization of 10 −6 .The mean-squared-error was employed as the objective function during all training.The coefficients for various quantities when included in the loss function were as follows: 1.0 for the energy per atom, 0.1 for the atomic forces, ρ is the total ionic number density, T is the temperature, N is the number of particles in the simulation cell, Δt is the ionic timestep, and N c is the total number of configurations.

Journal of Chemical Theory and Computation
and 10 −4 for the pressure.Once the DNNPs were trained, they were utilized in LAMMPS 35 to perform the large-scale MD simulations.
Table 1 provides details of the different kinds of simulations performed for the two systems considered in the study.In all cases, the simulation samples with N particles were enclosed inside a periodic cubic cell with a size determined by the experimental ionic number density.After thermal equilibration, the number of equilibrium configurations, N c , given in Table 1 was subsequently used for the evaluation of the static and dynamic properties of both systems.

RESULTS AND DISCUSSION
3.1.Static Properties.3.1.1.Liquid La. Figure 1 depicts the calculated static structure factor, S(q) as obtained from the AIMD calculation along with the ML result and the XD data of Waseda et al. 7 Notice that the AIMD and ML results are practically identical, and the comparison with the XD data shows agreement for the position of the main peak, q p = 2.1 Å −1 , although the calculated amplitudes are more marked; moreover, both simulations predict an asymmetric shape of the second peak of the S(q) which is not found as pronounced in the XD data.Nevertheless, we notice that this feature has been experimentally observed in several liquid metals 36−38 and has been related to the presence of a substantial amount of icosahedral order in the liquid.
Concerning the evaluation of q-dependent magnitudes, the use of ML interatomic potentials allows us to reach considerably smaller q-vectors, allowing for a more detailed investigation into the long-wavelength (q → 0) region.The smallest attainable value by the ML method was q = 0.099 Å −1 , whereas in the case of the AIMD calculation, the smallest value was q = 0.400 Å −1 .We have exploited this feature of the ML method to estimate the isothermal compressibility, κ T , by resorting to the relation S(q → 0) = ρk B Tκ T , where k B is Boltzmann's constant.First, the ML-calculated S(q) has been extrapolated to q → 0 by a leastsquares fit, S(q) = s 0 + s 2 q 2 , of the q-values for q ≤ 1.0 Å −1 .It yields a value of S(q → 0) = 0.0160 ± 0.001, which gave κ T = 3.60 ± 0.20 (in units of 10 11 m 2 N −1 ).Other evaluations for this magnitude at the melting point have been proposed by some authors; thus, McAlister et al. 8 estimated a value of κ T = 4.24 from speed of sound measurements, while Blairs 39,40 suggested a value κ T = 4.30 obtained by using an empirical formula connecting isothermal compressibility, surface tension, and the Ornstein−Zernike correlation length.However, no direct experimental measurement is available to date.
Figure 2 shows the pair distribution function, g(r), which provides some insights into the short-range order in the liquid.
Both AIMD and ML results are depicted along with the available experimental data. 6,7There is a fair agreement concerning the position of the oscillations, but the calculated amplitudes are greater with both AIMD and ML.
The average number of nearest neighbors (also called the coordination number, CN) around any given ion is evaluated by integrating the radial distribution function, 4πρr 2 g(r), up to the position of its first minimum, r min ≈ 4.97 Å.The calculated CN with AIMD and ML is ≈13.0, similar to those obtained for other simple liquid metals near melting. 41he common neighbor analysis 42−44 (CNA) method provides additional insights into the short-range order as it gives three-dimensional information about the ions surrounding each pair of ions which are near neighbors.Each pair is characterized by four indices which allow us to discern among different local structures like fcc, hcp, bcc, and icosahedral environments.The present calculations (both AIMD and ML) show that the five-fold symmetry dominates in l-La because the sum of perfect and distorted icosahedral structures reaches a value of ≈53% of the pairs with the number of perfect ones being around 30%.The amount of local bcc-type pairs is also significant, ≈27%; however, there is virtually no vestige of fcc and hcp-type pairs.
3.1.2.Liquid Li 4 Pb Alloy. Figure 3 depicts the ML calculated partial pair distribution functions g ij (r) and partial static structure factors S ij (q) along with the corresponding AIMD results.Both sets of results are practically identical, with very minor discrepancies concerning the amplitudes of the g PbPb (r).
We have evaluated the isothermal compressibility in the alloy by using the expression 45 where ρ is the number density of the alloy, x i is the concentration of component i, and S ij (0) are the q → 0 values of the partials.These latter values were derived by extrapolating the MLcalculated S ij (q) in the same way that was performed for l-La.Now the smallest wavevector attained was q = 0.080 Å −1 .Thus, we obtained (with an uncertainty ≈±0.001) S LiLi (0) = 0.0715,  S LiPb (0) = −0.0275and S PbPb (0) = 0.0920.The application of eq 1 yielded a value of κ T = 8.18 ± 0.20 (in 10 11 m 2 N −1 units).
There are no experimental data for this magnitude, but we will later compare it with the data for the adiabatic compressibility.
Finally, we note that information about the ordering tendencies in the alloy is given by the long-wavelength limit of the Bhatia− Thornton (BT) concentration−concentration partial structure factor, S cc (q), which can be evaluated in terms of the previous partials. 45For this alloy, we have obtained S cc (q → 0)/(x 1 x 2 ) = 0.0975, which clearly indicates a strong heterocoordinating tendency in the liquid Li 4 Pb alloy.

Collective Dynamics. 3.2.1. Liquid La.
The collective dynamics of density fluctuations in a liquid can be described by the intermediate scattering function, F(q, t), which is the autocorrelation function of the microscopic q-dependent number density. 41Its time Fourier transform (FT) gives the frequency spectrum known as the dynamic structure factor, S(q, ω).
Figure 4 shows, for several q-values, the AIMD results for F(q, t) along with the corresponding ML results.The calculated F(q, t) show a structure qualitatively similar to other liquid metals: at small q's, there is an oscillatory form, suggesting wave propagation, which is combined with a decaying term that indicates the existence of relaxation modes.The interplay between these two terms evolves with increasing q-values with the decaying term smoothly and progressively overcoming the propagating mode, leading at q ≈ q p to a very slow monotonic decay of F(q, t).
Further insights into the physical causes behind the propagation and relaxation mechanisms can be obtained from theoretical models where the F(q, t) is expressed using the formalism of memory functions, 41 more specifically its secondorder memory function N(q, t), i.e.
where F ̃(q, z) and N ̃(q, z) are the Laplace transforms of F(q, t) and N(q, t), respectively, M 0 (q) = k B Tq 2 /(mS(q)), and m is the atomic mass.The N(q, t) can be modeled by an analytical function with two exponentially decaying terms (slow and fast ones), i.e.

= + N q t
A q A q ( , ) ( )e ( )e where τ s (q) and τ f (q) represent slow and fast relaxation times, respectively.−48 We used F(q, t) obtained by the AIMD and ML simulations in order to evaluate the associated functions N(q, t).Then, they were fitted to eq 3 and were analyzed to find out whether they are consistent with a generalized hydrodynamic model (fast viscoelastic mode and a slow thermal one) or a generalized viscoelastic model, where the fast term is the thermal one.Specifically, we evaluated the generalized heat capacity ratio, γ(q), which in the q → 0 limit leads to the thermophysical value γ 0 , i.e., the ratio between the specific heats at constant pressure and constant volume.
If the slow mode is associated to the thermal relaxation, then A s (q) = (γ(q) − 1)M 0 (q), but if it is connected to the viscoelastic relaxation, then A s (q) = ω L 2 (q) − γ(q)M 0 (q), 47,48 where ω L 2 (q) is the second frequency moment of the longitudinal current correlation function (see below), which is obtained directly from the simulations.
By exploring both possibilities, we have calculated the functions γ th (q) and γ v (q), which represent the values obtained for γ(q) when either the thermal or the viscoelastic relaxations proceed by the slow mode, respectively.
The results are depicted in Figure 5, which shows the obtained γ th (q) and γ v (q) values within the range q/q p ≤ 0.70 Å −1 .Both functions exhibit a very different conformation, with γ v (q) taking values which are always greater than those of γ th (q); moreover, when q → 0, the γ v (q) becomes unphysically large.As for the γ th (q), we observe that the AIMD results show a monotonous increasing behavior down to the smallest attainable q-value (q/q p = 0.19 Å −1 ); however, the ML values show that a noticeable change around that same q-value and the associated γ th (q) displays a slowly decreasing behavior toward q = 0, with a γ th (q → 0) ≡ γ 0 ≈ 1.25.
We are not aware of any experimental data for the γ(q) of l-La but we note that the qualitative structure of the ML calculated γ th (q) is very similar to that reported in the experimental data of Hosokawa et al. 49 for the γ(q) of l-Fe near melting.Specifically, the γ exp,Fe (q) begins with a value of ≈1.40 at q/q p,Fe ≈ 0.07 Å −1 and smoothly increases reaching a value of ≈1.55 at q/q p,Fe ≈ 0.21 Å −1 and then decreases toward a value of ≈1.1 at q/q p,Fe ≈ 0.67 Å −1 (where q p,Fe ≈ 2.98 Å −1 ).Although there are no experimental data for γ(q → 0) ≡ γ 0 of l-La near melting, we can compare it with two estimates, namely, γ 0 = 1.23 39,40 and 1.31, 50 which were obtained by combining the experimental data of the adiabatic compressibility and the calculated isothermal compressibility.
The previous results suggest that the slow mode is associated with the thermal relaxation, and therefore, the generalized hydrodynamic model seems to be the appropriate one for describing the microscopic dynamics of l-La; moreover, it predicts a behavior for γ(q) which is qualitatively consistent with the experimental data for l-Fe and, furthermore, it also provides an estimate for γ 0 that compares well with some semiempirical values.
The dynamic structure factors, S(q, ω), are plotted in Figure 6, and show side peaks up to q ≈ (3/5)q p , whereas for greater qvalues, we observe a monotonic decreasing behavior.
We have determined the frequency of the side peaks as a function of the wavevector, namely, the function ω m (q), which represents the associated dispersion relation of the density excitations; therefrom, a q-dependent adiabatic sound velocity c s (q) = ω m (q)/q is obtained and has been depicted in Figure 7. Again, the AIMD and ML simulation results are practically identical but only the latter ones are able to clearly indicate the existence of positive dispersion.In the c s (q → 0) limit, this magnitude reduces to the bulk adiabatic sound velocity, c s , and by extrapolating the previous ML simulation results, we have obtained an estimate c s = 2230 (±150) m/s which qualitatively agrees with the experimental data c s,exp = 2022 m/s. 8nother important magnitude in the collective dynamics is the current due to the overall motion of the particles, j(q, t) where N is the total number of particles and v j (t) is the velocity of particle j at time t. 41From its longitudinal, j L (q, t), and transverse, j T (q, t), components, we have evaluated the associated longitudinal and transverse current correlation functions T T (6)   and their respective time FT, the associated spectra C L (q, ω) and C T (q, ω).
For any fixed q-value, when the C L (q, ω) is plotted as a function of ω, we observe a maximum, and its associated frequency, namely, ω L (q), stands for the dispersion relation of the longitudinal modes.These are plotted in Figure 8 where it is observed that the AIMD and ML simulation results are very similar, although the ML approach again provides closer insights into the small q-region.
Figure 9 depicts, for a range of q-values, the calculated AIMD and ML results for C T (q, t).For both small and large q, this function must be monotonically decreasing with time, but for intermediate q values, it shows oscillations.
Correspondingly, its spectrum, C T (q, ω), when plotted as a function of ω, shows within some q-range clear peaks that are related to propagating shear waves.According to the present ML results, the peaks start at q ≈ 0.15 Å −1 and last until q ≈ 3.0q p    when they fade away.Moreover, we have found a second, higher frequency peak that shows up over a much smaller q-range.
Some examples of the behavior of C T (q, ω), along with the corresponding longitudinal spectrum, C L (q, ω), are shown in Figure 10, where the appearance and later disappearance of the high-frequency peak for increasing q is observed.Moreover, it is also seen that the peak in the longitudinal spectrum at a given q is mostly unrelated to the peaks in the transverse one at that same value of q.
The high-and low-frequency dispersion relations for the transverse modes are plotted in Figure 11.Some proximity is observed between the values of the high-frequency transverse peak, which is almost nondispersive, and the maximum in the longitudinal dispersion relation, as shown in Figure 8, which is located around q p /2.Such behavior has been previously observed in other systems 51−53 and explained in terms of mode-coupling theories.This type of theories describes an indirect coupling between the transverse current at the wavevector considered, q, and the longitudinal magnitudes at all wavevectors, k, each one affected by a q and k (and also |q − k|)-dependent weight. 41,52,53It turns out that the maximum weight is obtained for k around q p /2 and q roughly in the range between q p /2 and q p , leading to the behavior observed in Figure 11.
Recently, it has been observed 53−60 that some liquid metals exhibit the simultaneous appearance of such high-frequency transverse branch and the presence of a marked high-frequency peak/shoulder in the spectra of its velocity autocorrelation function, namely, Z(ω).The latter is defined as the FT into the frequency domain of 2 , where the index 1 refers to a tagged ion in the fluid and ⟨••••⟩ stands for the ensemble average.As shown in Figure 11, the present results for l-La confirm this connection between the structure of the Z(ω) and the existence of one or two transverse dispersion branches.
There is not a complete theory that explains the common absence or presence of a high-frequency peak in Z(ω) and in C T (q, ω) for a certain q-range, especially because in case the peak exists in Z(ω), it is naturally assigned to the influence of longitudinal modes. 61Nevertheless, the possible existence of such a correlation has witnessed increasing computational evidence, and the present results for l-La engross the list of systems, where both features appear together.The self-diffusion coefficient, D, can be determined from the time integral of Z(t), which has led to a value of D = 0.22 ± 0.01 Å 2 /ps.There are no experimental self-diffusion data for l-La; nevertheless, Iida and Guthrie 62 inferred a value D = 0.34 Å 2 /ps, derived from some semiempirical expressions based on modified Stokes−Einstein-type formulas which relate the self-diffusion coefficient to other thermophysical magnitudes such as the density and the viscosity.
From the calculated C T (q, t), we also estimated the shear viscosity coefficient η, as follows. 41The memory function representation of C T (q, t), namely where the tilde denotes the Laplace transform and introduces a generalized shear viscosity coefficient η(q, z).The area under the normalized C T (q, t) gives = mC q z k T ( , 0)/( ) T B , and the extrapolation q → 0 of η(q, z = 0) ≡ η(q) gives the shear viscosity coefficient η.This procedure has been performed by fitting the Figure 9. Normalized transverse current correlation function, C T (q, t)/ C T (q, t = 0), of l-La at T = 1250 K for several q values.Figure 10.ML calculated spectral functions of the normalized longitudinal and transverse current correlation function, C L (q, ω) and C T (q, ω), respectively, for several q values.ML simulation results with the expression, η(q) = η/(1 + a 2 q 2 ), which was first proposed to analyze the results for a dense hardsphere system 1 but has also been applied to more complex systems. 63igure 12 shows the AIMD and ML simulation results for η(q) with the ML ones delving deeper into the small q-region.
We have evaluated its long-wavelength limit, yielding a value for the shear viscosity, η = 3.40 ± 0.25 (GPa ps) which is close to one of the two reported experimental data η = 2.65, 3.27. 9.2.2.Liquid Li 4 Pb Alloy.For a two-component system, the AL partial intermediate scattering functions, F ij (q, t), are defined as a straightforward generalization of the definition previously given for a one-component system.Along with the BT partial intermediate scattering functions, they provide information about the collective dynamics of the fluctuations in the partial number densities, componentwise in the case of the AL functions, and topological and chemical in the case of the BT ones.Their time FT yields the respective AL and BT partial dynamic structure factors, S ij (q, ω) (see refs 64 and 65 for more details).
Another interesting magnitude is the i-type component particle current because by evaluating the autocorrelation functions of its longitudinal, j i L (q,t), and transverse j i T (q, t) components, we obtain the partial longitudinal C ij L (q, t), and transverse, C ij T (q, t), current correlation functions, respectively. 64,65Moreover, it can be shown that the partial dynamic structure factors are related to the spectra of the partial longitudinal current correlation functions (both AL and BT), C ij L (q, ω), through the equation C ij L (q, ω) = (ω 2 /q 2 )S ij (q, ω).All of the abovementioned partial correlation functions have been evaluated by both AIMD and ML simulations, and below we report the more interesting results as yielded by the ML simulations.
Figure 13 shows, for several q-values, the AIMD and ML results for the number−number intermediate scattering function, F NN (q, t).This function represents the autocorrelation function of the total number density in the alloy, and therefore, its behavior is qualitatively similar to the intermediate scattering function in a one-component system; namely, it has an oscillatory structure at small q-values which is depleted by the effect of a decaying contribution.
The existence of propagating density fluctuations can be revealed by the appearance, within some q-range, of side peaks in the partial dynamic structure factors S ij (q, ω). Figure 14 shows, for the set of smallest attainable q-values, the obtained ML results for S LiLi (q, ω) and S NN (q, ω).The peaks are clearly visible, and with increasing q-values, its amplitude diminishes, and its position moves to higher frequencies.
From the frequencies associated with the positions of the side peaks, we have derived the corresponding dispersion curves, ω LiLi (q), ω PbPb (q), and ω NN (q), which are plotted in Figure 15.Notice that the ω LiLi (q) and ω NN (q) show two branches (highand low-frequency branches), whereas the ω PbPb (q) displays Figure 12. q-dependent shear viscosity, η(q) of l-La.The dashed line represents the fitting to the expression η(q) = η/(1 + a 2 q 2 ).
Figure 13.Normalized number−number intermediate scattering function, F NN (q, t)/F NN (q, t = 0), in the liquid Li 4 Pb alloy for several q values.Figure 14.ML simulation results for the partial dynamic structure factors, S LiLi (q, ω), (full lines) and S NN (q, ω), (broken lines) in the liquid Li 4 Pb alloy for several small q-values.
one, low-frequency, branch.The low-frequency ω LiLi (q) branch appears for a range of q ≤ 0.45 Å −1 and basically coincides with the ω PbPb (q) and ω NN (q) branches.All of them, when q → 0 smoothly merge into a hydrodynamic sound.On the other hand, the high-frequency branch of ω LiLi (q) begins at around the same q-value where the low frequency one vanishes and, in principle, would be compatible with a phase velocity of ≈5000 m/s which is similar to the value associated to the high-frequency mode revealed in the INS data of Alvarez et al. 20 Our results for the high-frequency branch of ω LiLi (q) show its appearance at q ≈ 0.45 Å −1 but do not allow us to conclude anything about its existence for smaller q values because it does not show up due to a small weight and/or, additionally, it is shielded by other, diffusive, modes.
The adiabatic sound velocity in the alloy, c s , has been obtained from the slope of ω NN (q → 0) with a value c s ≈ 2030 ± 150 m/s which agrees very well with the experimental data of ≈2000 m/s. 66From the previous value for c s , along with the formula = c ( ) , where ρ m is the mass density of the alloy, we have evaluated the adiabatic compressibility κ ad in the alloy, obtaining a value κ ad = 7.50 ± 0.20 (in units of 10 11 m 2 N −1 ) which is similar to the estimates suggested by Ruppersberg and Speicher. 67Moreover, taking into account that κ T = γ 0 κ ad , where γ 0 stands for the ratio of specific heats, and using the previous value for κ T , we obtain an estimate γ 0 ≈ 1.09 (±0.05) which compares well with the theoretical value of γ 0 ≈ 1.10 obtained by Bryk and Mryglod 17 in their GCM study of this alloy.
In order to gather additional information about the previous collective excitation modes, we calculated the spectra of the longitudinal current correlation functions, C ij L (q, ω).These functions allow the exposure of longitudinal modes that are not visible in the S ij (q, ω) because the appearance of the factor ω 2 attenuates the low-frequency modes (namely, the diffusive ones) and enhances the high frequency ones.
Figure 16 shows the ML calculated C ij L (q, ω) values for several small-q values.For q ≤ 0.26 Å −1 , the C PbPb L (q, ω) displays just one peak, whereas the C LiLi L (q, ω) and C NN L (q, ω) exhibit two peaks.At greater q's, all just show one peak, which is the high frequency one in the case of the LiLi and NN functions.On the other hand, C LiPb L (q, ω) has extrema that can be either positive or negative.Interestingly, in the region q ≤ 0.26 Å −1 , the positive peak in C LiPb L (q, ω) coincides in the position with the peaks in C LiLi L (q, ω) and C PbPb L (q, ω), which suggests the existence of an acoustic mode, where all the atoms/ions vibrate in phase.
From the positions of the peaks in C ij L (q, ω), we have obtained the longitudinal dispersion relations, ω ij L (q), which are plotted in Figure 17.There appear two branches for ω LiLi L (q) and ω NN L (q) and one branch for ω PbPb L (q).From the ML results, we can get full insights into the small q behavior of the high-frequency branch of ω LiLi L (q) which clearly tends to a finite value when q → 0; this result clearly shows the kinetic (nonhydrodynamic) nature of this excitation and discards previous claims concerning the fast sound character of this mode.
In order to make connection with the previous simulations that interpreted a merging of fast and slow sound into the hydrodynamic sound mode as q → 0, 16 we have calculated the phase velocities corresponding to the dispersions of the longitudinal currents, which are plotted in Figure 18.
As in ref 16, we observe a merging of two phase velocities into the hydrodynamic c s , one from below (corresponding to Pb) and the other one from above (corresponding to Li).However, it is clear that the higher one cannot be identified with the fast sound mode.The fast sound phase velocity, when q → 0, does not bend down abruptly to match the one that goes toward c s from above as inferred in ref 16.On the contrary, it increases toward infinity as a consequence of the optical character of the mode with a nonzero frequency at q = 0.
We have also evaluated the partial transverse current correlation functions as they provide information about the shear modes in the alloy and its viscosity.An important magnitude in the following discussion is the total transverse current correlation function C tt T (q, t) = ⟨j t T (q, t)j t T *(q, 0)⟩, where ML simulation results for the dispersion relations ω LiLi (q) (blue circles), ω PbPb (q) (red squares), and ω NN (q) (gray triangles) of the maxima in the partials S LiLi (q, ω), S PbPb (q, ω), and S NN (q, ω) for the Li 4 Pb liquid alloy.The slope of the dashed line corresponds to the experimental adiabatic sound velocity, whereas that of the full line corresponds to a phase velocity of 5000 m/s.(b) Closer insights into the low wavevector region.
Figure 16.ML simulation results for the partial longitudinal current correlation functions, C ij L (q, ω), in the liquid Li 4 Pb alloy.Full blue line: C LiLi L (q, ω), red line: C PbPb L (q, ω), and green line: C LiPb L (q, ω).The PbPb and LiPb components for q = 0.554 and 0.835 Å −1 have been multiplied by a factor of 5 to ease comparison.x m j q t x m j q t ( , ) ( , ) ( , ) is the total transverse current, m i (i = 1, 2) is the atomic mass of component i, and j i T (q, t) is the transverse component of the current defined in (8).
The C PbPb T (q,ω) and C tt T (q,ω) are found to have just a lowfrequency peak, whereas C LiLi T (q,ω) has a high-frequency peak and, for a limited q range, also another low-frequency one.The corresponding transverse dispersion relations, ω ii T (q), are shown in Figure 19.
The low-frequency transverse excitations correspond to shear waves that appear after a finite propagation gap.Near the edge of the propagation gap, the frequencies of the shear waves can be reasonably described by a viscoelastic model, 68 with where q c is the minimum wavevector for the existence of shear waves which in the present calculation gave the value q c = 0.401 Å −1 .For q values somewhat larger than q c the low-frequency branch can also be described by a linear form, ω tt T (q) ∼ c T (q − q T ), where the slope, c T , yields an estimate of the velocity of propagation of the shear modes in the alloy, c T ≈ 1250 ± 150 m/s.
The high-frequency branch of ω LiLi T (q) tends to a finite nonzero frequency in the long wavelength limit.This value is observed to coincide with the value taken by the longitudinal high-frequency branch, also shown in Figure 19, in this limit.This is the expected behavior for the optical longitudinal and transverse modes of a non-ionic binary liquid.Should unscreened Coulombic interactions between ions exist, then a gap between the longitudinal optic and transverse optic modes at q → 0 would appear. 68e finally report the alloy shear viscosity, which has been evaluated by a method that is just an extension of the onecomponent formulation (see eq 6) to the case of binary systems, using the total transverse current correlation function, C tt T (q,t), defined previously (for more details, consult refs 69 and 70).By resorting to the memory function representation of C tt T (q, t), we obtain a generalized shear viscosity coefficient, η(q, z), and the extrapolation to q → 0 of η(q) ≡ η(q, z = 0) gives the alloy shear viscosity.
Figure 20 shows the AIMD and ML results for η(q).The q → 0 limit of the ML results yields a value for the shear viscosity of η = 0.89 ± 0.05 (GPa ps).We are unaware of any experimental data to compare with, but there is some experimental trend according to which if an alloy exhibits heterocoordinating tendencies, as is the case in this then its shear viscosity usually shows a positive deviation from linearity.From the experimental values of the pure components at T = 1075 K, 9 i.e., η(Li) = 0.27 ± 0.05 and η(Pb) = 1.22 ± 0.10 (GPa ps), a linear variation gives a value that is half the present result.Therefore, a positive deviation from linearity clearly occurs, as might be expected from the strong heterocoordinating tendency of Li 4 Pb.

CONCLUSIONS
We have applied the ML method to develop accurate interatomic potentials for both l-La and the Li 4 Pb alloy.These potentials were trained on AIMD simulations and tested by comparing a range of static, dynamics, and transport properties with those previously calculated by the AIMD simulations.
The results reported for l-La constitute the first AIMD study of its static and dynamic properties, with previous studies employing only semiempirical models.We analyzed the static structure through the pair distribution function, g(r), and the static structure factor, S(q).We identified an asymmetric shape in the second peak in S(q) with a marked shoulder.This feature has been linked to the existence of icosahedral short-range order in the liquid.A further calculation using the CNA method revealed a significant presence of perfect and distorted   icosahedral structures, consistent with previous results linking the icosahedral structure to an asymmetric second peak in S(q).In addition, we have determined the generalized specific heat ratio, γ(q) by analyzing the relaxation processes behind the intermediate scattering function.The ML simulation results have proven to be crucial in revealing the low-q shape of γ(q) which appears qualitatively similar to the experimental γ(q) of l-Fe, as found by Hosokawa. 49Although no direct experimental data are available, the extrapolated value γ(q → 0) ≡ γ 0 agrees with some semiempirical estimates.
The calculated dynamic structure factors, S(q, ω), show side peaks that are indicative of collective density excitations.According to ML results, the associated dispersion relation shows the existence of a positive dispersion that could not be detected within the AIMD calculations.The transverse dispersion relation l-La exhibit two branches consistent with the suggested connection between the structure of the spectra of the VACF and the existence of one/two transverse dispersion branches.Results have also been reported for several transport coefficients, namely, the self-diffusion, adiabatic sound velocity, and shear viscosity coefficients.Taking into account the scarcity of data for most of these transport coefficients, we expect that the present results will be of interest for future experimental and theoretical research.
Despite all the previous work on the liquid Li 4 Pb alloy, there remained some controversial points about its properties, in particular about the nature of the high-frequency excitation associated with the LiLi dispersion relation, that have been addressed in detail in the present study.The reported ML simulations have delivered results for its static and dynamic properties with a precision matching the AIMD simulation results.However, the capability of the ML simulations to provide accurate information about the q → 0 behavior of any qdependent physical magnitude has been key to perform a very detailed investigation of the high-frequency excitation branch related to the LiLi dispersion relation.The ML simulation results show that this excitation is a kinetic (nonhydrodynamic) optic mode, and given the near ab initio precision of the present ML simulations, we can confidently rule out the fast sound character of this excitation.We have moreover reconciled the different views on the behavior of this alloy, showing that indeed, despite the kinetic optical nature of the high-frequency "fast sound" mode, a merging does exist between the phase velocities associated with the low-frequency peaks of the longitudinal currents associated with the heavy and light particles.
Overall, we have shown the power of ML interatomic potentials to analyze with higher accuracy multiple properties while, at the same time, finally discern the origin of some special features discovered in recent years.

Figure 1 .
Figure 1.Static structure factor, S(q), of l-La.Lines: present AIMD (dashed red) and ML (continuous black) calculations; open circles: XD data from ref 7. The inset shows a closer view of the second maximum.

Figure 2 .
Figure 2. Pair distribution function, g(r), of l-La.Lines: present AIMD (dashed red) and ML (continuous black) calculations; open circles: XD data from ref 7; and full diamonds: ND data for l-La at 1250 K from ref 6.

Figure 3 .
Figure 3. Ashcroft−Langreth (AL) partial pair distribution functions and static structure factors S ij (q) for the liquid Li 4 Pb alloy at 1075 K. Full blue, red, and green lines are the ML results, whereas dashed lines are the AIMD results.

Figure 5 .
Figure 5. Generalized specific heat ratio, γ(q), as obtained from the generalized hydrodynamic model (circles) and the generalized viscoelastic model (lozenges).The red and gray symbols are the AIMD and ML simulation results, respectively.

Figure 6 .
Figure 6.ML results for the dynamic structure factors, S(q, ω)/S(q), of l-La at T = 1250 K and several q values.

Figure 7 .
Figure 7. q-dependent adiabatic sound velocity of l-La at T = 1250 K.The red and gray symbols are the AIMD and ML simulation results, respectively.

Figure 8 .
Figure 8. Longitudinal dispersion relation of l-La, as obtained from the maxima in the spectra of the longitudinal current, C L (q, ω).

Figure 11 .
Figure 11.Transverse dispersion relation of l-La, as obtained from the maxima in the spectra of the transverse current, C T (q, ω).The blue dashed line represents the corresponding Z(ω).

Figure 15 .
Figure 15.ML simulation results for the dispersion relations ω LiLi (q) (blue circles), ω PbPb (q) (red squares), and ω NN (q) (gray triangles) of the maxima in the partials S LiLi (q, ω), S PbPb (q, ω), and S NN (q, ω) for the Li 4 Pb liquid alloy.The slope of the dashed line corresponds to the experimental adiabatic sound velocity, whereas that of the full line corresponds to a phase velocity of 5000 m/s.(b) Closer insights into the low wavevector region.Figure16.ML simulation results for the partial longitudinal current correlation functions, C ij L (q, ω), in the liquid Li 4 Pb alloy.Full blue line: C LiLi L (q, ω), red line: C PbPb

Figure 17 .
Figure 17.(a) ML simulation results for the longitudinal dispersion relations ω LiLi L (q) (blue circles), ω PbPb L (q) (red squares), and ω NN L (q) (gray triangles) for the Li 4 Pb liquid alloy.(b) Closer insights of the low wavevector region.

Figure 18 .
Figure 18.Phase velocities, ω(q)/q corresponding to the dispersion relations of the longitudinal currents.Symbols are the same as in the previous figure, and lines correspond to the hydrodynamic sound velocity and the fast sound velocity deduced from the peak positions of S LiLi (q, ω).

Figure 19 .
Figure 19.ML simulation results for the transverse dispersion relations ω LiLi T (q) (full and shaded blue circles), ω PbPb T (q) (shaded red squares), and ω tt T (q) (open green triangles), for the Li 4 Pb liquid alloy.The magenta diamonds correspond to the high-frequency longitudinal ω LiLi L (q).The continuous line represents the data obtained from eq 10.

Figure 20 .
Figure 20.q-dependent shear viscosity, η(q) of the liquid Li 4 Pb alloy.The dashed line represents the fitting to the expression mentioned in the text.