Vibrational properties of bcc U and Mo at different temperatures

An accurate description of the vibrational density of states is required for calculation of thermodynamic properties of crystal lattices as well as defect migration rates in the Vineyard theory framework. In this work we use EAM and MEAM interatomic potentials and apply the zero temperature lattice dynamics calculations and the velocity autocorrelation function calculations at finite temperatures to study vibrational density of states in bcc molybdenum and bcc uranium lattices. The latter case is especially interesting since γ-U is thermodynamically not favorable at low temperatures and is not observed in experiments.


Introduction
An adequate description of lattice vibrations in solids is necessary for calculation of their macroscopic properties including the equation of state in a wide range of temperatures and pressures, e.g. [1][2][3][4]. Vibrational properties are also crucial for microscopic description of defects dynamics in solids. The development of the corresponding atomistic approach was started with the Vineyard theory [5] that connected frequencies of collective atomic motions in a lattice with jumping rates of single point defects.
Radiation defects play the key role in the processes occurring in construction materials of nuclear reactors and nuclear fuels. Dynamics of defects determines mechanical properties of these materials, e.g. causes hardening, embrittlement, swelling, fission product release rates etc. Recent advances in the development of interatomic potentials (e.g. see [6,7]) result in the possibility to study properties of defects not only qualitatively but to address specific nuclear materials at the quantitative level. The corresponding examples of multiscale description of nuclear fuels based on atomistic theory have already appeared [8,9].
With the focus on metallic nuclear fuels [10][11][12] the EAM interatomic potential for uranium was developed [13]. This model was extended to the ternary U-Mo-Xe system [14] and used for simulation of displacement cascades [15], xenon bubbles in U-Mo matrix [16] and sputtering [17]. Later the non-central symmetric ADP potential for U-Mo system was shown to give somewhat higher accuracy [18].
The first non-central symmetric MEAM model for uranium was presented in [19]. It was deployed for simulation of cascades in γ-U [20]. Later it was extended into the interatomic model for U-Zr alloy [21]. Another MEAM model for uranium was presented by Fernández and Pascuet [22] and later extended into the U-Al model [23].
The charge optimized many body (COMB) formalism was used for the development of the interatomic potential for uranium in [24] and was later re-formulated and re-parameterized in order to describe adequately both α and γ lattices of uranium [25].
Despite this considerable progress in the development and deployment of the effective interatomic potentials for uranium and its alloys, vibrational properties of corresponding models of crystal lattices have not been analyzed as far as we know.
In this work we analyze the vibrational density of states (vDoS) for γ-uranium using EAM [13] and MEAM [19] models. The analysis of γ-uranium is not a straightforward procedure since the uranium bcc lattice is unstable at zero temperature [26]. The recently developed temperature dependent effective potential method opens the way to analyze the lattices dynamically stabilized at high temperatures [27]. However in this work we use less sophisticated approach combining lattice dynamics, structural relaxation and velocity autocorrelation functions calculations. For comparison we use the EAM model [14] of bcc molybdenum lattice that is stable at T = 0.

Model and calculation method 2.1. Interatomic potentials
In classical molecular dynamic (MD) calculations one of the key elements is a model of atomic interactions. For simulation of metals one of the most widely used models is the semi-empirical Embedded Atom Model (EAM) [28]. In this model the potential energy of a system is given as a sum of a pairwise interaction and an "embedding energy" The embedding energy is a function of the effective electron density at the atom position where electron density is a sum of densities from all neighbor atoms The functions ρ(R ij ), F (ρ) and U pair (R ij ) are fitted to experimental data and/or quantum simulations [6,7]. This model is rather universal and appropriate to simulation of pure metals and alloys. Although EAM potentials are good in many cases, a pure EAM model is spherically symmetric and cannot capture anisotropic properties of atomic electron density. For the more accurate treatment a Modified Embedded Atom Model (MEAM) [29] was developed. This model adds anisotropic terms to the EAM electron density: This model also introduces an explicit screening function which damps the interaction between high-order neighbors.
In this work we use the EAM potentials for uranium [13] and molybdenum [14] and the MEAM potential for uranium [19].

Simulation details
All simulations are carried out with the LAMMPS package [30]. The system of 5 3 , 7 3 and 15 3 bcc unit cells is simulated to assure that effect of size does not change the results critically. To eliminate boundary effects, periodic boundary conditions are used. All calculations are carried out at constant volume to avoid changes of vDoS due to lattice constant variations.
For calculation of vDoS of a static lattice we have written a code which uses LAMMPS as an engine to calculate the potential energy of the system. The program builds the Hessian matrix displacing all possible pairs of atoms and calculating the potential energy. Then it solves the eigenproblem (12) for this matrix to get eigenfrequencies.

