Ab initio molecular dynamics study of collective dynamics in liquid Tl: Thermo-viscoelastic analysis

. We studied collective dynamics of pure liquid metal Tl using a com-bination of ab initio molecular dynamics (AIMD) simulations and a thermo-viscoelastic model applied to calculations of dynamic eigenmodes and dispersion of collective excitations in particular. We found that for liquid Tl at ambient pressure the transverse current spectral functions obtained directly in ab initio simulations for wave numbers larger than ﬁrst pseudo-Brillouin-zone boundary contain two low- and high-frequency peaks that is an evidence of emergence of the unusually high-frequency transverse modes as it was observed before in liquid Li at very high pressures. The thermo-viscoelastic dynamic model shows perfect reproduction of the simulation-derived longitudinal current auto-correlation functions, and the acoustic eigenmodes are in nice agreement with the peaks of the longitudinal current spectral functions up to the ﬁrst pseudo-Brillouin-zone boundary. The deviation of the dynamic eigenmodes from peak positions at higher wave numbers gives evidence of L-T coupling e ﬀ ects.


Introduction
Collective dynamics in disordered systems is one of the most fascinating topics of modern condensed matter physics. In particular the increasing amount of information from inelastic X-ray (IXS) scattering experiments and computer simulations on liquids opens new insight in understanding and new problems in description of collective dynamics in liquids on atomistic scale. The theory of collective model in liquids is working well in frames of linearized hydrodynamics, which is able to describe correctly dynamic processes in liquids, like sound propagation, thermal diffusivity, etc only on macroscopic scales where atomistic structure of liquid is not taken into account [1,2]. The simulations and scattering experiments, instead, allow probing the correlations on atomistic scale, and a challenge for generalized hydrodynamic theory is to describe emerging effects in liquid dynamics due to non-hydrodynamic effects on molecular scale like structural relaxation [3], cage effect, shear waves, or optic-like excitations in binary liquids [4] and molten salts [5] etc.
The standard generalized hydrodynamic approach is based on the memory function formalism [1,2] and for simple liquids it provides quite good description of time-dependent processes. Recently there appeared reports on IXS experiments on liquid metals Ga [6,7], Sn [8], Na [9], Zn, Cu and Fe [10] in which a contribution from transverse excitations to the experimental dynamic structure factors were discussed as a consequence of some longitudinaltransverse (L-T) coupling. From the simulation side there was nothing similar observed from classical simulations based on effective pair potentials between particle, however when the ab initio simulations were applied it was found that for liquid Sn the spectral functions can have two-peak structure [11]. Also, it was observed in AIMD that at very high pressure in liquid Li there emerged unusual high-frequency peaks on transverse current correlation functions [12], which give evidence of two low-and high-frequency contribution to the transverse dynamics at large wave numbers, which can be also a consequence of the suggested before L-T coupling in liquid metals. Later on similar findings were reported from AIMD study of collective excitations in liquid Fe [13]. At the date there is no clear explanation by theory of the suggested L-T coupling effects.
On the theoretical side there is practically absence of theoretical analysis of AIMDderived time correlation functions, while the dispersions of collective excitations are routinely estimated from the peak (or shoulder) positions of the L and T current spectral functions despite low statistics which is inherent in AIMD simulations. Therefore other schemes which would be based in particular on the analysis of time correlation functions in order to avoid numerical time Fourier-transformations are very timing. In particular those theoretical schemes of analysis should include temperature fluctuations in order to be consistent in the long-wavelength limit with the hydrodynamic theory as well as offer explanations based on emergence of non-hydrodynamic collective processes [14] outside the hydrodynamic region.
Recently an analysis of collective dynamics in liquid Na [15] implied contributions to the density-density time correlation functions coming from non-hydrodynamic heat waves [16], which exist on nanoscales and take part in heat transport in liquids and solids.
The aim of this study is to study features in longitudinal and transverse collective excitation spectra of liquid Tl (as a metallic system having similarly to Ga three valence electrons per ion) from ab initio simulations and rationalize them using thermo-viscoelastic approach of generalized hydrodynamics. The remaining paper has the following structure: in the next section we supply information on the ab initio simulations and analysis used in this study. Section III contains the spectra calculations and discussion on the obtained results. The conclusions of this study are given in the last section.