Arrhenius law and Vineyard theory
Since one should simulate the evolution of defects at time scales compared to the period of operation of the reactor materials, the native molecular dynamic approach should be supplemented by other methods for higher time and length scales. As we want to keep the detailed consideration of atomic processes, a reasonable compromise is a kinetic approach. In this paper we will discuss only vacancies, keeping in mind that the same arguments can be applied to interstitials too.
Diffusion of vacancies can be described via the Transitional State Theory (TST) [31]. For rare events, the system considered is believed to be in thermodynamic equilibrium almost always. Then the probability of the reaction can be divided into two independent parts: the probability that particle reaches the "transition state", i.e. the saddle point of the potential surface, and the probability that the system crosses the ridge between the minima. Assuming the potential barrier of migration and the vibrational spectrum to be temperature-independent, the vacancy jump rate Γ should satisfy the Arrhenius equation: where E a is the free energy of activation of the hop from one lattice site to another, and ν is the effective frequency of hop attempts. There are several ways to estimate the pre-exponential factor in solids, but the most general is a Vineyard theory [5]. It considers the system of N particles showed schematically on the figure 1a. Firstly, we make the transformation from ordinary atomic coordinates x 1 . . . x 3N to collective coordinates y 2 . . . y 3N in such way that y 1 represents the process studied. It means that the initial and final states should lie on the y 1 axis. There are two minima A and B of the potential energy and the ridge hyperplane S between them. With the assumption that the potential energy surface is harmonic near the minima and almost harmonic near the saddle point P (i.e. harmonic along all coordinates y 2 . . . y 3N which are perpendicular to the reaction path) Vineyard showed that Γ can be written as where U A is the total potential energy in point A, U P is the energy in saddle point P ,ν j are all 3N normal frequencies of the lattice in the minimum,ν ′ j are 3N −1 frequencies in the hyperplane S.
Although Vineyard theory is pretty simple mathematically, its implementation in practice involves difficulties in both finding the energy barrier and the calculation ofν j andν ′ j . In this paper we concentrate on vDoS calculations for of bcc lattices of Mo and U. We discuss certain peculiarities of the latter case since bcc γ-phase of uranium is unstable at low temperatures.

Lattice dynamics approach
To build Vineyard model we need to calculate normal vibrational modes of the crystal at the equilibrium state (ν j ) and at the saddle point (ν ′ j ). The straightforward techique is to obtain the matrix of second derivatives of the potential energy (also known as a Hessian matrix) of the studied system [32]: where x 1 . . . x 3N are the cartesian coordinates of all atoms.
Adding masses to the model, we get a dynamical matrix: Solving the eigenproblem we obtain the normal frequencies as the square roots of the eigenvalues. The eigenvectorsx show which collective movement of the lattice corresponds to each frequency. If we calculate normal modes for a sufficiently large system, say, a few thousands of atoms, we get almost a continuous vDoS. The normal spectra for the molybdenum and uranium with two different potentials are shown in figures 2a, 3a and 4a. Here the significant difference between the uranium and the molybdenum becomes visible. Bcc phase of molybdenum is stable at zero temperature, so its normal frequency spectrum is completely real. It also means that no additional minimization needed, because the bcc lattice sites correspond to the absolute minima of the potential energy.
In contrast, bcc-U with the EAM potential [13] is unstable, so we see imaginary modes on its spectrum. But due to the limited sizes of the simulation box the bcc lattice remains with slight distortions only, and it is possible to explore its properties. The potential surface of the MEAM U potential seems to have a minimum at the ideal bcc lattice, but this minimum is a local one. To prove this fact, we replace standard minimization procedure with "annealing": we run dynamics on the almost cold lattice (at about 50 K) and after some time the system starts to self-heat and change its structure. To reach the convergence we artificially decrease velocities of atoms. After several iterations of this process, the system finds a new minimum, and the Hessian matrix is calculated again. The results are shown in the lower part of figure 4a. One can see now that standard minimization procedure does not have such a profound effect on the system and cannot escape the local minimum of energy.
This uncertainty makes the study of uranium lattice properties difficult that is why we try to find a more robust method to get vDoS.

Velocity autocorrelation function
There is a way to obtain vibrational properties of a system from the finite temperature dynamic simulation. If we calculate velocity autocorrelation function (VACF) for a long enough period, its Fourier transform gives us the vDoS for system studied [33]. For each potential a series of dynamic simulations at different temperatures are carried out. Resulting density spectra are given in figures 2b, 3b and 4b. One can see that molybdenum shows the good agreement between static and dynamic spectra, and the rise of temperature causes the broadening of peaks only.
The uranium spectra show different behavior. For the EAM model there are two different cases: at low temperatures the dynamic spectrum is very close to the static one for the minimized lattice, but has nothing similar to the ideal bcc spectrum. But as the temperature rises vDoS becomes bimodal. One can expect that at higher temperatures only the ideal bcc mode remains, because the bcc phase becomes dominant. For the MEAM potential the situation is the same: the low temperature spectrum is similar to the completely annealed structure static spectrum, and higher temperatures show resemblance to the spectrum of the ideal lattice.   Figure 4. (a) The normal frequency spectra of bcc uranium with the MEAM potential [19], obtained from the static simulation. The MEAM potential does not show imaginary modes, but the simple minimization procedure data (blue histogram) is not able to reach consistency with MD. The dynamic annealing data (green histogram) gives stable state with harder spectrum. (b) The vibrational density of states of bcc uranium with the MEAM potential, obtained via the Fourier transform of the velocity autocorrelation functions at different temperatures.

Conclusions
The vibrational density of states of bcc U and Mo have been analyzed using accurate interatomic potential models based on ab initio calculations. The normal modes calculations have been used for vDoS calculations at zero temperature. At higher temperatures vDoS has been calculated as a Fourier transform of VACF.
In the case of bcc Mo both methods give the perfectly consistent results. Therefore for bcc Mo it is possible to calculate Vineyard frequencies from static simulations.
The vDoS results for uranium have been obtained in EAM and MEAM models. The EAM and MEAM models show essentially different vDoS shapes. For the both U models the normal modes calculations become consistent with the low-temperature MD results only if the properly relaxed lattice is used instead of the ideal bcc structure.