Simulations and methodology of thermo-viscoelastic analysis
We performed ab initio simulations of liquid Tl at T=577 K using a collection of 300 particles interacting with the electron subsystem via PAW-potentials [17,18]. The exchangecorrelation functional was taken in Perdew-Burke-Ernzerhof (PBE) form [19], and because of quite large size of simulated system only Γ-point of the Brillouin zone was taken for calculations of electron density. The initial configuration was taken from classical simulations of liquid Tl [20]. The time step in ab initio simulations was 2 fs, and production runs were over 39000 configurations. Each time step the positions, velocities and forces acting on particles were stored, as well as every 10 steps we stored the instantaneous configuration of the electron density. The ionic pair distribution functions g ii (r) and electron-ion distribution g ie (r) are shown in Fig.1. The obtained in ab initio simulations pair distribution function is very similar as the one from classical simulations [20] with the main peak located at 3.29Å. Although the temperature in the AIMD simulations was right at the melting point the pair distribution functions are typical as for the liquid systems, without any signs of supercooling on the second maximum or possible solid-like clusters. The advantage of AIMD simulations is in a possibility to study simultaneously with atomistic structure the distributions of electron density. The electron-ion correlations for valence electrons in liquid Tl are represented by the ion-electron distribution function g ie (r) (Fig.1), which in contrast to g ii (r) does not start from zero. The main maximum of the ion-electron distribution function is located at 0.93Å, that is also a consequence of the preserved nodal structure of electron wave functions in PAW formalism [17] in contrast to the nodeless pseudo-wave functions when the simple norm-conserving pseudopotentials are used in AIMD. The calculated static structure factor S (k) has the same location of the first sharp diffraction peak, however with larger amplitude, than the experimental static structure factor for liquid Tl [21]. This fact can be a consequence of overbinding in local structure due to the applied PBE exchange-correlation, similar as this is usually observed in AIMD simulations of water [22].
We calculated the density-density F nn (k, t) and current-current F L/T JJ (k, t) time correlation functions by averaging them over all possible directions of k-vectors having the same absolute value. The smallest wave number accessible from our AIMD simulations was k min = 0.301Å −1 . The simplest methodology to roughly estimate the frequencies of collective excitations is from the peak positions of the longitudinal (L) and transverse (T) current spectral functions C L/T (k, ω), which are obtained from numerical time-Fourier transformation of the corresponding simulation-derived current-current time correlation functions F L/T JJ (k, t). More precise analysis of collective dynamics and calculations of dispersion of collective excitations consists in application of generalized hydrodynamic models in order to estimate contributions from different relaxing and propagating modes to the time correlation functions (or spectral functions) of interest. In particular, in [23] a thermo-viscoelastic analysis based on generalized collective modes (GCM) approach in connection with the ab initio simulations was suggested. In the thermo-viscoelastic approach a five-variable dynamic model is considered: where n(k, t), J L (k, t) and e(k, t) are the spatial Fourier-components of particle density, density of longitudinal current and energy density, respectively. The overdots in 1 denote the time derivative of corresponding dynamic variable, which can be usually represented in analytical form [24]. The set of dynamic variables A (5) (k, t) is applied for solving the generalized Langevin equation [1] in matrix form in terms of dynamic eigenmodes for so-called generalized hydrodynamic matrix T(k) [24]. The corresponding eigenvectors allow to estimate contribution of particular dynamic mode to the time correlation function of interest, and the quality of the theory is usually verified by comparison of theoretical and simulation-derived time correlation functions [23,25]. An important issue is the number of exact sum rules fulfilled for theoretical time correlation functions (short-time behavior), and the thermo-viscoelastic dynamic model provides the same first time derivatives at the origin up to the fourth order for density-density time correlations functions from theory and simulations, fulfilling in general six sum rules within the GCM scheme. For the longitudinal current-current time correlation functions this GCM model provides six sum rules too [23]. The 5 × 5 generalized hydrodynamic matrix T(k) is expressed in frames of the GCM approach [24] in the following way via the 5 × 5 matrices of static correlation functions F(k, t = 0) and of Laplace-transformed time correlation functions in Markovian approximationF(k, z = 0). For each wave number Figure 1. Obtained in AIMD ion-ion pair distribution function g ii (r) (a), ion-electron distribution function g ie (r) (b) and static structure factor S(k) (c) for liquid Tl at 577 K. Experimental static structure factor [21] in (c) is shown by cross symbols.
one has to calculate from ab initio simulations the following matrix elements of matrices [23] F(k, where the indices in matrix elements correspond to dynamic variables from the set (1). For simplicity the purely imaginary matrix elements in (3) and (4) are shown explicitly. In ab initio simulations in contrast to classical molecular dynamics it is extremely difficult to sample dynamic variables of energy density and its first time derivative. Therefore the corresponding matrix elements in (3) and (4) which involve these two dynamic variables for calculation were treated as fitting parameters. All the matrix elements with n(k, t), J L (k, t) andJ L (k, t) were directly calculated from ab initio simulations. The fitting procedure was applied simultaneously to the AIMD-derived density-density and longitudinal current-current time correlation functions having just six fitting parameters in matrices (3) and (4). Having some trial values for the fitting parameters (matrix elements) at each iteration the eigenvalues and associated eigenvectors of the 5 × 5 generalized hydrodynamic matrix T(k) were obtained, the GCMtheoretical time correlation functions were constructed from the eigenvalues and associated eigenvectors [24] and compared to the AIMD-derived functions.

Results and discussion
The time correlation functions for liquid Tl obtained directly in AIMD simulations were analyzed by the thermo-viscoelastic dynamic model (1). The quality of the fit procedure is demonstrated in Fig. 2 for longitudinal current-current time correlation functions with two wave numbers. One can see perfect recovering by the thermo-viscoelastic theory the damping oscillations in the time dependence of current correlations. As it was stated in [23] the nice quality of the fit is provided by the exact sum rules for short-time behavior and time moments of the corresponding AIMD-derived function. The proposed fitting scheme in that sense has advantages over the fitting of spectral functions with spectral shapes of a single or two damped harmonic oscillators (DHO) [26]. We also performed the numerical time-Fourier transformation for AIMD-derived longitudinal and transverse time correlation functions in order to check their peak structure. The longitudinal current-current spectra C L (k, ω) showed a single-peak shape for all studied wave numbers. The transverse current spectral function had a well-defined single-peak shape for wave numbers less than the first pseudo-Brillouin zone boundary at ∼1.15Å −1 , however for larger wave numbers we observed a noisy two-peak structure (see Fig.3). Note that similar two-peak structure of the transverse current spectral function was observed before for liquid Li at very high pressure, while for ambient pressure in liquid Li the two-peak structure was absent [12]. Here we stress that we obtained for liquid Tl at the melting point and ambient pressure an evidence of the similar extra high-frequency excitations contributing to transverse dynamics.
The behavior of peak positions of L-and T-current spectral functions is shown in Fig.4. The dispersion of T-excitations shows a tendency very similar as was observed before for the case of liquid Li at high-pressure [12] when the quite standard dispersion of the shear waves splits into two low-and high-frequency branches upon wave number reaching the pseudo-Brillouin zone boundary. The low-frequency transverse branch practically overlaps with the peaks of the longitudinal current spectral function for wave numbers k > 1.7Å −1 that supports an idea of the existing L-T coupling on atomic scale [6]. The thermo-viscoelastic dynamic model results in dispersion for acoustic modes, which perfectly recovers the peak positions of the L-current spectral function but only up to the maximum of the dispersion curve. Further increase of wave numbers results in a high-frequency deviation of the acoustic eigenvalues from the peak positions. Namely in this region of wave numbers emerges the second, high-frequency, branch of transverse excitations. It is obvious that the thermoviscoelastic dynamic model does not contain any mechanism of possible L-T coupling, and perhaps therefore we observed the discrepancy between numerical and theoretical dispersion of longitudinal propagating modes. A further development of the thermo-viscoelastic dynamic model is needed in order to describe the L-T coupled spectra on atomic scale.

Conclusion
We performed ab initio simulations of liquid Tl at the melting point with the purpose to study dispersion of longitudinal and transverse collective excitations. We calculated the dispersion by two ways: from peak positions of the AIMD-derived longitudinal and transverse current spectral functions, and for the longitudinal case -using analysis of AIMD-derived time correlation functions by the five-variable thermo-viscoelastic dynamic model. We found a nice agreement between the two methodologies for the acoustic dispersion in the region of wave numbers up to the first pseudo-Brillouin zone boundary. For larger wave numbers the numerically obtained transverse current spectral functions gave evidence for existing two branches of collective excitations contributing to transverse dynamics. Namely in this range of wave numbers we observed a high-frequency deviation of the eigenvalues corresponding to acoustic modes from relevant peak positions of the Lcurrent spectral functions. This deviation is attributed to the existing L-T coupling effects, which were not taken into account in the thermo-viscoelastic model.
These results imply a need in elaborating dynamic models which would explicitly treat the coupling between longitudinal and transverse processes on atomic scales. However by date there is no clear understanding of what kind of microscopic processes are responsible for the L-T coupling, because the local L-T coupling is prohibited by symmetry reasons. These effects can emerge as a consequence of rotational motion of bonded units like dimers with short lifetime. However even in that case it is a challenge how to introduce the short-time bondings into generalized hydrodynamic theory. Further ab initio simulations should answer the question on the role of local bonding in liquid (not really simple) metals in collective dynamics on atomic scale.
The computing time allocation for the ab initio simulations in frames of the UA-GRIDproject is gratefully acknowledged. The calculations have been performed using the ab-initio total-energy and molecular dynamics program VASP (Vienna ab-initio simulation program) developed at the Institute für Materialphysik of the Universität Wien [27][28][29].